Langevin dynamics with space-time periodic nonequilibrium forcing

Introduction

Nonequilibrium transport has attracted a lot of attention in recent years both in the mathematical physics and physics literature. It is by now well understood that simple low dimensional systems that are driven away from equilibrium can exhibit quite complicated behavior such as stochastic resonance , directed transport , absolute negative mobility and giant enhancement of diffusion . In particular, simple stochastic differential equations (SDEs) with space-time periodic coefficients have been proposed in recent years as models for Brownian motors. A natural question is whether SDEs with space-time mean zero coefficients can give rise to a nonzero effective drift. This question was studied in detail in . The goal of this paper is to refine and extend the results obtained in this paper.

The long time dynamics of SDEs with space-time periodic coefficients is characterized by an effective drift, whereas fluctuations around the macroscopic directed transport are described by an effective Brownian motion with covariance matrix $D$. The rigorous mathematical analysis of such models is based on proving the law of large numbers (ergodic theorem) at the hyperbolic timescale and a functional central limit theorem at the diffusive timescale. These problems are closely related to the theory of homogenization for parabolic PDEs (and of the corresponding stochastic differential equations) with space-time dependent coefficients [1, Chapter 3]. Homogenization problems for Brownian motion in a space-time periodic potential were also studied in .

On the other hand, nonequilibrium perturbations of systems at equilibrium can be used for calculating transport coefficients using linear response theory and the Green-Kubo formalism . Linear response theory can be rigorously analyzed for certain stochastic systems , and higher order corrections can also be obtained . The rigorous analysis of linear response theory and of the Green-Kubo formalism is very closely related to homogenization theory, in particular the systematic use of Poisson equations for obtaining formulas for transport coefficients. As an example we mention where the linear response theory/Poisson equation formalism for the Langevin dyanamics is used to calculate the shear viscosity coefficient.

Stochastic dynamics in a bistable potential under the influence of an external time periodic forcing can exhibit stochastic resonance, whereby the stochastic trajectories are tuned in an optimal way to the deterministic forcing . Frequency resonance can also appear for periodic potentials and for the underdamped Langevin dynamics, and it can be used in order to optimize the effective drift.

We will consider the nonequilibrium dynamics obtained by adding a space-time periodic driving force :

In these equations, $\gamma>0$ is the friction coefficient, $W_{t}$ is a standard $d$-dimensional Brownian motion, the mass matrix $M$ is a positive definite $d\times d$ matrix, and $\beta$ is inversely proportional to the temperature. Throughout this work, we will assume that

The dynamics (2) has been studied in detail for constant external forcings and for particular types of space-dependent forcings . In this paper we will focus on spatiotemporal periodic external forcings. It can be proved that the dynamics (2) has a well-defined steady state (see Propositions 1 and 3). When $\eta=0$ (so that the dynamics is (1)), this steady-state is given by the canonical measure

When the external force is constant, the systematic driving manifests itself through some average nonzero velocity in the system. When the external force is non-constant, and in particular when its space-time average vanishes, it is unclear whether a nontrivial average velocity can be observed. The surprising result by Collet and Martinez is that, in fact, there is in general a nonzero average velocity even if the space-time average of the external force is zero. Our aim in this article is to refine and complete the results obtained in where the case $V=0$ was studied. In particular, we present a more detailed analysis of the problem of convergence to equilibrium for the dynamics (2), we show that this system can exhibit the phenomenon of absolute negative mobility, we analyze the phenomenon of mobility resonance for (2) and we prove a functional central limit theorem (homogenization theorem) for the particle position, under the diffusive rescaling. Our theoretical results are supported by numerical simulations.

The main results of this work are the following.

Proposition 4 and the reformulation (21) show that the linear response of the average velocity is generically nontrivial even if the external force vanishes in average (provided its time-average is non-gradient, see Remark 2). The results we obtain extend the ones presented in , by taking into account a nonzero potential $V$ and giving explicitly the expression of the average time dependent velocity in the system.

Upon adding an appropriate constant force we are also able to find situations in which the average velocity and the average external force experienced by the system are in oppposite directions (see Section 3.2). This phenomenon therefore corresponds to a situation of negative mobility. We emphasize the fact that we observe negative mobility at the level of linear response. In the existing literature negative mobility has been observed at the nonlinear response level, for time-independent forcings; in particular, a subtle interplay between periodic and static forcings at low temperatures were needed for this effect to be observed .

Resonance effects for the average drift in (2) are studied in Section 3.3. We give mathematical properties of the amplitude of the time dependent response as a function of the period of the forcing, and present numerical results illustrating the phenomenon of resonance. To the best of our knowledge, these are the first results on such resonance effects for Langevin dynamics, while there are plenty of studies on the resonance of the average drift for overdamped Langevin dynamics, see for instance .

The effective diffusion obtained in a diffusive space-time scaling is studied in Section 4, for arbitrary forcing magnitudes $\eta$. We show in particular that the effective diffusion matrix varies at second order in $\eta$ when the time average of the forcing is 0.

The rest of the paper is organized as follows. After studying the convergence of the dynamics to its (time dependent) stationary state in Section 2 for arbitrary perturbations $\eta$, we prove in Section 3 various results on the linear response of the average velocity, and finally consider the diffusive regime in Section 4. The proofs of the results presented in Sections 3 to 4 are gathered in Section 5. Numerical simulations for a simple two-dimensional potential, used to illustrate our findings, are presented throughout the paper.

Convergence towards the nonequilibrium steady state

We first introduce some notation. We denote by

the extended phase-space. The (time dependent) generator of the process (2) is $\mathcal{A}_{0}+\eta\mathcal{A}_{1}$, with

We denote by $\mathcal{A}_{0}^{\dagger}$ and $\mathcal{A}_{1}^{\dagger}$ the adjoints of the generators (i.e. Fokker-Planck operators) on the space $L^{2}(\mathcal{E})$. We also introduce the family of Lyapunov functions for $n\geqslant 1$,

and the corresponding $L^{\infty}$ norms on functions $f(t,q,p)$ defined on $\mathcal{E}$:

The first convergence result shows that the process stabilizes around a limit cycle described by a time dependent (periodic) invariant measure.

The invariant distribution is smooth, positive ($\psi_{\eta}(t,q,p)>0$ for all $(t,q,p)\in\mathcal{E}$) and satisfies the Fokker-Planck equation

Finally, it has has finite moments of order $2n$ uniformly in the time variable

and has uniform marginals in the time variable:

Upon averaging in time, standard convergence results can be recovered, such as the following Law of Large Numbers, which will prove useful to study the effective diffusive behavior.

The invariant measure can be fully characterized in the linear response regime as a perturbation around the equilibrium measure $\mu$ defined in (3), for forcings sufficiently small. Let us emphasize that this result is perturbative, in contrast to the convergence statement given by Proposition 1. Similar results have been obtained for different stochastic systems in .

There exists $C,r>0$ such that, for $|\eta|<r$, the invariant measure is given by the following series expansion in $\eta$:

The functions $\varrho_{m}$ are not explicitely known, but are defined as solutions of appropriate Poisson equations (see Section 5.2). The leading order correction $\varrho_{1}$ is particularly important since it governs the linear response.

Linear response of the velocity

For a given perturbation strength, define the time dependent spatially averaged velocity

for any $t\in[0,T]$, and the associated linear response

and first decompose the real-valued external force as

Note that $F_{-n}=\overline{F_{n}}$ (the bar indicating here complex conjugation).

The following result (proved in Section 5.3) shows that each time harmonic of the linear response of the time dependent velocity is directly proportional to the corresponding harmonic of the external force. We will use the notation

for the marginal density of the canonical measure in the position variables ($\widetilde{Z}$ denotes the normalization constant), and, for a given operator $A$, consider the element $Ap$ as the vector with components $Ap_{i}$.

The linear response of the time dependent spatially averaged velocity can be related to the external force as

where the position-dependent diffusion matrix reads

the expectation in the first equality being with respect to canonically distributed initial momenta $p_{0}$, and for all realizations of the equilibrium Langevin dynamics (1) starting from $(q,p_{0})$. In particular, the average (time-independent) velocity depends only on the component $F_{0}$ of the external force:

Proposition 4 refines the results of in two ways: (i) it gives an expression of $\mathscr{V}(t)$ and not only of its time average, and (ii) it highlights the fact that $F_{0}$ solely determines whether the average velocity vanishes or not.

When $F(t,q)=F_{0}(q)=-\nabla W(q)$, the process has an invariant measure whose explicit expression is known: it is the canonical measure associated with the potential energy function $V+\eta W$. In this case, the average velocity should vanish. In fact, the average velocity $\overline{\mathscr{V}}$ is zero as soon as the time-averaged external force is given by the gradient of a scalar function: $F_{0}(q)=-\nabla W(q)$, as can be seen from (13). Indeed, with expectations taken for all initial conditions distributed according to the equilibrium steady state $\mu$ defined in (3) and for all realizations of the equilibrium Langevin dynamics (1),

The second expectation vanishes. For the first expectation we use the fact that, by ergodicity, the law of $q_{\tau}$ has some limiting behavior whatever the choice of $p_{0}$, so that

We illustrate the results from Proposition 4 with some numerical results. We consider a single particle of mass $1$ in the two-dimensional potential

where $q=(x,\,y)$. The numerical scheme is obtained by a Strang splitting between the Hamiltonian part and the fluctuation-dissipation part (including the nonequilibrium forcing). More precisely, denoting by $(q^{n},p^{n})$ approximations of $(q_{n\Delta t},p_{n\Delta t})$ (to simplify the notation, we do not explicitly denote the dependence on $\eta$),

The average, time-independent linear response $\overline{\mathscr{V}}$ is then approximated by fitting the linear dependence of the average velocity

as a function of $\eta$ using a least-square fit.

In the numerical experiments reported below, the numerical parameters are set to $\Delta t=0.01$, $\beta=1$, $\gamma=1$, and the dynamics was integrated over $N=4.5\times 10^{9}$ time-steps. The maximum value of the forcing strength is $\eta_{\rm max}=1$. We performed $R=100$ independent simulations with equally spaced intermediate values $\eta_{\rm max}/R,2\eta_{\rm max}/R,\dots,\eta_{\rm max}$.

We consider non-gradient external forcings of the form

2 Negative mobility

The physical interpretation of (13) is that a nontrivial velocity can be observed when the spatial modes of the position-dependent diffusion matrix $D_{0}(q)$ are excited by the external forcing. The aim of this section is to make this observation precise, and use it to construct situations in which a negative mobility is observed.

Note that the function $\Phi_{0}$ is time-independent. A priori we should consider the above equation as posed on $\mathcal{E}$ and look for $\Phi_{0}\in(L^{2}(\mathcal{E};\mu))^{d}$, but the existence of a time-independent solution (as given by Lemma 5) and the uniqueness of solutions, as given by Lemma 7, enable us to conclude that the time-dependence can be removed.

Since the matrix $D_{0}$ is real, $\mathcal{M}$-periodic and symmetric (by the same time-reversal argument as in Remark 2), it can be written as

where $\mathcal{L}^{*}$ is the reciprocal lattice associated with the lattice $\mathcal{L}$ whose unit cell is $\mathcal{M}$, and $a_{0,L},b_{0,L}$ are real, symmetric $d\times d$ matrices. The following result (see Section 5.4 for the proof) shows that there should be in general a non-trivial $q$-dependence in (19) since, in view of (17),

When the average on the right-hand side is non-zero, it indeed holds $D_{0}(q)\neq D_{0,0}$.

Assume that the following non-degeneracy condition holds:

Then, the smooth function $\Phi_{0}=(\Phi_{0,1},\dots,\Phi_{0,d})$ has a genuine spatial dependence, i.e. $\nabla_{q}\Phi_{0,i}\neq 0$ for all $i\in\{1,\dots,d\}$.

The non-degeneracy condition (20) ensures that the conservative force is sufficiently rich to obtain a coupling between the potential and the driving force. It prevents the use of trivial potential such as $V(q)=v(q_{1})$.

Upon decomposing $F_{0}\widetilde{\mu}$ in Fourier series, the time-averaged external force $F_{0}$ can be written in general as

The normalization is chosen such that $F_{0,0}$ is the canonical average of $F_{0}$. Mean zero forces are characterized by the condition

The space-time averaged linear response of the velocity reads, in view of (13),

This decomposition highlights the two mechanisms separately contributing to the mean velocity: (i) the response to a constant forcing with average $F_{0,0}$, in which case the relevant diffusion matrix is a spatial average of the position dependent diffusion matrix $D_{0}(q)$; (ii) an additional contribution arising from a spatial resonance effect between two terms whose spatial averages vanish, namely the centered diffusion matrix $D_{0}(q)-D_{0,0}$ and the centered (with respect to the canonical measure) force $F_{0}(q)-F_{0,0}\widetilde{\mu}(q)^{-1}$.

Proposition 5 shows that it is possible to obtain a nontrivial velocity even when the external force vanishes on average, upon choosing an element $K_{0}\in\mathcal{L}^{*}\backslash\{0\}$ such that $D_{0,K_{0}}\neq 0$, and choosing for instance a forcing

In fact, the decomposition (21) then shows how to obtain a situation with negative mobility: pick a function $F_{0}$ such that the second term on the righthand side of (21) (denoted by $\overline{\mathscr{V}}_{\rm res}$) is nonzero, and add then a not too large constant force $F_{0,0}$ leading to a contribution $\overline{\mathscr{V}}_{\rm cst}=D_{0,0}F_{0,0}$ with

We now use the numerical algorithm presented in the previous section to study numerically a situation in which negative mobility at the linear response level appears. We consider the force

The parameters are chosen so that the linear response in the $x$ direction is positive, while the average force in the same direction is negative. The average external force actually experienced by the system is estimated over trajectories $(q^{n},p^{n})_{n=1,\dots,N}$ of (15) as

and we retain only the $x$ component (the $y$ component vanishes at first order in $\eta$). The numerical results presented in Figure 2 indeed confirm that the average experienced force and the associated observed velocities are of opposite signs.

3 Resonance of the frequency-dependent mobility

We consider now the dependence of the effective velocity $\mathscr{V}(t)$ on the period $T$ of the external forcing. We restrict ourselves to monochromatic forcings at frequency $\omega=2\pi/T$:

In the sequel, we fix the force $F_{1}(q)=\overline{F_{-1}}(q)$ and only vary the frequency $\omega$. From (11), it is clear that the linear response $\mathscr{V}(t)$ involves a linear response at a single frequency $\omega$:

Frequency resonance corresponds to a nonmonotonic behavior of the magnitude $\left|\widehat{\mathscr{V}}(\omega)\right|$ of the linear response, which has some (local) maximal value for some frequency $0<\omega_{\rm res}<+\infty$. The existence of at least one maximizer in $[0,+\infty)$ is related to the fact that $\widehat{\mathscr{V}}$ has a definite value at $\omega=0$, and vanishes as $\omega\to+\infty$, as stated by the following result (proved in Section 5.5).

In particular, the vectorial amplitude $\widehat{\mathscr{V}}(\omega)$ vanishes as $\omega\to+\infty$.

The interpretation of this result is that, when the forcing period is very small (i.e. $\omega$ is very large), the system does not have time to follow the excitation provided by the external forcing, so that the observed response is very small. The above asymptotic result does not however give any information on the existence of a local maximum of $\left|\widehat{\mathscr{V}}(\omega)\right|$ for $\omega>0$. We therefore have to rely on numerical simulations to study this effect.

We use the notation introduced in the previous section. To obtain the amplitude of the linear response, we consider $R=100$ independent simulations with equally spaced intermediate values $\eta_{\rm max}/R,2\eta_{\rm max}/R,\dots,\eta_{\rm max}$. For each simulation, we compute approximations of the time dependent average velocity $\overline{v}_{\eta}(t)$ in a given direction ($x$ or $y$) as given by (16), and perform a fast Fourier transform to extract the amplitude of the mode at frequency $\omega$. The linear response of this quantity (computed using a least-square fit) gives the desired approximation of $\widehat{\mathscr{V}}(\omega)_{x}$ or $\widehat{\mathscr{V}}(\omega)_{y}$.

Figure 3 presents results on the resonance of the mobility obtained in the $x$ direction for forcings of the form

with a friction $\gamma=0.1$ and a maximal forcing strength $\eta_{\rm max}=0.5$, the other numerical parameters being the same as in Section 3.1. The time-step was refined to $\Delta t=0.001$ instead of $\Delta t=0.01$ for periods $T\leqslant 1$, and extended to $\Delta t=0.025$ for $T\geqslant 20$. The computations show the existence of a resonance frequency around $\omega/(2\pi)\simeq 0.45$. By Proposition 6, the decrease of the amplitude is at least of order $\omega^{-2}$ as $\omega\to+\infty$ since $\nu_{1}=0$. Such a fast decay is difficult to observe numerically.

To observe more easily the predicted decay of the amplitude, we consider the force

for which $\nu_{1}\neq 0$. The numerical results presented in Figure 4 show two interesting features: the resonance peak corresponds to a global maximum of the amplitude, and a decay of order $\omega^{-1}$ can be observed at large frequencies, as predicted by Proposition 6.

Longtime diffusive behavior

In the section we study the longtime diffusive behavior of the nonequilibrium dynamics (2); in particular we show that the diffusively rescaled particle position converges weakly to a Brownian motion with an appropriate diffusion matrix. This result is based on the well developed techniques for proving functional central limit theorems for additive functionals of Markov processes , in particular on the study of an appropriate Poisson equation and on the use of the martingale central limit theorem.

In order to properly define the diffusion matrix, the fundamental ingredient is estimates on the solution of appropriate Poisson equations. One possible result is the following, which gives a polynomial control on derivatives of arbitrary order.

Consider a smooth function $f$ with derivatives growing at most polynomially in $p$. Then, there exists a unique solution $\Phi_{\eta}$ to the Poisson equation (recall that the phase-space $\mathcal{E}$ is introduced in (4))

Moreover, for any $k\geqslant 1$ and $\eta_{*}>0$, there exists a real constant $C>0$ and integers $n,m,N\geqslant 1$ such that, for all $\eta\in[-\eta_{*},\eta_{*}]$,

This result is proved in Section 5.6, as an extension of a technique first employed by Talay in , which was recently carefully rewritten and extended in [24, Appendix A]. We emphasize the fact that the equation (25) is equipped with periodic boundary conditions in both $t$ and $q$. The boundary condition for $p$ is that $\Phi_{\eta}\in L^{2}(\psi_{\eta})$. We also remark that we do not need precise information on the invariant measure $\psi_{\eta}$. In particular, it need not be a perturbation of the invariant measure of the equilibrium dynamics $\mu$.

2 Definition of the effective diffusion matrix

The behavior of the process depends in the long time limit on the space-time scaling. In the hyperbolic scaling where time and space are renormalized by the same factor, the behavior is determined by the average velocity

In general, as discussed in Section 3, this average velocity is not zero, and is of order $\eta\overline{\mathscr{V}}$ when $\eta$ is small (see (59) below). The behavior of the process over longer times, in a diffusive space-time scaling, describes the deviations around the average velocity. We introduce for this study the process

The only difference between $q_{t}^{\eta}$ and $\mathscr{Q}_{t}^{\eta}$ is that $\mathscr{Q}_{t}$ is not reprojected in $\mathcal{M}$ by the periodization procedure, and hence can diverge as time passes. In the long time limit, $\mathscr{Q}_{t}^{\eta}$ drifts linearly as $t\mathcal{V}_{\eta}$. The diffusive behavior is captured by the centered process

considered in a diffusive space-time scaling: for any $\varepsilon>0$,

We assume that the process starts at stationary initial conditions, namely $(q^{\eta}_{0},p^{\eta}_{0})\sim\psi_{\eta}(0,q,p)\,dq\,dp$. The stationarity assumption can be removed at the expense of additional technical difficulties. Introduce finally the following Poisson equation (well defined by Proposition 7)

We can then state the following convergence result.

As $\varepsilon\to 0$, the rescaled process $Q_{t}^{\eta,\varepsilon}$ converges weakly on any finite time interval to the effective $d$-dimensional Brownian motion

with symmetric, positive definite covariance matrix $\mathscr{D}_{\eta}$ defined by its action on test vectors:

and initial conditions $\overline{Q}_{0}\sim\widetilde{\psi}_{\eta}(q)\,dq$, where

The proof of Theorem 3, based on the decomposition technique presented in (see also for an up-to-date account and further references) is quite standard. It is nonetheless provided in Section 5.7 for completeness. Note that a little more work would allow to obtain convergence rates, by slightly extending the approach from . Related homogenization results for the overdamped dynamics with space-time periodic coefficients can be found in .

Formal asymptotic expansions for the corresponding Fokker-Planck equation would lead to generally nonsymmetric homogenized diffusion matrix, whose symmetric part is the covariance matrix $\mathscr{D}_{\eta}$. Indeed, when the microscopic dynamics is reversible, then the homogenized diffusion matrix is symmetric . The homogenized diffusion matrix can be symmetric also for nonreversible dynamics, when additional symmetries are present . The antisymmetric part can affect the homogenized dynamics when the homogenized diffusion matrix is space dependent. This would be the case if we considered locally periodic coefficients, as in .

A simple computation (based on (58) below) shows that, when $\eta=0$, the effective diffusion matrix $\mathscr{D}_{0}$ is related to the position dependent diffusion matrices $D_{0}(q)$ introduced in Proposition 4 as

The homogenization problem for the equilibrium dynamics is well studied and the properties of the diffusion tensor $\mathscr{D}_{0}$ (i.e. its scaling with respect to the friction coefficient) are understood, at least in one dimension .

3 Properties of the diffusion matrix for small forcings

where $\mathscr{D}_{0}$ is the effective diffusion of the equilibrium dynamics corresponding to $\eta=0$ and $\widetilde{\mathpzc{D}}_{\eta,\xi}$ is uniformly bounded for $|\eta|\leqslant r$ and $|\xi|\leqslant 1$. When the external force has time average 0 for all configurations, namely when

then the first order correction vanishes: $\mathpzc{D}_{1}=0$.

We now present numerical experiments illustrating the previous theoretical results. We solve numerically the Langevin dynamics (2), discretized with (15) for the potential (14), so that we expect a coupling between the directions $x$ and $y$. We calculate the drift vector and the diffusion matrix using empirical averages of trajectories generated by (15). One could also, in principle, solve the stationary Fokker-Planck equation and the Poisson equation using a spectral method, as was done in for the equilibrium dynamics and for constant external forcings. Here we rely on Monte Carlo simulations. The drift vector and the diffusion matrix are respectively given by

We are interested both in the linear response regime $\eta\to 0$ (for which we expect $\mathscr{D}_{\eta}\to\mathscr{D}_{0}$) and the behavior of $\mathscr{D}_{\eta}$ for large $\eta$. To estimate the drift coefficients and the diffusion tensor, we proceed by a multiple replica strategy. The system is initialized by positions and momenta $(q^{0},p^{0})$ distributed according to the nonequilibrium measure $\psi_{\eta}(0,q,p)\,dq\,dp$. In practice we run the dynamics (15) for a simulation time $\tau_{\rm neq}$ sufficiently large for the system to converge towards the nonequilibrium steady-state.

For simplicity, we consider nonequilibrium forcings with general form

In the example we considered, we picked $F_{0}$ as given by (22) and chose either $\omega=2\pi$ (time dependent forcing) or $\omega=0$ (time-independent forcing). To perform the computations, we fix a (large) number $M$ of replicas and a simulation time $\tau_{\rm sim}$ sufficiently large for the asymptotic diffusive regime to be attained. The drift coefficient and diffusion matrix are evaluated at time $\tau_{\rm sim}$ by empirical averages. More precisely,

where the approximation of the unprojected position of the $k$th replica at time $n\Delta t$ is

In the simulations reported below, we used $M=2.5\times 10^{7}$, $\tau_{\rm sim}=1500$, $\tau_{\rm neq}=500$, $\Delta t=0.01$, $\gamma=1$, $\beta=1$. Since we use an independent replica strategy, error bars are deduced from the empirical variance. They are not reported here since statistical errors are in all cases below 10% in relative accuracy.

Figure 5 depicts the component of the diffusion matrix as a function of the forcing strength $\eta$. Figure 6 depicts the spectrum of the diffusion matrix as a function of $\eta$. As can be seen, the $\mathscr{D}_{\eta,xx}$ component increases much more than the other components as $\eta$ becomes large. The $\mathscr{D}_{\eta,yy}$ component on the other hand remains almost stationary. A zoom on the small $\eta$ variations of $\mathscr{D}_{\eta,xx}$ is presented in Figure 7, together with a least-square fit of the form

The leading coefficient $a$ is found to be more or less equal to 0 for time-dependent forcings, as predicted by Proposition 8. It is however also found to be more or less equal to 0 for the time-independent forcing under consideration.

Proofs of our results

The idea of the proof is the following. We first start by establishing the convergence of the distribution of a sampled Markov chain obtained by considering the time-inhomogeneous process (2) at times $nT$ (a standard idea, considered in for instance). An invariant measure is then obtained by evolving the invariant measure of the sampled process over a period.

We introduce the Markov chain $(Q^{\eta}_{n},P^{\eta}_{n})=(q^{\eta}_{nT},p^{\eta}_{nT})$, where the time-inhomogeneous process $(q^{\eta}_{t},p^{\eta}_{t})$ is started at time $t=0$ from $(Q^{\eta}_{0},P^{\eta}_{0})=(q_{0},p_{0})$. We have indicated explicitly the dependence on the forcing magnitude $\eta$ although all the estimates below will hold uniformly with respect to this parameter as long as it remains the bounded interval $[-\eta_{*},\eta_{*}]$. The generator $U_{T,\eta}$ of the Markov chain $(Q_{n}^{\eta},P_{n}^{\eta})$ is defined as

The convergence result for the sampled Markov chain, with rates uniform in $\eta$, is based on the following two lemmas (proved at the end of this section).

For any $n\geqslant 1$ and $\eta_{*}>0$, there exist $b>0$ and $a\in[0,1)$ such that, for all $\eta\in[-\eta_{*},\eta_{*}]$,

where $\mathscr{B}(X)$ are the Borel sets of $X$.

From the results of , we can then state the following uniform convergence result for the sampled chain.

Fix $\eta^{*}>0$ and $n\geqslant 1$. There exist $\lambda,C>0$ and probability measures $m_{\eta}(q,p)\,dq\,dp$ such that, for any $f\in L^{\infty}_{\mathcal{K}_{n}}$ and any $\eta\in[-\eta^{*},\eta^{*}]$,

Moreover, the integration of the inequality (33) with respect to $m_{\eta}$ gives the moment estimate

To obtain the Law of Large Numbers for all initial conditions, we use the following property (proved at the end of this section).

The generator $U_{T,\eta}$ has a transition kernel which is absolutely continuous with respect to Lebesgue measure and positive.

This property implies that the chain is irreducible with respect to the Lebesgues measure, and also gives the positivity of $m_{\eta}(q,p)$. In view of the Lyapunov condition (33) and relying on [31, Theorem 17.0.1], we can then conclude that, for any $f\in L^{\infty}_{\mathcal{K}_{n}}$,

for almost all initial conditions $(Q^{\eta}_{0},P^{\eta}_{0})$. In fact, convergence occurs for all initial conditions and not simply for almost all initial conditions. This is a consequence of the smoothness of the transition probability, ensured by Lemma 3 (the chain is Harris recurrent, see [45, Corollary 1], based on ).

Now that we have proved the exponential convergence to the steady-state, we characterize the invariant measure as the weak solution of an appropriate Fokker-Planck equation. Note that, for a smooth function $f$,

which gives (7) in the sense of distributions. The fact that $\psi_{\eta}$ is smooth and that (7) actually holds pointwise is a consequence of the following lemma.

The operators $\partial_{t}+\mathcal{A}_{0}+\eta\mathcal{A}_{1}$ and $-\partial_{t}+\mathcal{A}_{0}^{\dagger}+\eta\mathcal{A}_{1}^{\dagger}$ (considered on $L^{2}(\mathcal{E})$) are hypoelliptic.

We use Hörmander’s criterion, and present the proof for $\partial_{t}+\mathcal{A}_{0}+\eta\mathcal{A}_{1}$. The proof for the adjoint of this operator is similar. Define $X_{0}=\partial_{t}+M^{-1}p^{T}(\nabla_{q}-\gamma\nabla_{p})-\nabla V^{T}\nabla_{p}+\eta\mathcal{A}_{1}$ and

In addition, denoting by $[L_{1},L_{2}]=L_{1}L_{2}-L_{2}L_{1}$ the commutator between two operators $L_{1}$ and $L_{2}$, a simple computation shows that

from which derivatives in the directions $q_{j}$ are recovered. Linear combinations with $X_{0}$ and $X_{1,i}$ finally allow to recover $\partial_{t}$, so that the Lie algebra generated by $X_{0}$, $X_{1,i}$ and $[X_{1,i},X_{0}]$ is full. ∎

We conclude this section by giving the proofs of the technical lemmas used above.

Note that the martingale $\mathscr{G}_{k}$ is distributed according to a Gaussian with mean 0 and covariance $M(1-\exp(-2\gamma TM^{-1}))/\beta$, hence has finite moments of all orders; while a simple uniform bound on $\mathscr{F}_{k}$ is for instance

for any $\varepsilon>0$, chosen so that $\alpha^{2}(1+\varepsilon)<1$. Denoting by $\mathcal{F}_{k}$ the filtration induced by the Markov chain up to the $k$th step, it holds

for some constant $C_{\varepsilon}>0$ which can be chosen to be independent on $\eta\in[\eta_{*},\eta_{*}]$ (but depending on $\eta_{*}$ of course). This gives the result in the case $n=2$.

For a general index $n\geqslant 2$, we take the $n$th power of the bound (39). The leading term is $\alpha^{2n}(1+\varepsilon)^{n}$, which is still strictly smaller than 1 with the above choice of $\varepsilon$. Terms with odd powers of $\mathscr{G}_{k}$ vanish when taking the expectation. The expectation of the remainder is a linear combination of terms $(P_{k}^{\eta})^{2k}$ with $k\leqslant n-1$. Therefore,

where $\mathcal{P}$ is a polynomial of degree at most $2(n-1)$, with bounded, positive coefficients. As $|P_{k}^{\eta}|\to+\infty$,

This shows that, for any $\delta>0$, there exists $C_{\delta}>0$ such that $\mathcal{P}(|P_{k}^{\eta}|)\leqslant\delta|P_{k}^{\eta}|^{2n}+C_{\delta}$ (consider a radius $R_{\delta}>0$ for which $\mathcal{P}(r)/r^{2n}\leqslant\delta$ when $r\in[0,R_{\delta}]$ and define $C_{\delta}$ as the supremum of $\mathcal{P}$ over $[0,R_{\delta}]$). The conclusion easily follows upon choosing $\delta>0$ sufficiently small. ∎

is uniformly bounded by $mT\mathfrak{F}_{\eta}$ (where $\mathfrak{F}_{\eta}$ and $m$ are defined in (37) and (38), respectively). The random variables $\widetilde{\mathscr{G}}_{k},\mathscr{G}_{k}$ are correlated centered Gaussian random variables, independent of $\widetilde{\mathscr{G}}_{l},\mathscr{G}_{l}$ for $l\neq k$. Their covariance reads

where $\mathscr{Q}_{k}=Q_{k}^{\eta}+(1-\alpha_{T})P^{\eta}_{k}/\gamma+\widetilde{\mathscr{F}}_{k}$ and $\mathscr{P}_{k}=\alpha_{T}P_{k}^{\eta}+\mathscr{F}_{k}$ are uniformly bounded in norm by a constant $R>0$ which depends only on $p_{\rm max}$ and $\eta_{*}$. The proof is therefore concluded by defining the probability measure

To prove the smoothness of the transition kernel, we use the same argument as in the proof of [44, Lemma 2.2] and rewrite, thanks to Girsanov’s formula, the nonequilibrium evolution $(q_{t}^{\eta},p_{t}^{\eta})$ solution of (2) in terms of the equilibrium evolution $(q_{t}^{0},p_{t}^{0})$ given by (1). Starting from an initial condition $(q_{0},p_{0})$ at time $t=0$ and integrating over one period $T$, the Girsanov transform (see e.g. ) shows that the law of $(q_{t}^{\eta},p_{t}^{\eta})$ is absolutely continuous with respect to the law of $(q_{t}^{0},p_{t}^{0})$ and that, for any bounded, measurable function $f$,

The positivity of the transition kernel can be proved using a standard control argument, see for instance [41, Section 6] or . ∎

2 Proof of Proposition 3

The idea of the proof is to show, by a perturbative construction, the existence of a solution to the extended Fokker-Planck equation (7). Proposition 1 then gives the uniqueness. We follow the strategy summarized in and already used in . See also .

Instead of the flat space $L^{2}(\mathcal{E})$, we consider as a reference space the Hilbert space $\mathcal{H}=L^{2}(\mathcal{E},\mu)\cap\{\mathbf{1}\}^{\perp}$, the space of functions $f(t,q,p)$ such that

The Fokker-Planck equation (7) can be rewritten as

where adjoints are now taken on $\mathcal{H}$. Formally,

To make this argument rigorous, we use two technical results (proved at the end of this section).

The operator $\partial_{t}+\mathcal{A}_{0}$ is invertible on $\mathcal{H}$.

The operator $\mathcal{A}_{1}$ is $(\partial_{t}+\mathcal{A}_{0})$-bounded on $L^{2}(\mathcal{E},\mu)$: There exists $a,b>0$ such that, for all smooth functions $f$,

Therefore, the operator $(-\partial_{t}+\mathcal{A}_{0}^{*})^{-1}\mathcal{A}_{1}^{*}$ is bounded on $\mathcal{H}$. This proves the validity of (41) for $|\eta|$ smaller than the spectral radius of $(-\partial_{t}+\mathcal{A}_{0}^{*})^{-1}\mathcal{A}_{1}^{*}$.

We conclude this section with the proofs of Lemmas 5 and 6. We first recall the following result, which is a simple consequence of the results in .

Consider, for a given function $g\in\mathcal{H}$, the equation $(\partial_{t}+\mathcal{A}_{0})f=g$. We decompose $f$ and $g$ in Fourier series as

Therefore, $(\partial_{t}+\mathcal{A}_{0})f=g$ can be rewritten as

from which the conclusion easily follows. ∎

3 Proof of Proposition 4

We first make precise the equation satisfied by the linear response term $\varrho_{1}$. From (41),

Decomposing both sides of the equation in Fourier series, we obtain

This equation is well posed for $n\neq 0$ in view of Lemma 7. For $n=0$, it is also well posed since $p^{T}M^{-1}F_{0}$ has a vanishing average with respect to $\mu(q,p)\,dq\,dp$.

Using the expression (43) of the first order correction in $\eta$ of the invariant measure, the linear response of the average velocity is

Using this result, and since $M^{-1}p$ has all his components in $\widetilde{H}^{1}(\mu)$, we rewrite the linear response of the velocity as

4 Proof of Proposition 5

Recall first that, by hypoellipticity, $\Phi_{0}$ is smooth. Assume, by contradiction, that $\Phi_{0}$ does not depend on $q$ so that $\Phi_{0}=\Phi_{0}(p)=(\Phi_{0,1}(p),\dots,\Phi_{0,d}(p))$. In this case,

5 Proof of Proposition 6

The proof is an immediate consequence of the following lemma (upon replacing $\psi$ by the components $p_{i}$ of $p$, and changing the sign of $\omega$).

Consider $\psi\in L^{2}(\mu)$ such that $\mathcal{A}_{0}^{m}\psi\in L^{2}(\mu)$ for any $0\leqslant m\leqslant n$. Introduce, for a frequency $\omega\neq 0$, the unique solution $\varphi_{\omega}$ (well defined by Lemma 7) of

Then, there exist a constant $C_{n}>0$ and functions $\phi_{2},\dots,\phi_{n-1}$ such that, for all $\omega\geqslant 1$,

6 Proof of Proposition 7

The first remark is that, by linearity, $\Phi_{\eta}$ is the sum of the solution of (25) with right-hand side $f-\overline{f}_{\eta}(t)$ (see (6) for the definition of $\overline{f}_{\eta}(t)$) and

Estimates on solutions of the above equation are obtained by generalizing a technique first employed by Talay in , and recently carefully rewritten and extended in [24, Appendix A]. We briefly recall the beginning of the argument here in order to show that it can be extended to processes with time-periodic drivings. The crucial estimate is the one provided by Lemma 10 below.

where the process $(t_{s},\widetilde{q}_{s},\widetilde{p}_{s})$ has spatial initial conditions $(\widetilde{q}^{\eta}_{0},\widetilde{p}^{\eta}_{0})=(q_{0},p_{0})$ and evolves according to the time homogeneous Markovian dynamics with generator

The solutions of (48) for a given initial condition $(t_{0},q_{0},p_{0})$ can be seen as the solutions of (2) started at time $t_{0}$ rather than 0. The exponential convergences (5) and (34) (and their immediate generalizations for $t_{0}\neq 0$) imply that, for any $n\geqslant 1$ and all initial conditions $(t_{0},q_{0},p_{0})$,

An integration of this inequality from $s=0$ to $+\infty$ shows that the function

is a well defined, admissible solution of the Poisson equation (46). In fact it is also the only one. This already gives, for some constant $C>0$ depending only on $n$ and $\eta_{*}$, the following pointwise estimate:

Upon choosing $l>2n+d/2$, we also obtain the following weighted $L^{2}$ exponential convergence from (49):

where we have introduced the polynomial weight

The idea of Talay in [44, Section 3.4] is to use appropriate mixed derivatives to obtain pointwise decay estimates for the derivatives of $\Phi_{\eta}$, by first obtaining decay estimates in weighted $L^{2}$ spaces and then concluding by appropriate Sobolev embeddings. A key estimate to perform the computation is the following lemma (compare [24, Lemma A.6]).

The proof is based on the following computations:

In addition, denoting by $T_{\eta}^{\dagger}$ the adjoint of $T_{\eta}$ on the flat space $L^{2}(\mathcal{E})$,

A simple computation shows that there exists $A_{k,\eta_{*}}>0$ such that, for any $\eta\in[-\eta_{*},\eta_{*}]$ (see the expression of $T_{0}^{\dagger}\Pi_{k}$ in the proof of [28, Lemma 25])

In the remainder of the proof, the constants may vary from line to line. To use Lemma 10, we will need the following commutators:

is uniformly bounded for $s\geqslant 0$. We next consider $L_{i}=\partial_{q_{i}}-R\partial_{p_{i}}$ for some real parameter $R>0$, which mixes derivatives in $q$ and $p$ in order to retrieve some dissipation in the $q$ direction (this was one of the main ideas of hypocoercivity, as mentioned in the introduction of ). Note that

so that, using a Cauchy-Schwarz inequality and the boundedness of $\nabla(\partial_{q_{i}}V)+\eta\partial_{q_{i}}F$,

Lemma 10 together with (55) then shows that it is possible to choose $R>0$ sufficiently large so that, upon increasing the value of $K$ and choosing $l>2n^{\prime}+d/2$ (in view of the integral on the initial condition $u(0;t,q,p)=f(t,q,p)-\overline{f}_{\eta}(t)$),

Therefore, by integration (and upon changing $K$ again),

uniformly in $s\geqslant 0$. The combination of (55) and (56) finally gives a control on $\nabla_{q}u$ as

We can now use Lemma 10 with $L=\partial_{p_{i}}$ or $\partial_{q_{i}}$ (with a Cauchy-Schwarz inequality to treat the last term on the right-hand side of (53), namely the one involving the commutator) to conclude that

Sobolev embeddings then imply pointwise estimates of the form

from which (26) is obtained upon integrating in $s$ (see (50)).

7 Proof of Theorem 3

where $\mathcal{V}_{\eta,i}$ denotes the $i$th component of the vector $\mathcal{V}_{\eta}$. By Itô calculus (note that, in view of Proposition 7, the derivatives of $\Phi_{\eta}$ are well defined elements of $L^{\infty}(L^{\infty}_{{\mathcal{K}_{n}}})$ for some integer $n$ sufficiently large, hence $\mathcal{A}_{0}\Phi_{\eta}$, etc, make sense in this functional space),

The quadratic variation of $\mathscr{M}^{\eta,\xi}_{s}$ is

The proof of the weak convergence to a Brownian motion is very standard: we prove the convergence of the finite dimensional laws using the martingale central limit theorem , as well as the tightness of the stochastic process, and conclude with [3, Theorem 7.1]. Alternatively, we could also follow the more quantitative approach from (based on the use of an additional Poisson equation to get better control between $\xi^{T}Q^{\eta,\epsilon}_{t}$ and the limiting process, projected onto $\xi$) in view of the good resolvent estimates provided by Proposition 7.

Let us first prove the tightness by using Prohorov’s theorem. It is sufficient to prove that, for any $\alpha>0$ and any $\tau>0$,

Note that, using for instance the $L^{\infty}(L^{\infty}_{\mathcal{K}_{2}})$ bound, i.e. (51) with $n=2$,

uniformly in $\delta>0$ (using the finiteness of the moments (8) and the smoothness of $\psi_{\eta}$). On the other hand, by Doob’s inequality,

where the final bound relies on $\lceil x\rceil=[x]+1\leqslant x+1$ and $\lfloor x\rfloor=[x]-1\geqslant x-1$. To go from the second to the third line, we have used the trivial bound

This owes to the fact that the law of $(q^{\eta}_{\theta},p^{\eta}_{\theta})$ is $\psi_{\eta}([\theta],q,p)$ since the dynamics is started at time 0 from $\psi_{\eta}(0,q,p)\,dq\,dp$. Combining the estimates we have shown, (57) easily follows.

Define $n_{\varepsilon}(t)=\lfloor t/(T\varepsilon^{2})\rfloor$. Then,

with $\varepsilon\sqrt{n_{\varepsilon}(t)}\to\sqrt{t/T}$ and

the first condition of [2, Theorem 35.12] follows from the Law of Large Numbers (9):

by dominated convergence, which implies the second condition of [2, Theorem 35.12]. In conclusion, \xi^{T}\Big{(}Q_{t}^{\eta,\varepsilon}-Q_{0}^{\eta,\varepsilon}\Big{)} converges in law to a standard Gaussian distribution, with variance $2t\,\xi^{T}\mathscr{D}_{\eta}\xi$. The extension to a general finite dimensional distributions is done using an iterative procedure, as carefully documented in [1, Section 3.4] for instance.

It is clear that $\mathscr{D}_{\eta}$ is symmetric positive. To prove that it is definite, assume that there exists $\xi\neq 0$ such that $\xi^{T}\mathscr{D}_{\eta}\xi=0$, in which case

This shows that $\nabla_{p}\left(\xi^{T}\Phi_{\eta}\right)=0$ almost everywhere since $\psi_{\eta}(t,q,p)>0$ by Proposition 1, hence $\xi^{T}\Phi_{\eta}$ does not depend on $p$. But then, the equality

requires $\nabla_{q}\left(\xi^{T}\Phi_{\eta}\right)=\xi$, which cannot be satisfied for a smooth, periodic function of the $q$ variable. The contradiction shows that $\mathscr{D}_{\eta}$ is definite.

8 Proof of Proposition 8

Since the operator $T_{\eta}$ defined in (47) satisfies (54), and using, by definition of the invariant measure,

where $\overline{\mathscr{V}}$ is defined in (13) and $\widetilde{\mathcal{V}}_{\eta}^{2}$ collects the higher order terms and is uniformly bounded for $|\eta|$ sufficiently small. The expansion of the solution of the Poisson equation is given by the following lemma.

The solution $\Phi_{\eta}$ of the Poisson equation (27) can be expanded as

where $\Phi_{0}$ is defined in (18) and $\widetilde{\Phi}^{1}$ satisfies the Poisson equation

while there exists some $m\geqslant 1$ such that the remainder $\widetilde{\Phi}^{2}_{\eta}$ is uniformly bounded in the $L^{\infty}(L^{\infty}_{\mathcal{K}_{m}})$ norm for $|\eta|$ sufficiently small.

Gathering all these expansions, we see that

and where $\widetilde{\mathpzc{D}}_{\eta,\xi}$ is uniformly bounded for $|\eta|\leqslant r$ and $|\xi|\leqslant 1$.

We now turn to the specific case when the condition (28) holds. Note that there is only one function depending on time in each term on the right-hand side of the definition of $\xi^{T}\mathpzc{D}_{1}\xi$. The idea is to show that the time average of these functions is 0, which will then give the claimed result. We define the time-average of a function on $\mathcal{E}$ as

When (28) holds, an integration of (42) with respect to the time variable gives

which also gives $\widehat{\widetilde{\Phi}^{1}}=0$. Finally, it indeed holds $\xi^{T}\mathpzc{D}_{1}\xi=0$ when (28) holds.

The function $\widetilde{\Phi}^{1}$ introduced in (60) is indeed well defined in view of Lemmas 5. A simple computation shows that

Indeed, the linear term disappears since (see (42), (18) and the discussion after this equation)

which is the constant term on the right-hand side of the Poisson equation defining $\widetilde{\Phi}^{1}$ in (60). We then conclude from (61) and the estimates provided for instance by Proposition 7 in the case $\eta=0$. These estimates indeed show that $\Phi_{0}$ (and hence $\mathcal{A}_{1}\Phi_{0}$) is smooth with derivatives growing at most polynomially. Next, from (60), $\widetilde{\Phi}^{1}$ is smooth and this function and its derivatives grow at most polynomially. ∎

Acknowledgements

This work was initiated while GP was visiting the INRIA team MICMAC (now MATHERIALS) at CERMICS. The hospitality and financial support from INRIA are greatly acknowledged. RJ’s research is supported by the EPSRC through grant EP/J009636/1. GP’s research is partially supported by the EPSRC through grants EP/J009636/1 and EP/H034587/1. GS’s research is partially supported by the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement number 614492. The authors benefited from discussions with Stefano Olla and Stephan De Bièvre.

References