Reduced Markovian Descriptions of Brownian Dynamics: Toward an Exact Theory

1 Introduction

The derivation of a contracted description of a Brownian particle subject to a confining potential is a long-standing problem in statistical mechanics, which dates back to an old question originally posed by Uhlenbeck and Ornstein [1]. For more background details and a general outline of the methods, we refer the reader to the seminal paper by van Kampen [2] as well as to the more recent reviews [3, 4]. The Smoluchowski equation stands as a prominent example of a reduced description in the high friction regime, where the momentum variable rapidly thermalizes and the statistics of the particle is determined only by the distribution in the configuration space. In the last decades, a large research endeavour [5–7] has pointed towards the derivation of corrections to the Smoluchowski formula for finite values of the friction constant. A classical iterative scheme, described e.g. in [8], derives solutions of the Kramers equation in terms of matrix continued fractions, via an expansion in powers of the inverse friction coefficient. Further guidelines on the derivation of the Smoluchowski equation from the Kramers equation can also be found in [9–11]. A different approach, developed by Titulaer [12], implements a Chapman-Enskog reduction scheme on the Fokker-Planck equation. A systematic use of the same procedure applied for the adiabatic elimination of fast variables was also considered in [13]. More recently, a non-local version of the Smoluchowski equation was also obtained from the Kramers equation through the Chapman-Enskog procedure in [14]. In the set-up of kinetic theory of gases, the Chapman-Enskog method has proved successful in the derivation of the Euler and the Navier-Stokes equations of fluid dynamics from the Boltzmann equation. However, as it was first demonstrated by Bobylev for Maxwell’s molecules [15], the Chapman-Enskog expansion is prone to small wavelength instabilities: namely, sufficiently short acoustic waves increase with time instead of decaying. This creates difficulties for an extension of hydrodynamics, as derived from a kinetic description, in the regime of finite Knudsen numbers, where the Navier-Stokes approximation is inapplicable. The study of various kinetic models [16–19] has later revealed that such instabilities can be cured by taking into account also the remote terms of the expansion. The resulting hydrodynamic equations, obtained from an exact summation of the Chapman-Enskog scheme, are indeed stable for all wavelengths, at variance with the finite-order approximations. The structure of the slow invariant manifold for systems displaying an intrinsic separation of time scales can be determined through a non-perturbative reduction procedure known as the invariant manifold method [20]. A key ingredient, in the method, is an equation of dynamic invariance which, for a class of kinetic models known as linearized Grad’s moment systems, was shown to lead to the same result as the exact summation of the Chapman-Enskog expansion [21, 22]. Thus, the plan of this work is to adopt, in the context of Brownian dynamics, the invariant manifold set-up to derive meaningful corrections to the Smoluchowski equation, while preserving the Markovian structure of the original process. Unlike previous reduction schemes, the proposed procedure operates on the deterministic component of the dynamics, whereas the stochastic terms are properly handled via the Fluctuation-Dissipation theorem. Explicit calculations can be carried out for the Brownian oscillator model, in which the equivalence between the condition of dynamic invariance and the exact summation of the Chapman-Enskog expansion is also established. The special case with a harmonic potential can hence be used as a test-bed to illustrate the general formalism. With other models, instead, quasi-equilibrium manifolds typically constitute the starting point towards an iterative method of solution of the equation of invariance.

We envisage further developments of similar reduction schemes in the direction of coarse-graining of interacting particle systems as well as of partial differential equations with randomly fluctuating coefficients. Such research line might possibly connect this work with periodic and/or random homogenization questions; see [23–28] for recent applications of reduction schemes to statistical physics and epidemiological models. We also refer the reader to [29–32] for related matters, as well as to [33, 34] for applications of similar methods to the reduction of complex dynamics connected to climate change topics, where the need of developing innovative reduction techniques is growing. In this line of thinking, a rigorous characterization of slow invariant manifolds of random dynamical systems can be found in [35]. Other relevant questions related to the procedure of model reduction are also addressed in this work. One, for instance, concerns the derivation of a quantitative estimate of the reduction error stemming from the application of the method. The role of the initial data and of the defect of invariance will be properly highlighted. Furthermore, as the contracted description retains just some of the observables of the original model, it is a non-trivial task to quantify to which extent a certain reduced dynamics yields a response to (small) perturbations comparable to that expressed by the original dynamics. The study of the linear response formalism will enable us to answer also this question.

This paper is organized as follows. In Section 2, we introduce the model and also illustrate the general set-up of the invariant manifold method and the related Chapman-Enskog scheme. In Section 3 we present some iterative schemes, based on the Maximum Entropy principle, which can be used to solve the invariance equation for a class of Brownian dynamics. We focus, in particular, on the case of a harmonically bound Brownian particle, where the reduction procedure can be outlined in detail. We also provide a quantitative pointwise representation of the reduction error, and derive response functions due to the original and the reduced dynamics of the Brownian oscillator model. Finally, we draw our conclusions in Section 4.

2 The Model

In this work we exploit the Dynamic Invariance principle to derive reduced descriptions for a class of Brownian dynamics. Specifically, we consider the Brownian dynamics of particles subject to a power law potential V(x) = xⁿ, n ≥ 2 an integer, described by a system of stochastic differential equations (SDEs) written in the Itô form:

$\begin{array}{cc}\hfill dx\left(t\right)& =v\left(t\right)dt\hfill \\ \hfill dv\left(t\right)& =-\frac{1}{m}{\partial }_{x}V\left(x\right)dt-\gamma v\left(t\right)dt+\sqrt{2D{\gamma }^2}dW\left(t\right),\hfill \end{array}(1)$

where W(t) is a one-dimensional Wiener process, m is the mass, γ is the friction constant, D = (βmγ)⁻¹ is the diffusion coefficient, and β is the inverse temperature of the system. We recall that the leading high friction approximation of the set of Eq. 1 is commonly written in the Smoluchowski form:

$dx\left(t\right)=-\frac{1}{m\gamma }{\partial }_{x}V\left(x\right)dt+\sqrt{2D}dW\left(t\right).(2)$

We hence seek for a reduced description resembling the structure of Eq. 2 and based on the following SDE:

$dx\left(t\right)=-\chi {\partial }_{x}V\left(x\right)dt+\sqrt{2{\mathcal{D}}_r}dW\left(t\right),(3)$

where χ is the mobility and ${\mathcal{D}}_r$ the diffusion coefficient of the reduced dynamics. The two coefficients χ and ${\mathcal{D}}_r$ are expected to reduce to (mγ)⁻¹ and D, respectively, in the overdamped limit expressed by Eq. 2. The derivation of the mobility χ, in Eq. 3, shall be pursued in Section 3 by exploiting the Dynamic Invariance principle, which is shortly reviewed next.

2.1 The Dynamic Invariance Principle

The invariant manifold method is a procedure of model reduction that was originally introduced as a special analytical perturbation technique in the KAM theory of integrable Hamiltonian systems [36–38]. The method was later exploited in the kinetic theory of gases to derive the evolution equations of the hydrodynamic fields from the Boltzmann equation or related kinetic models [18, 21, 22]. The basic picture underlying the invariant manifold method can be shortly summarized as follows, see Refs. [20, 39]. There exists a manifold of slow motions, in the phase space of the system, parameterized by a set of distinguished macroscopic variables, which is positively invariant: if a trajectory starts on the manifold at time t₀, it will remain on the manifold for all times t 〉 t₀. Trajectories starting from arbitrary initial conditions quickly reach a neighborhood of the manifold, and then evolve along such slow manifold, until the equilibrium state is eventually attained. More explicitly, for kinetic equations we let f(r, v, t) denote the (single-particle) distribution function, whose evolution in a domain U is described by the equation:

${\partial }_{t}f=J\left(f\right).(4)$

Let m: f → M be a surjective linear map, with M denoting the macroscopic variables (moment fields) and also let f_M denote a manifold parametrized by a set of macroscopic fields M. We look for an invariant manifold f_M obeying the self-consistency condition m(f_M) = M. The “microscopic” and “macroscopic” time derivatives of f on the manifold f_M are defined as:

${\partial }_t^{\left(micro\right)}{f}_M≔J\left(f_M\right),(5)$

${\partial }_t^{\left(macro\right)}{f}_M≔\left(D_M{f}_M\right)m\left(J\left(f_M\right)\right)\left.\right),(6)$

where the differential D_Mf_M, in Eq. 6, is evaluated at the point M = m(f_M). While J(f_M) in Eq. 5, corresponds to a value of the vector field J evaluated on the manifold f_M, Eq. 6 codifies a chain rule: one computes the time derivative of the moment M via the map m, as $\dot{M}=m(J(f_M))$, and then exploits the time dependence of f which is expressed through the time dependence of M. The invariance equation is written in the form:

${\partial }_t^{\left(micro\right)}{f}_M={\partial }_t^{\left(macro\right)}{f}_M,(7)$

or, alternatively, as:

${\mathrm{\Delta }}_M≔{\partial }_t^{\left(micro\right)}{f}_M-{\partial }_t^{\left(macro\right)}{f}_M=0,(8)$

where Δ_M is called defect of invariance. The Dynamic Invariance principle requires that the equality in Eq. 8 is fulfilled for any values of the macroscopic variables M. Solutions to Eq. 8 have been obtained from the study of various kinetic models, see e.g., [16, 17, 19]. Determining the structure of f_M constitutes an instance of the closure problem in kinetic theory [40]. Note, indeed, that if a solution of Eq. 8 is found, then the moments M obey the following closed system of evolution equations:

$\dot{M}=m\left(J\left(f_M\right)\right).(9)$

We also recall that if a functional E(f) is conserved for the dynamics (Eq. 4), then E(f_M) remains constant along the trajectories of the reduced system described by Eq. 9. Moreover, if the time derivative of a functional H(f) is nonpositive due to the dynamics (Eq. 4), then so is also the time derivative of the functional H(f_M) due to the reduced dynamics (Eq. 9).

Appropriate iterative schemes have been developed to solve the invariance Eq. 8. One of the first systematic procedures of constructing invariant manifolds was the celebrated Chapman-Enskog method for the Boltzmann equation [41], whose main steps are also shortly recalled here. One starts with the singularly perturbed dynamics:

${\partial }_{t}f+A\left(f\right)=\frac{1}{\epsilon }Q\left(f\right),(10)$

with ɛ 〉 0 a small parameter. One requires that m(Q(f)) = 0 and that, for each M ∈ m(U), the equation Q(f) = 0, with m(f) = M, has a unique solution denoted by $f_M^{eq}$ (corresponding to the Maxwellian distributions, in Boltzmann’s theory), which is also asymptotically stable and globally attracting for the fast dynamics

${\partial }_{t}f=\frac{1}{\epsilon }Q\left(f\right).(11)$

The invariance Eq. 8 can be adapted to the singularly perturbed system (Eq. 10) in the form:

$\frac{1}{\epsilon }Q\left(f_M\right)=A\left(f_M\right)+\left(D_M{f}_M\right)\left(m\left(A\left(f_M\right)\right)\right).(12)$

In the Chapman-Enskog scheme, a solution to Eq. 12 is sought in the form of a series in powers of a small parameter ɛ:

$f_M=f_M^{eq}+\sum _{i=1}^{\infty }{\epsilon }^i{f}_M^{\left(i\right)}.(13)$

In the set-up of Boltzmann’s theory, the zero-order approximation $f_M\simeq f_M^{eq}$ leads, upon integrating the Boltzmann equation over the velocity space, to the inviscid Euler equations of fluid dynamics, whereas the first-order correction $f_M\simeq f_M^{(eq)}+\epsilon f_M^{(1)}$ gives rise to the compressible Navier-Stokes equations supplied with transport coefficients which depend on the underlying collision model.

3 Results

The aim of this Section is to exploit the method of the invariant manifold to derive suitable expressions for the mobility χ introduced in Eq. 2. In Section 3.1 we thus address the general case described by Eq. (1), and outline useful iterative schemes based on the Maximum Entropy principle. The Brownian oscillator model, corresponding to the case n = 2, is studied in detail in Section 3.2. Next, the analysis of the reduction error and of correlation functions is deferred to Sections 3.3, 3.4.

3.1 Quasi-Equilibrium Closures

We start with the Fokker-Planck equation associated to the SDE (Eq. 3), which reads:

$\frac{\partial \rho \left(x,t\right)}{\partial t}={\mathcal{D}}_r\frac{\partial }{\partial x}\left(\rho \left(x,t\right){\partial }_{x}U\left(x\right)+\frac{\partial }{\partial x}\rho \left(x,t\right)\right),(14)$

where we enforced the Einstein relation $\chi =\beta {\mathcal{D}}_r$ and set U(x) = βV(x). Letting x(t) = (x(t), v(t)), we denote by 〈…〉 the conditional average referring to a deterministic initial datum x(0) = x = (x, v). We then write M(t) ≡ 〈∂_xU(x)〉 and introduce the closure:

$⟨v\left(t\right)⟩=-\chi M\left(t\right),(15)$

where M(t) plays the role of a slow variable, driving the evolution of the fast variable 〈v(t)〉. Application of the Dynamic Invariance principle, introduced in Section 2.1, requires the evaluation of the two operators ${\partial }_t^{(micro)}⟨v(t)⟩$ and ${\partial }_t^{(macro)}⟨v(t)⟩$ for the dynamics obtained upon averaging Eq. 1 over noise. We thus find:

${\partial }_t^{\left(micro\right)}\langle v\left(t\right)\rangle ≔-\frac{1}{m}M\left(t\right)+\gamma \chi M\left(t\right)(16)$

${\partial }_t^{\left(macro\right)}\langle v\left(t\right)\rangle ≔-\chi \dot{M}\left(t\right)(17)$

Upon equalizing the two expressions in Eqs 16, 17 one establishes the invariance equation

$-\frac{1}{m}M\left(t\right)+\gamma \chi M\left(t\right)+\chi \dot{M}\left(t\right)=0,(18)$

to be solved for the mobility χ. Note that the classical form 1/(mγ) of the mobility is recovered from Eq. 18 as γ → ∞, m → 0 with γm finite (provided χ remains bounded as well). In order to solve Eq. 18, and find corrections to the Smoluchowski equation, it is thus necessary to determine an explicit expression for $\dot{M}(t)$. This task can be pursued iteratively via the Maximum Entropy method (MaxEnt), often invoked in statistical mechanics; see e.g., [20, 42]. We start by defining the entropy:

$S\left[\rho \right]=-\int \rho \left(x,t\right)\ln \frac{\rho \left(x,t\right)}{{\rho }_0\left(x\right)}dx,(19)$

which is a monotonically growing function attaining its maximum at equilibrium, i.e., when ρ = ρ₀. Let M = {M₀, …, M_k} be a set of linearly independent moments of ρ defined as:

$\int m_i\left(x\right)\rho \left(x,t\right)dx=M_i\left(t\right),i=0,\dots ,k,(20)$

where the m_i’s are the microscopic densities of the moments, with m₀ = 1. The so-called quasi-equilibrium density ρ*(x, M) is obtained by maximising the entropy S[ρ] under the constraints of fixed M, which yields:

${\rho }^{*}\left(x,M\right)={\rho }_0\left(x\right)\exp \left[\sum _{i=0}^k{\mathrm{\Lambda }}_i{m}_i\left(x\right)\right](21)$

where Λ = {Λ₁, …, Λ_k} denotes the set of Lagrange multipliers which depend on the set of moments M. The quasi-equilibrium projector is then defined as:

${\mathrm{\Pi }}^{*}•=\sum _{i=0}^k\frac{\partial {\rho }^{*}}{\partial M_i}\int m_i\left(x\right)•dx(22)$

By acting with the projector Π* on both sides of Eq. 14 one obtains the following moment equations in the quasi-equilibrium approximation:

${\dot{M}}_i=-\sum _{j=0}^k{\mathrm{\Lambda }}_j\int {\partial }_x{m}_i\left(x\right){\rho }^{*}\left(x,t\right){\partial }_x{m}_j\left(x\right)dx.(23)$

A step forward can be made by splitting the set of moments as M = M_I ∪ M_II, with M_I = {M₀, …, M_ℓ} and M_II = {M_ℓ+1, …, M_k}. We assume that the (first) quasi-equilibrium distribution ρ*(x) can be derived explicitly for the set M_I, i.e., ρ* = ρ*(x, M_I), and we hence seek for the second quasi-equilibrium closure in the form ρ = ρ*(1 + φ). An expansion of the functional (Eq. 19) in a neighbourhood of ρ* to second order in φ yields:

$\mathrm{\Delta }S\left[\phi \right]=-\int {\rho }^{*}\phi \left[1+\ln \frac{{\rho }^{*}}{{\rho }_0}\right]dx-\frac{1}{2}\int {\rho }^{*}{\phi }^{2}dx.(24)$

The deviation φ is determined from the following maximization problem, called “Triangle MaxEnt approximation” [43]:

$\mathrm{\Delta }S\to \mathrm{m}\mathrm{a}\mathrm{x},\int m_I\left(x\right){\rho }^{*}\phi dx=0,\int m_{II}\left(x\right){\rho }^{*}\phi dx=\mathrm{\Delta }{M}_{II},(25)$

where ΔM_II ≡ M_II − M_II(M_I) represents the deviations of the moments M_II from the values attained in the first quasi-equilibrium state.

To solve the invariance Eq. 18, we restrict to the one-moment quasi-equilibrium closures, i.e., we set M_I = M₀ (i.e., ρ*(x) = ρ₀(x)) and M_II = M, with ∫m(x)ρ dx = M. We thus find:

$\frac{\partial \phi }{\partial t}=\mathcal{L}\phi ,,\mathcal{L}•={\rho }_0^{-1}{\mathcal{D}}_r\frac{\partial }{\partial x}{\rho }_0\frac{\partial }{\partial x}•.(26)$

The triangle one-moment quasi-equilibrium distribution is found in the form:

${\phi }^{\left(0\right)}\left(x,t\right)=m^{\left(0\right)}\left(x\right)\left(M\left(t\right)-{⟨M⟩}_0\right),(27)$

with

$m^{\left(0\right)}\left(x\right)=\frac{m\left(x\right)-{⟨m⟩}_0}{{⟨m^2⟩}_0-{⟨m⟩}_0^2}.(28)$

We can thus construct a refinement of the quasi-equilibrium dynamics (Eq. 23), by defining the new projection operator:

${\mathrm{\Pi }}^{\left(0\right)}•={\rho }_0\left(x\right)\frac{m^{\left(0\right)}\left(x\right)}{{⟨m^{\left(0\right)}{m}^{\left(0\right)}⟩}_0}\int m^{\left(0\right)}\left(x\right)•dx.(29)$

After inserting the expression (Eq. 27) in (Eq. 26) and upon acting on both sides with Π⁽⁰⁾, one finally obtains:

$\dot{M}=-{\ell }_0\left(M\left(t\right)-{⟨M⟩}_0\right),(30)$

with

${\ell }_0=-\frac{{⟨{\partial }_x{m}^{\left(0\right)}{\mathcal{D}}_r{\partial }_x{m}^{\left(0\right)}⟩}_0}{{⟨m^{\left(0\right)}{m}^{\left(0\right)}⟩}_0},(31)$

which defines the relaxation dynamics of the moment M to the corresponding value attained in the equilibrium state. We observe that the strategy exploited so far can be prosecuted further, by using the iterative scheme discussed in Ref. [43]:

$\left(1-{\mathrm{\Pi }}^{\left(j\right)}\right)\mathcal{L}{m}^{\left(j+1\right)}=0,(32)$

with j = 0, 1, …, where $m^{(j+1)}=m^{\left(\right.(j}+{\mu }^{(j+1)}$, and in which the orthogonality condition ${⟨{\mu }^{(j+1)}{m}^{\left(\right.(j}⟩}_0=0$ is exploited. At the k-th iteration step, a refinement of the inverse relaxation time (Eq. 31) attains the form:

${\ell }_k=-\frac{{⟨{\partial }_x{m}^{\left(k\right)}{\mathcal{D}}_r{\partial }_x{m}^{\left(k\right)}⟩}_0}{{⟨m^{\left(k\right)}{m}^{\left(k\right)}⟩}_0}.(33)$

The sequence ${({\ell }_k)}_{k\in \mathbb{N}}$ is found to converge to the eigenvalue of the operator $\mathcal{L}$ with the minimal non-zero absolute value.

In the next Section we will address the case with n = 2, corresponding to a harmonically bound particle, which can be solved explicitly without adopting the foregoing iterative scheme based on the MaxEnt principle.

3.2 The Brownian Oscillator

We now turn to study in detail the case of a Brownian particle bounded in a harmonic potential $V(x)=1/2m{\omega }_0^2{x}^2$, which constitutes one of the classical exactly solvable models of nonequilibrium statistical mechanics [44]. We refer also the reader to Ref. [45] for some recent inspection of the ergodic properties of Ornstein-Uhlenbeck systems. The dynamics of the Brownian oscillator is described by a system of linear SDEs:

$\begin{array}{cc}\hfill dx\left(t\right)& =v\left(t\right)dt\hfill \\ \hfill dv\left(t\right)& =-{\omega }_0^{2}x\left(t\right)dt-\gamma v\left(t\right)dt+\sqrt{2D{\gamma }^2}dW\left(t\right),\hfill \end{array}(34)$

where ${\omega }_0=\sqrt{k/m}$ is the natural frequency of the oscillator with mass m and elastic constant k, see Supplementary Appendix S1 for details. We shall refer to Eq. 34 as the original dynamics of the Brownian oscillator. We introduce the drift matrix M, defined as:

$\mathbf{M}=\left(\begin{array}{cc}\hfill 0\hfill & \hfill -1\hfill \\ \hfill {\omega }_0^2\hfill & \hfill \gamma \hfill \end{array}\right),(35)$

whose eigenvalues read:

${\lambda }_{\pm }=\frac{\gamma \pm {\gamma }_s}{2},(36)$

with ${\gamma }_s=\sqrt{{\gamma }^2-4{\omega }_0^2}={\lambda }_{+}-\lambda -$. Henceforth we shall restrict our discussion to the overdamped regime, namely the region in the parameter space in which γ_s is real and larger than zero.

An exact reduced description, not requiring a separation of time scales, is available for the Brownian oscillator model [46]. This is obtained by integrating over time the second of Eq. 34 and by then inserting the obtained expression in the first equation. The resulting reduced dynamics, expressed in terms of the configuration variable x(t), turns out to be non-Markovian. Nevertheless, in the regime of high friction, and for times much longer than γ⁻¹, the Markovian structure of the reduced dynamics can be restored [7, 47].

Another contracted description of the model can instead be derived by considering the overdamped limit of Eq. 34 [48], that is worth briefly recalling. By letting x^ɛ(t) = x(ɛ⁻¹t), with ɛ = γ⁻¹, the original dynamics can be rescaled as follows:

$d{x}^{\epsilon }\left(t\right)={\epsilon }^{-1}vdt(37)$

$d{v}^{\epsilon }\left(t\right)=-{\epsilon }^{-1}{\omega }_0^2{x}^{\epsilon }\left(t\right)dt-{\epsilon }^{-2}{v}^{\epsilon }\left(t\right)dt+{\epsilon }^{-1}\sqrt{2{\left(\beta m\right)}^{-1}}dW\left(t\right),(38)$

where we exploited the scaling dW(ϵ⁻¹t) = ϵ^−1/2dW(t). It thus holds:

${\epsilon }^{-1}{v}^{\epsilon }\left(t\right)dt=-{\omega }_0^2{x}^{\epsilon }\left(t\right)dt+\sqrt{2{\left(\beta m\right)}^{-1}}dW\left(t\right)+\mathcal{O}\left(\epsilon \right),(39)$

and hence,

$d{x}^{\epsilon }\left(t\right)=-{\omega }_0^2{x}^{\epsilon }\left(t\right)dt+\sqrt{2{\left(\beta m\right)}^{-1}}dW\left(t\right)+\mathcal{O}\left(\epsilon \right).(40)$

As ϵ → 0, Eq. 40 leads to the Smoluchowski equation for the Brownian oscillator, which, after turning back to the original variables, attains the well-known structure:

$dx\left(t\right)=-\frac{{\omega }_0^2}{\gamma }x\left(t\right)dt+\sqrt{2D}dW\left(t\right).(41)$

Let us now turn to illustrate our reduction scheme, based on the Dynamic Invariance principle. We aim at setting up a reduced description which formally resembles the structure of Eq. 41, and is based on the linear SDE:

$dx\left(t\right)=-\alpha x\left(t\right)dt+\sqrt{2{\mathcal{D}}_r}dW\left(t\right),(42)$

where α and ${\mathcal{D}}_r$ denote the drift and diffusion coefficients of the reduced dynamics, respectively. These coefficients will properly recover the two corresponding expressions ${\omega }_0^2/\gamma$ and D attained in the overdamped limit, cf. Eq. 41. We remark that unlike the inverse friction expansion method, see e.g., [8], which typically starts from the Kramers equation associated to Eq. 34 and yields Eq. 41 as the leading-order term of an expansion in powers of ɛ = γ⁻¹, our scheme applies to the deterministic component of the Brownian dynamics. The drift coefficient α will be obtained from the solution of an invariance equation, which is derived either from an exact summation of the Chapman-Enskog expansion or, equivalently, by the application of the Dynamic Invariance principle. The coefficient ${\mathcal{D}}_r$ is then found by exploiting the Fluctuation-Dissipation theorem. In the sequel we shall address, separately, the two distinct steps of the reduction procedure.

3.2.1 Exact Summation of the Chapman-Enskog Expansion

The Chapman-Enskog scheme, introduced in Section 2.1, can be adapted to the reduction of the Brownian oscillator model as follows. The procedure starts from averaging Eq. 34 over noise,

$\langle \dot{\mathbf{x}}\left(t\right)\rangle =-\mathbf{M}\langle \mathbf{x}\left(t\right)\rangle .(43)$

One regards 〈x(t)〉 as the variable characterizing the reduced description, and assumes that the evolution of the fast variable 〈v(t)〉, after the initial layer, reaches a neighborhood of the slow manifold parameterized by 〈x(t)〉. Next, the variable 〈v(t)〉 is expanded in powers of ɛ = γ⁻¹, viz.

$⟨v\left(t\right)⟩=\sum _{j=0}^{\infty }{\epsilon }^j{v}^{\left(j\right)}\left(t\right).(44)$

The coefficients v^(j)(t) are found from the recurrence procedure

$v^{\left(j+1\right)}=-\sum _{k=0}^j{D}_{CE}^{\left(k\right)}{v}^{\left(j-k\right)},j\ge 1,(45)$

where the Chapman-Enskog operators $D_{CE}^{(k)}$ act on the coefficients v^(j) as follows:

$D_{CE}^{\left(k\right)}{v}^{\left(j\right)}≔\frac{\partial v^{\left(j\right)}}{\partial ⟨x⟩}{v}^{\left(k\right)}.(46)$

The recurrence Eq. 45 starts with v⁽⁰⁾ = 0 and $v^{(1)}=-{\omega }_0^2⟨x⟩$. A direct computation shows that the coefficients v^(j) have the following structure to an arbitrary order j ≥ 0:

$v^{\left(j\right)}\left(t\right)=-{\bar{\alpha }}_j⟨x\left(t\right)⟩,(47)$

with ${\bar{\alpha }}_{2j+1}\rangle 0$ and ${\bar{\alpha }}_{2j}=0$. Upon inserting the relation (Eq. 47) into the recurrence Eq. 45, the Chapman-Enskog method results in the nonlinear recurrence procedure for the coefficients ${\bar{\alpha }}_j$:

${\bar{\alpha }}_{j+1}=\sum _{k=0}^j{\bar{\alpha }}_k{\bar{\alpha }}_{j-k},j\ge 1,(48)$

with the initial conditions ${\bar{\alpha }}_0=0$ and ${\bar{\alpha }}_1={\omega }_0^2$. We observe that the term

${\alpha }_1=\frac{{\omega }_0^2}{\gamma },(49)$

recovers the drift coefficient in the Smoluchowski Eq. 41. The corresponding closure

$\langle v\left(t\right)\rangle =-{\alpha }_1\langle x\left(t\right)\rangle (50)$

does not allow to accurately reconstruct the behavior of the trajectories of the dynamics (Eq. 43) in presence of moderate damping effects, as visible in the left and central panels of Figure 1.

FIGURE 1

FIGURE 1. Behavior of the solutions of the ODE system (Eq. 43), for different initial data, with ω₀ = 1 and with γ_s = 0.04 (A), γ_s = 1.5 (B) and γ_s = 4 (C). The tiny solid lines correspond to individual trajectories, the thick solid lines denote the eigenvector u_M and the dashed lines represent the solution obtained with the closure given in Eq. 50.

We now aim at showing that the series

$\alpha =\sum _{j=0}^{\infty }{\bar{\alpha }}_j{\epsilon }^j=\sum _{j=0}^{\infty }{\alpha }_j(51)$

can be summed up in closed form: this procedure will single out an algebraic invariant manifold for the linear ODE system (Eq. 43). We start by multiplying both sides of (Eq. 48) by ɛ^j+1 and then sum in j from 1 to ∞. We obtain:

${\epsilon }^{-1}\left[\sum _{j=0}^{\infty }{\bar{\alpha }}_j{\epsilon }^j-{\bar{\alpha }}_0-{\bar{\alpha }}_1\epsilon \right]=\sum _{j=0}^{\infty }{\epsilon }^j\left(\sum _{k=0}^j{\bar{\alpha }}_k{\bar{\alpha }}_{j-k}\right)-{\bar{\alpha }}_0^2,(52)$

which, using (Eq. 51), yields the invariance equation:

${\alpha }^2-\gamma \alpha +{\omega }_0^2=0.(53)$

Note that Eq. 53 is readily established from Eq. 18 by setting $\chi =\alpha /(m{\omega }_0^2)$ and $M=\beta m{\omega }_0^2⟨x⟩$. The roots of the quadratic Eq. 53 coincide with the two real-valued eigenvalues λ_± of the drift matrix M. Note that the eigenvalue λ₊ diverges in the limit ɛ → 0. Hence, since we look for bounded solutions to the invariance equation, we set

$\alpha ={\lambda }_{-}.(54)$

We remark that the parameter α, in Eq. 54, corresponds to the exact summation of the Chapman-Enskog series (Eq. 51), and it thus yields the desired correction of the drift term (Eq. 49) of the Smoluchowski equation up to an arbitrary order of ɛ.

The same Eq. 53 can also be derived, in a non-perturbative fashion, via the principle of Dynamic Invariance. To this aim, we express the variable 〈v(t)〉 in terms of 〈x(t)〉 via the closure $\mathrm{\Phi }:\mathbb{R}\to \mathbb{R}$, which is endowed with the linear structure:

$⟨v\left(t\right)⟩=\mathrm{\Phi }\left[⟨x\left(t\right)⟩\right]=-\alpha ⟨x\left(t\right)⟩,(55)$

where the parameter α 〉 0 depends on γ and ω₀. The relation (Eq. 55) highlights a key aspect of the invariant manifold method: the variable 〈v(t)〉 depends on time only through the time dependence of the variable 〈x(t)〉. Next, upon inserting the closure (Eq. 55) in the ODE system (Eq. 43), one obtains the so-called “microscopic” time derivative of 〈v(t)〉:

${\partial }_t^{\left(micro\right)}\langle v\left(t\right)\rangle ≔\langle \dot{v}\left(t\right)\rangle =-{\omega }_0^2\langle x\left(t\right)\rangle +\gamma \alpha \langle x\left(t\right)\rangle .(56)$

We then define a projection operator ${\mathcal{P}}_x$, such that ${\mathcal{P}}_x⟨\dot{v}(t)⟩{|}_{⟨v(t)⟩=\mathrm{\Phi }[⟨x(t)⟩]}$ yields the evolution of the fast variable along the slow manifold parametrized by 〈x(t)〉. The action of ${\mathcal{P}}_x$ on $⟨\dot{v}(t)⟩$ is expressed, in this case, via the chain rule:

${\mathcal{P}}_x\langle \dot{v}\left(t\right)\rangle {|}_{\langle v\left(t\right)\rangle =\mathrm{\Phi }\left[\langle x\left(t\right)\rangle \right]}=D_{\langle x\left(t\right)\rangle }\mathrm{\Phi }\left[\langle x\left(t\right)\rangle \right]\langle \dot{\mathbf{x}}\left(t\right)\rangle .(57)$

In the sequel, to ease the notation, we shall denote the projected dynamics of 〈v(t) by ${\mathcal{P}}_x⟨\dot{v}(t)⟩$, where the closure (Eq. 55) is implicitly assumed. The “macroscopic” time derivative of 〈v(t)〉 is thus defined with the aid of the projection operator ${\mathcal{P}}_x$ as follows:

${\partial }_t^{\left(macro\right)}⟨v\left(t\right)⟩≔{\mathcal{P}}_x⟨\dot{v}\left(t\right)⟩={\alpha }^2⟨x\left(t\right)⟩.(58)$

The Dynamic Invariance principle, recall Eq. 8, states that the two “microscopic” and “macroscopic” time derivatives Eqs 56, 58 coincide, and the equality should hold independently of the value of the variable 〈x(t)〉: this leads again to the invariance Eq. 53. It is worth pointing out that by reintroducing the expansion (Eq. 51) in Eq. 53, one may reconstruct “backward” the recurrence relation (Eq. 48) with the corresponding initial conditions. This observation clarifies that the invariance Eq. 53 stands as the central result of the invariant manifold method, whereas the Chapman-Enskog expansion can be interpreted an iterative procedure for solving the invariance equation via a power series representation. Relying on approximate solutions is, in fact, the only feasible approach when the invariance equation can not be solved analytically. Alternative iterative methods (based e.g., on the Newton’s method), which may help circumvent some well-known instabilities appearing in low-order truncations of the Chapman-Enskog expansion, were considered in the framework of kinetic theory of gases [16, 17].

We conclude this Section by remarking that the invariant manifold method neglects, by construction, the fast relaxation dynamics in the initial layer, ruled by the eigenvalue λ₊, whereas it accurately captures the evolution along the slow manifold, encoded by the eigenvalue λ₋. A meaningful application of the method thus requires that the parameter γ_s be large enough to guarantee the existence of an appropriate separation of time scales [49]. This, in turn, allows to retain the Markovian approximation also in the contracted description, as commonly done in the context of the Mori-Zwanzig projection operator approach [46].

3.2.2 The Fluctuation-Dissipation Theorem

Let us characterize, then, the fluctuations in Eq. 42, by properly embedding the diffusion coefficient ${\mathcal{D}}_r$ in the framework of the Fluctuation-Dissipation theorem. By integrating Eq. 42 with a deterministic initial datum x(0) = x, one obtains:

$x\left(t\right)=e^{-\alpha t}x+{\int }_0^t{e}^{-\alpha \left(t-s\right)}dW\left(s\right).(59)$

The two-time correlation function of the position variable can be then calculated explicitly. It reads:

$\begin{array}{ll}\hfill ⟨x\left(s\right)x\left(t\right)⟩& =e^{-\alpha \left(t+s\right)}{x}^2+2{\mathcal{D}}_r{\int }_0^{\min \left(s,t\right)}{e}^{-\alpha \left(t+s-2\tau \right)}d\tau \hfill \\ \hfill & =\left(x^2-\frac{{\mathcal{D}}_r}{\alpha }\right)e^{-\alpha \left(t+s\right)}+\frac{{\mathcal{D}}_r}{\alpha }{e}^{-\alpha |t-s|}.\hfill \end{array}(60)$

We set s = t and require that the stationary value of ⟨x(t)²⟩ fulfills the Equipartition Theorem, namely:

$\underset{t\to \infty }{\lim }⟨x{\left(t\right)}^2⟩={\left(\beta m{\omega }_0^2\right)}^{-1}.(61)$

As a direct consequence, we obtain a relation establishing a connection between the exact drift coefficient α and the reduced diffusion coefficient ${\mathcal{D}}_r$:

$\alpha =\beta m{\omega }_0^2{\mathcal{D}}_r.(62)$

Eq. 62 is an instance of the Fluctuation-Dissipation theorem of the II kind [50] for the reduced dynamics (Eq. 42). In fact, since for the harmonically bound particle the mobility takes the form $\chi =\alpha /(m{\omega }_0^2)$, Eq. 62 establishes the Einstein relation $\chi =\beta {\mathcal{D}}_r$ connecting the mobility to the diffusion coefficient of the reduced dynamics.

We also note that using Eqs 49, 62, it is possible to relate ${\mathcal{D}}_r$ to the diffusion coefficient D of the original dynamics:

${\mathcal{D}}_r={\left({\alpha }_1\right)}^{-1}\alpha D,(63)$

which offers a multi-level characterization of the fluctuations in the Brownian oscillator model. One may expand D_r in a power series in ɛ, viz.:

${\mathcal{D}}_r=\sum _{j=0}^{\infty }{\bar{\mathcal{D}}}_j{\epsilon }^{j-1}=\sum _{j=0}^{\infty }{\mathcal{D}}_j.(64)$

Upon inserting (Eqs. 51, 64) into (Eq. 63), one obtains a hierarchy of equations relating, for each j ≥ 0, the coefficients ${\bar{\mathcal{D}}}_j$ to the coefficients ${\bar{\alpha }}_j$ in (Eq. 51):

${\bar{\mathcal{D}}}_j={\left({\bar{\alpha }}_1\right)}^{-1}{\bar{\alpha }}_{j}D,j\ge 0.(65)$

Thus, the leading-order term in Eq. 65 corresponds to

${\mathcal{D}}_1=D,(66)$

which recovers the diffusion coefficient in the Smoluchowski Eq.41.

3.3 Quantitative Control of the Reduction Error

We provide, here, some quantitative estimate of the error introduced by the application of the reduction method to the Brownian oscillator model. There are two relevant sources of error coming with the proposed scheme. The first source, which is somehow intrinsic in the procedure, traces back to the moment parameterization introduced in Section 2.1 and is connected to the existence of an invariant manifold parameterized by the values of the slow variable. A proper choice of the initial data allows one to control such first contribution. A second source, instead, is related to the defect of invariance, and keeps track of the approximation introduced in solving the invariance equation. To see this, we denote by 〈y(t)〉 = (〈y(t)〉, 〈w(t)〉) the solutions of the ODE system (Eq. 43) supplied with the closure (Eq. 55), viz.:

$\langle w\left(t\right)\rangle =-\alpha \langle y\left(t\right)\rangle ,(67)$

with deterministic initial datum y(0) = y. Our purpose here is to compare in a quantitative way 〈y(t)〉 with 〈x(t)⟩, the latter being the solution of the same ODE system when no closure is invoked. We shall finally give a pointwise in time a priori representation of the reduction error.

To be specific, we introduce the error terms e₁(t)≔〈x(t)〉 − 〈y(t)〉, e₂(t) ≔〈v(t)〉 − 〈w(t)〉, and $e_3(t)≔⟨\dot{v}(t)⟩-⟨\dot{w}(t)⟩$, as well as the defect of invariance:

${\mathrm{\Delta }}_y≔\left(1-P_y\right)⟨\dot{w}\left(t\right)⟩,(68)$

which, by virtue of the closure (Eq. 67), takes the form:

${\mathrm{\Delta }}_y=-\left({\omega }_0^2-\alpha \gamma +{\alpha }^2\right)⟨y\left(t\right)⟩.(69)$

Since Eq. 67 gives $⟨\dot{w}(t)⟩={\alpha }^2⟨y(t)⟩$, using Eq. 69, we arrive at:

$\begin{array}{ll}\hfill e_3\left(t\right)& =-{\omega }_0^2⟨x\left(t\right)⟩-\gamma ⟨v\left(t\right)⟩-{\alpha }^2⟨y\left(t\right)⟩\hfill \\ \hfill & =-{\omega }_0^2\left(⟨x\left(t\right)⟩-⟨y\left(t\right)⟩\right)-\gamma \left(⟨v\left(t\right)⟩-⟨w\left(t\right)⟩\right)+{\mathrm{\Delta }}_y,\hfill \end{array}(70)$

which can thus be rewritten as:

$e_3\left(t\right)=-{\omega }_0^2{e}_1\left(t\right)-\gamma e_2\left(t\right)+{\mathrm{\Delta }}_y.(71)$

Then, noticing that $e_3(t)={\dot{e}}_2(t)={\ddot{e}}_1(t)$, we may rewrite Eq. 71 as:

${\ddot{e}}_1\left(t\right)+\gamma {\dot{e}}_1\left(t\right)+{\omega }_0^2{e}_1\left(t\right)={\mathrm{\Delta }}_y,(72)$

in which Δ_y stands as the source term in the second-order linear ODE describing the dynamics of e₁(t).

A double integration over time of Eq. 72 yields:

$e_1\left(t\right)=e_1\left(0\right)+e_2\left(0\right)t-{\int }_0^{t}ds{\int }_0^{s}d\tau \left({\omega }_0^2{e}_1\left(\tau \right)+\gamma e_2\left(\tau \right)\right)+{\int }_0^{t}ds{\int }_0^{s}d\tau {\mathrm{\Delta }}_y,(73)$

which can be cast in the form:

$e_1\left(t\right)+{\omega }_0^2{\int }_0^{t}ds{\int }_0^{s}d\tau \left[e_1\left(\tau \right)\right]+\gamma {\int }_0^{t}ds{\int }_0^{s}d\tau \left[e_2\left(\tau \right)\right]=\mathcal{R}\left(t\right),(74)$

where the residual

$\mathcal{R}\left(t\right)≔e_1\left(0\right)+e_2\left(0\right)t+{\int }_0^{t}ds{\int }_0^{s}d\tau {\mathrm{\Delta }}_y(75)$

quantifies the quality of the reduction method. Note that the definition of the residual in Eq. 75 includes the defect of invariance Δ_y expressed by Eq. 69. If Δ_y vanishes, then controlling the error of the reduction method amounts to guessing the initial value y such that the first two terms on the r.h.s. of (Eq. 75) are small. The behavior of e₁(t) for different values of γ_s and for different choices of the closure (Eq. 67) is displayed in Figure 2. In each panel, solid, dotted and dashed lines correspond to the solution of the invariance equation Δ_y = 0 (see Eq. 54), the first-order term of Chapman-Enskog expansion (cf. Eq. 49) and the third-order approximation of the same expansion, respectively.

FIGURE 2

FIGURE 2. Solutions of the ODE (Eq. 72), with x = 2, v = −1, y = 1 and with different choices of the closure (Eq. 67). Shown are the reduction errors corresponding to the first-order term of the Chapman-Enskog expansion (dotted curve), the corresponding third-order approximation of the same expansion (dashed curve) and the solution of the invariance equation (solid curve). The values of ω₀ and γ_s in the various panels correspond to those considered in Figure 1.

3.4 Response and Correlation Functions

We now turn to the study of correlation functions, which constitute a useful test-bed to assess the range of applicability of the proposed reduced description of the Brownian oscillator model. Following the basic tenets of linear response theory, correlation functions are connected to the response of the system to an external stimulus; we refer the reader to Ref. [51] for an exhaustive review on this subject and also to the concise theoretical guidelines provided in Supplementary Appendix S2. An extension of linear response theory to chaotic non-Hamiltonian systems has also been discussed in [52]. We assume that the system described by Eq. 42 is initially in equilibrium with a heat bath at inverse temperature β. The stationary distribution of the reduced dynamics (Eq. 42) takes the form:

${\rho }_0\left(x\right)=\sqrt{\frac{\beta m{\omega }_0^2}{2\pi }}\exp \left\{-\frac{1}{2}\beta m{\omega }_0^2{x}^2\right\}.(76)$

We then probe the dynamics (Eq. 42) by adding on the right hand side, at time t = 0, a small, purely time-dependent, perturbation F(t).

The perturbation induces the following structure of the Fokker-Planck equation:

$\begin{array}{ccc}\hfill \frac{\partial \rho \left(x,t\right)}{\partial t}& =\hfill & \left({\mathcal{L}}^{*}+{\mathcal{L}}_{ext}^{*}\right)\rho \left(x,t\right)\hfill \\ \hfill \rho \left(x,0\right)& =\hfill & {\rho }_0\left(x\right),\hfill \end{array}(77)$

where the operator ${\mathcal{L}}^{*}$ and ${\mathcal{L}}_{ext}^{*}$ acts on probability densities as follows:

${\mathcal{L}}^{*}\rho \left(x,t\right)=\left(\alpha \frac{\partial }{\partial x}x+{\mathcal{D}}_r\frac{{\partial }^2}{\partial x^2}\right)\rho \left(x,t\right)(78)$

${\mathcal{L}}_{ext}^{*}\rho \left(x,t\right)=-F\left(t\right)\frac{\partial }{\partial x}\rho \left(x,t\right).(79)$

To write the response formula, we introduce the observable

$B\left(x\right)=-{\rho }_0{\left(x\right)}^{-1}{\partial }_x{\rho }_0\left(x\right),$

which takes here the form:

$B\left(x\right)=\beta m{\omega }_0^{2}x.(80)$

We then look at the response of the observable A(x) = x. To this aim, we shall denote by 〈…〉₀ the complete average taken with initial density ρ₀. According to the basic guidelines of linear response theory, the response function R_x,x(t), for t 〉 0, attains the structure:

$\begin{array}{ll}\hfill R_{x,x}\left(t\right)& =\int A\left(x\left(t\right)\right)B\left(x\right){\rho }_0\left(x\right)dx\hfill \\ \hfill & =C{\left(0\right)}^{-1}C\left(t\right),\hfill \end{array}(81)$

where C(t) = 〈x(t)x〉₀ is the autocorrelation function of the position variable, and

$C\left(0\right)={⟨x^2⟩}_0=\frac{{\mathcal{D}}_r}{\alpha }={\left(\beta m{\omega }_0^2\right)}^{-1}.(82)$

We also note in passing that the drift coefficient α is connected to the autocorrelation time τ_c, defined as:

${\tau }_c=C{\left(0\right)}^{-1}{\int }_0^{\infty }C\left(t\right)dt={\alpha }^{-1}.(83)$

Starting from (Eq. 42), C(t) is found to obey for any t 〉 0 the equation:

$\frac{d}{dt}C\left(t\right)+\alpha C\left(t\right)=0,(84)$

with the initial condition fixed by Eq. 82.

The connection between correlation and response functions can be further unveiled as follows. By the Wiener-Khinchin Theorem [8], the spectral density S(ω) of a stationary random process x(t) is equal to the Fourier transform of its autocorrelation function, i.e.:

$S\left(\omega \right)=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }{e}^{-i\omega t}C\left(t\right)dt.(85)$

The dynamics described by Eq. 42, admits a dynamical mobility [50] (or generalized susceptibility) μ(ω) of the form:

$\mu \left(\omega \right)=\frac{1}{\alpha +i\omega }.(86)$

By multiplying both sides of Eq. 84 by the factor e^−iωt, and by integrating over time from 0 to + ∞, an integration by parts gives:

$\left(\alpha +i\omega \right){\int }_0^{\infty }{e}^{-i\omega t}C\left(t\right)dt-C\left(0\right)=0,(87)$

which, using Eq. 86, leads to the following remarkable expression of the mobility:

$\mu \left(\omega \right)=C{\left(0\right)}^{-1}{\int }_0^{+\infty }{e}^{-i\omega t}C\left(t\right)dt.(88)$

Owing to the fact that C(t) = C( − t) is a real, symmetric function of time, we use the relation (Eq. 88) to reshape Eq. 85 in the form:

$S\left(\omega \right)=\frac{C\left(0\right)}{\pi }\mathfrak{R}\left[\mu \left(\omega \right)\right],(89)$

where $\mathfrak{R}[\mu (\omega )]$ denotes the real part of the mobility μ(ω). The relation (Eq. 89) is a classical version of the Fluctuation-Dissipation theorem of the I kind [50], as it connects the response to an external stimulus, represented by the dynamical mobility, to the fluctuations spontaneously produced in the system described by Eq. 42, encoded by the spectral density. By now putting the explicit expressions (Eqs 82, 86) in Eq. 89, we recover the standard Lorentzian form of the spectral density of the reduced dynamics (Eq. 42):

$S\left(\omega \right)=\frac{1}{\pi }\frac{{\mathcal{D}}_r}{{\omega }^2+{\alpha }^2}.(90)$

One may analogously repeat the foregoing derivation for the original dynamics of the Brownian oscillator, described by Eq. 34, which constitutes an exactly solvable example [8, 48]. The invariant density for the unperturbed dynamics has the explicit representation:

${\rho }_0\left(x,v\right)=\frac{m{\omega }_0\beta }{2\pi }\exp \left\{-\frac{m\beta }{2}\left(v^2+{\omega }_0^2{x}^2\right)\right\}.(91)$

We probe Eq. 34 by adding a time-dependent term F(t) in the dynamics of the position variable and check the response in the variable x itself, as above. The perturbed operator ${\mathcal{L}}_{ext}$ in the Fokker-Planck equation reads:

${\mathcal{L}}_{ext}\rho \left(x,v,t\right)=-F\left(t\right)\frac{\partial }{\partial x}\rho \left(x,v,t\right),(92)$

and it holds:

$A\left(x\right)=x,B\left(x\right)=\beta m{\omega }_0^{2}x.(93)$

Using the explicit expression of the element G_xx of the Green’s matrix, see Supplementary Appendix S1, we end up with the following response formula for the original dynamics:

${\mathcal{R}}_{x,x}\left(t\right)=\frac{{\lambda }_{+}{e}^{-{\lambda }_{-}t}-{\lambda }_{-}{e}^{-{\lambda }_{+}t}}{{\gamma }_s},(94)$

which inherits the contributions of both the “fast” and the “slow” time scales of the system, encoded by the eigenvalues λ_± of the matrix M in Eq. 35. For any finite time t 〉 0, it makes thus sense to compare the response formulae computed with both the reduced and the original dynamics, Eqs 81, 94, respectively. Thus, from Eqs 81, 84 one finds the following structure for the response function of the reduced dynamics:

$R_{x,x}\left(t\right)=e^{-\alpha t},t>0.(95)$

Recalling (Eq. 54), one thus finds:

$\left|{\mathcal{R}}_{x,x}\left(t\right)-R_{x,x}\left(t\right)\right|\le \frac{{\lambda }_{-}}{{\gamma }_s}.(96)$

Using the explicit dependence of λ₋ and γ_s on γ, one finds that λ₋/γ_s → 0 as γ → +∞, which implies the uniform-in-time convergence of R_x,x(t) to ${\mathcal{R}}_{x,x}(t)$.

We emphasize that the upper bound (Eq. 96) guarantees that the response function of the reduced dynamics converges to the response function of the original dynamics in the high friction limit, namely when the time scale separation, controlled by γ_s, grows. In this context, it is instructive to observe in Figure 1 the first plot in comparison with the third plot.

4 Conclusions

In this work we considered a classical problem of statistical mechanics, concerning the extraction of a reduced description for Brownian dynamics in a confining potential. Adiabatic elimination techniques were already introduced in the 1970s to derive the Smoluchowski equation from the Kramers equation in the high friction limit. Application of the Chapman-Enskog scheme to the Fokker-Planck equation paved the way to a systematic derivation of the Smoluchowski formula and its higher order corrections via an expansion in power of the inverse friction coefficient. The same procedure is traditionally exploited, in kinetic theory of gases, to obtain the Euler and the Navier-Stokes equations from the Boltzmann equation. Nevertheless, the method is also known to suffer from the onset of short wavelength instabilities, which violate the H-theorem. As evidenced by the study of different kinetic models, the failure of the Chapman-Enskog expansion does not lie in the scheme itself, but in its truncation to lower order levels. The invariant manifold method is a non-perturbative reduction technique leading to hydrodynamic equations that are, instead, stable at all wavelengths. The method stipulates a condition of dynamic invariance which, for the class of kinetic models known as linearized Grad’s moment systems, was shown to yield the same result as the exact summation of the Chapman-Enskog expansion. The aim of the present work is, hence, to outline the use of the invariant manifold set-up to the reduction of a class of Brownian dynamics.

The main results can be summarized as follows:

(1) By exploiting the Dynamic Invariance principle, we derived an equation of invariance whose solutions generalize the structure of the mobility coefficient beyond the overdamped limit, thus making the reduced description prone to real world applications. While analytical solutions to the invariance equation are available in the presence of power law potentials, in the simplest case of the Brownian oscillator model we also succeed to sum up exactly the Chapman-Enskog expansion.

(2) We obtained a quantitative estimate of the error encoded in the reduction procedure, which can be controlled through a suitable choice of the initial data and by minimizing the defect of invariance.

(3) We used linear response theory to shed light on the response functions due to the original and the reduced dynamics. We proved the convergence of the two response functions in the high friction limit.

We believe that the procedure outlined in this work can provide useful insights on a wider class of Brownian dynamics. Even when exact solutions of the invariance equation are not available, the combined use of an iterative scheme of solution of the equation and the Fluctuation-Dissipation relation may help unravel meaningful reduced descriptions.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

MC thanks L. Rondoni (Turin Polytechnic, Italy) for many useful discussions. AM thanks H. Duong (Birmingham, United Kingdom) for his KAAS seminar on a related topic as well as for his constructive ideas related to Section 3.3.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2022.903030/full#supplementary-material

References

1. Uhlenbeck GE, Ornstein LS. On the Theory of the Brownian Motion. Phys. Rev. (1930) 36:823–41. doi:10.1103/physrev.36.823

CrossRef Full Text | Google Scholar

2. van Kampen N. Elimination of Fast Variables. Phys Rep (1985) 124:69–160. doi:10.1016/0370-1573(85)90002-x

CrossRef Full Text | Google Scholar

3. Givon D, Kupferman R, Stuart A. Extracting Macroscopic Dynamics: Model Problems and Algorithms. Nonlinearity (2004) 17:R55–R127. doi:10.1088/0951-7715/17/6/r01

CrossRef Full Text | Google Scholar

4. Frédéric L, Lelièvre T, Sharma U. Effective Dynamics for Non-reversible Stochastic Differential Equations: a Quantitative Study. Nonlinearity (2019) 32:4779–816.

Google Scholar

5. Wilemski G. On the Derivation of Smoluchowski Equations with Corrections in the Classical Theory of Brownian Motion. J. Stat. Phys. (1976) 14:153–69. doi:10.1007/bf01011764

CrossRef Full Text | Google Scholar

6. Chaturvedi S, Shibata F. Time-convolutionless Projection Operator Formalism for Elimination of Fast Variables. Applications to Brownian Motion. Z. Phys B (1979) 35:297–308. doi:10.1007/bf01319852

CrossRef Full Text | Google Scholar

7. Miguel MS, Sancho JM. A Colored-Noise Approach to Brownian Motion in Position Space. Corrections to the Smoluchowski Equation. J. Stat. Phys. (1980) 22:605–24. doi:10.1007/bf01011341

CrossRef Full Text | Google Scholar

8. Risken H. The Fokker-Planck Equation. Berlin, Germany: Springer-Verlag (1996).

Google Scholar

9. Bocquet L. High Friction Limit of the Kramers Equation: The Multiple Time-Scale Approach. Am J Phys (1997) 65:140–4. doi:10.1119/1.18507

CrossRef Full Text | Google Scholar

10. Cerrai S, Freidlin M. On the Smoluchowski-Kramers Approximation for a System with an Infinite Number of Degrees of Freedom. Probab. Theory Relat. Fields (2006) 135:363–94. doi:10.1007/s00440-005-0465-0

CrossRef Full Text | Google Scholar

11. Goddard BD, Pavliotis GA, Kalliadasis S. The Overdamped Limit of Dynamic Density Functional Theory: Rigorous Results. Multiscale Model. Simul. (2012) 10:633–63. doi:10.1137/110844659

CrossRef Full Text | Google Scholar

12. Titulaer UM. A Systematic Solution Procedure for the Fokker-Planck Equation of a Brownian Particle in the High-Friction Case. Phys A Stat Mech its Appl (1978) 91:321–44. doi:10.1016/0378-4371(78)90182-6

CrossRef Full Text | Google Scholar

13. Theiss W, Titulaer UM. The Systematic Adiabatic Elimination of Fast Variables from a Many-Dimensional Fokker-Planck Equation. Phys A Stat Mech its Appl (1985) 130:123–42. doi:10.1016/0378-4371(85)90100-1

CrossRef Full Text | Google Scholar

14. Chavanis P-H, Laurençot P, Lemou M. Chapman-Enskog Derivation of the Generalized Smoluchowski Equation. Phys A Stat Mech its Appl (2004) 341:145–64. doi:10.1016/j.physa.2004.04.102

CrossRef Full Text | Google Scholar

15. Bobylev AV. Instabilities in the Chapman-Enskog Expansion and Hyperbolic Burnett Equations. J Stat Phys (2006) 124:371–99. doi:10.1007/s10955-005-8087-6

CrossRef Full Text | Google Scholar

16. Colangeli M, Karlin IV, Kröger M. From Hyperbolic Regularization to Exact Hydrodynamics for Linearized Grad's Equations. Phys Rev E Stat Nonlin Soft Matter Phys (2007) 75:051204. doi:10.1103/PhysRevE.75.051204

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Colangeli M, Karlin IV, Kröger M. Hyperbolicity of Exact Hydrodynamics for Three-Dimensional Linearized Grad's Equations. Phys Rev E Stat Nonlin Soft Matter Phys (2007) 76:022201. doi:10.1103/PhysRevE.76.022201

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Karlin IV, Colangeli M, Kröger M. Exact Linear Hydrodynamics from the Boltzmann Equation. Phys. Rev. Lett. (2008) 100:214503. doi:10.1103/physrevlett.100.214503

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Colangeli M, Kröger M, Ottinger HC. Boltzmann Equation and Hydrodynamic Fluctuations. Phys Rev E Stat Nonlin Soft Matter Phys (2009) 80:051202. doi:10.1103/PhysRevE.80.051202

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Gorban AN, Karlin IV. Invariant Manifolds for Physical and Chemical Kinetics, Vol. 660 of Lect. Notes Phys.. Berlin, Germany: Springer-Verlag (2005).

Google Scholar

21. Karlin IV. Exact Summation of the Chapman-Enskog Expansion from Moment Equations. J. Phys. A. Math. Gen. (2000) 33:837–46. doi:10.1088/0305-4470/33/45/303

CrossRef Full Text | Google Scholar

22. Karlin IV, Gorban AN. Hydrodynamics from Grad's Equations: What Can We Learn from Exact Solutions? Ann. Phys. (2002) 11:783–833. doi:10.1002/1521-3889(200211)11:10/11〈783::aid-andp783〉3.0.co;2-v

CrossRef Full Text | Google Scholar

23. Bensoussan A, Lions J-L, Papanicolaou G. Asymptotic Analysis for Periodic Structures. Amsterdam: North-Holland (1978).

Google Scholar

24. Duong MH, Lamacz A, Peletier MA, Schlichting A, Sharma U. Quantification of Coarse-Graining Error in Langevin and Overdamped Langevin Dynamics. Nonlinearity (2018) 31:4517–66. doi:10.1088/1361-6544/aaced5

CrossRef Full Text | Google Scholar

25. Lucarini V, Pavliotis GA, Zagli N. Response Theory and Phase Transitions for the Thermodynamic Limit of Interacting Identical Systems. Proc. R. Soc. A (2020) 476:20200688. doi:10.1098/rspa.2020.0688

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Ohtsuka H. On the Derivation of the Mean Field Equation of the Gibbs Distribution Function for Equilibrium Vortices in an External Field. RIMS Kokyuroku Bessatsu (2020) B82:67–85.

Google Scholar

27. Pavliotis GA, Stuart AM. Multiscale Methods: Averaging and Homogenization. Berlin, Germany: Springer-Verlag (2008).

Google Scholar

28. Colangeli M, Muntean A. Towards a Quantitative Reduction of the SIR Epidemiological Model. In: Libelli G,, and Bellomo N, editors. Crowd Dynamics, Vol. 3: Theory, Models and Safety Problems. Modeling and Simulation in Science, Engineering and Technology. Boston: Birkhäuser (2021).

Google Scholar

29. MacKay RS, Robinson JD. Aggregation of Markov Flows I: Theory. Phil. Trans. R. Soc. A (2018) 376:20170232. doi:10.1098/rsta.2017.0232

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Mohammed S-EA, Scheutzow MKR. The Stable Manifold Theorem for Stochastic Differential Equations. Ann Probab (1999) 27:615–52. doi:10.1214/aop/1022677380

CrossRef Full Text | Google Scholar

31. Roberts AJ. Resolving the Multitude of Microscale Interactions Accurately Models Stochastic Partial Differential Equations. LMS J. Comput. Math. (2006) 9:193–221. doi:10.1112/s146115700000125x

CrossRef Full Text | Google Scholar

32. Ruelle D. Characteristic Exponents and Invariant Manifolds in Hilbert Space. Ann Math (1982) 115:243–90. doi:10.2307/1971392

CrossRef Full Text | Google Scholar

33. Gritsun A, Lucarini V. Fluctuations, Response, and Resonances in a Simple Atmospheric Model. Phys D Nonlinear Phenom (2017) 349:62–76. doi:10.1016/j.physd.2017.02.015

CrossRef Full Text | Google Scholar

34. Majda AJ, Franzke C, Crommelin D. Normal Forms for Reduced Stochastic Climate Models. Proc. Natl. Acad. Sci. U.S.A. (2009) 106:3649–53. doi:10.1073/pnas.0900173106

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Arnold L. Random Dynamical Systems. Berlin, Germany: Springer-Verlag Berlin Heidelberg (1998).

Google Scholar

36. Kolmogorov AN. On Conservation of Conditionally Periodic Motions under Small Perturbations of the Hamiltonian. Dokl. Akad. Nauk SSSR (1954) 98:527–30.

Google Scholar

37. Arnold VI. Proof of a Theorem of A.N. Kolmogorov on the Invariance of Quasi-Periodic Motions under Small Perturbations of the Hamiltonian. Russ Math Surv (1963) 18:9–36.

Google Scholar

38. Moser JR. Convergent Series Expansions for Quasi-Periodic Motions. Math. Ann. (1967) 169:136–76. doi:10.1007/bf01399536

CrossRef Full Text | Google Scholar

39. Gorban AN, Karlin I. Hilbert's 6th Problem: Exact and Approximate Hydrodynamic Manifolds for Kinetic Equations. Bull. Amer. Math. Soc. (2013) 51:187–246. doi:10.1090/s0273-0979-2013-01439-3

CrossRef Full Text | Google Scholar

40. Cercignani C. The Boltzmann Equation and its Applications. New York: Springer-Verlag (1988).

Google Scholar

41. Chapman S, Cowling TG. The Mathematical Theory of Nonuniform Gases. New York: Cambridge University Press (1970).

Google Scholar

42. Auletta G, Rondoni L, Vulpiani A. On the Relevance of the Maximum Entropy Principle in Non-equilibrium Statistical Mechanics. Eur. Phys. J. Spec. Top. (2017) 226:2327–43. doi:10.1140/epjst/e2017-70064-x

CrossRef Full Text | Google Scholar

43. Karlin IV. Invariance Principle and Model Reduction for the Fokker-Planck Equation. Phil. Trans. R. Soc. A (2016) 374:20160142. doi:10.1098/rsta.2016.0142

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Chandrasekhar S. Stochastic Problems in Physics and Astronomy. Rev. Mod. Phys. (1943) 15:1–89. doi:10.1103/revmodphys.15.1

CrossRef Full Text | Google Scholar

45. Cherstvy AG, Thapa S, Mardoukhi Y, Chechkin AV, Metzler R. Time Averages and Their Statistical Variation for the Ornstein-Uhlenbeck Process: Role of Initial Particle Distributions and Relaxation to Stationarity. Phys. Rev. E (2018) 98:022134. doi:10.1103/PhysRevE.98.022134

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Zwanzig R. Nonequilibrium Statistical Mechanics. Oxford, UK: Oxford University Press (2001).

Google Scholar

47. Łuczka J. Non-markovian Stochastic Processes: Colored Noise. Chaos (2005) 15:026107.

Google Scholar

48. Pavliotis GA. Stochastic Processes and Applications. Diffusion Processes, the Fokker-Planck and Langevin Equations. Berlin, Germany: Springer-Verlag (2014).

Google Scholar

49. Roberts AJ. Normal Form Transforms Separate Slow and Fast Modes in Stochastic Dynamical Systems. Phys A Stat Mech its Appl (2008) 387:12–38. doi:10.1016/j.physa.2007.08.023

CrossRef Full Text | Google Scholar

50. Kubo R, Toda M, Hashitsume N. Statistical Physics II. Nonequilibrium Statistical Mechanics. Berlin, Germany: Springer-Verlag (1985).

Google Scholar

51. Marconi U, Puglisi A, Rondoni L, Vulpiani A. Fluctuation-Dissipation: Response Theory in Statistical Physics. Phys Rep (2008) 461:111–95. doi:10.1016/j.physrep.2008.02.002

CrossRef Full Text | Google Scholar

52. Colangeli M, Rondoni L, Vulpiani A. Fluctuation-Dissipation Relation for Chaotic Non-Hamiltonian Systems. J Stat Mech (2012) 2012 4:L04002. doi:10.1088/1742-5468/2012/04/L04002

CrossRef Full Text | Google Scholar