Homogenization for Generalized Langevin Equations with Applications to Anomalous Diffusion

Abstract

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

1 Introduction

1.1 Motivation

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

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

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

1.2 Definitions and Models

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

$$\begin{aligned} m \mathrm{d}\varvec{v}_t = \varvec{F}_0 \left( t,\varvec{x}_t,\varvec{v}_t,\varvec{\eta }_t \right) \mathrm{d}t + \varvec{F}_1\left( t, \{\varvec{x}_s,\varvec{v}_s\}_{s \in [0,t]}, \varvec{\xi }_t \right) \mathrm{d}t + \varvec{F}_e(t, \varvec{x}_t)\mathrm{d}t. \end{aligned}$$

(1.1)

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

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

• (i) $\varvec{F}_e$ is an external force field, which may be conservative (potential) or not.

• (ii) $\varvec{F}_0$ is a Markovian force of the form

  $$\begin{aligned} \varvec{F}_0\left( t, \varvec{x}_t,\varvec{v}_t, \varvec{\eta }_t\right) \mathrm{d}t = -\varvec{\gamma }_0(t, \varvec{x}_t)\varvec{v}_t \mathrm{d}t + \varvec{\sigma }_0(t, \varvec{x}_t)\mathrm{d}\varvec{W}^{(k)}_t, \end{aligned}$$

  (1.2)

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

• (iii) $\varvec{F}_1$ is a non-Markovian force of the form

  $$\begin{aligned}&\varvec{F}_1\left( t, \{\varvec{x}_s,\varvec{v}_s\}_{s \in [0,t]},\varvec{\xi }_t\right) \nonumber \\&\quad = - \varvec{g}(t, \varvec{x}_t) \left( \int _{0}^{t} \varvec{\kappa }(t-s) \varvec{h}(s, \varvec{x}_s) \varvec{v}_s \mathrm{d}s \right) + \varvec{\sigma }(t, \varvec{x}_t) \varvec{\xi }_t, \end{aligned}$$

  (1.3)

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

  $$\begin{aligned} \varvec{R}(t):= E [\varvec{\xi }_t \varvec{\xi }^{*}_0] = \varvec{R}^{*}(-t) \in {\mathbb {R}}^{r \times r}, \end{aligned}$$

  (1.4)

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

  $$\begin{aligned} \varvec{{\mathcal {S}}}(\omega ) = \int _{-\infty }^{\infty } \varvec{R}(t) e^{-i\omega t} \mathrm{d}t. \end{aligned}$$

  (1.5)

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

As the name itself suggests, GLEs are generalized versions of the Markovian Langevin equations, frequently employed to model physical systems. A basic form of the GLEs was first introduced by Mori in [53] and subsequently used in numerous statistical physics models [36, 72, 76]. The studies of GLEs have attracted increasing interest in recent years. We refer to, for instance, [24, 27, 39, 46, 47, 50, 69, 70, 75] for various applications of GLEs and [21, 40, 49, 56] for their asymptotic analysis. The main merit of GLEs from modeling point of view is that they take into account the effects of memory and the colored nature of noise on the dynamics of the system.

Remark 1.1

In general, there need not be any relation between $\varvec{\kappa }(t)$ and $\varvec{R}(t)$, or any relation between the damping coefficients and the noise coefficients appearing in the formula for $\varvec{F}_0$ and $\varvec{F}_1$. A particular but important case that we will revisit often in this paper is the case when a fluctuation-dissipation relation holds. In this case, $\varvec{\gamma }_0$ is proportional to $\varvec{\sigma }_0 \varvec{\sigma }_0^*$, $\varvec{h} = \varvec{g}^*$, $\varvec{g}$ is proportional to $\varvec{\sigma }$ and (without loss of generality¹) $\varvec{R}(t) = \varvec{\kappa }(t)$. Studies of microscopic Hamiltonian models for open classical systems lead to GLEs of the form (1.1) satisfying the above fluctuation-dissipation relation (see, for instance, Appendix A of [43] or [11]). On another note, GLEs of the form (1.1) are extended versions of the ones studied in our previous work [43]—here the GLEs are generalized to include a Markovian force, in addition to the non-Markovian one, as well as explicit time dependence in the coefficients.

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

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

Example 1

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

$$\begin{aligned} \mathrm{d}x_t&= v_t \mathrm{d}t, \end{aligned}$$

(1.6)

$$\begin{aligned} m \mathrm{d}v_t&= - \left( \int _0^t \kappa (t-s) v_s \mathrm{d}s\right) \mathrm{d}t + \xi _t \mathrm{d}t, \end{aligned}$$

(1.7)

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

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

• (M1) memory function of the bi-exponential form:

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

  (1.8)

  where the parameters satisfy $\Gamma _2> \Gamma _1 > 0$, and $\xi _t$ has the covariance function $R(t)= \kappa (t)$ and thus the spectral density,

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

  (1.9)

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

• (M2) memory function of the form:

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

  (1.10)

  where $\Gamma _1 > 0$, and $\xi _t$ has the covariance function $R(t)= \kappa (t)$ and thus the spectral density,

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

  (1.11)

  This model can be obtained from the one in (M1) by sending $\Gamma _2 \rightarrow \infty$ in the formula for $\kappa (t)$ in (1.8).

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

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

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

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

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

One of our main goals in this paper is to devise and study new homogenization procedures that yield reduced models retaining essential features of a more general class of models. This program is of importance for identification, parameter inference and uncertainty quantification of stochastic systems [26, 39, 46, 62] arising in the studies of anomalous diffusion [50, 52], climate modeling [23, 48] and molecular systems [10], among others. In particular, classical homogenization and averaging were performed within the Mori–Zwanzig formalism in [23]. There is increasing amount of effort striving to implement this or related programs, starting from microscopic models [61], using various techniques [7, 19, 20, 27, 58], for different systems of interest in the literature. The derived effective SDE models will be of particular interest for modelers of anomalous diffusion (see [22] for deterministic homogenization of anomalous-diffusive systems).

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

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

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

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

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

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

2 Application to One-Dimensional GLE Models

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

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

$$\begin{aligned} \mathrm{d}x_t&= v_t \mathrm{d}t, \end{aligned}$$

(2.1)

$$\begin{aligned} m \mathrm{d}v_t&= -g(x_t)\left( \int _0^t \kappa (t-s) h(x_s) v_s \mathrm{d}s\right) \mathrm{d}t + \sigma (x_t) \xi _t \mathrm{d}t + F_e(t,x_t)\mathrm{d}t, \end{aligned}$$

(2.2)

where

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

(2.3)

where $\Gamma _1 > 0$, and $\xi _t$ is the mean-zero stationary Gaussian process with the covariance function $R(t)=\kappa (t)$ and spectral density,

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

(2.4)

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

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

Corollary 2.1

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

$$\begin{aligned} \mathrm{d}X_t&= \frac{2}{\beta ^2 g h} F_e(t,X_t) \mathrm{d}t - \frac{2}{\beta h} Y_t \mathrm{d}t + S_1(X_t) \mathrm{d}t + \frac{2 \sigma }{\beta g h} (Z_t \mathrm{d}t + \mathrm{d}W_t), \end{aligned}$$

(2.5)

$$\begin{aligned} \mathrm{d}Y_t&= -\frac{\Gamma _1}{\beta g} F_e(t,X_t) \mathrm{d}t + S_2(X_t) \mathrm{d}t - \frac{\Gamma _1 \sigma }{g} (\mathrm{d}W_t + Z_t \mathrm{d}t), \end{aligned}$$

(2.6)

$$\begin{aligned} \mathrm{d}Z_t&= -\Gamma _1 Z_t \mathrm{d}t - \Gamma _1 \mathrm{d}W_t, \end{aligned}$$

(2.7)

where

$$\begin{aligned} S_1(X)&= \frac{2}{\beta ^2} \frac{\partial }{\partial X}\left( \frac{1}{g h} \right) \frac{\sigma ^2}{g h}, \ \ \ \ S_2(X) = -\frac{\Gamma _1}{\beta } \frac{\partial }{\partial X}\left( \frac{1}{g} \right) \frac{\sigma ^2}{g h}. \end{aligned}$$

(2.8)

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

$$\begin{aligned} \mathrm{d}X_t&= \frac{2}{\beta ^2 \phi ^2 } \frac{\partial }{\partial X}\left( \frac{1}{\sigma h}\right) \frac{\sigma }{h} \mathrm{d}t + \frac{2}{ \phi \sigma h \beta ^2} F_e(t,X_t) \mathrm{d}t - \frac{2}{\beta \phi h}U_t^{\phi } \mathrm{d}t + \frac{2}{\beta \phi h} \mathrm{d}W_t, \end{aligned}$$

(2.9)

$$\begin{aligned} \mathrm{d}U_t^\phi&= -\frac{\Gamma _1}{\beta \phi ^2} \frac{\partial }{\partial X}\left( \frac{1}{\sigma }\right) \frac{\sigma }{h} \mathrm{d}t - \frac{\Gamma _1}{ \beta \sigma } F_e(t,X_t) \mathrm{d}t, \end{aligned}$$

(2.10)

where $U_t^\phi = \phi Y_t-Z_t$.

The convergence is in the sense that for every $T>0$, $\sup _{t \in [0,T]} |x_t - X_t| \rightarrow 0$ in probability as $m \rightarrow 0$.

Proof

We apply Theorem 5.4 by setting $d=1, d_2 = d_4 = 2$, $\alpha _1 = \alpha _3 = 0$, $\alpha _2 = \alpha _4 = 1$, $\varvec{\gamma }_0 = \beta ^2 g h/2$, $\varvec{\sigma }_0 = \beta \sigma$, $\varvec{h} = h$, $\varvec{g} = g$, $\varvec{\sigma } = \sigma$, $\varvec{C}_2 = \varvec{C}_4 = \beta$, $\varvec{\Gamma }_2 = \Gamma _1$, $\varvec{M}_2 \varvec{C}_2^* = -\Gamma _1 \beta /2$, $\varvec{\Gamma }_4 = \Gamma _1$, $\varvec{\Sigma }_4 = -\Gamma _1$, and $\varvec{F}_e = F_e$. The assumptions of Theorem 5.4 can be verified in a straightforward way and so the results of the corollary follow. $\square$

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

$$\begin{aligned} \mathrm{d}x_t&= v_t \mathrm{d}t, \end{aligned}$$

(2.11)

$$\begin{aligned} m \mathrm{d}v_t&= -g(x_t)\left( \int _0^t \kappa (t-s) h(x_s) v_s \mathrm{d}s\right) \mathrm{d}t + \sigma (x_t) \xi _t \mathrm{d}t + F_e(t,x_t)\mathrm{d}t, \end{aligned}$$

(2.12)

where

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

(2.13)

with $\Gamma _2> \Gamma _1 > 0$, and $\xi _t$ is the mean-zero stationary Gaussian process with the covariance function $R(t)=\kappa (t)$ and spectral density,

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

(2.14)

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

For $\epsilon > 0$, we set $m = m_0 \epsilon$ and $\Gamma _2 = \gamma _2/\epsilon$ in (2.11) and (2.12), where $m_0$ and $\gamma _2$ are positive constants. This gives the family of equations:

$$\begin{aligned} \mathrm{d}x^\epsilon _t&= v^\epsilon _t \mathrm{d}t, \end{aligned}$$

(2.15)

$$\begin{aligned} m_0 \epsilon \mathrm{d}v^\epsilon _t&= -g(x^\epsilon _t)\left( \int _0^t \kappa ^\epsilon (t-s) h(x^\epsilon _s) v^\epsilon _s \mathrm{d}s\right) \mathrm{d}t + \sigma (x^\epsilon _t) \xi ^\epsilon _t \mathrm{d}t + F_e(t,x^\epsilon _t)\mathrm{d}t, \end{aligned}$$

(2.16)

where

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

(2.17)

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

Fig. 1

Plot of the memory function $\kappa (t)$ in (2.13) with $\Gamma _1 = 1$, $\beta = 1$ for different values of $\Gamma _2$

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

$$\begin{aligned} \kappa ^\epsilon (t) \rightarrow \kappa (t) := \frac{\beta ^2}{2} (\delta (t)-\Gamma _1 e^{-\Gamma _1|t|}) \end{aligned}$$

(2.18)

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

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

Corollary 2.2

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

$$\begin{aligned} \mathrm{d}X_t&= \frac{2}{\beta ^2 g h} F_e(t,X_t) \mathrm{d}t - \frac{2}{\beta h}Y_t \mathrm{d}t + S_1(X_t) \mathrm{d}t + \frac{2\sigma }{\beta g h } (\mathrm{d}W_t + Z_t \mathrm{d}t), \end{aligned}$$

(2.19)

$$\begin{aligned} \mathrm{d}Y_t&= -\frac{\Gamma _1}{\beta g} F_e(t,X_t) \mathrm{d}t + S_2(X_t) \mathrm{d}t - \frac{\Gamma _1 \sigma }{g}(\mathrm{d}W_t + Z_t \mathrm{d}t), \end{aligned}$$

(2.20)

$$\begin{aligned} \mathrm{d}Z_t&= -\Gamma _1 Z_t \mathrm{d}t - \Gamma _1 \mathrm{d}W_t, \end{aligned}$$

(2.21)

where $g=g(X_t)$, $h = h(X_t)$, $\sigma =\sigma (X_t)$, $W_t$ is a one-dimensional Wiener process, and

$$\begin{aligned} S_1&= \frac{2}{\beta ^2} \frac{\partial }{\partial X}\left( \frac{1}{gh}\right) \frac{\sigma ^2}{gh} - \frac{\partial }{\partial X}\left( \frac{1}{h}\right) \frac{4 \sigma ^2}{g(gh \beta ^2+4m_0\gamma _2)} \nonumber \\&\ \ \ \ + \frac{\partial }{\partial X}\left( \frac{\sigma }{gh}\right) \frac{4 \sigma }{ \beta ^2 g h +4m_0\gamma _2}, \end{aligned}$$

(2.22)

$$\begin{aligned} S_2&= -\frac{\Gamma _1}{\beta }\frac{\partial }{\partial X}\left( \frac{1}{g}\right) \frac{\sigma ^2}{g h} - \frac{\partial }{\partial X}\left( \frac{\sigma }{g}\right) \frac{2 \Gamma _1 \beta \sigma }{\beta ^2 g h +4m_0\gamma _2}. \end{aligned}$$

(2.23)

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

$$\begin{aligned} \mathrm{d}X_t&= \frac{2}{\beta ^2 \phi ^2 } \frac{\partial }{\partial X}\left( \frac{1}{\sigma h}\right) \frac{\sigma }{h} \mathrm{d}t + \frac{2}{ \phi \sigma h \beta ^2} F_e(t,X_t) \mathrm{d}t - \frac{2}{\beta \phi h}U_t^{\phi } \mathrm{d}t + \frac{2}{\beta \phi h} \mathrm{d}W_t, \end{aligned}$$

(2.24)

$$\begin{aligned} dU_t^\phi&= -\frac{\Gamma _1}{\beta \phi ^2} \frac{\partial }{\partial X}\left( \frac{1}{\sigma }\right) \frac{\sigma }{h} \mathrm{d}t - \frac{\Gamma _1}{ \beta \sigma } F_e(t,X_t) \mathrm{d}t, \end{aligned}$$

(2.25)

where $U_t^\phi = \phi Y_t-Z_t$.

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

Proof

Let $d=1$, $d_2 = d_4 = 2$ and denote the one-dimensional version of the variables, coefficients and parameters in Theorem 6.7 by non-bold letters (for instance, $x_t$, $B_2$, $\Gamma _{2,2}$ etc.). Furthermore, set $B_2 = B_4 = \beta > 0$, $\gamma _{2,2}=\gamma _{4,2}=\gamma _2 > 0$ and $\Gamma _{2,1}=\Gamma _{4,1}=\Gamma _1$. Then it can be verified that the assumptions of Theorem 6.7 hold and the results follow upon solving a Lyapunov equation.

$\square$

Remark 2.3

A few remarks on the contents of Corollary 2.2 follow.

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

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

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

$$\begin{aligned} \mathrm{d}X_t&= \frac{2}{\beta ^2 \sigma ^2} F_e(t,X_t) \mathrm{d}t - \frac{2}{\beta \sigma } U_t \mathrm{d}t + \frac{2}{\beta ^2}\frac{\partial }{\partial X_t}\left( \frac{1}{\sigma ^2} \right) \mathrm{d}t + \frac{2}{\beta \sigma } \mathrm{d}W_t, \end{aligned}$$

(2.26)

$$\begin{aligned} dU_t&= -\frac{\Gamma _1}{\beta \sigma } F_e(t,X_t) \mathrm{d}t - \frac{\Gamma _1}{\beta } \frac{\partial }{\partial X}\left( \frac{1}{\sigma } \right) \mathrm{d}t. \end{aligned}$$

(2.27)

To end this section, we remark that one could in principle repeat the above analysis for the case where the spectral density varies as $\omega ^{2l}$, for $l=2,4,\dots$ (i.e., the highly nonlinear case).

3 GLEs in Finite Dimensions

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

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

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

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

$$\begin{aligned} \varvec{\kappa }(t)= \alpha _1 \varvec{\kappa }_1(t) + \alpha _2\varvec{\kappa }_2(t) = \sum _{i=1}^2 \alpha _i \varvec{C}_i e^{-\varvec{\Gamma _i}|t|}\varvec{M}_i\varvec{C}_i^*, \end{aligned}$$

(3.1)

The noise process is defined as:

$$\begin{aligned} \varvec{\xi }_t = \alpha _3 \varvec{C}_3 \varvec{\beta }^3_t + \alpha _4 \varvec{C}_4 \varvec{\beta }^4_t, \end{aligned}$$

(3.2)

where the $\varvec{\beta }^{j}_t \in {\mathbb {R}}^{d_j}$ ($j=3,4$) are independent Ornstein–Uhlenbeck type processes, i.e., solutions of the SDEs:

$$\begin{aligned} \mathrm{d}\varvec{\beta }^j_t = -\varvec{\Gamma }_j \varvec{\beta }^j_t \mathrm{d}t + \varvec{\Sigma }_j \mathrm{d}\varvec{W}^{(q_j)}_t, \end{aligned}$$

(3.3)

with the initial conditions, $\varvec{\beta }^j_0$, normally distributed with mean-zero and covariance $\varvec{M}_j$. Here, $\varvec{W}^{(q_j)}_t$ denotes a $q_j$-dimensional Wiener process, independent of $\varvec{\beta }^j_0$. Also, the Wiener processes $\varvec{W}_t^{(q_3)}$ and $\varvec{W}_t^{(q_4)}$ are independent.

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

$$\begin{aligned} \varvec{\Gamma }_i \varvec{M}_i+\varvec{M}_i \varvec{\Gamma }_i^*=\varvec{\Sigma }_i \varvec{\Sigma }_i^*. \end{aligned}$$

(3.4)

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

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

$$\begin{aligned} \varvec{R}(t)=\alpha _3 \varvec{R}_3(t)+ \alpha _4 \varvec{R}_4(t) = \sum _{i=3}^4 \alpha _i \varvec{C}_i e^{-\varvec{\Gamma _i}|t|}\varvec{M}_i\varvec{C}_i^*, \end{aligned}$$

(3.5)

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

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

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

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

$$\begin{aligned} (\varvec{\Gamma }'_i=\varvec{T}_i \varvec{\Gamma }_i \varvec{T}^{-1}_i, \varvec{M}_i' = \varvec{T}_i \varvec{M}_i \varvec{T}_i^{*}, \varvec{C}'_i = \varvec{C}_i \varvec{T}_i^{-1}), \end{aligned}$$

(3.6)

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

Remark 3.1

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

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

Assumption 3.2

The memory function $\varvec{\kappa }(t)$ in the GLE (1.1) is a real-valued Bohl function defined by (3.1) and the noise process, $\{\varvec{\xi }_t, t \in {\mathbb {R}}^+ \}$, is a mean-zero, mean-square continuous, stationary Gaussian process with Bohl covariance function (hence, with rational spectral density), admitting a stochastic realization given by (3.2) and (3.3).

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

Definition 3.3

For $n \in {\mathbb {Z}}$, the nth order effective damping constant is defined as the constant matrix, parametrized by $\alpha _1, \alpha _2 \in \{0,1\}$:

$$\begin{aligned} \varvec{K}^{(n)}(\alpha _1,\alpha _2) := \alpha _1 \varvec{K}_1^{(n)} + \alpha _2 \varvec{K}_2^{(n)} \in {\mathbb {R}}^{q \times q}, \end{aligned}$$

(3.7)

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

$$\begin{aligned} \varvec{L}^{(n)}(\alpha _3,\alpha _4) := \alpha _3 \varvec{L}_3^{(n)} + \alpha _4 \varvec{L}_4^{(n)} \in {\mathbb {R}}^{r \times r}, \end{aligned}$$

(3.8)

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

Note that the first-order effective damping constant $\varvec{K}^{(1)}(\alpha _1,\alpha _2) = \int _0^{\infty } \varvec{\kappa }(t) \mathrm{d}t$ and the first-order effective diffusion constant $\varvec{L}^{(1)}(\alpha _3,\alpha _4) = \int _0^{\infty } \varvec{R}(t) \mathrm{d}t$ are simply the effective damping constant and effective diffusion constant introduced in [43]. The memory function and the covariance function of the noise process can be expressed in terms of these constants:

$$\begin{aligned} \varvec{\kappa }(t) = \sum _{i=1,2} \sum _{n=0}^{\infty } \alpha _i \frac{(-|t|)^n}{n!} \varvec{K}^{(-n)}_i, \ \ \ \varvec{R}(t) = \sum _{j=3,4} \sum _{n=0}^{\infty } \alpha _j \frac{(-|t|)^n}{n!} \varvec{L}^{(-n)}_j. \end{aligned}$$

(3.9)

Assumption 3.4

The matrix $\varvec{K}_1^{(1)}$ in the expression for first-order effective damping constant is invertible and the matrix $\varvec{K}_2^{(1)}$ equals zero. Similarly, in the expression for the first-order effective diffusion constant $\varvec{L}_3^{(1)}$, which is invertible, $\varvec{L}_4^{(1)} = \varvec{0}$.

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

Example 2

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

$$\begin{aligned} m \mathrm{d}\varvec{v}(t) = -\varvec{\gamma }_0 \varvec{v}(t) \mathrm{d}t -\varvec{g} \left( \int _0^t \varvec{\kappa }(t-s) \varvec{h} \varvec{v}(s) \mathrm{d}s \right) \mathrm{d}t + \varvec{\sigma }_0 \mathrm{d}\varvec{W}^{(k)}(t) + \varvec{\sigma } \varvec{\xi }(t) \mathrm{d}t. \end{aligned}$$

(3.10)

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

The following proposition gives the asymptotic behavior of the spectral densities (equivalently, covariance functions, or memory functions), the regularity² (in the mean-square sense) of the noise process, and, in the exactly solvable case of Example 2, the long-time mean-square displacement of the particle.

Proposition 3.5

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

• (i) We have $\varvec{{\mathcal {S}}}_3(\omega ) = O(1)$ as $\omega \rightarrow 0$. Also, let $k \ge 3$ be a positive odd integer and assume that $\varvec{L}_4^{(n)} = 0$ for $0< n < k$, where n is odd, and $\varvec{L}_4^{(k)} \ne 0$. Then $\varvec{{\mathcal {S}}}_{4}(\omega ) = O(\omega ^{k-1})$ as $\omega \rightarrow 0$. If there exists $h > 0$ such that the noise spectral density, $\varvec{{\mathcal {S}}}(\omega ) = O\left( \frac{1}{\omega ^{2h+1}}\right)$ as $\omega \rightarrow \infty$, then $\varvec{\xi }_t$ is n-times mean-square differentiable³ for $n < h$.

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

  $$\begin{aligned} E[\varvec{x}(t)\varvec{x}^*(t)]&= 2 \int _0^t \varvec{H}(s) \mathrm{d}s + 2m \left( \varvec{H}(t) {\mathcal {E}} \varvec{H}^*(t) - \int _0^t \varvec{H}(u) \dot{\varvec{H}^*}(u) \mathrm{d}u \right) \nonumber \\&\ \ \ \ + \int _0^t \varvec{H}(u) (\varvec{\sigma }_0 \varvec{\sigma }_0^* - 2 \varvec{\gamma }^*_0) \varvec{H}^*(u) \mathrm{d}u, \end{aligned}$$

  (3.11)

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

  $$\begin{aligned} \hat{\varvec{F}}(z) = (z^2(mz\varvec{I}+\varvec{\gamma }_0+ \varvec{g}\hat{\varvec{\kappa }}(z)\varvec{h}))^{-1}. \end{aligned}$$

  (3.12)

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

• (iii) Let $\alpha _1 = \alpha _3 = 1$ ($\alpha _i$, for $i=2,4$, can be 0 or 1 and $\varvec{F}_0$ can be zero or nonzero). Then $E[\varvec{x}(t)\varvec{x}^*(t)] = O(t)$ as $t \rightarrow \infty$, in which case we say that the particle diffuses normally.

• (iv) Let $\alpha _1 = 0$, $\alpha _2 = 1$ and $\varvec{F}_0 = \varvec{0}$ (the vanishing effective damping constant case). Then $E[\varvec{x}(t)\varvec{x}^*(t)] = O(t^{2})$ as $t \rightarrow \infty$, in which case we say that the particle exhibits a ballistic (super-diffusive) behavior.

Proof

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

  $$\begin{aligned} \varvec{{\mathcal {S}}}_i(\omega )&= \varvec{C}_i[(i\omega \varvec{I}+\varvec{\Gamma }_i)^{-1} + (-i\omega \varvec{I}+\varvec{\Gamma }_i)^{-1}]\varvec{M}_i \varvec{C}_i^* \end{aligned}$$

  (3.13)
  $$\begin{aligned}&= 2\varvec{C}_i [(i\omega \varvec{I}+\varvec{\Gamma }_i)^{-1} \varvec{\Gamma }_i (-i\omega \varvec{I}+\varvec{\Gamma }_i)^{-1}]\varvec{M}_i \varvec{C}_i^* \end{aligned}$$

  (3.14)
  $$\begin{aligned}&= 2\varvec{C}_i\varvec{\Gamma }_i^{-1}(\omega ^2 \varvec{\Gamma }_i^{-2} + \varvec{I})^{-1} \varvec{M}_i \varvec{C}_i^*, \end{aligned}$$

  (3.15)

  and so one has:

  $$\begin{aligned} \varvec{{\mathcal {S}}}_i(\omega )= & {} 2\varvec{C}_i\varvec{\Gamma }_i^{-1}\varvec{M}_i \varvec{C}_i^* - 2\varvec{C}_i\varvec{\Gamma }_i^{-3}\varvec{M}_i \varvec{C}_i^* \omega ^2 \nonumber \\&+ 2\varvec{C}_i\varvec{\Gamma }_i^{-5}\varvec{M}_i \varvec{C}_i^* \omega ^4 + \cdots , \end{aligned}$$

  (3.16)

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

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

  $$\begin{aligned} m \dot{\varvec{v}}(t)=-\varvec{\gamma }_0 \varvec{v}(t)+\varvec{\sigma }_0 \varvec{\eta }(t) -\varvec{g}\int _0^t \varvec{\kappa }(t-s) \varvec{h}\varvec{v}(s) \mathrm{d}s + \varvec{\sigma } \varvec{\xi }(t), \end{aligned}$$

  (3.17)

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

  $$\begin{aligned} z \hat{\varvec{x}}(z)&= \hat{\varvec{v}}(z), \end{aligned}$$

  (3.18)
  $$\begin{aligned} m (z \hat{\varvec{v}}(z) - \varvec{v}(0))&= -\varvec{g} \hat{\varvec{\kappa }}(z) \varvec{h} \hat{\varvec{v}}(z) - \varvec{\gamma }_0 \hat{\varvec{v}}(z) + \varvec{\sigma }_0 \hat{\varvec{\eta }}(z) + \varvec{\sigma } \hat{\varvec{\xi }}(z), \end{aligned}$$

  (3.19)

  and thus

  $$\begin{aligned} \hat{\varvec{x}}(z) = \hat{\varvec{H}}(z) (m \varvec{v}(0) + \varvec{\sigma }_0 \hat{\varvec{\eta }}(z) + \varvec{\sigma } \hat{\varvec{\xi }}(z)), \end{aligned}$$

  (3.20)

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

  $$\begin{aligned} \varvec{x}(t) = \varvec{H}(t) m \varvec{v} + \int _0^t \varvec{H}(t-s) (\varvec{\sigma }_0 \varvec{\eta }(s) + \varvec{\sigma } \varvec{\xi }(s)) \mathrm{d}s, \end{aligned}$$

  (3.21)

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

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

  $$\begin{aligned} E[\varvec{x}(t) \varvec{x}^T(t)]&= 2m \varvec{H}(t) {\mathcal {E}} \varvec{H}^*(t) + \int _0^t \varvec{H}(t-s) \varvec{\sigma }_0 \varvec{\sigma }_0^* \varvec{H}^*(t-s) \mathrm{d}s + \varvec{L}(t), \end{aligned}$$

  (3.22)

  where

  $$\begin{aligned} \varvec{L}(t)&= \int _0^t \mathrm{d}s \int _0^t \mathrm{d}u \ \varvec{H}(t-s) \varvec{\sigma } \varvec{\kappa }(|s-u|) \varvec{\sigma }^* \varvec{H}^*(t-u). \end{aligned}$$

  (3.23)

  To compute the double integral $\varvec{L}(t)$, we first rewrite it as $\varvec{L}(t) = \varvec{L}_1(t) + \varvec{L}_2(t)$, with

  $$\begin{aligned} \varvec{L}_1(t)&= \int _0^t \mathrm{d}s \ \varvec{H}(t-s) \int _s^t \mathrm{d}u \ \varvec{\sigma } \varvec{\kappa }(u-s) \varvec{\sigma }^* \varvec{H}^*(t-u), \end{aligned}$$

  (3.24)
  $$\begin{aligned} \varvec{L}_2(t)&= \int _0^t \mathrm{d}s \ \varvec{H}(t-s) \int _0^s \mathrm{d}u \ \varvec{\sigma } \varvec{\kappa }(s-u) \varvec{\sigma }^* \varvec{H}^*(t-u). \end{aligned}$$

  (3.25)

  We then compute:

  $$\begin{aligned} \varvec{L}_1(t)&= \int _0^t \mathrm{d}s \ \varvec{H}(t-s) \int _s^t d(t-u) \ \varvec{\sigma } \varvec{\kappa }(t-s-(t-u)) \cdot (-1) \varvec{\sigma }^* \varvec{H}^*(t-u), \end{aligned}$$

  (3.26)
  $$\begin{aligned}&= \int _0^t \mathrm{d}s \ \varvec{H}(t-s) \int _0^{t-s} d\tau \ \varvec{\sigma } \varvec{\kappa }(t-s-\tau ) \varvec{\sigma }^* \varvec{H}^*(\tau ), \end{aligned}$$

  (3.27)
  $$\begin{aligned}&= \int _0^t \mathrm{d}s \ \varvec{H}(t-s) (\varvec{\sigma } \varvec{\kappa } \varvec{\sigma }^* \star \varvec{H}^*)(t-s), \end{aligned}$$

  (3.28)
  $$\begin{aligned}&= \int _0^t \mathrm{d}u \ \varvec{H}(u) (\varvec{\sigma } \varvec{\kappa } \varvec{\sigma }^* \star \varvec{H}^*)(u), \end{aligned}$$

  (3.29)

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

  $$\begin{aligned} \varvec{L}_1(t) = \int _0^t \mathrm{d}u \ \varvec{H}(u) (\varvec{I} - m \dot{\varvec{H}^*}(u) - \varvec{\gamma }_0^* \varvec{H}^*(u)). \end{aligned}$$

  (3.30)

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

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

  $$\begin{aligned} E[\varvec{x}(t)\varvec{x}^*(t)]&= 2 \int _0^t \varvec{H}(s) \mathrm{d}s + 2m \left( \varvec{H}(t) {\mathcal {E}} \varvec{H}(t) - \int _0^t \varvec{H}(u) \dot{\varvec{H}}(u) \mathrm{d}u \right) . \end{aligned}$$

  (3.31)

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

  $$\begin{aligned} \hat{\varvec{H}}(z) = \frac{1}{z}\left[ mz\varvec{I} + \varvec{\gamma }_0 + \varvec{g}\sum _{i=1,2}\alpha _i \varvec{C}_i(z\varvec{I}+\varvec{\Gamma }_i)^{-1}\varvec{M}_i\varvec{C}_i^* \varvec{h} \right] ^{-1} \end{aligned}$$

  (3.32)

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

  $$\begin{aligned} \hat{\varvec{H}}(z)&\sim \frac{1}{z}\bigg [\varvec{\gamma }_0 + \alpha _1 \varvec{g} \varvec{K}_1^{(1)}\varvec{h} + \left( m\varvec{I}-\sum _{j=1,2}\alpha _j \varvec{g}\varvec{K}_j^{(2)}\varvec{h}\right) z \nonumber \\&\quad + \alpha _2 \varvec{g} \varvec{K}_2^{(3)}\varvec{h} z^2 + \alpha _2 \varvec{g} \varvec{K}_2^{(4)}\varvec{h} z^3 + \cdots \bigg ]^{-1}. \end{aligned}$$

  (3.33)

  Therefore, if $\varvec{\gamma }_0 = \varvec{\sigma }_0 \varvec{\sigma }_0^*/2$ is nonzero, then $\hat{\varvec{H}}(z) \sim 1/z$ as $z \rightarrow 0$. Otherwise, if in addition $\alpha _1 = 1$, then $\hat{\varvec{H}}(z) \sim 1/z$ as $z \rightarrow 0$, whereas if in addition $\alpha _1=0$, $\alpha _2=1$, then $\hat{\varvec{H}}(z) \sim 1/z^2$ as $z \rightarrow 0$. The results in (iii) and (iv) then follow by applying the Tauberian theorems [18], which say, in particular, that if $\hat{\varvec{H}}(z) \sim 1/z^\beta$ as $z \rightarrow 0$, then $\varvec{H}(t) \sim t^{\beta -1}$ as $t \rightarrow \infty$, for $\beta = 1, 2$ here. $\square$

Remark 3.6

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

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

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

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

$\varvec{z}_t := (\varvec{x}_{t}, \varvec{v}_{t}, \varvec{y}^1_{t}, \varvec{y}^2_t, \varvec{\beta }^3_{t}, \varvec{\beta }^4_{t}) \in {\mathbb {R}}^{d}\times {\mathbb {R}}^d \times {\mathbb {R}}^{d_1} \times {\mathbb {R}}^{d_2} \times {\mathbb {R}}^{d_3} \times {\mathbb {R}}^{d_4}$:

$$\begin{aligned} \mathrm{d}\varvec{x}_{t}&= \varvec{v}_{t} \mathrm{d}t, \end{aligned}$$

(3.34)

$$\begin{aligned} m \mathrm{d}\varvec{v}_{t}&= -\varvec{\gamma }_0(t, \varvec{x}_t) \varvec{v}_t \mathrm{d}t + \varvec{\sigma }_0(t, \varvec{x}_t) \mathrm{d}\varvec{W}_t^{(k)} - \varvec{g}(t, \varvec{x}_{t}) \sum _{i=1,2} \alpha _i \varvec{C}_i \varvec{y}^i_{t} \mathrm{d}t \nonumber \\&\ \ \ \ + \varvec{\sigma }(t, \varvec{x}_{t}) \sum _{j=3,4} \alpha _j \varvec{C}_j \varvec{\beta }^j_{t} \mathrm{d}t + \varvec{F}_e(t, \varvec{x}_{t})\mathrm{d}t, \end{aligned}$$

(3.35)

$$\begin{aligned} \mathrm{d}\varvec{y}^i_{t}&= -\varvec{\Gamma }_i \varvec{y}^i_{t} \mathrm{d}t + \varvec{M}_i \varvec{C}_i^* \varvec{h}(t,\varvec{x}_{t}) \varvec{v}_{t} \mathrm{d}t, \ \ i=1,2,\end{aligned}$$

(3.36)

$$\begin{aligned} \mathrm{d}\varvec{\beta }^j_{t}&= -\varvec{\Gamma }_j \varvec{\beta }^j_{t} \mathrm{d}t + \varvec{\Sigma }_j \mathrm{d}\varvec{W}^{(q_j)}_{t}, \ \ j=3,4, \end{aligned}$$

(3.37)

where we have defined the auxiliary memory processes:

$$\begin{aligned} \varvec{y}^i_{t} := \int _{0}^{t} e^{-\varvec{\Gamma }_i(t-s)} \varvec{M}_i \varvec{C}_i^* \varvec{h}(s,\varvec{x}_{s}) \varvec{v}_{s} \mathrm{d}s \in {\mathbb {R}}^{d_i}, \ \ i=1,2. \end{aligned}$$

(3.38)

Remark 3.7

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

4 On the Homogenization of Generalized Langevin Dynamics

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

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

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

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

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

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

$$\begin{aligned} (\varvec{f})_{\varvec{y}}(\varvec{y}) = (\varvec{\nabla }_{\varvec{y}}[f]_{jk}(\varvec{y}))_{j=1,\dots , n_2; k=1,\dots ,n_3}, \end{aligned}$$

(4.1)

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

5 Small Mass Limit of Generalized Langevin Dynamics

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

$$\begin{aligned} \mathrm{d}\varvec{x}_t^m&= \varvec{v}_t^m \mathrm{d}t, \end{aligned}$$

(5.1)

$$\begin{aligned} m \mathrm{d}\varvec{v}^m_t&= -\varvec{\gamma }_0(t, \varvec{x}_t^m) \varvec{v}_t^m \mathrm{d}t - \varvec{g}(t, \varvec{x}_t^m) \left( \int _0^t \varvec{\kappa }(t-s) \varvec{h}(s, \varvec{x}_s^m) \varvec{v}_s^m \mathrm{d}s \right) \mathrm{d}t \nonumber \\&\ \ \ \ + \varvec{\sigma }_0(t, \varvec{x}_t^m) \mathrm{d}\varvec{W}_t^{(k)} + \varvec{\sigma }(t, \varvec{x}_t^m) \varvec{\xi }_t \mathrm{d}t + \varvec{F}_e(t, \varvec{x}_t^m) \mathrm{d}t, \end{aligned}$$

(5.2)

where $\varvec{\kappa }(t)$ and $\varvec{\xi }_t$ are the memory function and noise process defined in (3.1) and (3.2), respectively, with each of the $\alpha _i$ ($i=1,2,3,4$) equal to zero or to one. Equations (5.1) and (5.2) are equivalent to the following system of SDEs for the Markov process $\varvec{z}^m_t := (\varvec{x}^m_{t}, \varvec{v}^m_{t}, \varvec{y}^{1,m}_{t}, \varvec{y}^{2,m}_t, \varvec{\beta }^{3,m}_{t}, \varvec{\beta }^{4,m}_{t}) \in {\mathbb {R}}^{d}\times {\mathbb {R}}^d \times {\mathbb {R}}^{d_1} \times {\mathbb {R}}^{d_2} \times {\mathbb {R}}^{d_3} \times {\mathbb {R}}^{d_4}$:

$$\begin{aligned} \mathrm{d}\varvec{x}^m_{t}&= \varvec{v}^m_{t} \mathrm{d}t, \end{aligned}$$

(5.3)

$$\begin{aligned} m \mathrm{d}\varvec{v}^m_{t}&= -\varvec{\gamma }_0(t, \varvec{x}^m_t) \varvec{v}^m_t \mathrm{d}t + \varvec{\sigma }_0(t, \varvec{x}^m_t) \mathrm{d}\varvec{W}_t^{(k)} - \varvec{g}(t, \varvec{x}^m_{t}) \sum _{i=1,2} \alpha _i \varvec{C}_i \varvec{y}^{i,m}_{t} \mathrm{d}t \nonumber \\&\ \ \ \ + \varvec{\sigma }(t, \varvec{x}^m_{t}) \sum _{j=3,4} \alpha _j \varvec{C}_j \varvec{\beta }^{j,m}_{t} \mathrm{d}t + \varvec{F}_e(t, \varvec{x}^m_{t})\mathrm{d}t, \end{aligned}$$

(5.4)

$$\begin{aligned} \mathrm{d}\varvec{y}^{i,m}_{t}&= -\varvec{\Gamma }_i \varvec{y}^{i,m}_{t} \mathrm{d}t + \varvec{M}_i \varvec{C}_i^* \varvec{h}(t, \varvec{x}^m_{t}) \varvec{v}^m_{t} \mathrm{d}t, \ \ i=1,2,\end{aligned}$$

(5.5)

$$\begin{aligned} \mathrm{d}\varvec{\beta }^{j,m}_{t}&= -\varvec{\Gamma }_j \varvec{\beta }^{j,m}_{t} \mathrm{d}t + \varvec{\Sigma }_j \mathrm{d}\varvec{W}^{(q_j)}_{t}, \ \ j=3,4, \end{aligned}$$

(5.6)

where we have defined the auxiliary memory processes:

$$\begin{aligned} \varvec{y}^{i,m}_{t} := \int _{0}^{t} e^{-\varvec{\Gamma }_i(t-s)} \varvec{M}_i \varvec{C}_i^* \varvec{h}(s, \varvec{x}^m_{s}) \varvec{v}^m_{s} \mathrm{d}s \in {\mathbb {R}}^{d_i}, \ \ i=1,2. \end{aligned}$$

(5.7)

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

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

Assumption 5.1

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

Assumption 5.2

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

Assumption 5.3

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

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

Theorem 5.4

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

$$\begin{aligned} \mathrm{d}\varvec{X}_t&= \varvec{\gamma }_0^{-1}(t, \varvec{X}_t)\bigg [ - \varvec{g}(t, \varvec{X}_t) \sum _{i=1}^2 \alpha _i \varvec{C}_i \varvec{Y}_t^i + \varvec{\sigma }(t, \varvec{X}_t) \sum _{j=3}^4 \alpha _j \varvec{C}_j \varvec{\beta }_t^j \nonumber \\&\ \ \ \ + \varvec{F}_e(t, \varvec{X}_t) \bigg ] \mathrm{d}t + \varvec{\gamma }_0^{-1}(t, \varvec{X}_t)\varvec{\sigma }_0(t, \varvec{X}_t)\mathrm{d}\varvec{W}_t^{(k)} +\varvec{S}^{(0)}(t, \varvec{X}_t)\mathrm{d}t, \end{aligned}$$

(5.8)

$$\begin{aligned} \mathrm{d}\varvec{Y}_t^k&= -\varvec{\Gamma }_k \varvec{Y}_t^k \mathrm{d}t + \varvec{M}_k \varvec{C}_k^* \varvec{h}(t, \varvec{X}_t)\varvec{\gamma }_0^{-1}(t, \varvec{X}_t) \bigg [ - \varvec{g}(t, \varvec{X}_t) \sum _{i=1}^2 \alpha _i \varvec{C}_i \varvec{Y}_t^i \nonumber \\&\ \ \ \ + \varvec{\sigma }(t, \varvec{X}_t) \sum _{j=3}^4 \alpha _j \varvec{C}_j \varvec{\beta }_t^j + \varvec{F}_e(t,\varvec{X}_t) \bigg ] \mathrm{d}t + \varvec{S}^{(k)}(t, \varvec{X}_t) \mathrm{d}t \nonumber \\&\ \ \ \ + \varvec{M}_k \varvec{C}_k^* \varvec{h}(t, \varvec{X}_t)\varvec{\gamma }_0^{-1}(t, \varvec{X}_t)\varvec{\sigma }_0(t, \varvec{X}_t)\mathrm{d}\varvec{W}_t^{(k)}, \ \ \text {for } k=1,2, \end{aligned}$$

(5.9)

$$\begin{aligned} \mathrm{d}\varvec{\beta }_t^l&= -\varvec{\Gamma }_l \varvec{\beta }_t^l \mathrm{d}t + \varvec{\Sigma }_l \mathrm{d}\varvec{W}_t^{(q_l)}, \ \ \text { for } l=3,4, \end{aligned}$$

(5.10)

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

$$\begin{aligned} S_i^{(0)}(t, \varvec{X})&= \frac{\partial }{\partial X_l}\left( (\varvec{\gamma }_0^{-1})_{ij}(t, \varvec{X}) \right) J_{lj}, \ \ j,l=1,\dots ,d, \end{aligned}$$

(5.11)

and for $k=1,2$,

$$\begin{aligned} S_i^{(k)}(t, \varvec{X})&= \frac{\partial }{\partial X_l}\left( (\varvec{M}_k \varvec{C}_k^* \varvec{h}(t, \varvec{X}) \varvec{\gamma }_0^{-1}(t, \varvec{X}))_{ij} \right) J_{lj}, \ \ j,l=1,\dots ,d, \end{aligned}$$

(5.12)

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

Proof

We prove the theorem by applying Theorem A.6. Using the notation in the statement of Theorem A.6, let $\epsilon = m$, $n_1 =d+d_1+d_2+d_3+d_4$, $n_2 = d$, $k_1 = q_3 + q_4$, $k_2 = k$, $\varvec{x}^\epsilon (t) = (\varvec{x}_t^m, \varvec{y}_t^{1,m}, \varvec{y}_t^{2,m}, \varvec{\beta }_t^{3,m}, \varvec{\beta }_t^{4,m})$, $\varvec{v}^\epsilon (t) = \varvec{v}_t^m$,

$$\begin{aligned} \varvec{a}_1&= [\varvec{I} \ \ \varvec{M}_1 \varvec{C}_1^* \varvec{h}(t, \varvec{x}_t^m) \ \ \varvec{M}_2 \varvec{C}_2^* \varvec{h}(t, \varvec{x}_t^m) \ \ \varvec{0} \ \ \varvec{0}], \end{aligned}$$

(5.13)

$$\begin{aligned} \varvec{a}_2&= -\varvec{\gamma }_0(t, \varvec{x}_t^m), \end{aligned}$$

(5.14)

$$\begin{aligned} \varvec{b}_1&= -(\varvec{0},\varvec{\Gamma }_1 \varvec{y}_t^{1,m}, \varvec{\Gamma }_2 \varvec{y}_t^{2,m}, \varvec{\Gamma }_3 \varvec{\beta }_t^{3,m}, \varvec{\Gamma }_4 \varvec{\beta }_t^{4,m}), \end{aligned}$$

(5.15)

$$\begin{aligned} \varvec{b}_2&= \varvec{F}_e(t, \varvec{x}_t^m) - \varvec{g}(t, \varvec{x}_t^m) \sum _{i=1,2} \alpha _i \varvec{C}_i \varvec{y}_t^{i,m} + \varvec{\sigma }(t, \varvec{x}_t^m) \sum _{j=3,4} \alpha _j \varvec{C}_j \varvec{\beta }_t^{j,m}, \end{aligned}$$

(5.16)

$$\begin{aligned} \varvec{\sigma }_1&= \begin{bmatrix} \varvec{0} &{}\quad \varvec{0} \\ \varvec{0} &{}\quad \varvec{0} \\ \varvec{0} &{}\quad \varvec{0} \\ \varvec{\Sigma }_3 &{}\quad \varvec{0} \\ \varvec{0} &{}\quad \varvec{\Sigma }_4 \end{bmatrix},\end{aligned}$$

(5.17)

$$\begin{aligned} \varvec{\sigma }_2&= \varvec{\sigma }_0(t, \varvec{x}_t^m), \end{aligned}$$

(5.18)

$\varvec{W}^{(k_1)}(t) = (\varvec{W}_t^{(q_3)}, \varvec{W}_t^{(q_4)})$ and $\varvec{W}^{(k_2)}(t) = \varvec{W}_t^{(k)}$. The initial conditions are $\varvec{x}(0) = (\varvec{x}, \varvec{0}, \varvec{0}, \varvec{\beta }_0^3, \varvec{\beta }_0^4)$ and $\varvec{v}(0) = \varvec{v}$, where $\varvec{\beta }_0^j$$(j=3,4$) are normally distributed with mean-zero and covariance $\varvec{M}_j$. They are independent of m.

Observe that in the above formula, $\varvec{a}_i$, $\varvec{b}_i$, $\varvec{\sigma }_i$ ($i=1,2$) do not depend explicitly on $\epsilon = m$, so by the convention adopted in Appendix A, we denote them $\varvec{A}_i$, $\varvec{B}_i$, $\varvec{\Sigma }_i$, respectively, and we put $a_i = b_i = c_i = d_i = \infty$, where $a_i, b_i, c_i, d_i$ are the rates in Assumption A.5.

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

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

Remark 5.5

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

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

• (ii) in the case when $\alpha _1$ and/or $\alpha _2$ are/is one, the limiting equation for the (slow) auxiliary memory variables contains correction drift terms ($\varvec{S}^{(1)}$ and/or $\varvec{S}^{(2)}$)—which could be called the memory-induced drifts. Interestingly, the memory-induced drifts disappear when $\varvec{h}$ is proportional to $\varvec{\gamma }_0$, a phenomenon that can be attributed to the interaction between the forces $\varvec{F}_0$ and $\varvec{F}_1$.

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

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

In this section, we consider the GLE (1.1), with $\varvec{F}_0 = \varvec{0}$, $\alpha _1=\alpha _3=0$, and $\alpha _2 = \alpha _4 = 1$. We explore a class of homogenization schemes, aiming to:

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

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

Remark 6.1

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

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

We assume that the $\varvec{\Gamma }_i$$(i=1,2,3,4)$ are already in the Jordan normal form and work in Jordan basis. Such form will reveal the slow-fast time scale structure of the system and so give us a rubric to develop homogenization schemes.

Assumption 6.2

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

$$\begin{aligned} \varvec{\Gamma }_i = \mathrm{diag}(\varvec{\Gamma }_{i,1},\dots ,\varvec{\Gamma }_{i,N_i}), \end{aligned}$$

(6.1)

where $N_i < d_i$, $\varvec{\Gamma }_{i,k} \in {\mathbb {R}}^{\nu (\lambda _{i,k}) \times \nu (\lambda _{i,k})}$ ($k=1, \dots , N_i$) is the Jordan block associated with the (controllable and observable) eigenvalue $\lambda _{i,k}$ (or time scale $\tau _{i,k}=1/\lambda _{i,k}$) and corresponds to the invariant subspace ${\mathcal {X}}_{i,k} = Ker(\lambda _{i,k}\varvec{I}-\varvec{\Gamma }_{i,k})^{\nu (\lambda _{i,k})}$, where $\nu (\lambda _{i,k})$ is the index of $\lambda _{i,k}$, i.e., the size of the largest Jordan block corresponding to the eigenvalue $\lambda _{i,k}$. Let $1 \le M_i < N_i$ and the eigenvalues be ordered as $0< \lambda _{i,1} \le \cdots \le \lambda _{i,M_i} < \lambda _{i,M_{i}+1} \le \cdots \le \lambda _{i,N_i}$, so that we have the invariant subspace decomposition, ${\mathbb {R}}^{d_i} = \bigoplus _{j=1}^{N_i} {\mathcal {X}}_{i,j}$, with $d_i = \sum _{k=1}^{N_i} \nu (\lambda _{i,k})$.

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

Algorithm 6.3

Procedure to study a class of homogenization problems.

• (1) Let $\alpha _1=\alpha _3 =0$, $\alpha _2=\alpha _4 =1$ and $\varvec{F}_0 = \varvec{0}$ in the GLE (1.1). Suppose that Assumption 6.2 holds and there exists $M_i$ such that $l_i = \sum _{k=1}^{M_i} \nu (\lambda _{i,k})$. Take this $M_i$.

• (2) For $i=2,4$, set $m=m' \epsilon$ and $\lambda _{i,k} = \lambda '_{i,k}/\epsilon$, for $k=M_i+1,\dots ,N_i$ (i.e., we scale the $(d_2-l_2)$ smallest memory time scales and the $(d_4-l_4)$ smallest noise correlation time scales with $\epsilon$), where $m'$ and the $\lambda '_{i,k}$ are positive constants.

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

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

We remark that while one has the above procedure to study homogenization schemes that achieve (P1) and (P2), the derivations and formulae for the limiting equations could become tedious and complicated as the $l_i$ and $d_i$ become large. Therefore, we consider a simple yet still sufficiently general instance of Algorithm 6.3 in the following.

Assumption 6.4

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

$$\begin{aligned} \varvec{\Phi }_i(z) = \varvec{Q}^{-1}_i(z) \varvec{P}_i(z), \end{aligned}$$

(6.2)

where the $\varvec{P}_i(z) \in {\mathbb {R}}^{p_i \times m_i}$ are matrix-valued monomials with degree $l_i$ :

$$\begin{aligned} \varvec{P}_i(z) = \varvec{B}_{l_i}z^{l_i} \end{aligned}$$

(6.3)

and the $\varvec{Q}_i(z) \in {\mathbb {R}}^{p_i \times p_i}$ are matrix-valued polynomials of degree $d_i$, i.e.,

$$\begin{aligned} \varvec{Q}_i(z) = \prod _{k=1}^{d_i} (z\varvec{I}+\varvec{\Gamma }_{i,k}). \end{aligned}$$

(6.4)

Here $p_2=q$, $p_4=r$, the $m_i$$(i=2,4$) are positive integers, the $\varvec{B}_{l_i} \in {\mathbb {R}}^{p_i \times m_i}$ are constant matrices, $\varvec{\Gamma }_{i,k}\in {\mathbb {R}}^{p_i \times p_i}$ are diagonal matrices with positive entries, and $\varvec{I}$ denotes identity matrix of appropriate dimension.

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

We assume $d_2$ and $d_4$ are even integers and consider in detail the case when $l_2 = l_4 = l = 1$, $d_2 = d_4 = h = 2$,

$$\begin{aligned} \varvec{\Gamma }_{2,1}&= \mathrm{diag}(\lambda _{2,1}, \dots , \lambda _{2,d_2/2}), \ \ \ \varvec{\Gamma }_{2,2}= \mathrm{diag}(\lambda _{2,d_2/2+1},\dots ,\lambda _{2,d_2}), \end{aligned}$$

(6.5)

$$\begin{aligned} \varvec{\Gamma }_{4,1}&=\mathrm{diag}(\lambda _{4,1},\dots ,\lambda _{4,d_4/2}), \ \ \ \varvec{\Gamma }_{4,2}= \mathrm{diag}(\lambda _{4,d_4/2+1},\dots ,\lambda _{4,d_4}), \end{aligned}$$

(6.6)

with $\lambda _{2,d_2} \ge \cdots \ge \lambda _{2,d_2/2+1}>\lambda _{2,d_2/2}\ge \cdots \ge \lambda _{2,1}>0$ and $\lambda _{4,d_4} \ge \cdots \ge \lambda _{4,d_4/2+1}>\lambda _{4,d_4/2}\ge \cdots \ge \lambda _{4,1}>0$ in Assumption 6.4, so that for $i=2,4$,

$$\begin{aligned} \varvec{\Gamma }_i = \mathrm{diag}(\varvec{\Gamma }_{i,1},\varvec{\Gamma }_{i,2}) \in {\mathbb {R}}^{d_i \times d_i}. \end{aligned}$$

(6.7)

We consider:

$$\begin{aligned} \varvec{C}_i&= [\varvec{B}_i \ \ \varvec{B}_i] \in {\mathbb {R}}^{p_i \times d_i}, \end{aligned}$$

(6.8)

$$\begin{aligned} \varvec{\Sigma }_i&= \left[ -\varvec{\Gamma }_{i,1}\varvec{\Gamma }_{i,2}(\varvec{\Gamma }_{i,2}-\varvec{\Gamma }_{i,1})^{-1} \ \ \ \varvec{\Gamma }_{i,2}^2(\varvec{\Gamma }_{i,2}-\varvec{\Gamma }_{i,1})^{-1} \right] ^* \in {\mathbb {R}}^{d_i \times d_i/2}, \end{aligned}$$

(6.9)

$$\begin{aligned} \text { so that } \nonumber \\ \varvec{M}_i&= \left[ \begin{array}{cc} \varvec{M}_i^{11} &{}\quad \varvec{M}_i^{12} \\ \varvec{M}_i^{21} &{} \quad \varvec{M}_i^{22} \end{array} \right] \in {\mathbb {R}}^{d_i \times d_i}, \end{aligned}$$

(6.10)

where

$$\begin{aligned} \varvec{M}_i^{11}&= \frac{1}{2}\varvec{\Gamma }_{i,1} \varvec{\Gamma }_{i,2}^2 (\varvec{\Gamma }_{i,1}-\varvec{\Gamma }_{i,2})^{-2}, \end{aligned}$$

(6.11)

$$\begin{aligned} \varvec{M}_i^{12}&= \varvec{M}_i^{21} = -\varvec{\Gamma }_{i,1} \varvec{\Gamma }^3_{i,2} (\varvec{\Gamma }_{i,1}+\varvec{\Gamma }_{i,2})^{-1} (\varvec{\Gamma }_{i,1}-\varvec{\Gamma }_{i,2})^{-2}, \end{aligned}$$

(6.12)

$$\begin{aligned} \varvec{M}_i^{22}&= \frac{1}{2}\varvec{\Gamma }^3_{i,2} (\varvec{\Gamma }_{i,1}-\varvec{\Gamma }_{i,2})^{-2}, \end{aligned}$$

(6.13)

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

$$\begin{aligned} \varvec{C}_ie^{-\varvec{\Gamma }_i|t|}\varvec{M}_i \varvec{C}_i^* = \frac{1}{2}\varvec{B}_i \varvec{\Gamma }_{i,2}^2(\varvec{\Gamma }_{i,2}^2-\varvec{\Gamma }_{i,1}^2)^{-1} \left( \varvec{\Gamma }_{i,2} e^{-\varvec{\Gamma }_{i,2} |t|} - \varvec{\Gamma }_{i,1} e^{-\varvec{\Gamma }_{i,1} |t|} \right) \varvec{B}_i^* \end{aligned}$$

(6.14)

and their Fourier transforms are:

$$\begin{aligned} \varvec{{\mathcal {S}}}_i(\omega )=\varvec{B}_i \varvec{\Gamma }_{i,2}^2 \varvec{B}_i^* \omega ^2 ((\omega ^2\varvec{I}+\varvec{\Gamma }_{i,1}^2)(\omega ^2\varvec{I}+\varvec{\Gamma }_{i,2}^2))^{-1}, \end{aligned}$$

(6.15)

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

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

$$\begin{aligned} m_0 \epsilon \mathrm{d}\varvec{v}_t^\epsilon&= -\varvec{g}(t, \varvec{x}_t^\epsilon ) \left( \int _0^t \varvec{\kappa }_2^\epsilon (t-s) \varvec{h}(s, \varvec{x}_s^\epsilon ) \varvec{v}_s^\epsilon \mathrm{d}s \right) \mathrm{d}t + \varvec{\sigma }(t, \varvec{x}_t^\epsilon ) \varvec{C}_4 \varvec{\beta }_t^{4,\epsilon } \mathrm{d}t \nonumber \\&\ \ \ \ + \varvec{F}_e(t, \varvec{x}_t^\epsilon ) \mathrm{d}t, \end{aligned}$$

(6.16)

$$\begin{aligned} \epsilon \mathrm{d}\varvec{\beta }_t^{4,\epsilon }&= -\varvec{\Gamma }_4 \varvec{\beta }_t^{4,\epsilon } \mathrm{d}t + \varvec{\Sigma }_4 \mathrm{d}\varvec{W}_t^{(q_4)}, \end{aligned}$$

(6.17)

where

$$\begin{aligned} \varvec{\kappa }_2^\epsilon (t) = \frac{1}{2}\varvec{B}_2 \varvec{B}_2^* \varvec{\gamma }_{2,2}^2(\varvec{\gamma }_{2,2}^2 - \epsilon ^2 \varvec{\Gamma }_{2,1}^2)^{-1} \left( \frac{\varvec{\gamma }_{2,2}}{\epsilon } e^{-\frac{\varvec{\gamma }_{2,2}}{\epsilon } |t|} - \varvec{\Gamma }_{2,1} e^{-\varvec{\Gamma }_{2,1} |t|} \right) \end{aligned}$$

(6.18)

and the covariance function of the noise process $\varvec{\xi }_t^\epsilon = \varvec{C}_4 \varvec{\beta }_t^{4,\epsilon }$ is given by

$$\begin{aligned} \varvec{R}_4^\epsilon (t) = \frac{1}{2}\varvec{B}_4 \varvec{B}_4^* \varvec{\gamma }_{4,2}^2(\varvec{\gamma }_{4,2}^2 - \epsilon ^2 \varvec{\Gamma }_{4,1}^2)^{-1} \left( \frac{\varvec{\gamma }_{4,2}}{\epsilon } e^{-\frac{\varvec{\gamma }_{4,2}}{\epsilon } |t|} - \varvec{\Gamma }_{4,1} e^{-\varvec{\Gamma }_{4,1} |t|} \right) . \end{aligned}$$

(6.19)

Note that $\varvec{\kappa }_2^\epsilon (t)$ and $\varvec{R}_4^\epsilon (t)$ converge (in the sense of distribution), as $\epsilon \rightarrow 0$, to

$$\begin{aligned} \frac{1}{2} \varvec{B}_i \varvec{B}_i^* (\delta (t)\varvec{I}-\varvec{\Gamma }_{i,1} e^{-\varvec{\Gamma }_{i,1} |t|}), \end{aligned}$$

(6.20)

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

$$\begin{aligned} \varvec{{\mathcal {S}}}_i(\omega ) = \varvec{B}_i\varvec{B}_i^*\omega ^2 (\omega ^2 \varvec{I}+\varvec{\Gamma }_{i,1}^2)^{-1}, \end{aligned}$$

(6.21)

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

Together with the equation for the particle’s position, Eqs. (6.16) and (6.17) form the SDE system:

$$\begin{aligned} \mathrm{d}\varvec{x}^\epsilon _t&= \varvec{v}^\epsilon _t \mathrm{d}t, \end{aligned}$$

(6.22)

$$\begin{aligned} \epsilon m_0 \mathrm{d}\varvec{v}^\epsilon _t&= -\varvec{g}(t, \varvec{x}^\epsilon _t)\varvec{B}_2(\varvec{y}_t^{2,1,\epsilon }+\varvec{y}_t^{2,2,\epsilon }) \mathrm{d}t + \varvec{\sigma }(t, \varvec{x}^\epsilon _t) \varvec{B}_4 (\varvec{\beta }_t^{4,1,\epsilon }+\varvec{\beta }_t^{4,2,\epsilon })\mathrm{d}t \nonumber \\&\ \ \ \ + \varvec{F}_e(t, \varvec{x}^\epsilon _t) \mathrm{d}t, \end{aligned}$$

(6.23)

$$\begin{aligned} \mathrm{d}\varvec{y}_t^{2,1,\epsilon }&= -\varvec{\Gamma }_{2,1}\varvec{y}_t^{2,1,\epsilon } \mathrm{d}t + {\mathcal {M}}_1^\epsilon \varvec{h}(t, \varvec{x}^\epsilon _t)\varvec{v}^\epsilon _t \mathrm{d}t, \end{aligned}$$

(6.24)

$$\begin{aligned} \epsilon \mathrm{d}\varvec{y}_t^{2,2,\epsilon }&= -\varvec{\gamma }_{2,2}\varvec{y}_t^{2,2,\epsilon }\mathrm{d}t+{\mathcal {M}}_2^\epsilon \varvec{h}(t, \varvec{x}^\epsilon _t) \varvec{v}^\epsilon _t \mathrm{d}t, \end{aligned}$$

(6.25)

$$\begin{aligned} \mathrm{d}\varvec{\beta }_t^{4,1,\epsilon }&= -\varvec{\Gamma }_{4,1} \varvec{\beta }_t^{4,1,\epsilon } \mathrm{d}t + \varvec{\sigma }_1^\epsilon \mathrm{d}\varvec{W}_t^{(q_4/2)}, \end{aligned}$$

(6.26)

$$\begin{aligned} \epsilon \mathrm{d}\varvec{\beta }_t^{4,2,\epsilon }&= -\varvec{\gamma }_{4,2} \varvec{\beta }_t^{4,2,\epsilon } \mathrm{d}t + \varvec{\sigma }_2^\epsilon \mathrm{d}\varvec{W}_t^{(q_4/2)}, \end{aligned}$$

(6.27)

where

$$\begin{aligned} {\mathcal {M}}_1^\epsilon&= \bigg ( (2(\epsilon \varvec{\Gamma }_{2,1}-\varvec{\gamma }_{2,2})^2)^{-1} \varvec{\Gamma }_{2,1}\varvec{\gamma }_{2,2}^2 \nonumber \\&\quad -((\epsilon \varvec{\Gamma }_{2,1}-\varvec{\gamma }_{2,2})^2(\epsilon \varvec{\Gamma }_{2,1} + \varvec{\gamma }_{2,2} ) )^{-1} \varvec{\Gamma }_{2,1} \varvec{\gamma }_{2,2}^3 \bigg ) \varvec{B}_2^* , \end{aligned}$$

(6.28)

$$\begin{aligned} {\mathcal {M}}_2^\epsilon&= \bigg ((2(\epsilon \varvec{\Gamma }_{2,1}-\varvec{\gamma }_{2,2})^2)^{-1} \varvec{\gamma }_{2,2}^3 \nonumber \\&\quad - \epsilon ((\epsilon \varvec{\Gamma }_{2,1}-\varvec{\gamma }_{2,2})^2(\epsilon \varvec{\Gamma }_{2,1} + \varvec{\gamma }_{2,2} ) )^{-1} \varvec{\Gamma }_{2,1} \varvec{\gamma }_{2,2}^3 \bigg )\varvec{B}_2^*, \end{aligned}$$

(6.29)

$$\begin{aligned} \varvec{\sigma }_1^\epsilon&= -(\varvec{\gamma }_{4,2}-\varvec{\Gamma }_{4,1}\epsilon )^{-1} \varvec{\Gamma }_{4,1} \varvec{\gamma }_{4,2}, \end{aligned}$$

(6.30)

$$\begin{aligned} \varvec{\sigma }_2^\epsilon&= (\varvec{\gamma }_{4,2}-\varvec{\Gamma }_{4,1} \epsilon )^{-1} \varvec{\gamma }_{4,2}^2. \end{aligned}$$

(6.31)

In the following, we take $\epsilon \in \mathcal {E} := (0,\epsilon _0], \epsilon _0 > 0$, to be small. We make the following assumptions, similar to those made in Theorem 5.4.

Assumption 6.5

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

Assumption 6.6

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

The following theorem describes the homogenized dynamics of the family of the GLEs (6.16) and (6.17) [or equivalently, of the SDEs (6.22)–(6.27)] in the limit $\epsilon \rightarrow 0$, i.e., when the inertial time scale, one half of the memory time scales and one half of the noise correlation time scales in the original generalized Langevin system tend to zero at the same rate.

Theorem 6.7

Consider the family of the GLEs (6.16) and (6.17) [or equivalently, of the SDEs (6.22)–(6.27)]. Suppose that Assumption 5.2 and Assumptions 6.4–6.6 hold, with the $\varvec{C}_i$, $\varvec{\Sigma }_i$, $\varvec{M}_i$ and $\varvec{\Gamma }_i$ ($i=2,4$) given in (6.7)–(6.13).

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

$$\begin{aligned} \varvec{I} + \varvec{g}(t, \varvec{x}) \tilde{\varvec{\kappa }}_\epsilon (\lambda ) \varvec{h}(t, \varvec{x})/\lambda m_0 \ \text { and } \ \varvec{I} + \varvec{g}(t, \varvec{x}) \tilde{\varvec{\kappa }}(\lambda ) \varvec{h}(t, \varvec{x})/\lambda m_0 \end{aligned}$$

(6.32)

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

$$\begin{aligned} \tilde{\varvec{\kappa }}_\epsilon (z) = \varvec{B}_2(z \varvec{I} + \varvec{\gamma }_{2,2})^{-1} {\mathcal {M}}_2^\epsilon \ \text { and } \ \tilde{\varvec{\kappa }}(z) = \frac{1}{2} \varvec{B}_2 (z\varvec{I} + \varvec{\gamma }_{2,2})^{-1} \varvec{\gamma }_{2,2} \varvec{B}_2^*. \end{aligned}$$

(6.33)

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

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

$$\begin{aligned} \mathrm{d}\varvec{\theta }_t&= \varvec{P}(t,\varvec{\theta }_t) \mathrm{d}t + \varvec{Q}(t, \varvec{\theta }_t) \mathrm{d}t + \varvec{R}(t, \varvec{\theta }_t) \mathrm{d}\varvec{W}_t^{(d_4/2)}, \end{aligned}$$

(6.34)

where

$$\begin{aligned} \varvec{P}(t,\varvec{\theta })= & {} \begin{bmatrix} \varvec{\nu }^{-1}(\varvec{F}_e-\varvec{g}\varvec{B}_2\varvec{Y}_t+\varvec{\sigma }\varvec{B}_4\varvec{Z}_t) \\ -\frac{1}{2} \varvec{\Gamma }_{2,1} \varvec{B}_2^* \varvec{h} \varvec{\nu }^{-1}(\varvec{F}_e-\varvec{g}\varvec{B}_2\varvec{Y}_t+\varvec{\sigma }\varvec{B}_4\varvec{Z}_t) - \varvec{\Gamma }_{2,1} \varvec{Y}_t \\ -\varvec{\Gamma }_{4,1} \varvec{Z}_t \end{bmatrix}, \qquad \end{aligned}$$

(6.35)

$$\begin{aligned} \varvec{R}(t, \theta )= & {} \begin{bmatrix} \varvec{\nu }^{-1} \varvec{\sigma } \varvec{B}_4 \\ -\frac{1}{2} \varvec{\Gamma }_{2,1} \varvec{B}_2^* \varvec{h} \varvec{\nu }^{-1} \varvec{\sigma } \varvec{B}_4 \\ -\varvec{\Gamma }_{4,1} \end{bmatrix}, \end{aligned}$$

(6.36)

and the ith component of $\varvec{Q}$, $i=1,\dots ,d+d_2/2+d_4/2$, is given by:

$$\begin{aligned} Q_i = \frac{\partial }{\partial X_l}\left[ H_{i,j}(t, \varvec{X}) \right] J_{j,l}, \ \ l=1,\dots ,d; \ j=1,\dots ,d+d_2/2+d_4/2, \end{aligned}$$

(6.37)

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

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

$$\begin{aligned} \varvec{T} = \begin{bmatrix} \varvec{I} &{}\quad \varvec{0} &{}\quad \varvec{0} \\ -\frac{1}{2}\varvec{\Gamma }_{2,1} \varvec{B}_2^* \varvec{h} &{}\quad \varvec{0} &{}\quad \varvec{0} \\ \varvec{0} &{}\quad \varvec{0} &{}\quad \varvec{0} \end{bmatrix}, \ \ \varvec{U} = \begin{bmatrix} \varvec{0} &{}\quad \varvec{g} \varvec{B}_2/m_0 &{}\quad -\varvec{\sigma } \varvec{B}_4/m_0 \\ -\frac{1}{2}\varvec{\gamma }_{2,2}\varvec{B}_2^* \varvec{h} &{}\quad \varvec{\gamma }_{2,2} &{}\quad \varvec{0} \\ \varvec{0} &{}\quad \varvec{0} &{}\quad \varvec{\gamma }_{4,2} \end{bmatrix}.\nonumber \\ \end{aligned}$$

(6.38)

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

Proof

We apply Theorem A.6 to the SDEs (6.22)–(6.27). To this end, we set, in Theorem A.6, $n_1 = n_2 = d+d_2/2+d_4/2$, $k_1 = k_2 = d_4/2$ and

$$\begin{aligned} \varvec{x}^\epsilon (t)= & {} (\varvec{x}^\epsilon _t, \varvec{y}_t^{2,1,\epsilon },\varvec{\beta }_t^{4,1,\epsilon }), \nonumber \\ \varvec{v}^\epsilon (t)= & {} (\varvec{v}^\epsilon _t, \varvec{y}_t^{2,2,\epsilon }, \varvec{\beta }_t^{4,2,\epsilon }) \in {\mathbb {R}}^{d+d_2/2+d_4/2}, \end{aligned}$$

(6.39)

$$\begin{aligned} \varvec{a}_1(t, \varvec{x}^\epsilon (t),\epsilon )= & {} \begin{bmatrix} \varvec{I} &{}\quad \varvec{0} &{}\quad \varvec{0} \\ {\mathcal {M}}_1^\epsilon \varvec{h}(t, \varvec{x}^\epsilon _t) &{}\quad \varvec{0} &{}\quad \varvec{0} \\ \varvec{0} &{}\quad \varvec{0} &{}\quad \varvec{0} \end{bmatrix} \in {\mathbb {R}}^{(d+d_2/2+d_4/2) \times (d+d_2/2+d_4/2)}, \nonumber \\\end{aligned}$$

(6.40)

$$\begin{aligned} \varvec{a}_2(t, \varvec{x}^\epsilon (t),\epsilon )= & {} \begin{bmatrix} \varvec{0} &{} \quad -\varvec{g}(t, \varvec{x}^\epsilon _t) \varvec{B}_2/m_0 &{}\quad \varvec{\sigma }(t, \varvec{x}^\epsilon _t)\varvec{B}_4/m_0 \\ {\mathcal {M}}_2^\epsilon \varvec{h}(t, \varvec{x}^\epsilon _t) &{}\quad -\varvec{\gamma }_{2,2} &{}\quad \varvec{0} \\ \varvec{0} &{}\quad \varvec{0} &{}\quad -\varvec{\gamma }_{4,2} \\ \end{bmatrix} \end{aligned}$$

(6.41)

$$\begin{aligned}&\ \ \ \ \in {\mathbb {R}}^{(d+d_2/2+d_4/2) \times (d+d_2/2+d_4/2)}, \nonumber \\ \varvec{b}_1(t, \varvec{x}^\epsilon (t),\epsilon )= & {} (\varvec{0}, -\varvec{\Gamma }_{2,1} \varvec{y}_t^{2,1,\epsilon }, -\varvec{\Gamma }_{4,1} \varvec{\beta }_t^{4,1,\epsilon }) \in {\mathbb {R}}^{d+d_2/2+d_4/2}, \end{aligned}$$

(6.42)

$$\begin{aligned} \varvec{b}_2(t,\varvec{x}^\epsilon (t),\epsilon )= & {} ((-\varvec{g}(t, \varvec{x}^\epsilon _t)\varvec{B}_2\varvec{y}_t^{2,1,\epsilon } + \varvec{\sigma }(t, \varvec{x}^\epsilon _t) \varvec{B}_4 \varvec{\beta }_t^{4,1,\epsilon }+\varvec{F}_e(t,\varvec{x}^\epsilon _t))/m_0, \nonumber \\&\ \ \ \ \ \ \varvec{0},\varvec{0}) \in {\mathbb {R}}^{d+d_2/2+d_4/2}, \end{aligned}$$

(6.43)

$$\begin{aligned} \varvec{\sigma }_1(t,\varvec{x}^\epsilon (t),\epsilon )= & {} [\varvec{0} \ \ \varvec{0} \ \ \varvec{\sigma }_1^\epsilon ]^* \in {\mathbb {R}}^{(d+d_2/2+d_4/2)\times d_4/2}, \end{aligned}$$

(6.44)

$$\begin{aligned} \varvec{\sigma }_2(t, \varvec{x}^\epsilon (t),\epsilon )= & {} [\varvec{0} \ \ \varvec{0} \ \ \varvec{\sigma }_2^\epsilon ]^* \in {\mathbb {R}}^{(d+d_2/2+d_4/2)\times d_4/2}. \end{aligned}$$

(6.45)

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

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

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

$$\begin{aligned} \varvec{M}_4^\epsilon = \begin{bmatrix} {\mathbb {E}}[|\varvec{\beta }_0^{4,1,\epsilon }|^2] &{}\quad {\mathbb {E}}[ \varvec{\beta }_0^{4,1,\epsilon } (\varvec{\beta }_0^{4,2,\epsilon })^*] \\ {\mathbb {E}}[\varvec{\beta }_0^{4,2,\epsilon } (\varvec{\beta }_0^{4,1,\epsilon })^*] &{}\quad {\mathbb {E}}[|\varvec{\beta }_0^{4,2,\epsilon }|^2] \end{bmatrix}, \end{aligned}$$

(6.46)

where

$$\begin{aligned} {\mathbb {E}}[|\varvec{\beta }_0^{4,1,\epsilon }|^2]= & {} \frac{1}{2}\varvec{\Gamma }_{4,1} \varvec{\gamma }_{4,2}^2 (\epsilon \varvec{\Gamma }_{4,1}-\varvec{\gamma }_{4,2})^{-2} = O(1), \end{aligned}$$

(6.47)

$$\begin{aligned} {\mathbb {E}}[\varvec{\beta }_0^{4,1,\epsilon } (\varvec{\beta }_0^{4,2,\epsilon })^*]= & {} {\mathbb {E}}[ \varvec{\beta }_0^{4,2,\epsilon } (\varvec{\beta }_0^{4,1,\epsilon })^*] \nonumber \\= & {} -\varvec{\Gamma }_{4,1} \varvec{\gamma }^3_{4,2} (\epsilon \varvec{\Gamma }_{4,1}+\varvec{\gamma }_{4,2})^{-1} (\epsilon \varvec{\Gamma }_{4,1}-\varvec{\gamma }_{4,2})^{-2}\nonumber \\= & {} O(1), \end{aligned}$$

(6.48)

$$\begin{aligned} {\mathbb {E}}[|\varvec{\beta }_0^{4,2,\epsilon }|^2]= & {} \frac{1}{2\epsilon }\varvec{\gamma }^3_{4,2} (\epsilon \varvec{\Gamma }_{4,1}-\varvec{\gamma }_{4,2})^{-2} = O\left( \frac{1}{\epsilon }\right) \end{aligned}$$

(6.49)

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

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

$$\begin{aligned} \varvec{\Sigma }_1&= [\varvec{0} \ \ \varvec{0} \ \ -\varvec{\Gamma }_{4,1} ]^*, \end{aligned}$$

(6.50)

$$\begin{aligned} \varvec{\Sigma }_2&= [\varvec{0} \ \ \varvec{0} \ \ \varvec{\gamma }_{4,2} ]^*, \end{aligned}$$

(6.51)

and $a_1 = a_2 = c_1 = c_2 = d_1 = d_2 = 1$, $b_1=b_2 = \infty$, where the $a_i$, $b_i$, $c_i$ and $d_i$ are from Assumption A.5 of Theorem A.6. Therefore, the first part of Assumption A.5 is satisfied.

It remains to verify the (uniform) Hurwitz stability of $\varvec{a}_2$ and $\varvec{A}_2$ (i.e., Assumption A.2 and the last part of Assumption A.5). This can be done using the methods of the proof of Theorem 2 in [43], and we omit the details here. The results then follow by applying Theorem A.6, and (6.34)–(6.38) follow from matrix algebraic calculations. $\square$

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

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

7 Conclusions and Final Remarks

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

The following summarizes the main conclusions of the paper:

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

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

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

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

Appendices

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

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

Let $n_1$, $n_2$, $k_1$, $k_2$ be positive integers. Let $\epsilon \in (0,\epsilon _0] =: {\mathcal {E}}$ be a small parameter and $\varvec{x}^{\epsilon }(t) \in {\mathbb {R}}^{n_1}$, $\varvec{v}^{\epsilon }(t) \in {\mathbb {R}}^{n_2}$ for $t \in [0,T]$, where $\epsilon _0>0$ and $T>0$ are finite constants. Let $\varvec{W}^{(k_1)}$ and $\varvec{W}^{(k_2)}$ denote independent Wiener processes, which are ${\mathbb {R}}^{k_1}$-valued and ${\mathbb {R}}^{k_2}$-valued, respectively, on a filtered probability space $(\Omega , {\mathcal {F}}, {\mathcal {F}}_t, {\mathbb {P}})$ satisfying the usual conditions [32].

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

$$\begin{aligned} \varvec{x}^{\epsilon }(t)&= ([x^{\epsilon }]_1(t),[x^{\epsilon }]_2(t),\dots , [x^{\epsilon }]_{n_1}(t)), \end{aligned}$$

(A.1)

$$\begin{aligned} \varvec{v}^{\epsilon }(t)&= ([v^{\epsilon }]_1(t),[v^{\epsilon }]_2(t),\dots , [v^{\epsilon }]_{n_2}(t)). \end{aligned}$$

(A.2)

We consider the following family of perturbed SDE systems⁴ for

$(\varvec{x}^\epsilon (t), \varvec{v}^\epsilon (t)) \in {\mathbb {R}}^{n_1+n_2}$:

$$\begin{aligned} \mathrm{d}\varvec{x}^{\epsilon }(t)&= \varvec{a}_{1}(t,\varvec{x}^{\epsilon }(t), \epsilon ) \varvec{v}^{\epsilon }(t) \mathrm{d}t + \varvec{b}_{1}(t,\varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}t + \varvec{\sigma }_{1}(t,\varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}\varvec{W}^{(k_1)}(t), \end{aligned}$$

(A.3)

$$\begin{aligned} \epsilon \mathrm{d}\varvec{v}^{\epsilon }(t)&= \varvec{a}_{2}(t,\varvec{x}^{\epsilon }(t),\epsilon ) \varvec{v}^{\epsilon }(t) \mathrm{d}t + \varvec{b}_{2}(t,\varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}t + \varvec{\sigma }_2(t,\varvec{x}^{\epsilon }(t), \epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(t), \end{aligned}$$

(A.4)

with the initial conditions, $\varvec{x}^{\epsilon }(0) = \varvec{x}^\epsilon$ and $\varvec{v}^{\epsilon }(0) = \varvec{v}^\epsilon$, where $\varvec{x}^\epsilon$ and $\varvec{v}^\epsilon$ are random variables that possibly depend on $\epsilon$. In the SDEs (A.3) and (A.4), the coefficients $\varvec{a}_1: {\mathbb {R}}^+ \times {\mathbb {R}}^{n_1} \times {\mathcal {E}} \rightarrow {\mathbb {R}}^{n_1 \times n_2}$, $\varvec{a}_2 : {\mathbb {R}}^+ \times {\mathbb {R}}^{n_1} \times {\mathcal {E}} \rightarrow {\mathbb {R}}^{n_2 \times n_2}$, $\varvec{\sigma }_2 : {\mathbb {R}}^+ \times {\mathbb {R}}^{n_1} \times {\mathcal {E}} \rightarrow {\mathbb {R}}^{n_2 \times k_2}$ are nonzero matrix-valued functions, whereas $\varvec{b}_1 : {\mathbb {R}}^+ \times {\mathbb {R}}^{n_1} \times {\mathcal {E}} \rightarrow {\mathbb {R}}^{n_1}$, $\varvec{b}_2 : {\mathbb {R}}^+ \times {\mathbb {R}}^{n_1} \times {\mathcal {E}} \rightarrow {\mathbb {R}}^{n_2}$, $\varvec{\sigma }_1 : {\mathbb {R}}^+ \times {\mathbb {R}}^{n_1} \times {\mathcal {E}} \rightarrow {\mathbb {R}}^{n_1 \times k_1}$ are matrix-valued or vector-valued functions, which may depend on $\varvec{x}^{\epsilon }$, as well as on t and $\epsilon$ explicitly, as indicated by the parentheses $(t, \varvec{x}^{\epsilon }(t), \epsilon )$. In the case where the coefficients do not depend on $\epsilon$ explicitly, we will denote them by the corresponding capital letters (for instance, if $\varvec{a}_i(t,\varvec{x},\epsilon )=\varvec{a}_i(t,\varvec{x})$, then $\varvec{a}_i(t,\varvec{x}) := \varvec{A}_i(t,\varvec{x})$ etc.).

We are interested in the limit as $\epsilon \rightarrow 0$ of the SDEs (A.3) and (A.4), in particular the limiting behavior of the process $\varvec{x}^{\epsilon }(t)$, under appropriate assumptions⁵ on the coefficients. In this appendix, we present a homogenization theorem that studies this limit and delay its proof to Appendix B.

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

Assumption A.1

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

Assumption A.2

The matrix-valued functions

$$\begin{aligned} \{ -\varvec{a}_2(t,\varvec{y}, \epsilon ); t \in [0,T], \varvec{y} \in {\mathbb {R}}^{n_1}, \epsilon \in {\mathcal {E}} \} \end{aligned}$$

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

Assumption A.3

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

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

$$\begin{aligned}&\max \bigg \{\Vert \varvec{a}_i(t, \varvec{x},\epsilon )-\varvec{a}_i(t,\varvec{y},\epsilon )\Vert ,\Vert (\varvec{a}_i)_{\varvec{x}}(t,\varvec{x},\epsilon )-(\varvec{a}_i)_{\varvec{x}}(t,\varvec{y},\epsilon )\Vert , \nonumber \\&\quad |\varvec{b}_i(t,\varvec{x},\epsilon )-\varvec{b}_i(t,\varvec{y},\epsilon )|, \Vert \varvec{\sigma }_i(t,\varvec{x},\epsilon )-\varvec{\sigma }_i(t,\varvec{y},\epsilon )\Vert ; \ i=1,2\bigg \} \nonumber \\&\le L(\epsilon )|\varvec{x}-\varvec{y}|. \end{aligned}$$

(A.5)

Assumption A.4

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

Assumption A.5

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

$$\begin{aligned} \Vert \varvec{a}_i(t,\varvec{x},\epsilon )-\varvec{A}_i(t,\varvec{x}) \Vert&\le \alpha _i(\epsilon ), \end{aligned}$$

(A.6)

$$\begin{aligned} |\varvec{b}_i(t,\varvec{x},\epsilon )-\varvec{B}_i(t,\varvec{x}) |&\le \beta _i(\epsilon ), \end{aligned}$$

(A.7)

$$\begin{aligned} \Vert \varvec{\sigma }_i(t,\varvec{x},\epsilon )-\varvec{\Sigma }_i(t,\varvec{x}) \Vert&\le \gamma _i(\epsilon ),\end{aligned}$$

(A.8)

$$\begin{aligned} \Vert (\varvec{a}_i)_{\varvec{x}}(t,\varvec{x},\epsilon )-(\varvec{A}_i )_{\varvec{x}}(t,\varvec{x}) \Vert&\le \theta _i(\epsilon ) \end{aligned}$$

(A.9)

where $\alpha _i(\epsilon ) = O(\epsilon ^{a_i})$, $\beta _i(\epsilon ) = O(\epsilon ^{b_i})$, $\gamma _i(\epsilon ) = O(\epsilon ^{c_i})$ and $\theta _i(\epsilon ) = O(\epsilon ^{d_i})$, as $\epsilon \rightarrow 0$, for some positive exponents $a_i$, $b_i$, $c_i$ and $d_i$. Moreover, we assume that $\varvec{A}_2(t,\varvec{x})$ is Hurwitz stable for every t and $\varvec{x}$.

Convention. In the case where the coefficients do not show explicit dependence on $\epsilon$ or in the case when any of the coefficients $\varvec{b}_1$, $\varvec{b}_2$ and $\varvec{\sigma }_1$ is zero, we set the exponent, describing the corresponding rate of convergence, to infinity. For instance, if $\varvec{a}_i(t,\varvec{x},\epsilon ) = \varvec{A}_i(t,\varvec{x})$, we set $a_i = \infty$. Meanwhile, if $\varvec{\sigma }_1 = \varvec{0}$, we set $c_1 = \infty$, etc.

We now state our homogenization theorem.

Theorem A.6

Suppose that the family of SDE systems (A.3) and (A.4) satisfies Assumptions A.1–A.5. Let $(\varvec{x}^{\epsilon }(t), \varvec{v}^{\epsilon }(t)) \in {\mathbb {R}}^{n_1} \times {\mathbb {R}}^{n_2}$ be their solutions, with the initial conditions $(\varvec{x}^\epsilon , \varvec{v}^\epsilon )$. Let $\varvec{X}(t) \in {\mathbb {R}}^{n_1}$ be the solution to the following Itô SDE with the initial position $\varvec{X}(0) = \varvec{x}$:

$$\begin{aligned} \mathrm{d}\varvec{X}(t)&= [\varvec{B}_1(t,\varvec{X}(t))-\varvec{A}_1(t,\varvec{X}(t))\varvec{A}_2^{-1}(t,\varvec{X}(t))\varvec{B}_2(t,\varvec{X}(t))] \mathrm{d}t \nonumber \\&\ \ \ \ + \varvec{S}(t,\varvec{X}(t)) \mathrm{d}t + \varvec{\Sigma }_1(t,\varvec{X}(t)) \mathrm{d}\varvec{W}^{(k_1)}(t) \nonumber \\&\ \ \ \ - \varvec{A}_1(t,\varvec{X}(t)) \varvec{A}_2^{-1}(t,\varvec{X}(t))\varvec{\Sigma }_2(t,\varvec{X}(t)) \mathrm{d}\varvec{W}^{(k_2)}(t), \end{aligned}$$

(A.10)

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

$$\begin{aligned}{}[S]_{i}(t,\varvec{X}) = -\frac{\partial }{\partial X_{l}} \bigg ([A_1 A_2^{-1}]_{i,j}(t,\varvec{X}) \bigg ) \cdot [A_1]_{l,k}(t,\varvec{X}) \cdot [J]_{j,k}(t,\varvec{X}), \end{aligned}$$

(A.11)

where $i,l=1,\dots ,n_1, \ j,k=1,\dots ,n_2$, or in index-free notation,

$$\begin{aligned} \varvec{S} = \varvec{A}_1 \varvec{A}_2^{-1} \varvec{\nabla }\cdot (\varvec{J}\varvec{A}_1^*) -\varvec{\nabla } \cdot (\varvec{A}_1 \varvec{A}_2^{-1} \varvec{J} \varvec{A}_1^*) , \end{aligned}$$

(A.12)

and $\varvec{J} \in {\mathbb {R}}^{n_2 \times n_2}$ is the unique solution to the Lyapunov equation:

$$\begin{aligned} \varvec{J} \varvec{A}_2^{*} + \varvec{A}_2 \varvec{J} = -\varvec{\Sigma }_2 \varvec{\Sigma }_2^{*}. \end{aligned}$$

(A.13)

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

$$\begin{aligned} {\mathbb {E}}\left[ \sup _{t \in [0,T]} |\varvec{x}^{\epsilon }(t) - \varvec{X}(t)|^p; \epsilon \le \epsilon _1 \right] = O(\epsilon ^{r}), \end{aligned}$$

(A.14)

in the limit as $\epsilon \rightarrow 0$, with $r>0$ is defined as:

$$\begin{aligned} r= {\left\{ \begin{array}{ll} \beta \ \text { for all } 0< \beta < \frac{p}{2}, &{} \text { if}\ a_i, b_i, c_i, d_i \ge \frac{1}{2} \text { for } i=1,2, \\ p \cdot \min (a_i, b_i, c_i, d_i; i=1,2) , &{} \text { otherwise}, \end{array}\right. } \end{aligned}$$

(A.15)

where the $a_i$, $b_i$, $c_i$, $d_i$ ($i=1,2$) are the positive constants from Assumption A.5. In particular, for all finite $T>0$,

$$\begin{aligned} \sup _{t \in [0,T]} |\varvec{x}^\epsilon (t) - \varvec{X}(t)| \rightarrow 0, \end{aligned}$$

(A.16)

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

Remark A.7

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

Appendix B: Proof of Theorem A.6

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

We start from an elementary calculus result.

Lemma B.1

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

$$\begin{aligned} \Vert \varvec{f}_i(\varvec{y},\epsilon )-\varvec{f}_i(\varvec{z},\epsilon )\Vert \le M_i(\epsilon )|\varvec{y}-\varvec{z}|, \end{aligned}$$

(B.1)

where $M_i(\epsilon )=O(1)$ as $\epsilon \rightarrow 0$.

• (i) Suppose that for each i and $\varvec{y} \in {\mathbb {R}}^{n}$, there exists a unique bounded $\varvec{F}_i(\varvec{y}):{\mathbb {R}}^n \rightarrow {\mathbb {R}}^{m_i \times n}$ and a constant $C_i>0$ such that $\Vert \varvec{f}_i(\varvec{y},\epsilon ) - \varvec{F}_i(\varvec{y})\Vert \le C_i \epsilon ^{r_i}$, for some positive constant $r_i$, as $\epsilon \rightarrow 0$ (i.e., the left-hand side is of order $O(\epsilon ^{r_i})$ as $\epsilon \rightarrow 0$). Then there exist constants D, $K_1, \dots , K_N >0$, such that

  $$\begin{aligned} \bigg \Vert \prod _{i=1}^{N} \varvec{f}_i(\varvec{y},\epsilon )-\prod _{i=1}^N \varvec{F}_i(\varvec{y})\bigg \Vert&\le K_1 \epsilon ^{r_1} + \cdots + K_N \epsilon ^{r_N} \le D \epsilon ^{\min (r_1, \dots , r_N)} \end{aligned}$$

  (B.2)
  $$\begin{aligned}&= O(\epsilon ^{\min (r_1, \dots , r_N)}), \end{aligned}$$

  (B.3)

  as $\epsilon \rightarrow 0$. If, in addition, $n=m_1$, $\varvec{f}_1(\varvec{y},\epsilon )$ and $\varvec{F}_1(\varvec{y})$ are invertible for every $\varvec{y} \in {\mathbb {R}}^n$ and $\epsilon > 0$, then $\Vert \varvec{f}_1^{-1}(\varvec{y},\epsilon )-\varvec{F}_1^{-1}(\varvec{y})\Vert = O(\epsilon ^{r_1})$ as $\epsilon \rightarrow 0$.

• (ii) Let $c_i \in {\mathbb {R}}$, $i=1,\dots ,N$. For every $\epsilon > 0$ and $\varvec{y} \in {\mathbb {R}}^n$, $\sum _{i=1}^{N} c_i \varvec{f}_i(\varvec{y},\epsilon )$ and $\prod _{i=1}^N c_i \varvec{f}_i(\varvec{y},\epsilon )$ are globally Lipschitz with a Lipschitz constant that is O(1) as $\epsilon \rightarrow 0$. Moreover, if $m_1=n$ and for every $\epsilon >0$, $\varvec{y} \in {\mathbb {R}}^n$, $\varvec{f}_1(\varvec{y},\epsilon )$ is invertible, then for every $\epsilon > 0$, $\varvec{y} \in {\mathbb {R}}^n$, $\varvec{f}^{-1}_1(\varvec{y},\epsilon )$ is globally Lipschitz in $\varvec{y}$ with a Lipschitz constant that is O(1) as $\epsilon \rightarrow 0$.

Proof

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

  $$\begin{aligned}&\bigg \Vert \prod _{i=1}^{k+1} \varvec{f}_i(\varvec{y},\epsilon )-\prod _{i=1}^{k+1} \varvec{F}_i(\varvec{y})\bigg \Vert \nonumber \\&\quad = \bigg \Vert \varvec{f}_{k+1}(\varvec{y},\epsilon )\cdot \prod _{i=1}^{k} \varvec{f}_i(\varvec{y},\epsilon )-\varvec{F}_{k+1}(\varvec{y})\cdot \prod _{i=1}^{k} \varvec{F}_i(\varvec{y})\bigg \Vert \end{aligned}$$

  (B.4)
  $$\begin{aligned}&\quad \le \Vert \varvec{f}_{k+1}(\varvec{y},\epsilon )\Vert \cdot \left\| \prod _{i=1}^{k} \varvec{f}_i(\varvec{y},\epsilon )-\prod _{i=1}^{k} \varvec{F}_i(\varvec{y})\right\| \nonumber \\&\qquad + \Vert \varvec{f}_{k+1}(\varvec{y},\epsilon )- \varvec{F}_{k+1}(\varvec{y})\Vert \cdot \left\| \prod _{i=1}^{k} \varvec{F}_i(\varvec{y})\right\| \end{aligned}$$

  (B.5)
  $$\begin{aligned}&\quad \le C ( D_k \epsilon ^{\min (r_1,\dots ,r_k)} + C_{k+1}\epsilon ^{r_{k+1}} ) \end{aligned}$$

  (B.6)
  $$\begin{aligned}&\quad \le C \max \{D_k, C_{k+1} \}(\epsilon ^{\min (r_1,\dots ,r_k)} + \epsilon ^{r_{k+1}}) \le D_{k+1} \epsilon ^{\min (r_1,\dots ,r_{k+1})}, \end{aligned}$$

  (B.7)

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

  $$\begin{aligned} \Vert \varvec{f}_1^{-1}(\varvec{y},\epsilon )-\varvec{F}_1^{-1}(\varvec{y})\Vert&= \Vert \varvec{f}_1^{-1}(\varvec{y},\epsilon ) (\varvec{F}_1(\varvec{y})-\varvec{f}_1(\varvec{y},\epsilon )) \varvec{F}_1^{-1}(\varvec{y}) \Vert \end{aligned}$$

  (B.8)
  $$\begin{aligned}&\le \Vert \varvec{f}_1^{-1}(\varvec{y},\epsilon )\Vert \cdot \Vert \varvec{F}_1(\varvec{y})-\varvec{f}_1(\varvec{y},\epsilon ) \Vert \cdot \Vert \varvec{F}_1^{-1}(\varvec{y}) \Vert \end{aligned}$$

  (B.9)
  $$\begin{aligned}&\le C \epsilon ^{r_1}, \end{aligned}$$

  (B.10)

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

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

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

$$\begin{aligned} \mathrm{d}\varvec{p}^{\epsilon }(t)&= \frac{\varvec{a}_{2}(t, \varvec{x}^{\epsilon }(t),\epsilon )}{\epsilon } \varvec{p}^{\epsilon }(t) \mathrm{d}t + \varvec{b}_{2}(t,\varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}t + \varvec{\sigma }_2(t,\varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(t). \end{aligned}$$

(B.11)

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

We need the following lemma, adapted from Proposition A.2.3 of [31], to obtain an exponential bound on fundamental matrix solutions of a linear equation.

Lemma B.2

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

$$\begin{aligned} \frac{\partial \varvec{\Phi }^\epsilon (t,s,\omega )}{\partial t} = \frac{\varvec{B}^\epsilon (t,\omega )}{\epsilon } \varvec{\Phi }^\epsilon (t,s,\omega ), \ \ \varvec{\Phi }^\epsilon (s,s,\omega )=\varvec{I}, \ \ 0 \le s \le t \le T. \end{aligned}$$

(B.12)

Then there exists a constant $C > 0$ and an (in general random⁶) $\epsilon _{1}=\epsilon _1(\omega )$ such that

$$\begin{aligned} \Vert \varvec{\Phi }^\epsilon (t,s,\omega )\Vert \le C e^{-\kappa (t-s)/\epsilon } \end{aligned}$$

(B.13)

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

Proof

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

$$\begin{aligned} \frac{\partial \varvec{\Phi }^\epsilon (t,s,\omega )}{\partial t} = \frac{\varvec{B}^\epsilon (u,\omega )}{\epsilon } \varvec{\Phi }^\epsilon (t,s,\omega ) + \frac{\varvec{B}^\epsilon (t,\omega )-\varvec{B}^\epsilon (u,\omega )}{\epsilon } \varvec{\Phi }^\epsilon (t,s,\omega ), \end{aligned}$$

(B.14)

and represent the solution to the IVP as:

$$\begin{aligned} \varvec{\Phi }^\epsilon (t,s,\omega ) = e^{(t-s)\frac{\varvec{B}^\epsilon (u,\omega )}{\epsilon }} + \frac{1}{\epsilon } \int _s^t e^{(t-r) \frac{\varvec{B}^\epsilon (u,\omega )}{\epsilon }} (\varvec{B}^\epsilon (r,\omega )-\varvec{B}^\epsilon (u,\omega )) \varvec{\Phi }^\epsilon (r,s,\omega ) dr. \end{aligned}$$

(B.15)

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

$$\begin{aligned}&\varvec{W}^\epsilon (t,s,\omega ) \nonumber \\&\quad = e^{\kappa (t-s)/\epsilon } e^{(t-s)\varvec{B}^\epsilon (t,\omega )/\epsilon } + \frac{1}{\epsilon } \int _s^t e^{\kappa (t-s)/\epsilon } e^{(t-r) \varvec{B}^\epsilon (t,\omega )/\epsilon } (\varvec{B}^\epsilon (r,\omega )-\varvec{B}^\epsilon (t,\omega )) \nonumber \\&\qquad \cdot \varvec{\Phi }^\epsilon (r,s,\omega ) dr \end{aligned}$$

(B.16)

$$\begin{aligned}&\quad = e^{\kappa (t-s)/\epsilon } e^{(t-s)\varvec{B}^\epsilon (t,\omega )/\epsilon } + \frac{1}{\epsilon } \int _s^t e^{\kappa (t-s)/\epsilon } e^{(t-r) \varvec{B}^\epsilon (t,\omega )/\epsilon } e^{-\kappa (r-s)/\epsilon } \nonumber \\&\qquad \cdot (\varvec{B}^\epsilon (r,\omega )-\varvec{B}^\epsilon (t,\omega )) \varvec{W}^\epsilon (r,s,\omega ) dr. \end{aligned}$$

(B.17)

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

$$\begin{aligned} \Vert e^{s \varvec{B}^\epsilon (t,\omega )/\epsilon }\Vert \le C e^{-2\kappa s/\epsilon } \end{aligned}$$

(B.18)

Using this, we obtain:

$$\begin{aligned}&\Vert \varvec{W}^\epsilon (t,s,\omega )\Vert \le C e^{-\kappa (t-s)/\epsilon } \nonumber \\&\qquad + \frac{C}{\epsilon } \int _s^t e^{-2 \kappa (t-r)/\epsilon } e^{-\kappa (r-s)/\epsilon } e^{\kappa (t-s)/\epsilon } \Vert \varvec{W}^\epsilon (r,s,\omega )\Vert \nonumber \\&\qquad \cdot \Vert \varvec{B}^\epsilon (r,\omega )-\varvec{B}^\epsilon (t,\omega )\Vert dr. \end{aligned}$$

(B.19)

This leads to the estimate:

$$\begin{aligned}&\sup _{s,t \in [0,T]} \Vert \varvec{W}^\epsilon (t,s,\omega )\Vert \le C + \sup _{r,s \in [0,T]}\Vert \varvec{W}^\epsilon (r,s,\omega )\Vert \cdot A_\epsilon (\omega ), \end{aligned}$$

(B.20)

where

$$\begin{aligned} A_\epsilon (\omega ) = \frac{C}{\epsilon } \sup _{t \in [0,T]} \int _0^t e^{- \frac{\kappa (t-r)}{\epsilon }} \left\| \varvec{B}^\epsilon (r,\omega )-\varvec{B}^\epsilon (t,\omega )\right\| dr. \end{aligned}$$

(B.21)

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

$$\begin{aligned} \sup _{s,t \in [0,T]} \Vert \varvec{W}^\epsilon (t,s,\omega )\Vert \le C + \frac{1}{2} \sup _{s,t \in [0,T]} \Vert \varvec{W}^\epsilon (t,s,\omega )\Vert \end{aligned}$$

(B.22)

for all $\epsilon \le \epsilon _1$. This implies that $\sup _{s,t \in [0,T]} \Vert \varvec{W}^\epsilon (t,s,\omega )\Vert \le 2C$, which is the claimed bound. $\square$

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

Lemma B.3

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

$$\begin{aligned}&\sup _{t \in [0,T]} \left| \varvec{\Phi }^\epsilon (t) \int _{0}^{t} (\varvec{\Phi }^\epsilon )^{-1}(s) \varvec{V}_{s} \mathrm{d}\varvec{H}_{s} \right| \nonumber \\&\quad \le C\left( 1+\frac{4}{\kappa } \sup _{s \in [0,T]} \Vert \varvec{B}^\epsilon (s)\Vert \right) \bigg (e^{-\kappa \delta /\epsilon } \sup _{t \in [0,T]} \left| \int _{0}^{t} \varvec{V}_{r} \mathrm{d}\varvec{H}_{r} \right| \nonumber \\&\qquad + \max _{k=0,1,\dots ,N-1} \sup _{t \in [k\delta , (k+2)\delta ]} \left| \int _{k\delta }^{t} \varvec{V}_{r} \mathrm{d}\varvec{H}_{r} \right| \bigg ) , \end{aligned}$$

(B.23)

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

Proof

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

$$\begin{aligned} \Vert \varvec{\Phi }^\epsilon (t) (\varvec{\Phi }^\epsilon )^{-1}(s)\Vert = \Vert \varvec{\Phi }^\epsilon (t,0) \varvec{\Phi }^\epsilon (0,s)\Vert = \Vert \varvec{\Phi }^\epsilon (t,s)\Vert \le C e^{-\kappa (t-s)/\epsilon }. \end{aligned}$$

(B.24)

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

$\square$

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

Proposition B.4

Suppose that Assumptions A.1–A.5 hold. For all $p \ge 1$, $T>0$, $0<\beta <p/2$, there exists a positive random variable $\epsilon _1$ such that:

$$\begin{aligned} {\mathbb {E}}\left[ \sup _{t \in [0,T]}|\varvec{p}^\epsilon (t)|^p; \epsilon \le \epsilon _1 \right] = O(\epsilon ^\beta ), \end{aligned}$$

(B.25)

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

$$\begin{aligned} {\mathbb {E}}\left[ \sup _{t \in [0,T]} \Vert \epsilon \varvec{v}^\epsilon (t)\varvec{v}^\epsilon (t)^*\Vert _{F}^p; \epsilon \le \epsilon _1 \right] = O(\epsilon ^{-\beta }), \end{aligned}$$

(B.26)

as $\epsilon \rightarrow 0$, where $\Vert \cdot \Vert _F$ denotes the Frobenius norm.

Proof

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

$$\begin{aligned} \frac{\partial \varvec{\Phi }_\epsilon (t)}{\partial t} = \frac{\varvec{a}_2(t, \varvec{x}^\epsilon (t), \epsilon )}{\epsilon } \varvec{\Phi }_\epsilon (t), \ \ \varvec{\Phi }_\epsilon (0) = \varvec{I}. \end{aligned}$$

(B.27)

Then,

$$\begin{aligned} \varvec{p}^\epsilon (t)&= \varvec{\Phi }_\epsilon (t)\epsilon \varvec{v}^\epsilon + \varvec{\Phi }_\epsilon (t) \int _0^t \varvec{\Phi }^{-1}_\epsilon (s) \varvec{b}_2(s,\varvec{x}^\epsilon (s),\epsilon )\mathrm{d}s \nonumber \\&\quad + \varvec{\Phi }_\epsilon (t) \int _0^t \varvec{\Phi }^{-1}_\epsilon (s) \varvec{\sigma }_2(s,\varvec{x}^\epsilon (s),\epsilon )\mathrm{d}\varvec{W}^{(k_2)}(s) \end{aligned}$$

(B.28)

$$\begin{aligned}&= \varvec{\Phi }_\epsilon (t)\epsilon \varvec{v}^\epsilon + \varvec{\Phi }_\epsilon (t) \int _0^t \varvec{\Phi }^{-1}_\epsilon (s) \varvec{B}_2(s,\varvec{x}^\epsilon (s))\mathrm{d}s \nonumber \\&\quad +\varvec{\Phi }_\epsilon (t) \int _0^t \varvec{\Phi }^{-1}_\epsilon (s) \left[ \varvec{b}_2(s,\varvec{x}^\epsilon (s),\epsilon ) -\varvec{B}_2(s,\varvec{x}^\epsilon (s)) \right] \mathrm{d}s \nonumber \\&\quad + \varvec{\Phi }_\epsilon (t) \int _0^t \varvec{\Phi }^{-1}_\epsilon (s) \varvec{\sigma }_2(s,\varvec{x}^\epsilon (s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s). \end{aligned}$$

(B.29)

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

$$\begin{aligned} \left| \sum _{i=1}^N a_i \right| ^p \le N^{p-1} \sum _{i=1}^N |a_i|^p \end{aligned}$$

(B.30)

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

$$\begin{aligned}&\sup _{t \in [0,T]}|\varvec{p}^\epsilon (t)|^p \nonumber \\&\quad \le 4^{p-1} \sup _{t \in [0,T]} \bigg [C^p e^{-\frac{\kappa p}{\epsilon }t} \epsilon ^p |\varvec{v}^\epsilon |^p + C^p \left( \int _0^t e^{-\frac{\kappa }{\epsilon }(t-s)} |\varvec{B}_2(s,\varvec{x}^\epsilon (s))| \mathrm{d}s \right) ^p \nonumber \\&\qquad + C^p \left( \int _0^t e^{-\frac{\kappa }{\epsilon }(t-s)} \bigg |[\varvec{b}_2(s,\varvec{x}^\epsilon (s),\epsilon ) - \varvec{B}_2(s,\varvec{x}^\epsilon (s))] \bigg | \mathrm{d}s \right) ^p \nonumber \\&\qquad + \bigg | \varvec{\Phi }_\epsilon (t) \int _0^t \varvec{\Phi }_\epsilon ^{-1}(s) \varvec{\sigma }_2(s,\varvec{x}^\epsilon (s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \bigg |^p \bigg ] \end{aligned}$$

(B.31)

$$\begin{aligned}&\quad \le 4^{p-1} \bigg (C^p \epsilon ^p |\varvec{v}^\epsilon |^p + \frac{C^p \epsilon ^p}{\kappa ^p} \bigg (\sup _{s \in [0,T]}|\varvec{B}_2 (s,\varvec{x}^\epsilon (s))|^p \nonumber \\&\qquad + \sup _{s \in [0,T]}|\varvec{b}_2(s,\varvec{x}^\epsilon (s),\epsilon )-\varvec{B}_2(s,\varvec{x}^\epsilon (s))|^p\bigg ) \nonumber \\&\qquad + \sup _{t \in [0,T]} \bigg | \varvec{\Phi }_\epsilon (t) \int _0^t \varvec{\Phi }_\epsilon ^{-1}(s) \varvec{\sigma }_2(s,\varvec{x}^\epsilon (s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \bigg |^p \bigg ), \end{aligned}$$

(B.32)

for $\epsilon \le \epsilon _1$, where $C>0$, $\kappa >0$, and $\epsilon _1>0$ is the random variable whose existence was proven in Lemma B.2.

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

$$\begin{aligned} \sup _{s \in [0,T]}|\varvec{b}_2(s,\varvec{x}^\epsilon (s),\epsilon )-\varvec{B}_2(s,\varvec{x}^\epsilon (s))|^p \le |\beta _2(\epsilon )|^p, \end{aligned}$$

(B.33)

where $\beta _2(\epsilon ) \le K \epsilon ^{b_2}$.

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

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

$$\begin{aligned}&{\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} |\varvec{p}^\epsilon (t)|^p\right] \nonumber \\&\quad \le C_1(p)(\epsilon ^\alpha + \epsilon ^{b_2 p} + \epsilon ^p) \nonumber \\&\qquad + C_2(p) {\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} \bigg | \varvec{\Phi }_\epsilon (t) \int _0^t \varvec{\Phi }_\epsilon ^{-1}(s) \varvec{\sigma }_2(s,\varvec{x}^\epsilon (s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \bigg |^p \right] , \end{aligned}$$

(B.34)

where $C_1(p), C_2(p) > 0$ are constants.

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

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

$$\begin{aligned}&{\mathbb {E}}_1\left[ \sup _{t \in [0,T]} \left| \varvec{\Phi }_{\epsilon }(t) \int _0^t \varvec{\Phi }_{\epsilon }^{-1}(s) \varvec{\sigma }_2(s,\varvec{x}^\epsilon (s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \right| ^p \right] \nonumber \\&\quad \le 2^{p-1} C^p {\mathbb {E}}_1 \left[ \left( 1 + \frac{4}{\kappa } \sup _{s \in [0,T]} \Vert \varvec{a}_2(s, \varvec{x}^\epsilon (s),\epsilon )\Vert \right) ^p \cdot \Pi \right] , \end{aligned}$$

(B.35)

$$\begin{aligned}&\quad \le 2^{p-1} C^p \left( 1 + \frac{4}{\kappa } \Vert \varvec{a}_2(t,\varvec{x}^\epsilon (t),\epsilon )\Vert _{\infty } \right) ^p \cdot {\mathbb {E}}_1 [\Pi ], \end{aligned}$$

(B.36)

where $\Vert \varvec{a}_2(t,\varvec{x}^\epsilon (t),\epsilon )\Vert _{\infty } := \sup _{ t \in [0,T], \varvec{y} \in {\mathbb {R}}^{n_1},\epsilon \in {\mathcal {E}}} \Vert \varvec{a}_2(t,\varvec{y},\epsilon )\Vert$ and

$$\begin{aligned} \Pi&= e^{-p \delta \kappa /\epsilon } \sup _{t \in [0,T]} \bigg | \int _0^t \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \bigg |^p \nonumber \\&\quad + \max _{k=0,\dots ,N-1} \sup _{t \in [k \delta , (k+2)\delta ]} \bigg | \int _{k \delta }^t \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \bigg |^p. \end{aligned}$$

(B.37)

We estimate:

$$\begin{aligned} {\mathbb {E}}_1 [\Pi ]&= e^{-p \delta \kappa /\epsilon } {\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} \bigg | \int _0^t \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \bigg |^p\right] \nonumber \\&\quad + {\mathbb {E}}_1 \left[ \max _{k=0,\dots ,N-1} \sup _{t \in [k \delta , (k+2)\delta ]} \bigg | \int _{k \delta }^t \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \bigg |^p \right] \end{aligned}$$

(B.38)

$$\begin{aligned}&\le e^{-p \delta \kappa /\epsilon } {\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} \bigg | \int _0^t \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \bigg |^p \right] \nonumber \\&\quad + {\mathbb {E}}_1 \left[ \left( \sum _{k=0}^{N-1} \sup _{t \in [k \delta , (k+2)\delta ]} \left( \int _{k \delta }^t \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \right) ^{pq} \right) ^{1/q} \right] \end{aligned}$$

(B.39)

$$\begin{aligned}&\le e^{-p \delta \kappa /\epsilon } {\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} \bigg | \int _0^t \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \bigg |^p \right] \nonumber \\&\quad + \left( \sum _{k=0}^{N-1} {\mathbb {E}}_1 \left[ \sup _{t \in [k \delta , (k+2)\delta ]} \left( \int _{k \delta }^t \varvec{\sigma }_2(s, \varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \right) ^{pq} \right) ^{1/q} \right] , \end{aligned}$$

(B.40)

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

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

$$\begin{aligned} {\mathbb {E}}_1[\Pi ]&\le C_{p,q} \left[ e^{-p \kappa /\epsilon ^h} {\mathbb {E}}_1 \bigg [ \bigg ( \int _0^T \Vert \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s),\epsilon )\Vert _{F}^2 \mathrm{d}s \bigg )^{\frac{pq}{2}} \bigg ]^{1/q} \right. \nonumber \\&\quad \left. + \left( \sum _{k=0}^{N-1} {\mathbb {E}}_1 \left( \int _{k \delta }^{(k+2)\delta } \Vert \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s),\epsilon ) \Vert _F^2 \mathrm{d}s \right) ^{\frac{pq}{2}} \right) ^{1/q} \right] \end{aligned}$$

(B.41)

$$\begin{aligned}&\le C_{p,q} \Vert \varvec{\sigma }_2(s,\varvec{x}^\epsilon (s),\epsilon )\Vert ^p_{F,\infty } (e^{-p \kappa /\epsilon ^h} T^{p/2} + 2^{p/2} (N \delta ^{\frac{pq}{2}})^{1/q}), \end{aligned}$$

(B.42)

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

$$\begin{aligned} \Vert \varvec{\sigma }_2(s,\varvec{x}^\epsilon (s),\epsilon )\Vert _{F,\infty } := \sup _{t \in [0,T], \varvec{y} \in {\mathbb {R}}^{n_{1}}, \epsilon \in {\mathcal {E}}} \Vert \varvec{\sigma }_2(t, \varvec{y},\epsilon )\Vert _F < \infty . \end{aligned}$$

(B.43)

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

Therefore, we have

$$\begin{aligned} {\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} \left| \varvec{\Phi }_{\epsilon }(t) \int _0^t \varvec{\Phi }_{\epsilon }^{-1}(s) \varvec{\sigma }_2(s,\varvec{x}^\epsilon (s),\epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \right| ^p \right] = O(\epsilon ^\beta ) \end{aligned}$$

(B.44)

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

Combining all the estimates obtained, one has:

$$\begin{aligned} {\mathbb {E}}_1\left[ \sup _{t \in [0,T]} |\varvec{p}^\epsilon (t)|^p \right] \le C_1 \epsilon ^\alpha + C_2 \epsilon ^{p} + C_3 \epsilon ^{p b_2} + C_4 \epsilon ^\beta \end{aligned}$$

(B.45)

where the $C_i$ are positive constants, $\alpha \ge p/2$ is some constant, and $b_2 > 0$ is the constant from Assumption A.5. The statement of the proposition follows.

$\square$

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

Proposition B.5

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

Then for any $p \ge 1$, $T>0$, $0< \beta < p/2$, $i,j = 1, \dots ,n_2$, in the limit as $\epsilon \rightarrow 0$ we have

$$\begin{aligned} {\mathbb {E}}\left[ \sup _{t \in [0,T]} \left| \int _{0}^{t} h^\epsilon (s,\varvec{x}^{\epsilon }(s)) d( [\varvec{p}^{\epsilon }]_{i}(s)\cdot [\varvec{p}^{\epsilon }]_{j}(s)) \right| ^{p}; \epsilon \le \epsilon _1 \right] = O(\epsilon ^\beta ), \end{aligned}$$

(B.46)

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

Proof

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

$$\begin{aligned}&\int _{0}^{t} h^\epsilon (s,\varvec{x}^{\epsilon }(s)) d( [\varvec{p}^{\epsilon }]_{i}(s)\cdot [\varvec{p}^{\epsilon }]_{j}(s)) \nonumber \\&\quad = h^\epsilon (t, \varvec{x}^\epsilon (t)) [\varvec{p}^{\epsilon }]_{i}(t) [\varvec{p}^{\epsilon }]_{j}(t) - h^\epsilon (t, \varvec{x}^\epsilon ) [\varvec{p}^{\epsilon }]_{i} [\varvec{p}^{\epsilon }]_{j} \nonumber \\&\qquad - \int _0^t [\varvec{p}^{\epsilon }]_{i}(s) [\varvec{p}^{\epsilon }]_{j}(s) \left( \varvec{\nabla }_{\varvec{x}^\epsilon } h^\epsilon (s, \varvec{x}^\epsilon (s)) \cdot \frac{\varvec{p}^\epsilon (s)}{\epsilon } + \frac{\partial }{\partial s} h^\epsilon (s, \varvec{x}^\epsilon (s)) \right) \mathrm{d}s. \end{aligned}$$

(B.47)

Using the notation ${\mathbb {E}}_1[\cdot ] = {\mathbb {E}}[\cdot ; \epsilon \le \epsilon _1]$, we estimate, for $p \ge 1$,

$$\begin{aligned}&{\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} \left| \int _{0}^{t} h^\epsilon (s,\varvec{x}^{\epsilon }(s)) d( [\varvec{p}^{\epsilon }]_{i}(s)\cdot [\varvec{p}^{\epsilon }]_{j}(s)) \right| ^{p}\right] \nonumber \\&\quad \le 4^{p-1}\bigg ( {\mathbb {E}}_1 \sup _{t \in [0,T]} \left| h^\epsilon (t, \varvec{x}^\epsilon (t)) [\varvec{p}^{\epsilon }]_{i}(t) [\varvec{p}^{\epsilon }]_{j}(t) \right| ^p \nonumber \\&\qquad + {\mathbb {E}}_1 \sup _{t \in [0,T]} \left| h^\epsilon (t, \varvec{x}^\epsilon ) [\varvec{p}^{\epsilon }]_{i} [\varvec{p}^{\epsilon }]_{j} \right| ^p \nonumber \\&\qquad + {\mathbb {E}}_1 \sup _{t \in [0,T]} \left| \int _0^t [\varvec{p}^{\epsilon }]_{i}(s) [\varvec{p}^{\epsilon }]_{j}(s) \varvec{\nabla }_{\varvec{x}^\epsilon } h^\epsilon (s, \varvec{x}^\epsilon (s)) \cdot \frac{\varvec{p}^\epsilon (s)}{\epsilon } \mathrm{d}s \right| ^p \nonumber \\&\qquad + {\mathbb {E}}_1 \sup _{t \in [0,T]} \left| \int _0^t [\varvec{p}^{\epsilon }]_{i}(s) [\varvec{p}^{\epsilon }]_{j}(s) \frac{\partial }{\partial s} h^\epsilon (s, \varvec{x}^\epsilon (s)) \mathrm{d}s \right| ^p \bigg )\end{aligned}$$

(B.48)

$$\begin{aligned}&\quad \le C(p,T) \bigg [ \Vert h^\epsilon \Vert ^p_{\infty } \left( {\mathbb {E}}_1 \sup _{t \in [0,T]} |\varvec{p}^\epsilon (t)|^{2p} + {\mathbb {E}}_1 |\varvec{p}^\epsilon |^{2p} \right) \nonumber \\&\qquad + \frac{1}{\epsilon ^p} {\mathbb {E}}_1 \sup _{t \in [0,T]} \left| \int _0^t [\varvec{p}^{\epsilon }]_{i}(s) [\varvec{p}^{\epsilon }]_{j}(s) [\varvec{\nabla }_{\varvec{x}^\epsilon } h^\epsilon ]_k (s, \varvec{x}^\epsilon (s)) [\varvec{p}^\epsilon ]_k(s) \mathrm{d}s \right| ^p \nonumber \\&\qquad + \left\| \frac{\partial }{\partial s} h^\epsilon \right\| _{\infty }^p \cdot {\mathbb {E}}_1 \sup _{t \in [0,T]} |\varvec{p}^\epsilon (t)|^{2p} \bigg ], \end{aligned}$$

(B.49)

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

Now, estimating as before, we obtain:

$$\begin{aligned}&{\mathbb {E}}_1 \sup _{t \in [0,T]} \left| \int _0^t [\varvec{p}^{\epsilon }]_{i}(s) [\varvec{p}^{\epsilon }]_{j}(s) [\varvec{\nabla }_{\varvec{x}^\epsilon } h^\epsilon ]_k (s, \varvec{x}^\epsilon (s)) [\varvec{p}^\epsilon ]_k(s) \mathrm{d}s \right| ^p \nonumber \\&\quad \le D(p,T) \Vert \varvec{\nabla }_{\varvec{x}^\epsilon } h^\epsilon \Vert _{\infty } \cdot {\mathbb {E}}_1 \sup _{t \in [0,T]} |\varvec{p}^\epsilon (t)|^{3p}, \end{aligned}$$

(B.50)

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

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

$$\begin{aligned}&{\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} \left| \int _{0}^{t} h^\epsilon (s,\varvec{x}^{\epsilon }(s)) d( [\varvec{p}^{\epsilon }]_{i}(s)\cdot [\varvec{p}^{\epsilon }]_{j}(s)) \right| ^{p}\right] = O(\epsilon ^\beta ), \end{aligned}$$

(B.51)

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

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

$$\begin{aligned} {\mathbb {E}}\left[ \sup _{t \in [0,T]} |\varvec{x}_t^\epsilon - \varvec{X}_t|^p; \epsilon \le \epsilon _1 \right] \rightarrow 0, \end{aligned}$$

(B.52)

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

Proof of Theorem A.6

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

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

$$\begin{aligned} \varvec{v}^{\epsilon }(t) \mathrm{d}t&= \epsilon \varvec{a}_{2}^{-1}(t,\varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}\varvec{v}^{\epsilon }(t) - \varvec{a}_{2}^{-1}(t,\varvec{x}^{\epsilon }(t),\epsilon ) \varvec{b}_{2}(t,\varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}t \nonumber \\&\quad - \varvec{a}_{2}^{-1}(t,\varvec{x}^{\epsilon }(t),\epsilon ) \varvec{\sigma }_2(t,\varvec{x}^{\epsilon }(t), \epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(t). \end{aligned}$$

(B.53)

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

$$\begin{aligned} \mathrm{d}\varvec{x}^{\epsilon }(t)&= \epsilon \varvec{a}_{1}(t,\varvec{x}^{\epsilon }(t), \epsilon ) \varvec{a}_{2}^{-1}(t,\varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}\varvec{v}^{\epsilon }(t) \nonumber \\&\quad - \varvec{a}_{1}(t, \varvec{x}^{\epsilon }(t), \epsilon ) \varvec{a}_{2}^{-1}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \varvec{b}_{2}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}t \nonumber \\&\quad - \varvec{a}_{1}(t, \varvec{x}^{\epsilon }(t), \epsilon ) \varvec{a}_{2}^{-1}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \varvec{\sigma }_2(t, \varvec{x}^{\epsilon }(t), \epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(t) \nonumber \\&\quad + \varvec{b}_{1}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}t + \varvec{\sigma }_{1}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}\varvec{W}^{(k_1)}(t). \end{aligned}$$

(B.54)

In integral form, we have:

$$\begin{aligned} \varvec{x}^{\epsilon }(t)&= \varvec{x}^\epsilon + \epsilon \int _0^t \varvec{a}_{1}(s, \varvec{x}^{\epsilon }(s), \epsilon ) \varvec{a}_{2}^{-1}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}\varvec{v}^{\epsilon }(s) \nonumber \\&\quad + \int _0^t \{ \varvec{b}_{1}(s, \varvec{x}^{\epsilon }(s),\epsilon ) -\varvec{a}_{1}(s, \varvec{x}^{\epsilon }(s), \epsilon ) \varvec{a}_{2}^{-1}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \varvec{b}_{2}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \} \mathrm{d}s \nonumber \\&\quad - \int _0^t \varvec{a}_{1}(s, \varvec{x}^{\epsilon }(s), \epsilon ) \varvec{a}_{2}^{-1}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \varvec{\sigma }_2(s, \varvec{x}^{\epsilon }(s), \epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s) \nonumber \\&\quad + \int _0^t \varvec{\sigma }_{1}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}\varvec{W}^{(k_1)}(s). \end{aligned}$$

(B.55)

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

$$\begin{aligned}{}[\varvec{x}^{\epsilon }]_i(t)&= [\varvec{x}^\epsilon ]_i + \epsilon \int _0^t [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s), \epsilon ) \cdot \mathrm{d}[\varvec{v}^{\epsilon }]_j(s) \nonumber \\&\quad + \int _0^t \{ [\varvec{b}_{1}]_i(s, \varvec{x}^{\epsilon }(s),\epsilon ) -[\varvec{a}_{1}\varvec{a}_{2}^{-1}\varvec{b}_{2}]_{i}(s, \varvec{x}^{\epsilon }(s), \epsilon ) \} \mathrm{d}s \nonumber \\&\quad - \int _0^t [\varvec{a}_{1}\varvec{a}_{2}^{-1}\varvec{\sigma }_2]_{i,j}(s, \varvec{x}^{\epsilon }(s), \epsilon ) \cdot \mathrm{d}[\varvec{W}^{(k_2)}]_j(s) \nonumber \\&\quad + \int _0^t [\varvec{\sigma }_{1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \cdot \mathrm{d}[\varvec{W}^{(k_1)}]_j(s). \end{aligned}$$

(B.56)

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

$$\begin{aligned}&\int _0^t [S^{\epsilon }]_i(s, \varvec{x}^\epsilon (s),\varvec{v}^\epsilon (s),\epsilon )\mathrm{d}s :=\epsilon \int _0^t [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s), \epsilon ) \cdot \mathrm{d}[\varvec{v}^{\epsilon }]_j(s) \end{aligned}$$

(B.57)

$$\begin{aligned}&\quad = \epsilon [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \cdot [\varvec{v}^{\epsilon }]_j(t) - \epsilon [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(0, \varvec{x},\epsilon ) \cdot [\varvec{v}^\epsilon ]_j \nonumber \\&\qquad - \int _0^t \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \cdot \mathrm{d}[\varvec{x}^\epsilon ]_l(s) \cdot \epsilon [\varvec{v}^{\epsilon }]_j(s) \nonumber \\&\qquad - \int _0^t \frac{\partial }{\partial s}\left( [\varvec{a}_1 \varvec{a}_2^{-1}]_{i,j}(s, \varvec{x}^\epsilon (s),\epsilon ) \right) \cdot \epsilon [\varvec{v}^\epsilon ]_j(s) \mathrm{d}s. \end{aligned}$$

(B.58)

Substituting the following expression for $\mathrm{d}[\varvec{x}^\epsilon ]_l(s)$:

$$\begin{aligned} \mathrm{d}[\varvec{x}^\epsilon ]_l(s)&= [\varvec{a}_1]_{l,k}(s,\varvec{x}^\epsilon (s),\epsilon )[\varvec{v}^\epsilon ]_k(s) \mathrm{d}s + [\varvec{b}_1]_l(s,\varvec{x}^\epsilon (s),\epsilon )\mathrm{d}s \nonumber \\&\quad + [\varvec{\sigma }_1]_{l,k}(s,\varvec{x}^\epsilon (s),\epsilon )\mathrm{d}[\varvec{W}^{(k_1)}]_k(s) \end{aligned}$$

(B.59)

into (B.58), we obtain:

$$\begin{aligned}&\int _0^t [S^{\epsilon }]_i(s, \varvec{x}^\epsilon (s),\varvec{v}^\epsilon (s),\epsilon )\mathrm{d}s \nonumber \\&\quad = \epsilon [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \cdot [\varvec{v}^{\epsilon }]_j(t) - \epsilon [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(0, \varvec{x},\epsilon ) \cdot [\varvec{v}^\epsilon ]_j \nonumber \\&\qquad - \int _0^t \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \cdot [\varvec{b}_1]_l(s, \varvec{x}^\epsilon (s),\epsilon ) \cdot \epsilon [\varvec{v}^{\epsilon }]_j(s) \mathrm{d}s \nonumber \\&\qquad - \int _0^t \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \bigg )\nonumber \\&\qquad [\varvec{\sigma }_1]_{l,k}(s, \varvec{x}^\epsilon (s),\epsilon ) \epsilon [\varvec{v}^{\epsilon }]_j(s) \mathrm{d}[\varvec{W}^{(k_1)}]_{k}(s) \nonumber \\&\qquad - \int _0^t \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \nonumber \\&\qquad \times [\varvec{a}_1]_{l,k}(s, \varvec{x}^\epsilon (s),\epsilon ) \epsilon [\varvec{v}^\epsilon ]_k(s) [\varvec{v}^{\epsilon }]_j(s) \mathrm{d}s \nonumber \\&\qquad - \int _0^t \frac{\partial }{\partial s}\left( [\varvec{a}_1 \varvec{a}_2^{-1}]_{i,j}(s, \varvec{x}^\epsilon (s),\epsilon ) \right) \cdot \epsilon [\varvec{v}^\epsilon ]_j(s) \mathrm{d}s. \end{aligned}$$

(B.60)

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

$$\begin{aligned}&d[\epsilon \varvec{v}^{\epsilon }(t) (\epsilon \varvec{v}^{\epsilon }(t))^{*}] \nonumber \\&\quad = \epsilon \mathrm{d}\varvec{v}^{\epsilon }(t) \cdot \epsilon (\varvec{v}^{\epsilon }(t))^* + \epsilon \varvec{v}^{\epsilon }(t) \cdot \epsilon d(\varvec{v}^{\epsilon }(t))^{*} + d[\epsilon \varvec{v}^{\epsilon }(t)] \cdot d[ (\epsilon \varvec{v}^{\epsilon }(t))^{*}] \end{aligned}$$

(B.61)

$$\begin{aligned}&\quad = \left[ \varvec{a}_{2}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \varvec{v}^{\epsilon }(t) \mathrm{d}t + \varvec{b}_{2}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}t + \varvec{\sigma }_2(t, \varvec{x}^{\epsilon }(t), \epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(t) \right] \epsilon \varvec{v}^\epsilon (t)^{*} \nonumber \\&\qquad + \epsilon \varvec{v}^\epsilon (t)\left[ \varvec{a}_{2}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \varvec{v}^{\epsilon }(t) \mathrm{d}t + \varvec{b}_{2}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}t + \varvec{\sigma }_2(t, \varvec{x}^{\epsilon }(t), \epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(t) \right] ^{*}\nonumber \\&\qquad + \varvec{\sigma }_2(t, \varvec{x}^{\epsilon }(t), \epsilon )\varvec{\sigma }_2^*(t, \varvec{x}^{\epsilon }(t), \epsilon ) \mathrm{d}t. \end{aligned}$$

(B.62)

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

$$\begin{aligned} -\varvec{a}_{2}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \varvec{J}^\epsilon (t)\mathrm{d}t - \varvec{J}^\epsilon (t) \varvec{a}_{2}^*(t, \varvec{x}^{\epsilon }(t),\epsilon )\mathrm{d}t = \varvec{F}^{\epsilon }_1(t) \mathrm{d}t + \varvec{F}^{\epsilon }_2(t) \mathrm{d}t + \varvec{F}^{\epsilon }_3(t) \mathrm{d}t, \end{aligned}$$

(B.63)

where

$$\begin{aligned} \varvec{F}^{\epsilon }_1(t) \mathrm{d}t&= -d[\epsilon \varvec{v}^{\epsilon }(t) (\epsilon \varvec{v}^{\epsilon }(t))^{*}], \end{aligned}$$

(B.64)

$$\begin{aligned} \varvec{F}^{\epsilon }_2(t) \mathrm{d}t&= (\varvec{b}_{2}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}t + \varvec{\sigma }_2(t, \varvec{x}^{\epsilon }(t), \epsilon )\mathrm{d}\varvec{W}^{(k_2)}(t) )\epsilon (\varvec{v}^{\epsilon }(t))^{*} \nonumber \\&\qquad + \epsilon \varvec{v}^{\epsilon }(t)(\varvec{b}_{2}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \mathrm{d}t + \varvec{\sigma }_2(t,\varvec{x}^{\epsilon }(t), \epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(t))^{*}, \end{aligned}$$

(B.65)

$$\begin{aligned} \varvec{F}^{\epsilon }_3(t)&= \varvec{\sigma }_2(t,\varvec{x}^{\epsilon }(t), \epsilon ) \varvec{\sigma }_2(t,\varvec{x}^{\epsilon }(t), \epsilon )^{*}. \end{aligned}$$

(B.66)

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

$$\begin{aligned} \varvec{J}^\epsilon (t) = \varvec{J}_1^\epsilon (t) + \varvec{J}_2^\epsilon (t) + \varvec{J}_3^\epsilon (t), \end{aligned}$$

(B.67)

where

$$\begin{aligned} \varvec{J}_n^\epsilon (t)&= \int _0^\infty e^{\varvec{a}_2(t, \varvec{x}^\epsilon (t),\epsilon )y} \varvec{F}^{\epsilon }_n(t) e^{\varvec{a}^*_2(t, \varvec{x}^\epsilon (t),\epsilon )y} \mathrm{d}y \end{aligned}$$

(B.68)

for $n=1,2,3$.

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

$$\begin{aligned}&\epsilon [\varvec{v}^\epsilon ]_j(s) [\varvec{v}^\epsilon ]_k(s)\mathrm{d}s \nonumber \\&\quad =-\int _0^\infty \bigg [e^{\varvec{a}_2(s, \varvec{x}^\epsilon (s),\epsilon )y}\bigg ]_{j,p_1} \cdot \bigg [ d[\epsilon \varvec{v}^{\epsilon }(s) (\epsilon \varvec{v}^{\epsilon }(s))^{*}] \bigg ]_{p_1,p_2} \nonumber \\&\qquad \cdot \bigg [ e^{\varvec{a}^*_2(s, \varvec{x}^\epsilon (s),\epsilon )y}\bigg ]_{p_2,k} \mathrm{d}y \nonumber \\&\qquad + \int _0^\infty \bigg [ e^{\varvec{a}_2(s, \varvec{x}^\epsilon (s),\epsilon )y} \bigg ]_{j,p_1} \cdot \bigg [ (\varvec{b}_{2}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}s \nonumber \\&\qquad + \varvec{\sigma }_2(s, \varvec{x}^{\epsilon }(s), \epsilon )\mathrm{d}\varvec{W}^{(k_2)}(s) )\epsilon (\varvec{v}^{\epsilon }(s))^{*} \bigg ]_{p_1,p_2} \cdot \bigg [ e^{\varvec{a}^*_2(s, \varvec{x}^\epsilon (s),\epsilon )y} \bigg ]_{p_2,k} \mathrm{d}y \nonumber \\&\qquad + \int _0^\infty \bigg [ e^{\varvec{a}_2(s, \varvec{x}^\epsilon (s),\epsilon )y} \bigg ]_{j,p_1} \cdot \bigg [ \epsilon \varvec{v}^{\epsilon }(s)(\varvec{b}_{2}(s,\varvec{x}^{\epsilon }(s),\epsilon ) \mathrm{d}s \nonumber \\&\qquad + \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s), \epsilon ) \mathrm{d}\varvec{W}^{(k_2)}(s))^{*} \bigg ]_{p_1,p_2} \cdot \bigg [ e^{\varvec{a}^*_2(s,\varvec{x}^\epsilon (s),\epsilon )y} \bigg ]_{p_2,k} \mathrm{d}y \nonumber \\&\qquad + \int _0^\infty \bigg [ e^{\varvec{a}_2(s, \varvec{x}^\epsilon (s),\epsilon )y} \bigg ]_{j,p_1} \cdot \bigg [ \varvec{\sigma }_2(s, \varvec{x}^{\epsilon }(s), \epsilon ) \varvec{\sigma }_2(s, \varvec{x}^{\epsilon }(s), \epsilon )^{*} \mathrm{d}s \bigg ]_{p_1,p_2} \nonumber \\&\qquad \cdot \bigg [ e^{\varvec{a}^*_2(s, \varvec{x}^\epsilon (s),\epsilon )y} \bigg ]_{p_2,k} \mathrm{d}y. \end{aligned}$$

(B.69)

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

$$\begin{aligned} \varvec{X}(t)&= \varvec{x} + \int _0^t [\varvec{B}_1(s, \varvec{X}(s))-\varvec{A}_1(s, \varvec{X}(s))\varvec{A}_2^{-1}(s, \varvec{X}(s))\varvec{B}_2(s, \varvec{X}(s))] \mathrm{d}s\nonumber \\&\quad + \int _0^t \varvec{S}(s, \varvec{X}(s)) \mathrm{d}s + \int _0^t \varvec{\Sigma }_1(s, \varvec{X}(s)) \mathrm{d}\varvec{W}^{(k_1)}(s) \nonumber \\&\quad - \int _0^t \varvec{A}_1(s, \varvec{X}(s)) \varvec{A}_2^{-1}(s, \varvec{X}(s))\varvec{\Sigma }_2(s, \varvec{X}(s)) \mathrm{d}\varvec{W}^{(k_2)}(s). \end{aligned}$$

(B.70)

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

For any $p > 2$, $T>0$, $i=1,\dots ,n_1$ (recall that $[\varvec{b}]_i$ denotes the ith component of vector $\varvec{b}$), we estimate:

$$\begin{aligned}&{\mathbb {E}}_1\left[ \sup _{t \in [0,T]} |[\varvec{x}^{\epsilon }(t) - \varvec{X}(t)]_i|^p \right] \nonumber \\&\quad \le 6^{p-1}\bigg \{ {\mathbb {E}}_1\left[ |\varvec{x}^\epsilon - \varvec{x}|^p\right] \nonumber \\&\qquad + {\mathbb {E}}_1\left[ \sup _{t \in [0,T]} \bigg | \int _0^t \bigg [\varvec{S}_\epsilon (s, \varvec{x}^\epsilon (s),\varvec{v}^\epsilon (s),\epsilon )- \varvec{S}(s, \varvec{X}(s))\bigg ]_i \mathrm{d}s \bigg |^p \right] \nonumber \\&\qquad + {\mathbb {E}}_1 \bigg [ \sup _{t \in [0,T]} \bigg ( \int _0^t \bigg | \bigg [ \varvec{a}_{1}(s, \varvec{x}^{\epsilon }(s), \epsilon ) \varvec{a}_{2}^{-1}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \varvec{b}_{2}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \nonumber \\&\qquad - \varvec{A}_1(s, \varvec{X}(s))\varvec{A}_2^{-1}(s,\varvec{X}(s))\varvec{B}_2(s,\varvec{X}(s)) \bigg ]_i \bigg | \mathrm{d}s \bigg )^p \bigg ] \nonumber \\&\qquad + {\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} \bigg ( \int _0^t \bigg | \bigg [ \varvec{b}_{1}(s,\varvec{x}^{\epsilon }(s),\epsilon ) - \varvec{B}_1(s,\varvec{X}(s)) \bigg ]_i \bigg | \mathrm{d}s \bigg )^p \right] \nonumber \\&\qquad + {\mathbb {E}}_1 \bigg [ \sup _{t \in [0,T]} \bigg | \int _0^t \bigg [ \varvec{a}_{1}(s,\varvec{x}^{\epsilon }(s), \epsilon ) \varvec{a}_{2}^{-1}(s,\varvec{x}^{\epsilon }(s),\epsilon ) \varvec{\sigma }_2(s,\varvec{x}^{\epsilon }(s), \epsilon ) \nonumber \\&\qquad - \varvec{A}_1(s,\varvec{X}(s)) \varvec{A}_2^{-1}(s,\varvec{X}(s))\varvec{\Sigma }_2(s,\varvec{X}(s)) \bigg ]_{i,j} \mathrm{d}[\varvec{W}^{(k_2)}]_j(s) \bigg |^p \bigg ] \nonumber \\&\qquad + {\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} \bigg | \int _0^t \bigg [\varvec{\sigma }_{1}(s,\varvec{x}^{\epsilon }(s),\epsilon ) - \varvec{\Sigma }_1(s,\varvec{X}(s)) \bigg ]_{i,j} \mathrm{d}[\varvec{W}^{(k_1)}]_j(s) \bigg |^p \right] \ \bigg \} \end{aligned}$$

(B.71)

$$\begin{aligned}&\quad =: 6^{p-1}\left( \sum _{k=0}^5 R_k \right) . \end{aligned}$$

(B.72)

By Assumption A.4, $R_0 = {\mathbb {E}}_1 \left[ |\varvec{x}^\epsilon - \varvec{x}|^p\right] \le {\mathbb {E}}\left[ |\varvec{x}^\epsilon - \varvec{x}|^p \right] = O(\epsilon ^{ p r_0})$ as $\epsilon \rightarrow 0$, where $r_0 > 1/2$ is a constant. We now estimate each of the $R_k$, $k=1,\dots ,5$.

We have:

$$\begin{aligned} R_3&\le {\mathbb {E}}_1 \sup _{t \in [0,T]} \bigg ( \int _0^t | \varvec{b}_{1}(s,\varvec{x}^{\epsilon }(s),\epsilon ) - \varvec{B}_1(s,\varvec{X}(s)) | \mathrm{d}s \bigg )^p \end{aligned}$$

(B.73)

$$\begin{aligned}&= {\mathbb {E}}_1 \sup _{t \in [0,T]} \bigg ( \int _0^t |\varvec{b}_{1}(s,\varvec{x}^\epsilon (s),\epsilon ) - \varvec{b}_1(s,\varvec{X}(s),\epsilon ) + \varvec{b}_1(s,\varvec{X}(s),\epsilon ) \nonumber \\&\quad - \varvec{B}_1(s,\varvec{X}(s)) | \mathrm{d}s \bigg )^p \nonumber \\&\le 2^{p-1} \bigg [ {\mathbb {E}}_1 \sup _{t \in [0,T]} \left( \int _0^t |\varvec{b}_{1}(s,\varvec{x}^{\epsilon }(s),\epsilon ) - \varvec{b}_1(s,\varvec{X}(s),\epsilon )| \mathrm{d}s \right) ^p \nonumber \\&\quad + {\mathbb {E}}_1 \sup _{t \in [0,T]} \left( \int _0^t |\varvec{b}_1(s,\varvec{X}(s),\epsilon ) - \varvec{B}_1(s,\varvec{X}(s))| \mathrm{d}s \right) ^p \bigg ] \end{aligned}$$

(B.74)

$$\begin{aligned}&\le 2^{p-1}\left[ L^p(\epsilon ) {\mathbb {E}}_1 \sup _{t \in [0,T]} \int _0^t |\varvec{x}^\epsilon (s)-\varvec{X}(s) |^p \mathrm{d}s + T^p \beta _1(\epsilon )^p \mathbb {1}_{\{\varvec{b}_1 \ne \varvec{B}_1\}} \right] \end{aligned}$$

(B.75)

$$\begin{aligned}&\le L_3(\epsilon ,p,T) \int _0^T {\mathbb {E}}_1 \sup _{u \in [0,s]} |\varvec{x}^\epsilon (u)-\varvec{X}(u) |^p \mathrm{d}s + C_3(p,T) \beta _1(\epsilon )^p \mathbb {1}_{\{\varvec{b}_1 \ne \varvec{B}_1\}}, \end{aligned}$$

(B.76)

on the set $S_1 := \{\epsilon : \epsilon \le \epsilon _1\}$, where $\mathbb {1}_{A}$ denotes the indicator function of a set A, $L_3(\epsilon ,p,T) = O(1)$ as $\epsilon \rightarrow 0$ and $C_3(p,T)$ is a constant dependent on p and T. In the last two lines of the above estimate, we have used Assumption A.3, Assumption A.5, and the inequality:

$$\begin{aligned} {\mathbb {E}}_1 \sup _{t\in [0,T]} \left( \int _0^t |\varvec{u}(s)| \mathrm{d}s \right) ^p \le T^{p-1} {\mathbb {E}}_1 \int _0^T |\varvec{u}(s)|^p \mathrm{d}s, \end{aligned}$$

(B.77)

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

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

$$\begin{aligned} R_2&\le L_2(\epsilon ,p,T) \int _0^T {\mathbb {E}}_1 \sup _{u \in [0,s]} |\varvec{x}^\epsilon (u)-\varvec{X}(u) |^p \mathrm{d}s \nonumber \\&\qquad + C_2(p,T)\left[ \alpha _1(\epsilon )^p \mathbb {1}_{\{\varvec{a}_1 \ne \varvec{A}_1\}} + \alpha _2(\epsilon )^p \mathbb {1}_{\{\varvec{a}_2 \ne \varvec{A}_2\}} + \beta _2(\epsilon )^p \mathbb {1}_{\{\varvec{b}_2 \ne \varvec{B}_2\}} \right] , \end{aligned}$$

(B.78)

on $S_1$, where $\alpha _1(\epsilon )$, $\alpha _2(\epsilon )$, $\beta _2(\epsilon )$ are from Assumption A.3, $L_2(\epsilon , p,T) = O(1)$ as $\epsilon \rightarrow 0$ and $C_2(p,T)$ is a constant.

To estimate $R_5$, we use the Burkholder-Davis-Gundy inequality:

$$\begin{aligned} R_5&\le C'_p {\mathbb {E}}_1 \bigg ( \int _0^T \Vert \varvec{\sigma }_{1}(s,\varvec{x}^{\epsilon }(s),\epsilon ) - \varvec{\Sigma }_1(s,\varvec{X}(s))\Vert _{F}^2 \mathrm{d}s \bigg )^{p/2}, \end{aligned}$$

(B.79)

where $C'_p$ is a positive constant and $\Vert \cdot \Vert _F$ denotes the Frobenius norm. Using Hölder’s inequality, Assumption A.3, Assumption A.5, and the above techniques, we obtain:

$$\begin{aligned} R_5&\le C''_p {\mathbb {E}}_1 \bigg ( \int _0^T \Vert \varvec{\sigma }_1(s,\varvec{x}^\epsilon (s),\epsilon )-\varvec{\sigma }_1(s,\varvec{X}(s),\epsilon )\Vert _{F}^2 \mathrm{d}s \bigg )^{p/2} \nonumber \\&\quad + C''_p {\mathbb {E}}_1 \bigg ( \int _0^T \Vert \varvec{\sigma }_1(s,\varvec{X}(s),\epsilon ) - \varvec{\Sigma }_1(s,\varvec{X}(s)) \Vert _{F}^2 \mathrm{d}s \bigg )^{p/2} \end{aligned}$$

(B.80)

$$\begin{aligned}&\le C''_p T^{\frac{p}{2}-1} \int _0^T {\mathbb {E}}_1 \Vert \varvec{\sigma }_1(s,\varvec{x}^\epsilon (s),\epsilon )-\varvec{\sigma }_1(s,\varvec{X}(s),\epsilon )\Vert _{F}^p \mathrm{d}s \nonumber \\&\quad + C'''_p |\gamma _1(\epsilon )|^p T^{\frac{p}{2}} \mathbb {1}_{\{\varvec{\sigma }_1 \ne \varvec{\Sigma }_1 \}}\end{aligned}$$

(B.81)

$$\begin{aligned}&\le L_5(\epsilon ,p,T) \int _0^T {\mathbb {E}}_1 \sup _{u \in [0,s]} |\varvec{x}^\epsilon (u)-\varvec{X}(u) |^p \mathrm{d}s + C_5(p,T)\gamma _1(\epsilon )^p \mathbb {1}_{\{\varvec{\sigma }_1 \ne \varvec{\Sigma }_1\}}, \end{aligned}$$

(B.82)

on the set $S_1$, where $C_p''$ and $C_p'''$ are constants, $\gamma _1(\epsilon )$ is from Assumption A.3, $L_5(\epsilon ,p,T) = O(1)$ as $\epsilon \rightarrow 0$ and $C_5(p,T)$ is a constant.

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

$$\begin{aligned} R_4&\le L_4(\epsilon ,p,T) \int _0^T {\mathbb {E}}_1 \sup _{u \in [0,s]} |\varvec{x}^\epsilon (u)-\varvec{X}(u) |^p \mathrm{d}s \nonumber \\&\quad + C_4(p,T)\left[ \alpha _1(\epsilon )^p \mathbb {1}_{\{\varvec{a}_1 \ne \varvec{A}_1\}} + \alpha _2(\epsilon )^p \mathbb {1}_{\{\varvec{a}_2 \ne \varvec{A}_2\}} + \gamma _2(\epsilon )^p \mathbb {1}_{\{\varvec{\sigma }_2 \ne \varvec{\Sigma }_2\}} \right] , \end{aligned}$$

(B.83)

on $S_1$, where $\gamma _2(\epsilon )$ is from Assumption A.3, $L_4(\epsilon ,p,T) = O(1)$ as $\epsilon \rightarrow 0$ and $C_4(p,T)$ is a constant.

To obtain a bound for $R_1$, first we estimate:

$$\begin{aligned}&\bigg | \int _0^t \bigg [\varvec{S}^\epsilon (s,\varvec{x}^\epsilon (s),\varvec{v}^\epsilon (s),\epsilon )- \varvec{S}(s,\varvec{X}(s))\bigg ]_i \mathrm{d}s \bigg | \nonumber \\&\quad \le \bigg |\epsilon [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(t, \varvec{x}^{\epsilon }(t),\epsilon ) \cdot [\varvec{v}^{\epsilon }]_j(t) - \epsilon [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(0, \varvec{x},\epsilon ) \cdot [\varvec{v}]_j \bigg | \nonumber \\&\qquad + \left| \int _0^t \frac{\partial }{\partial s} \left( [\varvec{a}_1 \varvec{a}_2^{-1}]_{i,j}(s, \varvec{x}^\epsilon (s), \epsilon ) \right) \cdot \epsilon [\varvec{v}^\epsilon ]_j(s) \mathrm{d}s \right| \nonumber \\&\qquad + \int _0^t \bigg | \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \cdot [\varvec{b}_1]_l(s, \varvec{x}^\epsilon (s),\epsilon ) \cdot \epsilon [\varvec{v}^{\epsilon }]_j(s) \bigg | \mathrm{d}s \nonumber \\&\qquad + \bigg | \int _0^t \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \cdot [\varvec{\sigma }_1]_{l,k}(s, \varvec{x}^\epsilon (s),\epsilon ) \nonumber \\&\qquad \cdot \epsilon [\varvec{v}^{\epsilon }]_j(s) \mathrm{d}[\varvec{W}^{(k_1)}]_{k}(s) \bigg | \nonumber \\&\qquad + \bigg | \int _0^t \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \cdot [\varvec{a}_1]_{l,k}(s, \varvec{x}^\epsilon (s),\epsilon ) \cdot [\varvec{J}_1^\epsilon ]_{j,k}(s) \mathrm{d}s \bigg | \nonumber \\&\qquad + \bigg | \int _0^t \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \cdot [\varvec{a}_1]_{l,k}(s, \varvec{x}^\epsilon (s),\epsilon ) \cdot [\varvec{J}_2^\epsilon ]_{j,k}(s) \mathrm{d}s \bigg | \nonumber \\&\qquad + \bigg | \int _0^t \frac{\partial }{\partial [\varvec{X}]_l(s)}\left( [\varvec{A}_1\varvec{A}_2^{-1}]_{i,j}(s, \varvec{X}(s)) \right) \cdot [\varvec{A}_1]_{l,k}(s, \varvec{X}(s))\cdot [\varvec{J}]_{j,k}(s) \nonumber \\&\qquad - \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \cdot [\varvec{a}_1]_{l,k}(s, \varvec{x}^\epsilon (s),\epsilon ) \cdot [\varvec{J}_3^\epsilon ]_{j,k}(s) \mathrm{d}s \bigg |\end{aligned}$$

(B.84)

$$\begin{aligned}&\quad =: \sum _{k=0}^6 \Pi _k, \end{aligned}$$

(B.85)

and so $R_1 \le 6^{p-1} \sum _{k=0}^6 \left( {\mathbb {E}}_1 \sup _{t\in [0,T]} |\Pi _k|^p \right) =: 6^{p-1} \sum _{k=0}^6 M_k.$

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

$$\begin{aligned} M_k \le C_k(p,T) \cdot {\mathbb {E}}_1 \sup _{t\in [0,T]} |\epsilon \varvec{v}^\epsilon (t)|^p, \end{aligned}$$

(B.86)

where the $C_k$ are positive constants.

Applying Proposition B.5, we obtain:

$$\begin{aligned} M_4 := {\mathbb {E}}_1 \sup _{t \in [0,T]} |\Pi _4|^p \le C_4(p,T) \epsilon ^\beta , \end{aligned}$$

(B.87)

on $S_1$, for all $0< \beta < p/2$, as $\epsilon \rightarrow 0$, where $C_4(p,T)$ is a positive constant.

We now estimate $M_6$:

$$\begin{aligned}&M_6 \le {\mathbb {E}}_1 \sup _{t \in [0,T]}\bigg ( \int _0^t \bigg | \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s,\varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \cdot [\varvec{a}_1]_{l,k}(s,\varvec{x}^\epsilon (s),\epsilon ) \nonumber \\&\qquad \cdot [\varvec{J}_3^\epsilon ]_{j,k}(s) - \frac{\partial }{\partial [\varvec{X}]_l(s)}\left( [\varvec{A}_1\varvec{A}_2^{-1}]_{i,j}(s,\varvec{X}(s)) \right) \cdot [\varvec{A}_1]_{l,k}(s,\varvec{X}(s)) \nonumber \\&\qquad \cdot [\varvec{J}]_{j,k}(s) \bigg | \mathrm{d}s \bigg )^p \end{aligned}$$

(B.88)

$$\begin{aligned}&\quad \le C(p) {\mathbb {E}}_1 \sup _{t \in [0,T]}\bigg ( \int _0^t \bigg | \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s,\varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \cdot [\varvec{a}_1]_{l,k}(s, \varvec{x}^\epsilon (s),\epsilon ) \nonumber \\&\qquad - \frac{\partial }{\partial [\varvec{X}]_l(s)}\bigg ([\varvec{A}_1\varvec{A}_2^{-1}]_{i,j}(s,\varvec{X}(s)) \bigg ) \cdot [\varvec{A}_1]_{l,k}(s, \varvec{X}(s)) \bigg |^p \cdot |[\varvec{J}_3^\epsilon ]_{j,k}(s)|^p \mathrm{d}s \bigg ) \nonumber \\&\qquad + C(p) {\mathbb {E}}_1 \sup _{t\in [0,T]} \bigg (\int _0^t \bigg | \frac{\partial }{\partial [\varvec{X}]_l(s)}\bigg ([\varvec{A}_1\varvec{A}_2^{-1}]_{i,j}(s, \varvec{X}(s)) \bigg ) \cdot [\varvec{A}_1]_{l,k}(s, \varvec{X}(s)) \bigg |^p \nonumber \\&\qquad \cdot |[\varvec{J}_3^\epsilon -\varvec{J}]_{j,k}(s)|^p \mathrm{d}s\bigg ) \end{aligned}$$

(B.89)

$$\begin{aligned}&\quad \le C(p) {\mathbb {E}}_1 \sup _{t \in [0,T]}\bigg ( \int _0^t \bigg | \frac{\partial }{\partial [\varvec{x}^{\epsilon }]_l(s)}\bigg ( [\varvec{a}_{1}\varvec{a}_{2}^{-1}]_{i,j}(s, \varvec{x}^{\epsilon }(s),\epsilon ) \bigg ) \cdot [\varvec{a}_1]_{l,k}(s, \varvec{x}^\epsilon (s),\epsilon ) \nonumber \\&\qquad - \frac{\partial }{\partial [\varvec{X}]_l(s)}\bigg ([\varvec{A}_1\varvec{A}_2^{-1}]_{i,j}(s, \varvec{X}(s)) \bigg ) \cdot [\varvec{A}_1]_{l,k}(s, \varvec{X}(s)) \bigg |^p \cdot |[\varvec{J}_3^\epsilon ]_{j,k}(s)|^p \mathrm{d}s \bigg ) \nonumber \\&\qquad + C(p) {\mathbb {E}}_1 \sup _{t \in [0,T]} \int _0^t \Vert \varvec{J}_3^\epsilon (s) - \varvec{J}(s) \Vert _F^p \mathrm{d}s, \end{aligned}$$

(B.90)

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

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

$$\begin{aligned} \varvec{a}_2(s, \varvec{x}^\epsilon (s),\epsilon )\varvec{J}_3^\epsilon (s) +\varvec{J}_3^\epsilon (s) \varvec{a}_2^*(s, \varvec{x}^\epsilon (s),\epsilon ) = -(\varvec{\sigma }_2 \varvec{\sigma }^*_2)(s, \varvec{x}^\epsilon (s),\epsilon ) \end{aligned}$$

(B.91)

and

$$\begin{aligned} \varvec{A}_2(s, \varvec{X}(s)) \varvec{J}(s) +\varvec{J}(s) \varvec{A}_2^{*}(s, \varvec{X}(s)) = -(\varvec{\Sigma }_2 \varvec{\Sigma }_2^{*})(s, \varvec{X}(s)). \end{aligned}$$

(B.92)

respectively.

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

$$\begin{aligned}&\varvec{A}_2(s, \varvec{X}(s)) \varvec{H}^\epsilon (s) + \varvec{H}^\epsilon (s) \varvec{A}_2^*(s, \varvec{X}(s)) \nonumber \\&\quad = (\varvec{\Sigma }_2 \varvec{\Sigma }_2^*)(s, \varvec{X}(s)) - (\varvec{\sigma }_2 \varvec{\sigma }_2^*)(s, \varvec{x}^\epsilon (s),\epsilon ) - \varvec{G}^\epsilon (s)\varvec{J}_3^\epsilon (s) - \varvec{J}_3^\epsilon (s) (\varvec{G}^\epsilon )^*(s). \end{aligned}$$

(B.93)

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

$$\begin{aligned} \varvec{H}^\epsilon (s)&= \int _0^\infty e^{\varvec{A}_2(s, \varvec{X}(s)) y} \bigg ( -(\varvec{\Sigma }_2 \varvec{\Sigma }_2^*)(s, \varvec{X}(s)) + (\varvec{\sigma }_2 \varvec{\sigma }_2^*)(s, \varvec{x}^\epsilon (s),\epsilon ) \nonumber \\&\quad + \varvec{G}^\epsilon (s)\varvec{J}_3^\epsilon (s) + \varvec{J}_3^\epsilon (s) (\varvec{G}^\epsilon )^*(s)\bigg ) e^{\varvec{A}^*_2(s, \varvec{X}(s))y} \mathrm{d}y. \end{aligned}$$

(B.94)

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

$$\begin{aligned} {\mathbb {E}}_1 \sup _{t \in [0,T]} \int _0^t \Vert \varvec{J}^\epsilon _3(s) - \varvec{J}(s)\Vert ^p_F \mathrm{d}s&\le C(\epsilon , p, T) \int _0^T {\mathbb {E}}_1 \sup _{u \in [0,s]} |\varvec{x}^\epsilon (u) - \varvec{X}(u)|^p \mathrm{d}s \nonumber \\&\quad + D(p,T)[\alpha _2(\epsilon )^p \mathbb {1}_{\varvec{a}_2 \ne \varvec{A}_2} + \gamma _2(\epsilon )^p \mathbb {1}_{\varvec{\sigma }_2 \ne \varvec{\Sigma }_2}] \end{aligned}$$

(B.95)

on the set $S_1$, where $C(\epsilon , p, T) = O(1)$ as $\epsilon \rightarrow 0$ and D(p, T) is a positive constant, $\alpha _2(\epsilon )$ and $\gamma _2(\epsilon )$ are from Assumption A.5.

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

$$\begin{aligned} M_6&\le L_6(\epsilon ,p,T) \int _0^T {\mathbb {E}}_1 \sup _{u \in [0,s]} |\varvec{x}^\epsilon (u)-\varvec{X}(u) |^p \mathrm{d}s \nonumber \\&\quad + C_6(p,T)\bigg [\alpha _1(\epsilon )^p \mathbb {1}_{\{\varvec{a}_1 \ne \varvec{A}_1\}} + \alpha _2(\epsilon )^p \mathbb {1}_{\{\varvec{a}_2 \ne \varvec{A}_2\}} + \gamma _2(\epsilon )^p \mathbb {1}_{\{\varvec{\sigma }_2 \ne \varvec{\Sigma }_2\}} \nonumber \\&\quad + \theta _1(\epsilon )^p \mathbb {1}_{\{(\varvec{a}_1)_{\varvec{x}} \ne (\varvec{A}_1)_{\varvec{x}}\}} + \theta _2(\epsilon )^p \mathbb {1}_{\{(\varvec{a}_2)_{\varvec{x}} \ne (\varvec{A}_2)_{\varvec{x}}\}} \bigg ], \end{aligned}$$

(B.96)

on $S_1$, where $L_6(\epsilon ,p,T)=O(1)$ as $\epsilon \rightarrow 0$, $C_6(p,T)$ is a positive constant, and $\alpha _i(\epsilon )$, $\theta _i(\epsilon )$ ($i=1,2$) and $\gamma _2(\epsilon )$ are from Assumption A.5.

Collecting the above estimates for the $M_k$, we obtain:

$$\begin{aligned} R_1&\le C_1(p,T) \bigg ( {\mathbb {E}}_1 \sup _{t \in [0,T]} |\epsilon \varvec{v}^\epsilon (t)|^p \nonumber \\&\quad + \alpha _1(\epsilon )^p \mathbb {1}_{\{\varvec{a}_1 \ne \varvec{A}_1\}} + \alpha _2(\epsilon )^p \mathbb {1}_{\{\varvec{a}_2 \ne \varvec{A}_2\}} + \gamma _2(\epsilon )^p \mathbb {1}_{\{\varvec{\sigma }_2 \ne \varvec{\Sigma }_2\}} + \nonumber \\&\quad + \theta _1(\epsilon )^p \mathbb {1}_{\{(\varvec{a}_1)_{\varvec{x}} \ne (\varvec{A}_1)_{\varvec{x}}\}} + \theta _2(\epsilon )^p \mathbb {1}_{\{(\varvec{a}_2)_{\varvec{x}} \ne (\varvec{A}_2)_{\varvec{x}}\}} \bigg ) \nonumber \\&\quad + C_2(\epsilon , p,T) \int _0^T {\mathbb {E}}_1 \sup _{u \in [0,s]} |\varvec{x}^\epsilon (u) - \varvec{X}(u)|^p \mathrm{d}s + C_3(p,T) M_4 \end{aligned}$$

(B.97)

on $S_1$, where $C_1(p,T)$ and $C_3(p,T)$ are constants, $C_2(\epsilon , p, T) = O(1)$ as $\epsilon \rightarrow 0$, and $M_4$ satisfies the bound in (B.87).

Using all the estimates for the $R_i$, we have:

$$\begin{aligned}&{\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} |\varvec{x}^{\epsilon }(t) - \varvec{X}(t)|^p \right] = {\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} \sum _{k=1}^{n_1} |[\varvec{x}^{\epsilon }- \varvec{X}]_k(t)|^p \right] \end{aligned}$$

(B.98)

$$\begin{aligned}&\quad \le n_1 \max _{k=1,\dots ,n_1} \left\{ {\mathbb {E}}_1 \sup _{t \in [0,T]} |[\varvec{x}^{\epsilon }- \varvec{X}]_k(t)|^p \right\} \end{aligned}$$

(B.99)

$$\begin{aligned}&\quad \le L(\epsilon ,p,T,n_1) \int _0^T {\mathbb {E}}_1 \sup _{u \in [0,s]} |\varvec{x}^\epsilon (u)-\varvec{X}(u) |^p \mathrm{d}s \nonumber \\&\qquad + C(p,T,n_1) \bigg ( \epsilon ^{p r_0} + {\mathbb {E}}_1 \sup _{t \in [0,T]} |\epsilon \varvec{v}^\epsilon (t)|^p + M_4 \nonumber \\&\qquad + \alpha _1(\epsilon )^p \mathbb {1}_{\{\varvec{a}_1 \ne \varvec{A}_1\}} + \alpha _2(\epsilon )^p \mathbb {1}_{\{\varvec{a}_2 \ne \varvec{A}_2\}} + \gamma _1(\epsilon )^p \mathbb {1}_{\{\varvec{\sigma }_1 \ne \varvec{\Sigma }_1\}} \nonumber \\&\qquad + \gamma _2(\epsilon )^p \mathbb {1}_{\{\varvec{\sigma }_2 \ne \varvec{\Sigma }_2\}} + \beta _1(\epsilon )^p \mathbb {1}_{\{\varvec{b}_1 \ne \varvec{B}_1\}} + \beta _2(\epsilon )^p \mathbb {1}_{\{\varvec{B}_2 \ne \varvec{B}_2\}} \nonumber \\&\qquad +\theta _1(\epsilon )^p \mathbb {1}_{\{(\varvec{a}_1)_{\varvec{x}} \ne (\varvec{A}_1)_{\varvec{x}}\}} + \theta _2(\epsilon )^p \mathbb {1}_{\{(\varvec{a}_2)_{\varvec{x}} \ne (\varvec{A}_2)_{\varvec{x}}\}} \bigg ) \end{aligned}$$

(B.100)

$$\begin{aligned}&\quad \le L(\epsilon ,p,T,n_1) \int _0^T {\mathbb {E}}_1 \sup _{u \in [0,s]} |\varvec{x}^\epsilon (u)-\varvec{X}(u) |^p \mathrm{d}s \nonumber \\&\qquad + C(p,T,n_1)\epsilon ^{r}, \end{aligned}$$

(B.101)

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

Finally, applying the Gronwall lemma gives:

$$\begin{aligned}&{\mathbb {E}}_1 \left[ \sup _{t \in [0,T]} |\varvec{x}^{\epsilon }(t) - \varvec{X}(t)|^p \right] \le \epsilon ^{r} \cdot C(p,T,n_1) e^{L(\epsilon , p, T,n_1) T} \end{aligned}$$

(B.102)

on $S_1$.

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

$$\begin{aligned} {\mathbb {E}}_1 \left[ \sup _{t\in [0,T]} |\varvec{x}^\epsilon (t)-\varvec{X}(t)|^p \right]&\le \bigg [ {\mathbb {E}}_1 \bigg ( \sup _{t\in [0,T]} |\varvec{x}^\epsilon (t)-\varvec{X}(t)|^p \bigg )^{q/p} \bigg ]^{p/q} \end{aligned}$$

(B.103)

$$\begin{aligned}&= O(\epsilon ^\beta ), \end{aligned}$$

(B.104)

for all $0< \beta < p'$, as $\epsilon \rightarrow 0$. The statement on convergence in probability follows from Lemma 1 in [43]. $\square$

Appendix C: An Implementation of Algorithm 6.3 Under Assumption 6.4

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

$$\begin{aligned} \varvec{Q}_i(z) = z^{d_i}\varvec{I} + \varvec{a}_{i,d_i-1}z^{d_i-1}+\cdots +\varvec{a}_{i,1} z + \varvec{a}_{i,0}, \end{aligned}$$

(C.1)

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

$$\begin{aligned} \varvec{a}_{i,0}&= \prod _{k=1}^{d_i} \varvec{\Gamma }_{i,k}, \nonumber \\ \varvec{a}_{i,1}&= \sum _{k_1, \dots , k_{d_i-1}=1,\dots ,d_i: k_1> \cdots> k_{d_i-1}} \varvec{\Gamma }_{i,k_{1}} \varvec{\Gamma }_{i,k_2} \cdots \varvec{\Gamma }_{i,k_{d_i-1}}, \nonumber \\&\vdots \nonumber \\ \varvec{a}_{i,d_i-2}&= \sum _{k_1,k_2=1,\dots ,d_i: k_1>k_2} \varvec{\Gamma }_{i,k_1} \varvec{\Gamma }_{i,k_2}, \nonumber \\ \varvec{a}_{i,d_i-1}&= \sum _{k=1}^N \varvec{\Gamma }_{i,k}. \end{aligned}$$

(C.2)

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

$$\begin{aligned} \varvec{H}_i = [\varvec{0} \ \cdots \ \varvec{0} \ \ \varvec{B}_{l_i} \ \ \varvec{0} \ \cdots \ \varvec{0}] \in {\mathbb {R}}^{p_i \times p_i d_i}, \end{aligned}$$

(C.3)

where $\varvec{B}_{l_i}$ is in the $l_i$th slot,

$$\begin{aligned} \varvec{F}_i= & {} \begin{bmatrix} \varvec{0} &{}\quad -\varvec{I} &{}\quad \\ &{}\quad \varvec{0} &{}\quad -\varvec{I} &{}\quad \\ &{}\quad &{}\quad \ddots &{}\quad \ddots \\ &{}\quad &{}\quad &{}\quad \varvec{0} &{}\quad -\varvec{I} \\ \varvec{a}_{i,0} &{}\quad \varvec{a}_{i,1} &{}\quad \dots &{}\quad \varvec{a}_{i,d_{i-2}} &{}\quad \varvec{a}_{i,d_{i-1}} \end{bmatrix} \in {\mathbb {R}}^{p_i d_i \times p_i d_i}, \end{aligned}$$

(C.4)

$$\begin{aligned} \varvec{G}_i= & {} [\varvec{0} \ \cdots \ \varvec{0} \ \ \varvec{I}]^* \in {\mathbb {R}}^{p_i d_i}. \end{aligned}$$

(C.5)

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

$$\begin{aligned} \varvec{F}_i\varvec{M}_i + \varvec{M}_i \varvec{F}_i^* =: \varvec{\Sigma }_i \varvec{\Sigma }_i^* \ge 0, \ \ \varvec{M}_i \varvec{H}_i^* = \varvec{G}_i \end{aligned}$$

(C.6)

for $\varvec{M}_i = \varvec{M}_i^*$ [74].

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

$$\begin{aligned} d(\varvec{\beta }^i)'_t = -\varvec{\Gamma }_i' (\varvec{\beta }^i)'_t \mathrm{d}t + \varvec{\Sigma }_i' \mathrm{d}\varvec{W}_t^{(q_4)}, \end{aligned}$$

(C.7)

and one can compute the spectral density to be:

$$\begin{aligned} \varvec{{\mathcal {S}}}_i(\omega ) = \varvec{\Phi }_i(-i\omega )\varvec{\Phi }^*_i(i\omega ) =\varvec{B}_{l_i} \omega ^{2 l_i} ((\omega ^2\varvec{I}+\varvec{\Gamma }_{i,1})^2) \cdots (\omega ^2\varvec{I}+\varvec{\Gamma }_{i,d_i})^2) )^{-1} \varvec{B}_{l_i}^*. \end{aligned}$$

(C.8)

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

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