Thermodynamic uncertainty relations in the presence of non-linear friction and memory

We consider a Generalized Langevin equation with exponential memory kernel [29] for a particle with position θ and velocity ω, subjected to a nonlinear force $\mathcal{F}(\omega)$ and dragged by a constant force F_ext [30]:

$$\begin{eqnarray} \dot{\theta}& = &\omega, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \dot{\omega}\left(t\right)& = &-\int_{-\infty}^t\gamma\left(t-t^{^{\prime}}\right)\omega\left(t^{^{\prime}}\right) dt^{^{\prime}}+\eta_s\left(t\right) +F_{\text{ext}} +\mathcal{F}\left(\omega\right), \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \gamma\left(t\right)& = &\frac{2}{\tau}\delta\left(t\right)+\sum_{k = 1}^n a_k e^{-\frac{t}{\tau_k}}, \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} \langle \eta_s\left(t\right)\eta_s\left(t^{^{\prime}}\right) \rangle& = &\frac{2q}{\tau}\delta\left(|t-t^{^{\prime}}|\right)+ \sum_{k = 1}^n q_k a_k e^{-\frac{|t-t^{^{\prime}}|}{\tau_k}}. \end{eqnarray} \tag{ 5 }$$

Upon defining the n auxiliary variables:

$$\begin{equation} \Omega_k = -b_k\int_{-\infty}^t dt^{^{\prime}} e^{-\frac{t-t^{^{\prime}}}{\tau_k}}\left[a_k\omega\left(t^{^{\prime}}\right)-\sqrt{\frac{2q_k a_k}{\tau_k}}\xi_k\left(t^{^{\prime}}\right)\right] \quad k \in \left[1,n\right], \end{equation} \tag{ 6 }$$

the system can be written as a Markovian SDE for the n + 1 dimensional vector $\boldsymbol{X} = \{\omega,\Omega_1,\ldots,\Omega_n\}$:

$$\begin{equation} \dot{\omega} = -\frac{\omega}{\tau}+\sum_k\frac{\Omega_k}{b_k}+\sqrt{\frac{2q}{\tau}}\xi\left(t\right)+ F_{\text{ext}}+\mathcal{F}\left(\omega\right), \end{equation} \tag{ 7a }$$

$$\begin{equation} \dot{\Omega}_k = -\frac{\Omega_k}{\tau_k}-a_k b_k\omega+\sqrt{\frac{2q_k a_k b_k^2}{\tau_k}}\xi_k\left(t\right). \end{equation} \tag{ 7b }$$

where ξ_k s are independent Gaussian white noises with zero mean $\langle \xi_k(t) \rangle = 0$ and unitary variance $\langle \xi_k(t)\xi_j(t^{^{\prime}}) \rangle = \delta_{kj}\delta(t-t^{^{\prime}})$. We define accordingly the noise amplitudes B_i , such as $B_0 = \sqrt{2q/\tau}$ and $B_i = \sqrt{2q_i a_i b_i^2/\tau_i}$ for i > 0. Taking the average on both sides of equation (7b ) and imposing stationarity, we obtain the following relation between the average values of ω and $\Omega_k$s in the stationary state:

$$\begin{equation} \quad \langle \Omega_k\rangle = -\tau_k a_k b_k \langle \omega\rangle. \end{equation} \tag{ 8 }$$

The strategy to derive the new inequality is to introduce a perturbation in equations (7) by adding to the r.h.s. of each equation the corresponding component of the vector $h\boldsymbol{V}(\omega)$, which is defined as:

$$\begin{equation} V_0\left(\omega\right) = \langle \omega\rangle/\tau - \langle \omega\rangle\mathcal{F}^{^{\prime}}\left(\omega\right) -\sum_k u_k/b_k \quad \text{and} \quad V_k = a_k b_k\langle \omega\rangle + u_k/\tau_k, \end{equation} \tag{ 9 }$$

where u_k are arbitrary parameters, to be fixed later.

Considering that the h-perturbation can be arbitrarily small, we are allowed to take $\mathcal{F}(\omega)-h\langle\omega\rangle\mathcal{F}^{^{\prime}}(\omega)\simeq \mathcal{F}(\omega-h\langle\omega\rangle)$. Then, the perturbed dynamics is equivalent to the dynamics described by equations (7) for the vector $\boldsymbol{X}_h = \{\omega -h\langle \omega \rangle,\Omega_1 - hu_1,\ldots,\Omega_n-hu_n\}$. In other words, the perturbed $h-$dynamics is the same as the original one, but with shifted averaged values:

$$\begin{equation} \langle \omega \rangle_h = \left(1+h\right)\langle \omega \rangle, \quad \langle \Omega_k \rangle_h = \langle \Omega_k \rangle+hu_k, \end{equation} \tag{ 10 }$$

and for the perturbed steady state we have $P_h(\boldsymbol{X}) = P(\boldsymbol{X}_h)$, where P_h and P are the stationary distributions for the perturbed and unperturbed dynamics, respectively. Now we apply the Cramér–Rao inequality that, in our case, reads:

$$\begin{equation} \frac{\text{Var}_{h = 0}\left(\theta\left(t\right)\right)}{\left[\partial_{h = 0} \langle \theta\left(t\right) \rangle_{h}\right]^2} = \frac{\langle \Delta\theta\left(t\right)^2 \rangle}{\left(\langle \omega \rangle t\right)^2} \unicode{x2A7E} \frac{1}{\mathcal{I}_{\text{F}}\left(h = 0\right)}, \end{equation} \tag{ 11 }$$

where $\Delta \theta(t) = \theta(t) - \langle \theta\rangle$,

$$\begin{equation} \mathcal{I_{\text{F}}}\left(h = 0\right) = -\langle \partial^2_h \ln P_h\left(\boldsymbol{X}\right) \rangle_{h = 0} + \left\langle \int_0^t dt^{^{\prime}}\sum_i \left( \frac{V_i\left(\omega\right)}{B_{i}} \right)^2 \right\rangle_{h = 0} = \mathcal{I}/2+\mathcal{B}t, \end{equation} \tag{ 12 }$$

and

$$\begin{align} \mathcal{I}& = 2\int d\boldsymbol{X} \frac{\left[\partial_h P_h\left(\boldsymbol{X}\right)\right]^2|_{h = 0}}{P\left( \boldsymbol{X}\right)}, \quad \mathcal{B} = \left(\frac{\langle \omega \rangle}{\tau} - \langle \omega\rangle \langle\mathcal{F}^{^{\prime}}\rangle -\sum_j\frac{u_j}{b_j} \right)^2 B_{0}^{-2}\nonumber\\ &\quad +\sum_k\left(a_kb_k\langle \omega \rangle +\frac{u_k}{\tau_k}\right)^2B_{k}^{-2}. \end{align} \tag{ 13 }$$

Note that $\langle\mathcal{F}^{^{\prime}}\rangle$ requires the knowledge of the stationary distribution $P(\boldsymbol{X}_h)$ to be computed. The values of u_k as a function of the model parameters can be fixed by minimising $\mathcal{B}(\boldsymbol{u})$. In this way, one obtains a bound which is optimal in this set of possible perturbations. This is done by solving the system of equations defined by $\partial \mathcal{B}(\boldsymbol{u})/\partial u_k = 0 ~\forall k$. The system is formally solved by $\boldsymbol{u} = \hat{T}^{-1}\boldsymbol{z}$ where:

$$\begin{equation} T_{ii} = \frac{1}{\tau_i^2B_{i}^2}+\frac{1}{b_i^2B_{0}^2},\quad T_{i\neq j} = \frac{1}{b_ib_j B_{0}^2}, \quad z_i = \langle\omega\rangle\left[\frac{\tau^{-1}+\langle\mathcal{F}^{^{\prime}}\rangle}{b_i B_{0}^2}-\frac{a_i b_i}{\tau_iB_{i}^2} \right]. \end{equation} \tag{ 14 }$$

First, we note that with equal thermostats $q_i = q~\forall i$ and $\langle\mathcal{F}^{^{\prime}}\rangle = 0$, one has $z_i = 0$, which means that the optimal solution is $u_k = 0~\forall~k$. This is consistent with the analytical result for the linear system where the bound for the diffusion coefficient is saturated at equilibrium.

In order to find the analytical expression of the set of values $\{u_k\}$, an exact inversion of the matrix $\hat{T}$ is required, leading to very cumbersome formulae. However, from equation (14), one has that the solutions can be written as $u_k = \langle \omega\rangle\tilde{u}_k$. It is important to note that this is not an assumption of linearity in $\langle \omega\rangle$, but a simple rescaling that will allow us to obtain some simplifications in the final inequality.

Now we can define:

$$\begin{equation} \mathcal{I} = \langle \omega\rangle^2\tilde{\mathcal{I}}, \quad \mathcal{B} = \langle \omega\rangle^2\tilde{\mathcal{B}} \end{equation} \tag{ 15 }$$

and and simplify $\langle \omega\rangle^2$ in the bound for the msd equation (11) obtaining equation (1) showed at the beginning of the paper. We note that, from equation (13), it is easy to verify that $\mathcal I$ is proportional to $\langle \omega\rangle^2$. Note that from the equation above and equations (13), one has that the case without memory ($a_k = u_k = 0~\forall~k$) could lead to the pathological condition $1/\tau-\langle \mathcal{F}^{^{\prime}} \rangle = 0$ for which equation (1) predicts a non-physical asymptotic ballistic regime. However, for frictional forces, which constitute the main focus of this paper, $\mathcal{F}^{^{\prime}}$ is non-positive, thus preventing any irregularity of the bound. Anyway, even when $\mathcal{F}^{^{\prime}}\unicode{x2A7E} 0$, the system without memory admits a potential solution so that the quantity $1/\tau-\langle \mathcal{F}^{^{\prime}} \rangle$ can be evaluated analytically to check the validity of equation (1). Moreover, from equations (13) and (1), one can also check that, in the non-optimal case ($u_k = 0~\forall~k$) and without nonlinear friction ($\mathcal{F}(\omega) = 0$), the TUR derived in [15] for linear forces is recovered. In this condition, the specific contribution to the EPR due to the presence of external forcing appears in the denominator of the bound.

Both $\tilde{\mathcal{B}}$ and $\tilde{\mathcal{I}}$ are not known in general, since they require the knowledge of the stationary distribution $P(\boldsymbol{X})$. In the following, we will concentrate on cases where the l.h.s. of equation (1) is not known but the r.h.s. side can be found analytically. An alternative and purely data-driven application of the derived bound could be the following: one could measure the MSD and then use equation (1) to infer information about the nonlinear force $\mathcal{F}$ and/or the timescales involved in the memory effects, if they are not known a priori (e.g. if only positional data are available). However, the exploration of such strategies is beyond the scope of the present paper.

Equation (1) can be used as an analytical bound for the diffusion coefficient in the following classes of models: (i) models with no memory, featuring generic nonlinear friction; (ii) models with memory, including nonlinear friction with a bounded derivative, in the large time limit. In the first case, we have a gradient system for a single degree of freedom. We can then find the stationary probability distribution and in turn compute $\tilde{\mathcal{B}}$ and $\tilde{\mathcal{I}}$. In the second case, exploiting the maximum of the modulus of the derivative of the non-linear force $\max(|\mathcal{F}^{^{\prime}}|)$, we can use the following chain of inequalities

$$\begin{align} \left(\frac{\langle \omega \rangle}{\tau} - \langle \omega\rangle \langle\mathcal{F}^{^{\prime}}\rangle -\sum_j\frac{u_j}{b_j} \right)^2& \unicode{x2A7D} \left(\frac{\langle \omega \rangle}{\tau} + \langle \omega\rangle |\langle\mathcal{F}^{^{\prime}}\rangle| -\sum_j\frac{u_j}{b_j} \right)^2 \nonumber\\ & \unicode{x2A7D} \left(\frac{\langle \omega \rangle}{\tau} + \langle \omega\rangle \max\left(|\mathcal{F}^{^{\prime}}|\right) -\sum_j\frac{u_j}{b_j} \right)^2, \end{align} \tag{ 16 }$$

and define

$$\begin{equation} \mathcal{B}_{\text{max}} = \left(\frac{\langle \omega \rangle}{\tau} + \langle \omega\rangle \max\left(|\mathcal{F}^{^{\prime}}|\right) -\sum_j\frac{u_j}{b_j} \right)^2B_{0}^{-2}+\sum_k\left(a_kb_k\langle \omega \rangle +\frac{u_k}{\tau_k}\right)^2B_{k}^{-2}\unicode{x2A7E}\mathcal{B}, \end{equation} \tag{ 17 }$$

to have an approximate but analytical bound for the MSD at long times (i.e. for the diffusion coefficient):

$$\begin{equation} \lim_{t\to \infty} \frac{1}{t}\langle \Delta\theta\left(t\right)^2 \rangle \unicode{x2A7E} \frac{1}{\tilde{\mathcal{B}}_{\text{max}}}, \end{equation} \tag{ 18 }$$

where $\tilde{\mathcal{B}}_{\text{max}} = \mathcal{B}_{\text{max}}/\langle \omega \rangle^2$. This is an explicit bound with a closed analytical form. It is however weaker than the non-analytic one in which $\langle\mathcal{F}^{^{\prime}}\rangle$ explicitly appears. Nevertheless, it can still be improved by choosing the optimal u_k that minimize $\tilde{\mathcal{B}}_{\text{max}}(\boldsymbol{u})$ (note that the chain of inequalities (16) is written for generic u_k ).

For a large class of stochastic processes, which include the model equations (7) under study, a bound for asymptotic diffusion can be derived through the standard TUR involving the EPR $\langle\dot{S}\rangle$ in the steady state [1, 15]:

$$\begin{equation} \lim_{t\to \infty} \frac{1}{t}\langle \Delta\theta\left(t\right)^2 \rangle\unicode{x2A7E} \frac{2\langle \omega \rangle^2}{\langle\dot{S}\rangle}. \end{equation} \tag{ 19 }$$

Here, we compare this previous result with the bound equation (1) derived in this paper. The EPR in the steady state for our model can be easily computed by splitting the non-linear velocity-dependent force into a reversible and irreversible part: $\mathcal{F}(\omega) = \mathcal{F}_{\text{rev}}(\omega)+\mathcal{F}_{\text{irr}}(\omega)$, where $\mathcal{F}_{\text{rev}}(\omega) = 1/2 [\mathcal{F}(\omega)+\mathcal{F}(-\omega)]$ and $\mathcal{F}_{\text{irr}}(\omega) = 1/2 [\mathcal{F}(\omega)-\mathcal{F}(-\omega)]$ . Then, using the results in [30], we obtain:

$$\begin{equation} \begin{aligned} \langle\dot{S}\rangle = &\frac{1}{q}\langle\omega F_{\text{ext}}\rangle+\frac{1}{q}\langle\omega\mathcal{F}_{\text{rev}}\left(\omega\right)\rangle-\langle \omega \rangle^2\sum_ia_i\tau_i\left(\frac{1}{q} -\frac{1}{q_i}\right)+\sum_i\frac{1}{b_i}\left(\frac{1}{q} -\frac{1}{q_i}\right)\sigma_{0i} \\ & -\frac{\tau}{q}\sum_i\frac{1}{b_i}\langle \mathcal{F}_{\text{irr}}\left(\omega\right)\Omega_i \rangle-\frac{\tau}{q}\langle \mathcal{F}_{\text{irr}}\left(\omega\right)F_{\text{ext}} \rangle-\frac{\tau}{q}\langle \mathcal{F}_{\text{irr}}\left(\omega\right)\mathcal{F}_{\text{rev}}\left(\omega\right) \rangle, \end{aligned} \end{equation} \tag{ 20 }$$

where $\sigma_{0i} = \langle (\omega-\langle\omega\rangle)(\Omega_i-\langle\Omega_i\rangle) \rangle$.

In the cases of interest for this paper, the velocity-dependent force has only the irreversible part ($\mathcal{F}_{\text{rev}}(\omega) = 0$) and the expression substantially simplifies:

$$\begin{align} \langle\dot{S}\rangle& = \frac{1}{q}F_{\text{ext}}\langle\omega\rangle-\langle \omega \rangle^2\sum_ia_i\tau_i\left(\frac{1}{q} -\frac{1}{q_i}\right) +\sum_i\frac{1}{b_i}\left(\frac{1}{q} -\frac{1}{q_i}\right)\sigma_{0i}-\frac{\tau}{q}\sum_i\frac{1}{b_i}\langle \mathcal{F}_{\text{irr}}\left(\omega\right)\Omega_i \rangle\nonumber\\ & \quad - \frac{\tau}{q}F_{\text{ext}}\langle \mathcal{F}_{\text{irr}}\left(\omega\right) \rangle, \end{align} \tag{ 21 }$$

where $\mathcal{F}_{\text{irr}}(\omega)$ refers to any frictional force. The case without memory is simpler:

$$\begin{equation} \langle\dot{S}\rangle = \frac{1}{q} F_{\text{ext}}\langle\omega\rangle -\frac{\tau}{q}F_{\text{ext}} \langle \mathcal{F}_{\text{irr}}\left(\omega\right)\rangle. \end{equation} \tag{ 22 }$$

As we expect from [31], the dynamics is reversible when there is a single thermostat and $F_{\text{ext}} = 0$. Note that in the presence of a spatial potential, the system would have been out-of-equilibrium also in the case $F_{\text{ext}} = 0$ [31].

The bound predicted by the standard TUR coincides with equation (19) with $\langle\dot{S}\rangle$ given by equations (21) and (22). However, there are many crucial differences between this bound and equation (1):

• 1.  
  The standard bound is valid only in the long-time regime.
• 2.  
  It requires the computation of non-trivial averages on the steady state, that are more difficult to bound with respect to $\langle\mathcal{F}^{^{\prime}} \rangle$. Thus, making the case with memory analytical is non-trivial.
• 3.  
  In general, its denominator does not scale with $\langle \omega^2 \rangle$ and therefore the case of free diffusion (i.e. without any applied force) becomes more delicate to treat; in our result, on the contrary, the scaling with $\langle \omega^2 \rangle$ allows the free diffusion limit to have a well defined bound.

From the above remarks, it follows that a connection with our bound is not a direct task for the case with $F_ \textrm{ext} = 0$. A comparison can be done (see section 3.3) in the absence of memory: exploiting the potential form of the solution of the gradient system, one can explicitly compute the averages and check if a good scaling of $\langle \mathcal{F}_{\text{irr}}(\omega)\rangle$ with $\langle \omega\rangle$ occurs. In general, the bound derived from our procedure is more appropriate for the estimate of the diffusion coefficient.

In figure 1 we show the case where both variables are at temperature $q = q_1 = 1$, i.e. the system is at equilibrium unless the external force is switched on. Several curves in the main graph illustrate the behavior of the mean squared displacement (rescaled by time) for different values of the memory coupling a, where we remind that a = 0 corresponds to a lack of memory. In the inset (black dots), we consider the diffusion coefficient $2D = \lim_{t\to\infty}\langle [\theta(t)-\theta(0)]^2\rangle/t$, as a function of a. The first observation is the non-monotonic behaviour of the diffusion coefficient as a function of a. In the descending regime (the values of a for which D decreases with a) the curve of $MSD/t$ appears to be also non-monotonous. Apparently, the second time-scale introduced by the memory is responsible for a kind of dynamical slowing down which eventually influences the asymptotic diffusion, reducing D. In the inset, we also report the two bounds $1/B_{\text{max}}$ and $1/B_{\text{max,2}}$. The tightest bound fairly approximates the diffusion coefficient, including its non-monotonic behaviour. The non-optimised bound, on the contrary, gets close to D only in the limit a → 0 (no memory), while in the rest of the range it significantly underestimates D, showing a monotonic trend.

Zoom In Zoom Out Reset image size

In figure 2 we repeat the same analysis in a weakly out-of-equilibrium case q = 1.5 and q = 1 (the variable representing memory is slightly hotter than the main velocity variable). The scenario looks similar to the equilibrium case, with a non-monotonic behavior of both $MSD/t$ vs. t (for large a) and of D vs a. In the inset, we report the two bounds $1/B_{\text{max}}$ and $1/B_{\text{max,2}}$, finding the same behaviour observed in the previous case.

Zoom In Zoom Out Reset image size

Last, we consider a strongly out-of-equilibrium case with q = 10 and q₁ = 1, in figure 3. Here, while the curves of $MSD/t$ vs. t are still non-monotonous when memory is strong (large a), the scenario for D vs. a is different from the previous cases, as it is monotonously decreasing. The bound is confirmed and also monotonously decreasing, but it heavily underestimates D in the limit of weak memory.

Zoom In Zoom Out Reset image size

In this section, we study the PT frictional model. Figures 4–6 show the results of simulations for equilibrium, weakly out-of-equilibrium and strongly out-of-equilibrium cases, respectively. The scenario is quite similar to that observed in the Coulomb model. At large memory the $MSD/t$ vs t curve is non-monotonous, and D vs. a is also non-monotonous in the equilibrium and weakly out-of-equilibrium regimes. In these two regimes the optimised bound works very well, and much better than the non-optimised case. In the strongly out-of-equilibrium case, the non-monotonicy of D vs. a disappears and the two bounds become indistinguishable. Remarkably, for this model, the bound provides a nice approximation of D for a much larger range of memory coupling values a, with respect to the Coulomb case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The non-linear friction models considered here are, usually, not analytically solvable. Nevertheless, in the absence of memory, the single variable SDE becomes a gradient system and this gives easily an expression for the steady-state probability distribution:

$$\begin{equation} P\left(\omega\right) = \frac{\mathcal{N}\tau}{2q}\exp \left[\frac{1}{q} \int_0^\omega d\omega^{^{\prime}} \left(-\omega^{^{\prime}}+\tau F_{\text{ext}}+\tau\mathcal{F}_{\text{Coulomb}}\left(\omega^{^{\prime}}\right)\right) \right], \end{equation} \tag{ 31 }$$

where $\mathcal{N}$ is the normalization factor. The above formula is useful to get analytical results for the quantities $\mathcal{B}$ and $\mathcal{I}$ which are included in the most general bound expression, equation (1).

For the sake of simplicity, we consider only the Coulomb model without memory, that is equation (23) with a = 0 and with the force given by equation (25). We analyse the effect of τ on the diffusion coefficient. As a preliminary analysis, we verified that at small τ, where the force-velocity curve is almost linear, the steady state pdf is quasi-Gaussian, while at large τ, where the force non-linearity is important, the steady distribution is substantially non-Gaussian (see figure 7).

Zoom In Zoom Out Reset image size

We also take into account the effect of an external force F_ext, which is necessary in order to compare our results with the standard TUR, since its expression requires $F_{\text{ext}} \neq 0$. Indeed, contrary to what happens for our bound in equation (1), where we have simplified the terms $\langle \omega\rangle^2$, in the bound derived from the standard TUR equation (19) the limit $F_{\text{ext}}\to 0$ can be singular depending on the implicit dependence of $\langle \omega\rangle^2$ and $\langle \dot{S}\rangle$ on F_ext.

In figure 8 we show the behavior of the fully analytical bound as a function of τ, for a case with a vanishing external force (panel a) and a non-zero external force (panel b). For this model, the limit $F_{\text{ext} \to 0}$ of $\langle \omega\rangle^2/\langle \dot{S}\rangle$ is regular so we can do the comparison with the bound coming from the standard TUR in panel a. The blue curves correspond to the diffusion coefficient measured from numerical simulations, the black dashed lines to the bound obtained from the standard TUR, the red solid lines correspond to our bound discussed here and the green dash-dotted curves correspond to our bound where $\langle \mathcal{F^{^{\prime}}}\rangle$ is replaced with $\max(|\mathcal{F}^{^{\prime}}|)$. We first note that the non-monotonous behavior in τ is well captured by all the bounds. Then, we underline that when $F_{\text{ext}} \sim 0$ the standard TUR usually surpasses the TUR discussed here. The opposite occurs when an external force is present, particularly the TUR discussed here is largely dominating in the region, at large τ, where non-linear and non-Gaussian effects are important.

Zoom In Zoom Out Reset image size

In figure 9 we show the behavior of the fully analytical bound as functions of the external force F_ext. In panel (a), the effect discussed previously is enhanced: the TUR derived in this paper gives a tighter bound with respect to the one obtained from the standard TUR for $F_{\text{ext}}\gtrsim 0.6$. In panel (b), we show the results for regular Brownian diffusion. As we expect, our bound and the one coming from the standard TUR coincide in this case.

Zoom In Zoom Out Reset image size