Nonequilibrium Langevin dynamics: a demonstration study of shear flow fluctuations in a simple fluid

TL;DR

This paper demonstrates how second-order Langevin dynamics can accurately model shear flow fluctuations in a simple fluid, combining theoretical derivation with numerical simulation to validate the approach.

Contribution

It determines all stochastic parameters for a second-order Langevin model of shear flow and compares the model directly with molecular dynamics data.

Findings

The Langevin model accurately reproduces flow autocorrelations.

A linear relation between external force and flow fluctuation cumulants is identified.

Introduction of a new quantity: athermal cumulant ratio.

Abstract

The present study is based on a recent success of the second-order stochastic fluctuation theory in describing time autocorrelations of equilibrium and nonequilibrium physical systems. In particular, it was shown to yield values of the related deterministic parameters of the Langevin equation for a Couette flow in a microscopic Molecular Dynamics model of a simple fluid. In this paper we find all the remaining constants of the stochastic dynamics, which is then numerically simulated and directly compared with the original physical system. By using these data, we study in detail the accuracy and precision of a second-order Langevin model for nonequilibrium physical systems, theoretically and computationally. In addition, an intriguing relation is found between an applied external force and cumulants of the resulting flow fluctuations. This is characterized by a linear dependence of athermal cumulant ratio, a new quantity introduced here.

Figures

Tables

Table 1. Table 1: Comparison of statistics, in reduced units, between our MD and LD simulations for the shear current α ​ ( t ) = P y ​ x 𝛼 𝑡 subscript 𝑃 𝑦 𝑥 \alpha(t)=P_{yx} . We implemented separately two LD models, based on (ex) athermal white exponential noise, and (sh) athermal white exponential shot noise.

	LD(ex)	LD(sh)	MD
Mean, ${\kappa }_1\left(\alpha \right)$	$-1.515\pm 0.004$	$-1.518\pm 0.004$	$-1.521\pm 0.004$
Variance, ${\kappa }_2\left(\alpha \right)$	$0.185\pm 0.003$	$0.185\pm 0.003$	$0.184\pm 0.003$
Skewness, ${\kappa }_3\left(\alpha \right)/{\kappa }_2{\left(\alpha \right)}^{3/2}$	$-0.11\pm 0.03$	$-0.16\pm 0.03$	$-0.14\pm 0.03$
Excess kurtosis, ${\kappa }_4\left(\alpha \right)/{\kappa }_2{\left(\alpha \right)}^2$	$0.15\pm 0.09$	$0.19\pm 0.07$	$0.08\pm 0.07$

Equations

$$\ddot{\alpha }(t)+a\dot{\alpha }(t)+b^2\alpha (t)=r(t),$$

$$\ddot{\alpha }(t)+a\dot{\alpha }(t)+b^2\alpha (t)=f,$$

$$r^{\mathrm{NE}}(t)=A\omega (t)+Bϵ(t/\tau ),$$

$$\ddot{\alpha }(t)+a\dot{\alpha }(t)+b^2\alpha (t)=F+A\omega (t)+Bϵ(t/\tau ).$$

$$P_{yx}=V^{-1}\sum _{i=1}^N(p_{yi}{p}_{xi}/m+y_i{F}_{xi}),$$

$$R(t)={\int }_0^{t}d{t}^{\prime }[F+A\omega (t^{\prime })+Bϵ(t^{\prime })].$$

$${\kappa }_1(\alpha )$$

$${\kappa }_3(\alpha )$$

$${\kappa }_1(\alpha )$$

$${\kappa }_3(\alpha )$$

$${\kappa }_1(\alpha )$$

$${\kappa }_3(\alpha )$$

$${\kappa }_2(\alpha |\gamma =0)=\frac{A^2}{2a{b}^2};$$

$$\chi (\alpha )={\kappa }_2(\alpha )-\frac{A^2}{2a{b}^2}.$$

$$B_{\mathrm{ex}}$$

$${\tau }_{\mathrm{ex}}$$

$$k(P_{yx})=\frac{{\chi }^2(P_{yx})}{{\kappa }_1(P_{yx}){\kappa }_3(P_{yx})},$$

$$k\propto ({\mathrm{const}}_{1}F\tau /B+{\mathrm{const}}_2{)}^{-1},$$

$$\ddot{\alpha }(t)+a\dot{\alpha }(t)+b^2\alpha (t)\approx f+A\omega (t),$$

$$\Omega (t)={\int }_0^{t}d{t}^{\prime }A\omega (t^{\prime })\text{,}$$

$$\mu (\tilde{x})$$

$$\kappa (\tilde{x})$$

$$\kappa (\tilde{\xi })=\sum _{i=0}^n{\kappa }_i(\xi )\frac{{\tilde{\xi }}^i}{i!}={\kappa }_2(\xi ){\tilde{\xi }}^2/2+\mathcal{O}({\tilde{\xi }}^4),$$

$${\kappa }_G(\tilde{\xi })={\kappa }_2(\xi ){\tilde{\xi }}^2/2,$$

$$\kappa (\tilde{x})\approx n{\kappa }_2(\xi ){\tilde{x}}^2/2.$$

$$\kappa (\tilde{x})\approx \frac{{\kappa }_2(\xi )t{\tilde{x}}^2}{2\tau },$$

$${\mu }_{[1/1]}(\tilde{\xi })=\frac{a_{1,1}+a_{1,2}\tilde{\xi }}{a_{2,1}+a_{2,2}\tilde{\xi }},$$

$${\mu }_{[1/1]}(\tilde{\xi })=\frac{1+b_1\tilde{\xi }}{1+b_2\tilde{\xi }}=\frac{b_1}{b_2}+\frac{1-b_1/b_2}{1+b_2\tilde{\xi }}=q+\frac{1-q}{1-b\tilde{\xi }}.$$

$$p_{[1/1]}(\xi )=q\delta (\xi )+(1-q)H[\mathrm{sign}(B)\xi ]\exp (-\xi /B)/B,$$

$$E_q(t/\tau )={\int }_0^{t}d{t}^{\prime }B{ϵ}_q(t/\tau ),$$

$$\kappa ({\tilde{E}}_q)$$

$$p({\tilde{E}}_q)$$

$$\kappa ({\tilde{P}}_q)=\frac{[\mu (\tilde{\xi })-1]t}{{\tau }^{\prime }}{\text{vline}}_{\tilde{\xi }={\tilde{P}}_q},$$

Full text

keywords

I Introduction

The Langevin dynamics, originally inspired by the problem of Brownian motion (Coffey and Kalmykov, 2012, Chapters 1-2), has now become the fundamental stochastic model of fluctuations for equilibrium physical systems at mesoscopic scales Nyquist (1928); Onsager and Machlup (1953); Machlup and Onsager (1953). Its generalization to nonequilibrium steady states, though, is still actively developed, as suggested by a number of recent publications Kanazawa et al. (2015a, b); Morgado and Queirós (2016); Duarte Queirós (2016); Belousov et al. (2016a, b); Belousov and Cohen (2016). A common objective of these studies is to provide a statistical account of an externally applied force together with its spontaneous variations, which manifest themselves in fluctuations of the resulting conjugate current in a system of interest.

In our most recent paper Belousov and Cohen (2016) we showed, that a second-order Langevin equation, suggested in Ref. Machlup and Onsager (1953), provides excellent means for quantitative studies of fluctuations at mesoscopic scales in both, equilibrium and nonequilibrium, steady-state systems. In particular, it yields an accurate analytical model of the time autocorrelation function for currents, which was successfully applied to evaluate the Green-Kubo formula for a transport coefficient. For a general fluctuating quantity $\alpha(t)$, the Langevin equation of the second order in time, $t$, reads:

[TABLE]

where $a>0$, $b>0$ are constants, while $r(t)$ is a random noise.

In the macroscopic limit, Eq. (1) transforms into a deterministic equation (Onsager and Machlup, 1953, Sec. 2-3) and, therefore, $r(t)$ must then in general become a constant, $f\in\mathbb{R}$, so that we have

[TABLE]

which is a differential equation for a damped harmonic oscillator, subject to an externally applied macroscopic force $f$. Furthermore, its solution for $\alpha(t)$ converges with time to the steady-state ensemble average of Eq. (1), e.g. Belousov and Cohen (2016), with $\langle\alpha(t)\rangle=f/b^{2}$ and $\langle\dot{\alpha}(t)\rangle=0$. By comparing Eqs. (1) and (2), one can see that, at the mesoscopic scales, $r(t)$ represents the external force, $f=\langle r(t)\rangle$, as well as its spontaneous variations Onsager and Machlup (1953).

In the equilibrium regime Nyquist (1928); Onsager and Machlup (1953); Belousov et al. (2016b), $f=0$, the stochastic term $r(t)$ represents effects of the thermal fluctuations and is commonly given by $r(t)=A\omega(t)$, where $A>0$ is a constant proportional to the square root of the system’s temperature, while $\omega(t)$ is a Gaussian white noise with zero mean and unit variance parameters. In general, this model yields a dynamics, which is accurate up to the third-order statistics, due to symmetry considerations Belousov et al. (2016b); Belousov and Cohen (2016).

Experiments and molecular dynamics simulations, e.g. Refs. Belousov et al. (2016a, b), provide an evidence that, besides the thermal fluctuations, an additional source of spontaneous variations, $\epsilon$, is present in the nonequilibrium (NE) regime, $f\neq 0$,

[TABLE]

where $B\in\mathbb{R}$ and $\tau>0$ are constants, discussed later. Within the white noise approximation, essentially two models of this athermal contribution were suggested. One is formulated through a Poisson process Kanazawa et al. (2015a, b); Morgado and Queirós (2016), and another uses simple exponential noise Belousov et al. (2016b).

When the athermal fluctuations assume values in a continuous domain, rather than a discrete one Kanazawa et al. (2015a, b), the first of the above models is given by exponential shot noise; cf. Refs. Morgado and Queirós (2016); Belousov et al. (2016b). Physically, it can be interpreted as an external force, applied in the form of discrete impulses. Their magnitude is distributed exponentially with the scale parameter $B$, while their number, imparted per unit time, is determined by the Poisson law with the rate parameter $\tau^{-1}$; cf. Eq. (3). In contrast, the second model, based on simple exponential noise, describes a nonequilibrium force undergoing continuous variations in time; cf. LABEL:PRE2016II.

What seems unnoticed so far in the discussions of the Langevin dynamics with athermal noise, is that a nonequilibrium source of fluctuations in $r^{\mathrm{NE}}(t)$ does not exclude a possible presence of an additional deterministic constant term, $F$, on the right hand side of Eq. (1), which then reads

[TABLE]

For both models of athermal noise $\epsilon(t/\tau)$ mentioned above, the macroscopic limit of the right-hand side of Eq. (4) yields $f=F+B/\tau$; cf. Eq. (2) and Refs. Belousov et al. (2016b); Belousov and Cohen (2016). The necessity of the constant term $F$ in Eq. (4) is demonstrated in Sec. II, where we study in detail nonequilibrium aspects of the fluctuations of $\alpha(t)$.

The purpose of this paper is to demonstrate a complete Langevin representation of a nonequilibrium physical system. To do this, in Sec. II we first determine all the parameters of Eq. (4) for the microscopic model of a shear flow in a simple fluid, which was already studied in LABEL:PRE2016III. Then, in Sec. III we perform simulations of the Langevin dynamics and compare their results with the original system.

In our computational study, the two models of athermal noise, discussed earlier, are treated separately. They generate qualitatively very similar trajectories of $\alpha(t)$, as observed in stochastic simulations of Sec. III. Although exponential shot noise Morgado and Queirós (2016) has discrete singularities, which make its physical interpretation distinct from simple exponential noise Belousov et al. (2016b), in the second-order Langevin dynamics this affects only $\dot{\alpha}(t)$.

In addition, Appendix A reveals a mathematical analogy between the models of athermal noise, considered here. In particular, they can be regarded as Padé approximants for certain generating functions of the exact probability distribution, associated with $\epsilon(t/\tau)$. The mathematical structure, described in Appendix A, is applicable also to other families of approximations. This, in principle, allows formulation of alternative models for athermal noise.

Finally, for simulations of the Langevin dynamics, Appendix C proposes an algorithm, which is used to integrate Eq. (4) in Sec. III. Our numerical scheme minimizes sampling of random variables, which is usually a most intensive part of the computations. In this regard, the Langevin equation with simple exponential noise has an advantage over exponential shot noise, because it requires only one random number generation per step of numerical integration; cf. Appendix C. In the macroscopic limit our simulation method coincides with a second-order symplectic algorithm, which respects the time reversibility of Eq. (2).

II Molecular Dynamics Simulations

In this section we continue to study the same Molecular Dynamics (MD) model of a shear flow in a simple fluid, as we did in LABEL:PRE2016III. Here we briefly summarize, that we consider $N$ particles of equal mass $m$, interacting through the Weeks-Chandler-Anderson (WCA) potential Weeks et al. (1971) in three dimensions. A constant shear rate $\gamma$, which is the nonequilibrium force driving a current of linear momentum, is applied by the SLLOD equations of motion with Lees-Edwards boundary conditions (Evans and Morriss, 2007, Chapter 6). A constant temperature was maintained by the Nosé-Hoover (NH) thermostat (Frenkel and Smit, 2002, Chapter 6).

The current in our MD system is given by the $yx$-pressure tensor component, $P_{yx}$, which expresses a flow of $x$ momentum in the $y$-direction (Evans and Morriss, 2007, Sec. 3.8):

[TABLE]

where $V$ is the volume of the MD simulation cell, while, for the $i$-th particle, $p_{yi}$, $p_{xi}$, $y_{i}$, and $F_{xi}$ are, respectively, the $x$ and $y$ components of peculiar linear momentum, $y$-coordinate, and $x$-component of the force due to interactions with the other particles.

Data and results of our simulations are reported in the reduced units (r.u.), described in (Belousov and Cohen, 2016, Appendix C). Here we are interested in a nonequilibrium Langevin equation (4) for $\alpha(t)=P_{yx}$ accurate up to the third-order moment, which is most reliably estimated for the steady-state probability distribution of $P_{yx}$ observed in systems of a small size; cf. LABEL:PRE2016. Therefore, we conduct our simulations for $N=125$ particles, although at the same number density $0.8$ r.u. and temperature $1$ r.u., as in LABEL:PRE2016III.

As shown in LABEL:PRE2016III, the constants $a$ and $b$ in Eq. (4) do not depend on $\gamma$ and can be estimated from measurements of the time autocorrelation function of $P_{yx}$. Also, since the temperature is fixed in our simulations by the NH thermostat, the parameter $A$ can be determined from the variance for $P_{yx}$, measured in the equilibrium simulations, i.e. $\kappa_{2}(P_{yx}|\gamma=0)$; cf. Sec. I. Hence in this section we are mainly concerned with the remaining parameters of Eq. (4), $F$, $B$ and $\tau$, which can be obtained, by fitting the first three cumulants of $P_{yx}$, $\kappa_{i}(P_{yx}),i=1,2,3$, in a nonequilibrium steady-state $\gamma\neq 0$, as follows.

Formulas for the steady-state cumulants of $\alpha(t)$, which evolves according to Eq. (4), can be derived, by using the method of Ref. (Belousov and Cohen, 2016, Appendix A). For this we need to consider a random variable $R(t)$, given by the time integral of the force terms in the Langevin dynamics:

[TABLE]

The steady-state cumulants of $\alpha(t)$ can be then expressed as

[TABLE]

cf. Ref. (Belousov and Cohen, 2016, Appendix A).

Equation (II) can be further expanded, once the form of athermal noise in Eq. (4), $B\epsilon(t/\tau)$, is specified. We discussed in Sec. I, that in this paper we consider separately two models, simple exponential noise ($\mathrm{ex}$) and exponential shot noise ($\mathrm{sh}$). From now on, to distinguish them we will use subscripts in the parameters of Eq. (4), respectively, $F_{\mathrm{ex}}+B_{\mathrm{ex}}\epsilon_{\mathrm{ex}}(t/\tau_{\mathrm{ex}})$ and $F_{\mathrm{sh}}+B_{\mathrm{sh}}\epsilon_{\mathrm{sh}}(t/\tau_{\mathrm{sh}})$. In the case of simple exponential noise, Eq. (II) thus yields:

[TABLE]

And similarly, for exponential shot noise, we get

[TABLE]

cf. Belousov et al. (2016b); Belousov and Cohen (2016). The constants $A$, $a$ and $b$ do not depend on the shear rate, and we also have

[TABLE]

cf. Sec. I. Equations (II)-(9) suggest then to define a nonequilibrium part of the second cumulant as

[TABLE]

Formulas (II) and (II) form a system of equations, which can be solved for the parameters of Eq. (4),

[TABLE]

whereas $F$ is calculated by the residual principle, $F=b^{2}\kappa_{1}(\alpha)-B/\tau$, with $B$ and $\tau$ being the parameters of the chosen athermal noise model.

In Sec. I we promised to show, that the parameter $F$ is essential in Eq. (4). In particular, it allows our Langevin model for $P_{yx}$ to achieve the statistical accuracy of third order. Consider the following quantity, which relates the first three cumulants of $P_{yx}$,

[TABLE]

Below we refer to $k(P_{yx})$ as the athermal cumulant ratio, because it captures effects of the athermal noise in Eq. (4) and vanishes in equilibrium systems together with $\chi(P_{yx})$; cf. Eq. (9).

As can be deduced from Eqs. (II) and (II), the athermal cumulant ratio, for both models of athermal noise considered in Eq. (4), has a general form

[TABLE]

where $\mathrm{const}_{1}$ and $\mathrm{const}_{2}$ are some constants, which do not depend on the shear rate. By substituting $F=0$ into Eq. (11), we see, that without a force term $F$ on the right-hand side of Eq. (4) in both cases, simple exponential noise and exponential shot noise, $k$ would be independent of $\gamma$.

We plot $k(P_{yx})$ as a function of $\gamma$ for our MD simulations in Fig. 1, where a nearly linear trend can be observed. In particular, this implies, that $\chi(P_{yx})$, which essentially expresses the athermal variance of $P_{yx}$, vanishes faster than the mean and asymmetry of its steady-state probability distribution, respectively, $\kappa_{1}(P_{yx})$ and $\kappa_{3}(P_{yx})$. Such behavior of athermal statistics is consistent with the observation, that the fluctuations of the nonequilibrium force become negligible in the near-equilibrium systems Belousov et al. (2016a), i.e. $\chi(P_{yx})\approx 0$ and $\kappa_{3}(P_{yx})\approx 0$, and Eq. (4) thus can be approximated by

[TABLE]

which agrees in the macroscopic limit with Eq. (2), cf. Sec. I.

In Eq. (11), the linear dependence of the athermal cumulant ratio on the shear rate can only arise due to the term proportional to $F$, which therefore can not vanish. In fact, Eq. (4) admits an alternative form, in which one of its parameters is related linearly to the athermal cumulant ratio, as shown in Appendix B.

As anticipated in Sec. I, Eq. (4), with $\epsilon$ represented by white exponential noise or white exponential shot noise, is accurate up to the third-order statistics. For this, three parameters, $F$, $B$, and $\tau$, are required to fit exactly three cumulants $\kappa_{1}$, $\kappa_{2}$ and $\kappa_{3}$; cf. Eqs. (II)-(II).

III Langevin Dynamics Simulation

Once all the parameters of the stochastic differential equation (4) are found, it can be simulated numerically. In Appendix C we design an integration algorithm, which is consistent with the time reversibility of the Langevin Dynamics (LD) in the macroscopic limit, Eq. (2). Below we report the results of simulations, obtained with this numerical scheme.

As described in Sec. II, we determined the parameters of Eq. (4) for the shear current $P_{yx}$ and the two models of athermal noise, $\epsilon_{\mathrm{ex}}(t)$ and $\epsilon_{\mathrm{sh}}(t)$, from our MD calculations, obtained at $\gamma=1$ r.u. Table 1 summarizes results of our simulations. Statistics of the original MD data are reproduced very accurately and precisely by our numerical LD model up to the third-order statistics.

A notable discrepancy between the simulations is observed in Table 1 only for the excess kurtosis. This should be expected, because the parameters of Eq. (4) were obtained by fitting solely the first three cumulants of the MD data, and the two LD models, considered here, differ by their fourth cumulants in principle. Such level of precision, however, is beyond the statistical and numerical resolution of our simulations. Indeed, the magnitude of the excess kurtosis for the MD data is comparable to its standard deviation. Its theoretical value in our LD models is $17\cdot 10^{-5}$ for the case of simple exponential noise, and $12\cdot 10^{-5}$ for the case of exponential shot noise. These predictions of the excess kurtosis are also much less than the respective standard deviations in Table 1. Therefore, within the statistical uncertainties, the excess kurtosis of our MD and LD models can effectively be considered zero, i.e. Gaussian-like.

To integrate numerically Eq. (4), we used the same time step, $\Delta{t}=10^{-3}$, as in the original MD simulations; cf. LABEL:PRE2016III. Therefore, we can compare sample trajectories, which are traced pointwise in time for our LD models and the original MD system in Fig. 2. Both, simple exponential noise and exponential shot noise, produce a qualitatively similar behavior of $\alpha(t)$. In particular, the two stochastic models match closely the amplitude of $\alpha(t)$ fluctuations, as well as the time scale of their onset and decay, observed in the MD simulations.

The trajectory of the original MD system in Fig. 2 mostly resembles those of the LD simulations, perhaps, except for one subtle detail. Namely, on short time intervals, the stochastic dynamics of the time derivative $\dot{\alpha}(t)$ may undergo quick alterations, which makes the trajectory of $\alpha(t)$ look noisy. This is especially conspicuous for $\dot{\alpha}(t)\approx 0$, where this behavior leads to quick alterations of sign $\dot{\alpha}(t)\gtrless 0$ with a nearly constant value of $\alpha(t)$, like in Fig. 2(a) around $0.35<t<0.45$ r.u. or in Fig. 2(b) $0.50<t<0.55$. In contrast, the trajectory of our MD simulation appears more resistant to changes of its direction, $\dot{\alpha}(t)$, and thus slightly smoother.

The alterations of $\dot{\alpha}(t)$, mentioned above, might be an artifact of the white noise approximation. In principle, a colored noise, which introduces additional correlations into the stochastic forces on the right-hand side of Eq. (4), e.g. Duarte Queirós (2016), could make the dynamics of $\dot{\alpha}(t)$ more inertial. Possibly this would make the trajectory of $\alpha(t)$ smoother, by reducing the chance of quick sign alterations in its time derivative.

IV Conclusion

In this paper we have demonstrated, how efficiently the Langevin dynamics can represent certain nonequilibrium systems in a steady state. Section II describes a way to find the parameters of such representation from the statistics of the original system. The numerical simulations of Eq. (4), thus obtained, and the integration algorithm, proposed in Appendix C, yield results, which are both, statistically and qualitatively, very accurate and precise; cf. Sec. III.

Two models of athermal noise for Eq. (4), i.e. simple white exponential noise and white exponential shot noise, were compared in Sec. III. Within the third-order statistics they provide equivalent results. As explained in Appendix A, these two models belong to two, among many other, possible families of approximants for certain generating functions of the true probability distribution, associated with athermal noise $\epsilon(t)$. Therefore, from the statistical and mathematical perspectives, they are equivalent and interchangeable. The difference of the physical interpretations between the considered models of athermal noise, which might have distinct consequences for the first-order Langevin dynamics Machlup and Onsager (1953); Belousov et al. (2016b), in the second-order equation (4) applies mostly to $\dot{\alpha}(t)$, and does not affect significantly the behavior of $\alpha(t)$.

The above models of athermal noise differ by their values of the fourth cumulant. Within the statistical uncertainties of the original system considered in Sec. II, this discrepancy is, however, completely negligible; cf Sec. III. Statistically the excess kurtosis of our MD and LD data for $\alpha(t)$ is comparable to that of a Gaussian probability distribution, which equals zero.

To characterize the third-order statistics of athermal noise, we introduced a new quantity, athermal cumulant ratio, which has a nearly linear dependence on the external shear rate. This observation was analyzed in Sec. II, to justify the constant term $F$ in Eq. (4). Hence, the second-order Langevin representation of our original model needs, in total, specification of six constants, $a$, $b$, $A$, $B$, $\tau$, and $F$. Furthermore, Eq. (4) can be reparametrized so, that the athermal cumulant ratio is linearly related to one of the new parameters; see Appendix B.

Finally, at short-time scales we observed in Sec. III some features of the LD simulations, which appear to be artifacts of the white noise approximation. As discussed there, a correlated noise could, in principle, improve this aspect of the Langevin equation. A description of noise correlations may require to add another parameter to the model. This, however, seems far from being practical, since there would be more unknowns than statistically significant measurements to fit; cf. Secs. II-III. Instead, one could try to construct a first-order Langevin equation with a correlated noise, which would reproduce the time autocorrelation function for $\alpha(t)$ in Eq. (4). Indeed, the term $a\dot{\alpha}(t)$ of the second-order dynamics is required for the accurate description of time autocorrelations in the modeled system; cf. LABEL:PRE2016III. This introduces essentially a nuisance parameter $a$, if we are interested merely in the behavior of $\alpha(t)$. One may hope, however, that the time correlation function of a first-order Langevin dynamics, with a properly designed colored noise, would match that of $\alpha(t)$ in Eq. (4).

Appendix A Models of athermal noise

In LABEL:PRE2016II we have presented a derivation of exponential white noise for the Langevin dynamics of nonequilibrium physical systems, as an alternative to white exponential shot noise Morgado and Queirós (2016). As mentioned in Sec. I, these two models essentially account for the spontaneous variations of a nonequilibrium force, which are assumed independent from the thermal fluctuations and, therefore, called athermal Morgado and Queirós (2016).

The derivation of exponential white noise Belousov et al. (2016b) followed a procedure similar to that of the Gaussian white noise in Ref. (Chandrasekhar, 1943, Sec. I.1), used to model thermal fluctuations in both, equilibrium and nonequilibrium, systems. This procedure suggests to design a discrete random walk, which incorporates physically relevant assumptions and properties of the fluctuating force. Then, a definition of stochastic noise, which can be used in a Langevin equation, naturally arises in a properly chosen continuous limit.

In this section we revisit the above procedure and show, that it relies on a more general mathematical structure of asymptotic approximations, which allows to derive various stochastic noises with desired statistical properties. In fact both, simple exponential noise and exponential shot noise, belong to two families of random processes, which can be obtained in this manner.

Let us begin with a review of the Gaussian white noise $\omega(t)$, which is defined as a time derivative of the Gaussian random process $\Omega(t)$ with the zero mean and the variance $A^{2}t$, so that

[TABLE]

where $A>0$ is a constant.

In the derivation of Gaussian white noise we consider a random walk, where the walker moves in a discrete series of steps from an initial position $x_{0}=0$. The length of each step $\xi\approx A\omega(t)dt$ is an independent random variable with a given probability density $p(\xi)$. A physically relevant property, which we want to incorporate in this model, is a spatial symmetry for the forward ($\xi>0$) and backward ($\xi<0$) displacements. Statistically this means $p(\xi)=p(-\xi)$. The walker’s position $x$ after $n$ steps is then given by a sum of $n$ independent realizations of $\xi$. Therefore, the moment-generating function of the random variable $x$ equals the moment-generating function of $\xi$ raised to the power $n$, and thus

[TABLE]

where $\tilde{\xi}$ and $\mu(\tilde{\xi})$ are, respectively, the reciprocal dual variable111 The moment-generating function, which is the Laplace transform of a probability density, i.e. $\langle\exp(\xi\tilde{\xi})\rangle$, depends on a variable $\tilde{\xi}$, a reciprocal dual of $\xi$. of $\xi$ and its moment-generating function, and likewise for $x$, while $\kappa(\cdot)$ stands for the cumulant-generating function. From now on we consistently denote duals of random variables by the tilde.

By expanding the cumulant-generating function of $\xi$ in Taylor series at $\tilde{\xi}=0$,

[TABLE]

we use the fact that, for a symmetric random variable, all cumulants of odd orders $\kappa_{1}(\xi)$, $\kappa_{3}(\xi)$, and so on, must vanish, whereas the normalization of a probability density function always requires $\kappa_{0}(\cdot)=0$. The Taylor series Eq. (14), truncated at the second term is a cumulant-generating function of a Gaussian probability distribution

[TABLE]

which can be used to approximate $\kappa(\tilde{\xi})\approx\kappa_{G}(\tilde{\xi})$ up to the third order statistics; cf. Eq. (14). Hence, by virtue of Eq. (13) we further obtain,

[TABLE]

Not all orders of truncated power series are valid cumulant-generating functions. The second order approximation, however, always corresponds to a Gaussian probability distribution.

Notice that we have not invoked the Central Limit Theorem so far. Now we do this, in order to ensure that Eq. (15) holds asymptotically in the limit of large $n$, at least for fluctuations of size $\sqrt{n}$. At last, the continuous limit of the described random walk is taken, by requiring that, in an infinitesimal time $dt$, the walker makes a number of steps $dn$ with a rate $\tau=dt/dn$, so that Eq. (15) becomes

[TABLE]

which is the cumulant-generating function for the process $\Omega(t)$ with $A^{2}=\kappa_{2}(\xi)/\tau$.

Instead of simply reviewing exponential noise, however, we now derive a more general family of random processes, to which it pertains. To do this, we again consider a random walk problem $\xi\approx B\epsilon(t/\tau)dt$, similar to the one described above. Due to the physical assumptions discussed in LABEL:PRE2016II, this time we require $\kappa_{3}(\xi)>0$ and that the probability density $p(\xi)$ vanishes for $\xi<0$, which implies $\kappa_{1}(\xi)\geq 0$.

Earlier we used power series to approximate the cumulant-generating function for the symmetric random walk. This time we will use a Padé approximant Baker Jr. (1975) for the moment-generating function $\mu(\tilde{\xi})$.

Just like with power series, not all orders of the numerator and denominator, $i$ and $j$, respectively, in a Padé approximant $\mu_{[i/j]}(\tilde{\xi})\approx\mu(\tilde{\xi})$, are valid moment-generating functions. We set $i=j=1$, which is accurate at least up to the second-order statistics, for $i+j=2$, and has a general form (Baker Jr., 1975, Chapter 1)

[TABLE]

where $a_{1,1}$, $a_{1,2}$, etc. are constants. The normalization of the probability distribution, though, requires that $\mu_{[1/1]}(0)=1$ and, thus, $a_{2,1}=a_{1,1}$. With $b_{1}=a_{1,2}/a_{1,1}$ and $b_{2}=a_{2,1}/a_{1,1}$, Eq. (17) becomes

[TABLE]

where in the last equality we used another substitution $q=b_{1}/b_{2}$ and $B=-b_{2}$. The probability density, which corresponds to a moment-generating function $\mu_{[1/1]}(\tilde{\xi})$, is

[TABLE]

where $\delta(\cdot)$ and $H(\cdot)$ are, respectively, the Dirac delta and Heaviside step functions. The above equation is a valid probability density function $p_{[1/1]}(\xi)>0$ only for $0\leq q\leq 1$. Furthermore, our assumptions $\kappa_{1}(\xi)>0$ and $\kappa_{3}(\xi)>0$ impose another restriction $B>0$. The continuous limit of the described approximation can be obtained by requiring, that the walker makes steps with a rate $\tau=dt/dn$, just like we did earlier in the derivation of $\omega(t)$. The result is a defintion of random noise $\epsilon_{q}(t/\tau)$, which upon the time integration generates the following stochastic process

[TABLE]

with the cumulant-generating function and the probability density, respectively,

[TABLE]

where $p_{\Gamma}(\cdot;i,B)$ is a probability density of the Gamma distribution with the shape parameter $i$ and the scale $B$.

The stochastic noise $\epsilon_{q}(t/\tau)$ reduces to simple exponential noise when $q=0$, as it does under the conditions of random walk problem assumed in LABEL:PRE2016II. The general case of $q\neq 0$ introduces one more parameter into Eq. (4), which is unpractical within the statistical precision of our model; cf. Secs. II and III.

A shot noise approximation of the athermal noise Eq. (A) can be derived in two ways. First, we may alter slightly the interpretation of a continuous limit for the random walk. Above the walker attempts $dn=dt/\tau$ steps in an infinitesimal time $dt$ deterministically. Instead, we could require, that the walker attempts a step with a probability proportional to a time interval $dt$, so small that it can make at most one displacement forward with the probability $dt/\tau^{\prime}$, or to remain motionless with a probability $1-dt/\tau^{\prime}$. Given that the length of each step is the random variate $\xi$, this sheme generates a compound Poisson process, $P_{q}(t)$; cf. (Kardar, 2007, Sec. 2.3). Its cumulant-generating function is

[TABLE]

which together with the approximant given by Eq. (18) yields

[TABLE]

The dimensionless parameter $(1-q)$ can be further absorbed into the constant $\tau=\tau^{\prime}/(1-q)$, and we thus obtain the compound Poisson process $P(t)$ with exponentially distributed intensity

[TABLE]

cf. LABEL:PRE2016II.

Another way to derive Eq. (23) is to seek directly a Padé approximant for the cumulant generating function $\kappa(\tilde{\xi})\approx\kappa_{[1/1]}(\tilde{\xi})$. Like for $\mu_{[1/1]}$ in Eq. (17), we begin with a most general form

[TABLE]

where $a_{1,1}$, $a_{1,2}$, etc. are constants. Since a cumulant-generating function must vanish at $\tilde{\xi}=0$, as required by the normalization of probability densities, we then have $a_{1,1}=0$. With $B=-a_{2,2}/a_{2,1}$ and $a_{2}=-a_{1,2}/(Ba_{2,2})$, Eq. (24) becomes

[TABLE]

Now we can take a continuous limit in the manner we did for the Gaussian and simple exponential noise, to obtain a process $P^{\prime}(t)$

[TABLE]

where again the dimensionless constant $a_{2}$ can be absorbed into $\tau=\tau^{\prime}/a_{2}$ to yield Eq. (23).

Finally, we would like to remark, that the Central Limit Theorem ensures that the above approximations hold asymptotically. Indeed, all the random walk problems considered here tend to a Gaussian probability distribution for $n\to\infty$. Our Padé approximants have the same asymptotic limit, since they are accurate up to the second-order statistics.

Above we demonstrated, how the Padé approximation can be applied to formulate a random noise with presupposed statistical properties. In principle, however, other forms of approximants could be used in exactly the same way. For example, LABEL:DanchenkoChunaev2011 suggests a general approach, which allows to generate various mathematical series, which fit the Taylor expansion of the approximated function. In particular, one of these series can be generated from the simple fraction $(1-x)^{-1}$ Danchenko and Chunaev (2011). Under certain conditions, this yields a hyperexponential approximant for the moment-generating function $\mu(\tilde{\xi})$ Feldmann and Whitt (1998). In other words, the Central Limit Theorem, as used in this section, and the approach of LABEL:DanchenkoChunaev2011 give rise to a multitude of stochastic processes, which can fit various dynamical models.

Appendix B Athermal cumulant ratio revisited

In this section of Appendix we would like to show an alternative parameterization of Eq. (4), which provides a particularly simple expression for the athermal cumulant ratio introduced in Sec. II. We again consider the random walk problem, formulated in Appendix A for athermal noise. This time, as an approximant of the moment-generating function, we would like to use the following expression

[TABLE]

which is a moment-generating function of the exponentially distributed random variate, shifted by a location parameter $(1-q)\phi$, and corresponds to the probability density

[TABLE]

with $f\geq 0$ and $0\leq q\leq 1$. The Taylor expansion of the cumulant-generating function for the this distribution is

[TABLE]

in which we always can choose the two parameters $\phi=\kappa_{1}(\xi)$ and $q=\sqrt{\kappa_{2}(\xi)/\phi^{2}}$, so that

[TABLE]

This equation is inspired by the approach of LABEL:DanchenkoChunaev2011, which is based on matching the Taylor expansion coefficients of the approximated function and the approximant; cf. Appendix apx.

In the continuous limit, with the number of steps per unit time $\tau^{-1}=dn/dt$, Eq. (28) yields a cumulant-generating function of a process $E_{q,\phi}(t)$,

[TABLE]

If we put $F=(1-q)\phi/\tau$ and $B=q\phi$, the process $E_{q,\phi}(t)$ above can be expressed as

[TABLE]

where $\epsilon(t/\tau)$ is simple exponential noise; cf. Sec. II. By using Eq. (29), Eq. (4) can be written as

[TABLE]

In this new parameterization of the Langevin equation, which replaces constants $B$ and $F$ by $q$ and $\phi$, the athermal cumulant ratio, introduced in Sec. II, has a quite simple form

[TABLE]

An alternative parameterization, similar to Eq. (29), can be obtained also for the Langevin dynamics with exponential shot noise. For this, one needs merely to introduce a location parameter into the compound Poisson process, considered in Appendix A, like we did in Eq. (27). The constant of proportionality in Eq. (30), for simple exponential noise, $\mathrm{const}_{\mathrm{ex}}$, and exponential shot noise, $\mathrm{const}_{\mathrm{sh}}$, are

[TABLE]

Appendix C A macroscopically symplectic algorithm for simulations of Langevin Dynamics

The differential equation (4) is mathematically equivalent to the following system:

[TABLE]

The macroscopic limit of Eq. (32) is

[TABLE]

cf. Eqs. (1)-(4) in Sec. I. Equation (33) is time-reversible, like the equations of motion in our MD model of Sec. II.

In order to design an algorithm for the LD simulations, which respects the macroscopic time reversibility, consider first the deterministic equation (33). The phase space of this system, as well as of Eq. (32), is two-dimensional, $\bm{\Gamma}=(\alpha\,\beta)$. We are looking for a numerical scheme of the second order in time, which provides the accuracy level of our MD simulations in Sec. II and LABEL:PRE2016III, by using the operator splitting formalism Tuckerman et al. (1992); Martyna and Tuckerman (1995); Martyna et al. (1996); cf. (Belousov and Cohen, 2016, Appendix C). Let us rewrite Eq. (32) as

[TABLE]

with the Liouvillian (Evans and Morriss, 2007, Chapter 3)

[TABLE]

which is a sum of three mutually non-commuting terms

[TABLE]

A formal solution of Eq. (34) is

[TABLE]

If we neglect for a moment the dissipative part of the dynamics in Eq. (36), $\mathrm{i}\mathcal{L}_{\beta\beta}$, we are left with two possible choices of a symplectic algorithm, described in LABEL:TuckermanMartyna, namely, a widely used velocity Verlet scheme,

[TABLE]

and a perhaps less known position Verlet scheme

[TABLE]

where $\Delta{t}$ is a time step of numerical integration.

Now we recall, that $\mathrm{i}\mathcal{L}_{\alpha\beta}$ contains a force term $f$, which is the macroscopic limit of the right-hand side in Eq. (4). Therefore, to obtain a mesoscopic version of Eqs. (37) and (38), we make a substitution

[TABLE]

where we used Eq. (5).

In the Verlet algorithms above, Eq. (39) affects only an operator of the form $\exp(\mathrm{i}\mathcal{L}_{\alpha\beta}\Delta{t})$; cf Eq. (35). As explained shortly, its implementation involves sampling of random variables, which is usually the most computationally expansive part of stochastic simulations. Note that, on the right-hand side of Eq. (37), this operator is applied two times, in the form $\exp(\mathrm{i}\mathcal{L}_{\alpha\beta}\Delta{t}/2)$, while Eq. (38) uses it only once. Therefore we continue with the latter integration scheme, which is computationally more efficient.

By following the approach of Refs. Martyna and Tuckerman (1995); Martyna et al. (1996), the dissipative part of the dynamics, due to $\mathrm{i}\mathcal{L}_{\beta\beta}$, which was until now disregarded by Eq. (38), can be numerically implemented in two qualitatively equivalent ways. The one, which we chose for the simulations in Sec. III, is given by222 If, everywhere in Eq. (40), we exchange the order, in which the operators $\exp(\mathrm{i}\mathcal{L}_{\alpha}\Delta{t}/2)$ and $\exp(\mathrm{i}\mathcal{L}_{\alpha\beta}\Delta{t}/2)$ are applied, we would obtain the second implementation of the dissipative dynamics for the position Verlet scheme.

[TABLE]

The exact steps of the simulation algorithm can now be read from Eq. (40), as explained in detail in Ref. (Frenkel and Smit, 2002, Appendix E). On the right-hand side, we first operate on $\bm{\Gamma}(t)$ with $\exp(\mathrm{i}\mathcal{L}_{\alpha}\Delta{t}/2)$, which generates a translation of $\alpha$ by $\beta\Delta{t}/2$ and leaves intact $\beta$. Then, from right to left, $\exp(\mathrm{i}\mathcal{L}_{\beta\beta}\Delta{t}/2)$ scales $\beta$ by $\exp(-a\Delta{t}/2)$, while $\exp(\mathrm{i}\mathcal{L}_{\alpha\beta}\Delta{t})$ increments it by $b^{2}\alpha\Delta{t}+R(\Delta{t})$, and so on. In a compact form, a complete step of the simulation, $\Delta{t}$, can be written as

[TABLE]

Note, that the above algorithm requires generation of the random variable $R(\Delta{t})$ only in the second intermediate step, Eq. (42). When Eq. (4) contains exponential noise, the athermal part of $R(\Delta{t})$ is sampled once from a Gamma probability distribution Belousov et al. (2016b). In the case of exponential shot noise, first, one has to generate a random integer from the Poisson distribution. This number then determines how many exponentially distributed terms must be sampled to evaluate $R(\Delta{t})$. For the whole duration of a LD simulation, on average, exponential shot noise requires generation of $1+\Delta{t}/\tau$ athermal random variates per integration step $\Delta{t}$, while solely one is needed for simple exponential noise.