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 [57] 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 [911]. 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 [1619] 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 [2328] for recent applications of reduction schemes to statistical physics and epidemiological models. We also refer the reader to [2932] 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) = xn, n ≥ 2 an integer, described by a system of stochastic differential equations (SDEs) written in the Itô form:

dxt=vtdtdvt=1mxVxdtγvtdt+2Dγ2dWt,(1)

where W(t) is a one-dimensional Wiener process, m is the mass, γ is the friction constant, D = (βmγ)−1 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:

dxt=1mγxVxdt+2DdWt.(2)

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

dxt=χxVxdt+2DrdWt,(3)

where χ is the mobility and

Dr

the diffusion coefficient of the reduced dynamics. The two coefficients χ and

Dr

are expected to reduce to ()−1 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 [3638]. 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 t0, it will remain on the manifold for all times tt0. 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:

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

tmacrofMDMfMmJfM,(6)

where the differential DMfM, in Eq. 6, is evaluated at the point M = m(fM). While J(fM) in Eq. 5, corresponds to a value of the vector field J evaluated on the manifold fM, Eq. 6 codifies a chain rule: one computes the time derivative of the moment M via the map m, as

Ṁ=m(J(fM))

, 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:

tmicrofM=tmacrofM,(7)

or, alternatively, as:

ΔMtmicrofMtmacrofM=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 fM 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:

We also recall that if a functional E(f) is conserved for the dynamics (Eq. 4), then E(fM) 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(fM) 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:

with ɛ 〉 0 a small parameter. One requires that m(Q(f)) = 0 and that, for each Mm(U), the equation Q(f) = 0, with m(f) = M, has a unique solution denoted by

fMeq

(corresponding to the Maxwellian distributions, in Boltzmann’s theory), which is also asymptotically stable and globally attracting for the fast dynamics

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

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

fM=fMeq+i=1εifMi.(13)

In the set-up of Boltzmann’s theory, the zero-order approximation

fMfMeq

leads, upon integrating the Boltzmann equation over the velocity space, to the inviscid Euler equations of fluid dynamics, whereas the first-order correction

fMfM(eq)+εfM(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:

ρx,tt=Drxρx,txUx+xρx,t,(14)

where we enforced the Einstein relation

χ=βDr

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:

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

t(micro)v(t)

and

t(macro)v(t)

for the dynamics obtained upon averaging Eq. 1 over noise. We thus find:

tmicrovt1mMt+γχMt(16)
tmacrovtχṀt(17)

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

1mMt+γχMt+χṀt=0,(18)

to be solved for the mobility χ. Note that the classical form 1/() 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

Ṁ(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ρ=ρx,tlnρx,tρ0xdx,(19)

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

mixρx,tdx=Mit,i=0,,k,(20)

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

ρ*x,M=ρ0xexpi=0kΛimix(21)

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

Π*=i=0kρ*Mimixdx(22)

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

Ṁi=j=0kΛjxmixρ*x,txmjxdx.(23)

A step forward can be made by splitting the set of moments as M = MIMII, with MI = {M0, , M} and MII = {M+1, , Mk}. We assume that the (first) quasi-equilibrium distribution ρ*(x) can be derived explicitly for the set MI, i.e., ρ* = ρ*(x, MI), 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:

ΔSφ=ρ*φ1+lnρ*ρ0dx12ρ*φ2dx.(24)

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

ΔSmax,mIxρ*φdx=0,mIIxρ*φdx=ΔMII,(25)

where ΔMIIMIIMII(MI) represents the deviations of the moments MII 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 MI = M0 (i.e., ρ*(x) = ρ0(x)) and MII = M, with ∫m(x)ρ dx = M. We thus find:

φt=Lφ,,L=ρ01Drxρ0x.(26)

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

φ0x,t=m0xMtM0,(27)

with

m0x=mxm0m20m02.(28)

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

Π0=ρ0xm0xm0m00m0xdx.(29)

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

Ṁ=0MtM0,(30)

with

0=xm0Drxm00m0m00,(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]:

with j = 0, 1, …, where

m(j+1)=m(j+μ(j+1)

, and in which the orthogonality condition

μ(j+1)m(j0=0

is exploited. At the k-th iteration step, a refinement of the inverse relaxation time (Eq. 31) attains the form:

k=xmkDrxmk0mkmk0.(33)

The sequence

(k)kN

is found to converge to the eigenvalue of the operator

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ω02x2

, 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:

dxt=vtdtdvt=ω02xtdtγvtdt+2Dγ2dWt,(34)

where

ω0=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:

whose eigenvalues read:

with

γs=γ24ω02=λ+λ

. 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 γ−1, 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(ɛ−1t), with ɛ = γ−1, the original dynamics can be rescaled as follows:

dvεt=ε1ω02xεtdtε2vεtdt+ε12βm1dWt,(38)

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

ε1vεtdt=ω02xεtdt+2βm1dWt+Oε,(39)

and hence,

dxεt=ω02xεtdt+2βm1dWt+Oε.(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:

dxt=ω02γxtdt+2DdWt.(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:

dxt=αxtdt+2DrdWt,(42)

where α and

Dr

denote the drift and diffusion coefficients of the reduced dynamics, respectively. These coefficients will properly recover the two corresponding expressions

ω02/γ

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 ɛ = γ−1, 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

Dr

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,

ẋt=Mxt.(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 ɛ = γ−1, viz.

vt=j=0εjvjt.(44)

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

vj+1=k=0jDCEkvjk,j1,(45)

where the Chapman-Enskog operators

DCE(k)

act on the coefficients v(j) as follows:

DCEkvjvjxvk.(46)

The recurrence Eq. 45 starts with v(0) = 0 and

v(1)=ω02x

. A direct computation shows that the coefficients v(j) have the following structure to an arbitrary order j ≥ 0:

vjt=ᾱjxt,(47)

with

ᾱ2j+10

and

ᾱ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

ᾱj

:

ᾱj+1=k=0jᾱkᾱjk,j1,(48)

with the initial conditions

ᾱ0=0

and

ᾱ1=ω02

. We observe that the term

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

vt=α1xt(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.

www.frontiersin.org

FIGURE 1. Behavior of the solutions of the ODE system (Eq. 43), for different initial data, with ω0 = 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 uM and the dashed lines represent the solution obtained with the closure given in Eq. 50.

We now aim at showing that the series

α=j=0ᾱjεj=j=0α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:

ε1j=0ᾱjεjᾱ0ᾱ1ε=j=0εjk=0jᾱkᾱjkᾱ02,(52)

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

Note that Eq. 53 is readily established from Eq. 18 by setting

χ=α/(mω02)

and

M=βmω02x

. 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

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

Φ:RR

, which is endowed with the linear structure:

vt=Φxt=αxt,(55)

where the parameter α 〉 0 depends on γ and ω0. 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)〉:

tmicrovtv̇t=ω02xt+γαxt.(56)

We then define a projection operator

Px

, such that

Pxv̇(t)|v(t)=Φ[x(t)]

yields the evolution of the fast variable along the slow manifold parametrized by 〈x(t)〉. The action of

Px

on

v̇(t)

is expressed, in this case, via the chain rule:

Pxv̇t|vt=Φxt=DxtΦxtẋt.(57)

In the sequel, to ease the notation, we shall denote the projected dynamics of 〈v(t) by

Pxv̇(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

Px

as follows:

tmacrovtPxv̇t=α2xt.(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

Dr

in the framework of the Fluctuation-Dissipation theorem. By integrating Eq. 42 with a deterministic initial datum x(0) = x, one obtains:

xt=eαtx+0teαtsdWs.(59)

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

xsxt=eαt+sx2+2Dr0mins,teαt+s2τdτ=x2Drαeαt+s+Drαeα|ts|.(60)

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

limtxt2=βmω021.(61)

As a direct consequence, we obtain a relation establishing a connection between the exact drift coefficient α and the reduced diffusion coefficient

Dr

:

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

χ=α/(mω02)

, Eq. 62 establishes the Einstein relation

χ=βDr

connecting the mobility to the diffusion coefficient of the reduced dynamics.

We also note that using Eqs 49, 62, it is possible to relate

Dr

to the diffusion coefficient D of the original dynamics:

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

Dr=j=0D̄jεj1=j=0Dj.(64)

Upon inserting (Eqs. 51, 64) into (Eq. 63), one obtains a hierarchy of equations relating, for each j ≥ 0, the coefficients

D̄j

to the coefficients

ᾱj

in (Eq. 51):

D̄j=ᾱ11ᾱjD,j0.(65)

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

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.:

wt=αyt,(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 e1(t)≔〈x(t)〉 − 〈y(t)〉, e2(t) ≔〈v(t)〉 − 〈w(t)〉, and

e3(t)v̇(t)ẇ(t)

, as well as the defect of invariance:

Δy1Pyẇt,(68)

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

Δy=ω02αγ+α2yt.(69)

Since Eq. 67 gives

ẇ(t)=α2y(t)

, using Eq. 69, we arrive at:

e3t=ω02xtγvtα2yt=ω02xtytγvtwt+Δy,(70)

which can thus be rewritten as:

e3t=ω02e1tγe2t+Δy.(71)

Then, noticing that

e3(t)=ė2(t)=ë1(t)

, we may rewrite Eq. 71 as:

ë1t+γė1t+ω02e1t=Δy,(72)

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

A double integration over time of Eq. 72 yields:

e1t=e10+e20t0tds0sdτω02e1τ+γe2τ+0tds0sdτΔy,(73)

which can be cast in the form:

e1t+ω020tds0sdτe1τ+γ0tds0sdτe2τ=Rt,(74)

where the residual

Rte10+e20t+0tds0sdτΔ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 e1(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.

www.frontiersin.org

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 ω0 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:

ρ0x=βmω022πexp12βmω02x2.(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:

ρx,tt=L*+Lext*ρx,tρx,0=ρ0x,(77)

where the operator

L*

and

Lext*

acts on probability densities as follows:

L*ρx,t=αxx+Dr2x2ρx,t(78)
Lext*ρx,t=Ftxρx,t.(79)

To write the response formula, we introduce the observable

which takes here the form:

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

Rx,xt=AxtBxρ0xdx=C01Ct,(81)

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

C0=x20=Drα=βmω021.(82)

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

τc=C010Ctdt=α1.(83)

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

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ω=12π+eiωtCtdt.(85)

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

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

α+iω0eiωtCtdtC0=0,(87)

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

μω=C010+eiωtCtdt.(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:

where

R[μ(ω)]

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):

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:

ρ0x,v=mω0β2πexpmβ2v2+ω02x2.(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

Lext

in the Fokker-Planck equation reads:

Lextρx,v,t=Ftxρx,v,t,(92)

and it holds:

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

Rx,xt=λ+eλtλeλ+tγ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:

Recalling (Eq. 54), one thus finds:

Rx,xtRx,xtλγs.(96)

Using the explicit dependence of λ and γs on γ, one finds that λ/γs → 0 as γ → +, which implies the uniform-in-time convergence of Rx,x(t) to

Rx,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


Disclaimer:

This article is autogenerated using RSS feeds and has not been created or edited by OA JF.

Click here for Source link (https://www.frontiersin.org/)