---
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.'
author:
- Roman Belousov
- 'E.G.D. Cohen'
- Lamberto Rondoni
bibliography:
- 'References.bib'
title: 'Nonequilibrium Langevin dynamics: a demonstration study of shear flow fluctuations in a simple fluid'
---

\[sec:intro\]Introduction
=========================

The Langevin dynamics, originally inspired by the problem of Brownian motion [@Coffey2012Langevin Chapters 1-2], has now become the fundamental stochastic model of fluctuations for equilibrium physical systems at mesoscopic scales [@Nyquist1928; @OM1953; @MO1953]. Its generalization to nonequilibrium steady states, though, is still actively developed, as suggested by a number of recent publications [@KSSH2015I; @KSSH2015II; @MorgadoQ2016; @Queiros2016; @PRE2016; @PRE2016II; @PRE2016III]. 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 [@PRE2016III] we showed, that a second-order Langevin equation, suggested in Ref. [@MO1953], 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: $${\label{eq:lgv1}}
    \ddot\alpha(t) + a \dot\alpha(t) + b^2 \alpha(t) = r(t),$$ where $a > 0$, $b > 0$ are constants, while $r(t)$ is a random noise.

In the macroscopic limit, [Eq. (\[eq:lgv1\])]{} transforms into a deterministic equation [@OM1953 Sec. 2-3] and, therefore, $r(t)$ must then in general become a constant, $f \in \mathbb{R}$, so that we have $${\label{eq:mac1}}
    \ddot\alpha(t) + a \dot\alpha(t) + b^2 \alpha(t) = f,$$ 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. (\[eq:lgv1\])]{}, e.g. [@PRE2016III], with ${\langle\alpha(t)\rangle} = f/b^2$ and ${\langle\dot\alpha(t)\rangle} = 0$. By comparing Eqs. (\[eq:lgv1\]) and (\[eq:mac1\]), one can see that, at the mesoscopic scales, $r(t)$ represents the external force, $f = {\langler(t)\rangle}$, as well as its spontaneous variations [@OM1953].

In the equilibrium regime [@Nyquist1928; @OM1953; @PRE2016II], $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 [@PRE2016II; @PRE2016III].

Experiments and molecular dynamics simulations, e.g. Refs. [@PRE2016; @PRE2016II], provide an evidence that, besides the thermal fluctuations, an additional source of spontaneous variations, $\epsilon$, is present in the nonequilibrium (NE) regime, $f\ne0$, $${\label{eq:rNE}}r^\mathrm{NE}(t) = A \omega(t) + B \epsilon(t/\tau),$$ 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 [@KSSH2015I; @KSSH2015II; @MorgadoQ2016], and another uses simple exponential noise [@PRE2016II].

When the athermal fluctuations assume values in a continuous domain, rather than a discrete one [@KSSH2015I; @KSSH2015II], the first of the above models is given by exponential shot noise; cf. Refs. [@MorgadoQ2016; @PRE2016II]. 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. (\[eq:rNE\])]{}. In contrast, the second model, based on simple exponential noise, describes a nonequilibrium force undergoing continuous variations in time; cf. [Ref. [@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. (\[eq:lgv1\])]{}, which then reads $${\label{eq:lgv2}}
    \ddot\alpha(t) + a \dot\alpha(t) + b^2 \alpha(t) = F + A\omega(t) + B\epsilon(t/\tau).$$ For both models of athermal noise $\epsilon(t/\tau)$ mentioned above, the macroscopic limit of the right-hand side of [Eq. (\[eq:lgv2\])]{} yields $f = F + B/\tau$; cf. [Eq. (\[eq:mac1\])]{} and Refs. [@PRE2016II; @PRE2016III]. The necessity of the constant term $F$ in [Eq. (\[eq:lgv2\])]{} is demonstrated in [Sec. \[sec:mic\]]{}, 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. \[sec:mic\]]{} we first determine all the parameters of [Eq. (\[eq:lgv2\])]{} for the microscopic model of a shear flow in a simple fluid, which was already studied in [Ref. [@PRE2016III]]{}. Then, in [Sec. \[sec:mes\]]{} 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. \[sec:mes\]]{}. Although exponential shot noise [@MorgadoQ2016] has discrete singularities, which make its physical interpretation distinct from simple exponential noise [@PRE2016II], in the second-order Langevin dynamics this affects only $\dot\alpha(t)$.

In addition, Appendix \[sec:apx\] 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 \[sec:apx\], 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 \[sec:alg\] proposes an algorithm, which is used to integrate [Eq. (\[eq:lgv2\])]{} in [Sec. \[sec:mes\]]{}. 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 \[sec:alg\]. In the macroscopic limit our simulation method coincides with a second-order symplectic algorithm, which respects the time reversibility of [Eq. (\[eq:mac1\])]{}.

\[sec:mic\]Molecular Dynamics Simulations
=========================================

![\[fig:kfig\] Plot of the athermal cumulant ratio, $k(P_{yx})$, as a function of the shear rate, $\gamma$, for our MD simulations. Error bars are given by one standard deviation. A nearly linear trend can be fit to a line, $k(P_{yx}) = {\mathrm{const}}\gamma$. ](kfig){width="1\columnwidth"}

                                                                         LD(ex)               LD(sh)                   MD
  -------------------------------------------------------- -------------------- -------------------- --------------------
  Mean, $\kappa_1(\alpha)$                                   $-1.515 \pm 0.004$   $-1.518 \pm 0.004$   $-1.521 \pm 0.004$
  Variance, $\kappa_2(\alpha)$                               $ 0.185 \pm 0.003$   $ 0.185 \pm 0.003$   $ 0.184 \pm 0.003$
  Skewness, $\kappa_3(\alpha)/\kappa_2(\alpha)^{3/2}$        $ -0.11 \pm  0.03$   $ -0.16 \pm  0.03$   $ -0.14 \pm  0.03$
  Excess kurtosis, $\kappa_4(\alpha)/\kappa_2(\alpha)^2$     $  0.15 \pm  0.09$   $  0.19 \pm  0.07$   $  0.08 \pm  0.07$

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 [Ref. [@PRE2016III]]{}. Here we briefly summarize, that we consider $N$ particles of equal mass $m$, interacting through the Weeks-Chandler-Anderson (WCA) potential [@WCA] 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 [@EvansMorriss Chapter 6]. A constant temperature was maintained by the Nosé-Hoover (NH) thermostat [@FrenkelSmit 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 [@EvansMorriss Sec. 3.8]: $$P_{yx} = V^{-1} \sum_{i=1}^{N} (p_{yi} p_{xi} /m + y_i F_{xi}),$$ 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 [@PRE2016III Appendix C]. Here we are interested in a nonequilibrium Langevin equation (\[eq:lgv2\]) 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. [Ref. [@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 [Ref. [@PRE2016III]]{}.

As shown in [Ref. [@PRE2016III]]{}, the constants $a$ and $b$ in [Eq. (\[eq:lgv2\])]{} 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. \[sec:intro\]]{}. Hence in this section we are mainly concerned with the remaining parameters of [Eq. (\[eq:lgv2\])]{}, $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\ne0$, as follows.

Formulas for the steady-state cumulants of $\alpha(t)$, which evolves according to [Eq. (\[eq:lgv2\])]{}, can be derived, by using the method of Ref. [@PRE2016III 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: $${\label{eq:Rt}} R(t) = \int_0^{t} d{t^\prime} [F + A \omega({t^\prime}) + B \epsilon({t^\prime})].$$ The steady-state cumulants of $\alpha(t)$ can be then expressed as $$\begin{aligned}
{\label{eq:kappar}}
    \kappa_1(\alpha) &=& \frac{\kappa_1(R)}{b^2 t};\quad
    \kappa_2(\alpha) = \frac{\kappa_2(R)}{2 a b^2 t};\nonumber\\
    \kappa_3(\alpha) &=& \frac{2 \kappa_3(R)}{3 b^2 (2 a^2 + b^2) t}.\end{aligned}$$ cf. Ref. [@PRE2016III Appendix A].

Equation (\[eq:kappar\]) can be further expanded, once the form of athermal noise in [Eq. (\[eq:lgv2\])]{}, $B \epsilon(t/\tau)$, is specified. We discussed in [Sec. \[sec:intro\]]{}, 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. (\[eq:lgv2\])]{}, 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. (\[eq:kappar\])]{} thus yields: $$\begin{aligned}
{\label{eq:kappaex}}
    \kappa_1(\alpha) &=& \frac{F_\mathrm{ex} + B_\mathrm{ex}/\tau_\mathrm{ex}}{b^2};\quad
    \kappa_2(\alpha) = \frac{A^2 + B_\mathrm{ex}^2/\tau_\mathrm{ex}}{2 a b^2};\nonumber\\
    \kappa_3(\alpha) &=& \frac{4 B_\mathrm{ex}^3}{3 b^2 (2 a^2 + b^2) \tau_\mathrm{ex}}.\end{aligned}$$ And similarly, for exponential shot noise, we get $$\begin{aligned}
{\label{eq:kappash}}
    \kappa_1(\alpha) &=& \frac{F_\mathrm{sh} + B_\mathrm{sh}/\tau_\mathrm{sh}}{b^2};\quad
    \kappa_2(\alpha) = \frac{A^2 + 2 B_\mathrm{sh}^2/\tau_\mathrm{sh}}{2 a b^2};\nonumber\\
    \kappa_3(\alpha) &=& \frac{4 B_\mathrm{sh}^3}{b^2 (2 a^2 + b^2) \tau_\mathrm{sh}}.\end{aligned}$$ cf. [@PRE2016II; @PRE2016III]. The constants $A$, $a$ and $b$ do not depend on the shear rate, and we also have $${\label{eq:eqkappa2}} \kappa_2(\alpha|\gamma=0) = \frac{A^2}{2 a b^2};$$ cf. [Sec. \[sec:intro\]]{}. Equations (\[eq:kappaex\])-(\[eq:eqkappa2\]) suggest then to define a nonequilibrium part of the second cumulant as $$\chi(\alpha) = \kappa_2(\alpha) - \frac{A^2}{2 a b^2}.$$

![image](trj){width="2\columnwidth"}

Formulas (\[eq:kappaex\]) and (\[eq:kappash\]) form a system of equations, which can be solved for the parameters of [Eq. (\[eq:lgv2\])]{}, $$\begin{aligned}
{\label{eq:kappas}}
    B_\mathrm{ex} &=& \frac{3 (2 a^2 + b^2) \kappa_3(\alpha)}{8 a \chi(\alpha)} = 3 B_\mathrm{sh} / 2;
    \nonumber\\
    \tau_\mathrm{ex} &=& \frac{9 (2 a^2 + b^2)^2 \kappa_3^2(\alpha)}{128 a^3 b^2 \chi^3(\alpha)} = 9 \tau_\mathrm{sh} / 8;\end{aligned}$$ 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. \[sec:intro\]]{} we promised to show, that the parameter $F$ is essential in [Eq. (\[eq:lgv2\])]{}. 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}$, $$k(P_{yx}) = \frac{\chi^2(P_{yx})}{\kappa_1(P_{yx}) \kappa_3(P_{yx})},$$ Below we refer to $k(P_{yx})$ as [*the athermal cumulant ratio*]{}, because it captures effects of the athermal noise in [Eq. (\[eq:lgv2\])]{} and vanishes in equilibrium systems together with $\chi(P_{yx})$; cf. [Eq. (\[eq:eqkappa2\])]{}.

As can be deduced from Eqs. (\[eq:kappaex\]) and (\[eq:kappash\]), the athermal cumulant ratio, for both models of athermal noise considered in [Eq. (\[eq:lgv2\])]{}, has a general form $${\label{eq:kform}}
    k \propto ({\mathrm{const}}_1 F \tau / B + {\mathrm{const}}_2)^{-1}{,}$$ 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. (\[eq:kform\])]{}, we see, that without a force term $F$ on the right-hand side of [Eq. (\[eq:lgv2\])]{} 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. \[fig:kfig\]]{}, 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 [@PRE2016], i.e. $\chi(P_{yx})\approx0$ and $\kappa_3(P_{yx})\approx0$, and [Eq. (\[eq:lgv2\])]{} thus can be approximated by $$\ddot\alpha(t) + a \dot\alpha(t) + b^2 \alpha(t) \approx f + A\omega(t),$$ which agrees in the macroscopic limit with [Eq. (\[eq:mac1\])]{}, cf. [Sec. \[sec:intro\]]{}.

In [Eq. (\[eq:kform\])]{}, 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. (\[eq:lgv2\])]{} admits [*an alternative form*]{}, in which one of its parameters is related linearly to the athermal cumulant ratio, as shown in Appendix \[sec:acr\].

As anticipated in [Sec. \[sec:intro\]]{}, [Eq. (\[eq:lgv2\])]{}, 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. (\[eq:kappaex\])-(\[eq:kappas\]).

\[sec:mes\]Langevin Dynamics Simulation
=======================================

Once all the parameters of the stochastic differential equation (\[eq:lgv2\]) are found, it can be simulated numerically. In Appendix \[sec:alg\] we design an integration algorithm, which is consistent with the time reversibility of the Langevin Dynamics (LD) in the macroscopic limit, [Eq. (\[eq:mac1\])]{}. Below we report the results of simulations, obtained with this numerical scheme.

As described in [Sec. \[sec:mic\]]{}, we determined the parameters of [Eq. (\[eq:lgv2\])]{} 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 \[tab:stat\]]{} 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 \[tab:stat\]]{} only for the excess kurtosis. This should be expected, because the parameters of [Eq. (\[eq:lgv2\])]{} 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\cdot10^{-5}$ for the case of simple exponential noise, and $12\cdot10^{-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 \[tab:stat\]]{}. 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. (\[eq:lgv2\])]{}, we used the same time step, $\Delta{t} = 10^{-3}$, as in the original MD simulations; cf. [Ref. [@PRE2016III]]{}. Therefore, we can compare sample trajectories, which are traced pointwise in time for our LD models and the original MD system in [Fig. \[fig:trj\]]{}. 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. \[fig:trj\]]{} 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)\approx0$, where this behavior leads to quick alterations of sign $\dot\alpha(t)\gtrless0$ with a nearly constant value of $\alpha(t)$, like in [Fig. \[fig:trj\]]{}(a) around $0.35<t<0.45$ r.u. or in [Fig. \[fig:trj\]]{}(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. (\[eq:lgv2\])]{}, e.g. [@Queiros2016], 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.

\[sec:final\]Conclusion
=======================

In this paper we have demonstrated, how efficiently the Langevin dynamics can represent certain nonequilibrium systems in a steady state. Section \[sec:mic\] describes a way to find the parameters of such representation from the statistics of the original system. The numerical simulations of [Eq. (\[eq:lgv2\])]{}, thus obtained, and the integration algorithm, proposed in Appendix \[sec:alg\], yield results, which are both, statistically and qualitatively, very accurate and precise; cf. [Sec. \[sec:mes\]]{}.

Two models of athermal noise for [Eq. (\[eq:lgv2\])]{}, i.e. simple white exponential noise and white exponential shot noise, were compared in [Sec. \[sec:mes\]]{}. Within the third-order statistics they provide equivalent results. As explained in Appendix \[sec:apx\], 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 [@MO1953; @PRE2016II], in the second-order equation (\[eq:lgv2\]) 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. \[sec:mic\]]{}, this discrepancy is, however, completely negligible; cf [Sec. \[sec:mes\]]{}. 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. \[sec:mic\]]{}, to justify the constant term $F$ in [Eq. (\[eq:lgv2\])]{}. 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. (\[eq:lgv2\])]{} can be reparametrized so, that the athermal cumulant ratio is linearly related to one of the new parameters; see Appendix \[sec:acr\].

Finally, at short-time scales we observed in [Sec. \[sec:mes\]]{} 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. \[sec:mic\]-\[sec:mes\]. 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. (\[eq:lgv2\])]{}. 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. [Ref. [@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. (\[eq:lgv2\])]{}.

\[sec:apx\]Models of athermal noise
===================================

In [Ref. [@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 [@MorgadoQ2016]. As mentioned in [Sec. \[sec:intro\]]{}, 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*]{} [@MorgadoQ2016].

The derivation of exponential white noise [@PRE2016II] followed a procedure similar to that of the Gaussian white noise in Ref. [@Chandrasekhar1943 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 $$\Omega(t) = \int_0^{t} d{t^\prime} A \omega({t^\prime})\text{,}$$ 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 $$\begin{aligned}
    {\label{eq:mu}} \mu(\tilde{x}) &=& [\mu(\tilde\xi)]^n\vline_{\tilde\xi=\tilde{x}};\\
    {\label{eq:kappa}} \kappa(\tilde{x}) &=& \ln \mu(\tilde{x}) = n \kappa(\tilde\xi)\vline_{\tilde\xi=\tilde{x}}\text{,}\end{aligned}$$ where $\tilde\xi$ and $\mu(\tilde\xi)$ are, respectively, the reciprocal dual variable[^1] 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$, $${\label{eq:Taylor}}
    \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}),$$ 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. (\[eq:Taylor\])]{}, truncated at the second term is a cumulant-generating function of a Gaussian probability distribution $$\kappa_G(\tilde\xi) = \kappa_2(\xi) \tilde\xi^2/2,$$ which can be used to approximate $\kappa(\tilde\xi) \approx \kappa_G(\tilde\xi)$ up to the third order statistics; cf. [Eq. (\[eq:Taylor\])]{}. Hence, by virtue of [Eq. (\[eq:kappa\])]{} we further obtain, $${\label{eq:gaussn}}\kappa(\tilde{x}) \approx n \kappa_2(\xi) \tilde{x}^2/2.$$ 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. (\[eq:gaussn\])]{} 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. (\[eq:gaussn\])]{} becomes $${\label{eq:gausst}}\kappa(\tilde{x}) \approx \frac{\kappa_2(\xi)  t \tilde{x}^2}{2\tau},$$ 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 [Ref. [@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)\ge0$.

Earlier we used power series to approximate the cumulant-generating function for the symmetric random walk. This time we will use a Padé approximant [@BakerEssentials] 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 [@BakerEssentials Chapter 1] $${\label{eq:MuPade1}}
    \mu_{[1/1]}(\tilde\xi) = \frac{a_{1,1} + a_{1,2} \tilde\xi}{a_{2,1} + a_{2,2} \tilde\xi},$$ 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. (\[eq:MuPade1\])]{} becomes $${\label{eq:MuPade2}}
    \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}.$$ 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 $$p_{[1/1]}(\xi) = q \delta(\xi) + (1-q) H[\operatorname{sign}(B) \xi] \exp(-\xi/B)/B,$$ 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 \le q \le 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 $$E_q(t/\tau) = \int_0^t d{t^\prime} B \epsilon_q(t/\tau),$$ with the cumulant-generating function and the probability density, respectively, $$\begin{aligned}
{\label{eq:exp}}
    \kappa(\tilde{E}_q) &=& \ln[q + (1-q)/(1-B \tilde{E}_q)] t/\tau,\nonumber\\
    p(\tilde{E}_q) &=& q^{t/\tau} \delta(E_q/B)\nonumber\\
        &+& \sum_{i=0}^{\infty} \binom{t/\tau}{i} q^{t/\tau-i}(1-q)^{i} p_\Gamma(E_q;i,B)\text{,}
        \nonumber\\\end{aligned}$$ 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 [Ref. [@PRE2016II]]{}. The general case of $q\ne0$ introduces one more parameter into [Eq. (\[eq:lgv2\])]{}, which is unpractical within the statistical precision of our model; cf. Secs. \[sec:mic\] and \[sec:mes\].

A shot noise approximation of the athermal noise [Eq. (\[eq:exp\])]{} 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. [@KardarI Sec. 2.3]. Its cumulant-generating function is $$\kappa(\tilde{P}_q) = \frac{[\mu(\tilde\xi)-1] t}{{\tau^\prime}}\vline_{\tilde\xi=\tilde{P}_q},$$ which together with the approximant given by [Eq. (\[eq:MuPade2\])]{} yields $$\kappa(\tilde{P}_q) \approx \frac{(1-q) B \tilde{P}_q t}{\tau(1-B \tilde{P}_q)}.$$ 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 $${\label{eq:pois}}
    \kappa(\tilde{P}) = \frac{B \tilde{P} t}{(1-B x) \tau},$$ cf. [Ref. [@PRE2016II]]{}.

Another way to derive [Eq. (\[eq:pois\])]{} 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. (\[eq:MuPade1\])]{}, we begin with a most general form $${\label{eq:KaPade1}}
    \kappa_{[1/1]}(\tilde\xi) = \frac{a_{1,1} + a_{1,2} \tilde\xi}{a_{2,1} + a_{2,2} \tilde\xi},$$ 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} / (B a_{2,2})$, [Eq. (\[eq:KaPade1\])]{} becomes $${\label{eq:KaPade2}}
    \kappa_{[1/1]}(\tilde\xi) = \frac{a_2 B \tilde\xi}{1 - B \tilde\xi}.$$

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)$ $$\kappa({\tilde{P}^\prime}) = \frac{a_2 B {\tilde{P}^\prime}t}{(1-B \tilde{P}) {\tau^\prime}},$$ where again the dimensionless constant $a_2$ can be absorbed into $\tau = {\tau^\prime}/a_2$ to yield [Eq. (\[eq:pois\])]{}.

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, [Ref. [@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}$ [@DanchenkoChunaev2011]. Under certain conditions, this yields a hyperexponential approximant for the moment-generating function $\mu(\tilde\xi)$ [@Hyperexp]. In other words, the Central Limit Theorem, as used in this section, and the approach of [Ref. [@DanchenkoChunaev2011]]{} give rise to a multitude of stochastic processes, which can fit various dynamical models.

\[sec:acr\]Athermal cumulant ratio revisited
============================================

In this section of Appendix we would like to show an alternative parameterization of [Eq. (\[eq:lgv2\])]{}, which provides a particularly simple expression for the athermal cumulant ratio introduced in [Sec. \[sec:mic\]]{}. We again consider the random walk problem, formulated in Appendix \[sec:apx\] for athermal noise. This time, as an approximant of the moment-generating function, we would like to use the following expression $${\label{eq:muqf}}\mu_{q,\phi}(\tilde\xi) = \exp[(1-q) \phi \tilde\xi]/(1-q \phi \tilde\xi),$$ 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 $$p_{q,\phi}(\xi) = \frac{H[\xi - (1-q) \phi]}{q\phi} \exp\left(-\frac{\xi}{q\phi}\right),$$ with $f\ge0$ and $0\le q \le 1$. The Taylor expansion of the cumulant-generating function for the this distribution is $$\kappa_{q,\phi}(\tilde\xi) = \phi\tilde\xi + q^2 \phi^2 \tilde\xi^2/2 + q^3 \phi^3 \tilde\xi^3/6 + O(\xi^4),$$ in which we always can choose the two parameters $\phi=\kappa_1(\xi)$ and $q=\sqrt{\kappa_2(\xi)/\phi^2}$, so that $${\label{eq:cumf}} \kappa(\tilde\xi) = \kappa_{q,\phi}(\tilde\xi) + O(\tilde\xi^3).$$ This equation is inspired by the approach of [Ref. [@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. (\[eq:cumf\])]{} yields a cumulant-generating function of a process $E_{q,\phi}(t)$, $$\kappa(\tilde{E}_{q,\phi}) = \frac{t (1-q) \phi \tilde{E}_{q,\phi}}{\tau} - \frac{t\ln(1-q \phi \tilde{E}_{q,\phi})}{\tau}.$$ If we put $F = (1-q) \phi/\tau$ and $B = q \phi$, the process $E_{q,\phi}(t)$ above can be expressed as $${\label{eq:Efq}}
    E_{q,\phi}(t) = \int_0^{t} {dt^\prime}[F + \epsilon({t^\prime}/\tau)] \text{,}$$ where $\epsilon(t/\tau)$ is simple exponential noise; cf. Sec. \[sec:mic\]. By using [Eq. (\[eq:Efq\])]{}, [Eq. (\[eq:lgv2\])]{} can be written as $$\ddot\alpha(t) + a \dot\alpha(t) + b^2 \alpha(t) = A\omega(t) + dE_{q,\phi}(t)/dt.$$

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. \[sec:mic\]]{}, has a quite simple form $${\label{eq:propto}} k(\alpha) \propto q \text{.}$$ An alternative parameterization, similar to [Eq. (\[eq:Efq\])]{}, 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 \[sec:apx\], like we did in [Eq. (\[eq:muqf\])]{}. The constant of proportionality in [Eq. (\[eq:propto\])]{}, for simple exponential noise, ${\mathrm{const}}_\mathrm{ex}$, and exponential shot noise, ${\mathrm{const}}_\mathrm{sh}$, are $${\mathrm{const}}_\mathrm{ex}= \frac{3}{16} (2 + \frac{b^2}{a^2});\quad
    {\mathrm{const}}_\mathrm{sh}= \frac{1}{4} (2 + \frac{b^2}{a^2}).$$

\[sec:alg\]A macroscopically symplectic algorithm for simulations of Langevin Dynamics
======================================================================================

The differential equation (\[eq:lgv2\]) is mathematically equivalent to the following system: $${\label{eq:lgv}}
    \begin{cases}
        \dot\alpha(t) = \beta(t)\\
        \dot\beta(t) = - a \beta(t) - b^2 \alpha(t) + F + A\omega(t) + B\epsilon(t/\tau)\\
    \end{cases}.$$ The macroscopic limit of [Eq. (\[eq:lgv\])]{} is $${\label{eq:mac}}
    \begin{cases}
        \dot\alpha(t) = \beta(t)\\
        \dot\beta(t) = - a \beta(t) - b^2 \alpha(t) + f\\
    \end{cases};$$ cf. Eqs. (\[eq:lgv1\])-(\[eq:lgv2\]) in [Sec. \[sec:intro\]]{}. Equation (\[eq:mac\]) is time-reversible, like the equations of motion in our MD model of [Sec. \[sec:mic\]]{}.

In order to design an algorithm for the LD simulations, which respects the macroscopic time reversibility, consider first the deterministic equation (\[eq:mac\]). The phase space of this system, as well as of [Eq. (\[eq:lgv\])]{}, 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. \[sec:mic\]]{} and [Ref. [@PRE2016III]]{}, by using the operator splitting formalism [@TuckermanMartyna; @MartynaTuckermanI; @MartynaTuckermanII]; cf. [@PRE2016III Appendix C]. Let us rewrite [Eq. (\[eq:lgv\])]{} as $${\label{eq:dff}} \dot{{\bm{\Gamma}}}(t) = {\mathrm{i} \mathcal{L}}\Gamma(t),$$ with the Liouvillian [@EvansMorriss Chapter 3] $${\mathrm{i} \mathcal{L}}= \beta\partial_\alpha - a\beta\partial_\beta - b^2 \alpha\partial_\beta + f\partial_\beta,$$ which is a sum of three mutually non-commuting terms $${\label{eq:Lvls}}
    {\mathrm{i} \mathcal{L}}_\alpha = \beta\partial_\alpha;\quad
    {\mathrm{i} \mathcal{L}}_{\beta\beta} = -a \beta\partial_\beta;\quad
    {\mathrm{i} \mathcal{L}}_{\alpha\beta} = (f - b^2 \alpha)\partial_\beta.$$ A formal solution of [Eq. (\[eq:dff\])]{} is $${\label{eq:sln}} {\bm{\Gamma}}(t) = \exp({\mathrm{i} \mathcal{L}}t) \Gamma(0).$$

If we neglect for a moment the dissipative part of the dynamics in [Eq. (\[eq:sln\])]{}, ${\mathrm{i} \mathcal{L}}_{\beta\beta}$, we are left with two possible choices of a symplectic algorithm, described in [Ref. [@TuckermanMartyna]]{}, namely, a widely used velocity Verlet scheme, $$\begin{aligned}
{\label{eq:velV}}
    {\bm{\Gamma}}(t + \Delta{t}) &=& \exp({\mathrm{i} \mathcal{L}}_{\alpha\beta} \Delta{t}/2) \exp({\mathrm{i} \mathcal{L}}_\alpha \Delta{t})
        \exp({\mathrm{i} \mathcal{L}}_{\alpha} \Delta{t}/2)\nonumber\\
        &\times& {\bm{\Gamma}}(t) + \mathcal{{\bm{O}}}(\Delta{t}^3),\end{aligned}$$ and a perhaps less known position Verlet scheme $$\begin{aligned}
{\label{eq:posV}}
    {\bm{\Gamma}}(t + \Delta{t}) &=& \exp({\mathrm{i} \mathcal{L}}_\alpha \Delta{t}/2) \exp({\mathrm{i} \mathcal{L}}_{\alpha\beta} \Delta{t})
        \exp({\mathrm{i} \mathcal{L}}_{\alpha} \Delta{t}/2)\nonumber\\
        &\times& {\bm{\Gamma}}(t) + \mathcal{{\bm{O}}}(\Delta{t}^3),\end{aligned}$$ 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. (\[eq:lgv2\])]{}. Therefore, to obtain a mesoscopic version of Eqs. (\[eq:velV\]) and (\[eq:posV\]), we make a substitution $${\label{eq:macmes}}
f \partial_\beta \Delta{t} \to \int_0^{\Delta{t}}dt[F + A \omega(t) + B \epsilon(t/\tau)] \partial_\beta
    = R(\Delta{t})\partial_\beta,$$ where we used [Eq. (\[eq:Rt\])]{}.

In the Verlet algorithms above, [Eq. (\[eq:macmes\])]{} affects only an operator of the form $\exp({\mathrm{i} \mathcal{L}}_{\alpha\beta}\Delta{t})$; cf [Eq. (\[eq:Lvls\])]{}. 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. (\[eq:velV\])]{}, this operator is applied two times, in the form $\exp({\mathrm{i} \mathcal{L}}_{\alpha\beta}\Delta{t}/2)$, while [Eq. (\[eq:posV\])]{} uses it only once. Therefore we continue with the latter integration scheme, which is computationally more efficient.

By following the approach of Refs. [@MartynaTuckermanI; @MartynaTuckermanII], the dissipative part of the dynamics, due to ${\mathrm{i} \mathcal{L}}_{\beta\beta}$, which was until now disregarded by [Eq. (\[eq:posV\])]{}, can be numerically implemented in two qualitatively equivalent ways. The one, which we chose for the simulations in [Sec. \[sec:mes\]]{}, is given by[^2]

$$\begin{aligned}
{\label{eq:sim}}
    {\bm{\Gamma}}(t + \Delta{t}) =
        \exp\left(\frac{{\mathrm{i} \mathcal{L}}_\alpha \Delta{t}}{2}\right) \exp\left(\frac{{\mathrm{i} \mathcal{L}}_{\beta\beta} \Delta{t}}{2}\right)
        \exp({\mathrm{i} \mathcal{L}}_{\alpha\beta} \Delta{t})
    \exp\left(\frac{{\mathrm{i} \mathcal{L}}_{\beta\beta} \Delta{t}}{2}\right) \exp\left(\frac{{\mathrm{i} \mathcal{L}}_{\alpha} \Delta{t}}{2}\right)
    {\bm{\Gamma}}(t) + \mathcal{{\bm{O}}}(\Delta{t}^3).\end{aligned}$$

The exact steps of the simulation algorithm can now be read from [Eq. (\[eq:sim\])]{}, as explained in detail in Ref. [@FrenkelSmit 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

$$\begin{aligned}
{\label{eq:s1}}
    \alpha(t+\Delta{t}/2) &=& \alpha(t) + \beta(t)\Delta{t}/2;\\
{\label{eq:s2}}
    \beta(t+\Delta{t}) &=& \left[
        \exp\left(- \frac{a\Delta{t}}{2}\right)\beta(t) - b^2\alpha(t+\Delta{t}/2)\Delta{t} + R(\Delta{t})
    \right]\exp\left(- \frac{a\Delta{t}}{2}\right);\\
{\label{eq:s3}}
    \alpha(t+\Delta{t}) &=& \alpha(t+\Delta{t}/2) + \beta(t+\Delta{t})\Delta{t}/2.\end{aligned}$$

Note, that the above algorithm requires generation of the random variable $R(\Delta{t})$ only in the second intermediate step, [Eq. (\[eq:s2\])]{}. When [Eq. (\[eq:lgv2\])]{} contains exponential noise, the athermal part of $R(\Delta{t})$ is sampled once from a Gamma probability distribution [@PRE2016II]. 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.

[^1]: 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$.

[^2]: If, everywhere in [Eq. (\[eq:sim\])]{}, 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.
