# The OU2 process: Characterising dissipative confinement in noisy traps

Luca Cocconi1, Henry Alston2, Jacopo Romano1 and Thibault Bertrand2 1 Max
Planck Institute for Dynamics and Self-Organization (MPIDS), 37077 Göttingen,
Germany
2 Department of Mathematics, Imperial College London, South Kensington, London
SW7 2AZ, United Kingdom<EMAIL_ADDRESS>
###### Abstract

The Ornstein-Uhlenbeck (OU) process describes the dynamics of Brownian
particles in a confining harmonic potential, thereby constituting the
paradigmatic model of overdamped, mean-reverting Langevin dynamics. Despite
its widespread applicability, this model falls short when describing physical
systems where the confining potential is itself subjected to stochastic
fluctuations. However, such stochastic fluctuations generically emerge in
numerous situations, including in the context of colloidal manipulation by
optical tweezers, leading to inherently out-of-equilibrium trapped dynamics.
To explore the consequences of stochasticity at this level, we introduce a
natural extension of the OU process, in which the stiffness of the harmonic
potential is itself subjected to OU-like fluctuations. We call this model the
OU2 process. We examine its statistical, dynamic, and thermodynamic properties
through a combination of analytical and numerical methods. Importantly, we
show that the probability density for the particle position presents power-law
tails, in contrast to the Gaussian decay of the standard OU process. In turn,
this causes the trapping behavior, extreme value statistics, first passage
statistics, and entropy production of the OU2 process to differ qualitatively
from their standard OU counterpart. Due to the wide applicability of the
standard OU process and of the proposed OU2 generalisation, our study sheds
light on the peculiar properties of stochastic dynamics in random potentials
and lays the foundation for the refined analysis of the dynamics and
thermodynamics of numerous experimental systems.

  *

##  1 Introduction

The Ornstein-Uhlenbeck (OU) process is a continuous-time Gaussian stochastic
process with linear mean-reverting properties first introduced to describe the
fluctuating velocity of a Brownian particle immersed in a fluid [1, 2]. It has
found over the years countless applications across various subfields of
physics (as well as other disciplines), where it plays a similarly
paradigmatic role as that of the harmonic oscillator in classical and quantum
mechanics. It can be understood as the overdamped limit of any Langevin
dynamics exploring the local neighborhood of a differentiable minimum of an
arbitrary potential landscape. Its Langevin equation of motion generically
reads

$\frac{dx(t)}{dt}=-\bar{k}x(t)+\sqrt{2D_{x}}\zeta(t)\,.$ (1)

where, depending on the context, $\bar{k}>0$ might be interpreted as a
friction coefficient or potential stiffness coefficient, while $D_{x}$ denotes
the diffusivity and $\zeta(t)$ is a delta-correlated zero mean and unit
variance white noise. Key results including steady-state probability density
function, formal solution and Green’s function for the standard OU process are
reviewed for completeness in A. To give but one example of its wide
applicability for instance in spatially extended settings, a lattice of
elastically coupled OU processes formally defines the Gaussian free field
around which the perturbative expansion of non-conserved dynamical field
theories is typically constructed [3, 4].

While originally modelling the frictional contribution to Brownian motion, the
first term in the right-hand side of Eq. (1) is often interpreted in the case
of overdamped dynamics as a restoring force resulting from an effective
harmonic potential $V(x)=\bar{k}x^{2}/2$ acting on the coordinate $x$. This is
the case, for instance, in many physical models of micro-particle manipulation
by optical tweezers [5, 6, 7, 8, 9], where force gradients are established
through the inhomogeneous electric field within a highly focused laser beam
[10, 11, 12, 13]. However, the nature of the potential $V(x)$ might even be
more abstract, as exemplified by models of mean reverting portfolios in
finance [14] or continuous trait evolution in ecology [15].

In all such cases, it is reasonable to expect that the underlying processes
governing the potential are themselves subject to some degree of
stochasticity, implying that $V(x,t)$ may itself be a stochastic process [16,
17]. A case in point is that of optical tweezers controlled by real laser
systems, which are characterised by small fluctuations in power output around
its mean [12, 13]; these fluctuations in power lead in turn to fluctuations in
the stiffness of the potential experienced by the dielectric particle.

Inspired by this rather simple idea, we define here a generic model of
diffusion in a noisy trap. Namely, we introduce continuous, zero-mean
fluctuations in the potential stiffness of the original OU process [Eq. (1)].
More precisely, we model these fluctuations themselves by an OU process: we
characterise this second process by an effective stiffness $\mu$ and effective
diffusivity $D_{k}$ (see Fig. 1 for a schematic illustration and example
trajectories). Overall, the resulting coupled dynamics of the particle
position $x(t)$ and the fluctuations in the confining potential stiffness
$k(t)$ read

$\displaystyle\frac{dx(t)}{dt}$
$\displaystyle=-\left[\bar{k}+k(t)\right]x(t)+\sqrt{2D_{x}}\zeta_{x}(t)$ (2a)
$\displaystyle\frac{dk(t)}{dt}$ $\displaystyle=-\mu
k(t)+\sqrt{2D_{k}}\zeta_{k}(t)$ (2b)

where we fix the average stiffness $\bar{k}>0$, $k(t)$ is the zero-mean
fluctuating contribution and
$\langle\zeta_{i}(t)\zeta_{j}(t)\rangle=\delta_{ij}\delta(t-t^{\prime})$. The
particle dynamics reduce to a standard OU process upon setting $D_{k}=0$. We
call this generic composite stochastic process the OU2 process and dedicate
the rest of this paper work to its extensive characterisation111The OU2 model,
which we introduce here, is not to be confused with the squared-OU models, a
term sometimes used in the context of the modelling of interest rates by the
Cox-Ingersoll-Ross model [18].. To the best of our knowledge, this model has
been introduced for the first time in two recent works by the authors [16, 17]
in the context of nonequilibrium thermodynamics of diffusion in fluctuating
potentials. Alternative generalisations of the OU process have also been
investigated including models in which the location of the confining potential
minimum undergoes stochastic [19] or oscillatory [20] “sliding” dynamics at
fixed stiffness. The OU2 model is closely related to a number of other models,
which we mention briefly below.

Figure 1: The OU2 process as a minimal model of dissipative confinement. (a)
Gaussian fluctuations around a positive mean in the stiffness coefficient of
the harmonic potential acting on an overdamped Brownian particle lead to the
establishment of futile breathing cycles, with transients of
increasing/decreasing stiffness driving the particle closer to/further from
the origin. (b) Example trajectories for the particle position $x(t)$ and the
total confining potential stiffness $\bar{k}+k(t)$, with $\bar{k}=1$,
$D_{x}=1$, $\mu=0.01$ and $D_{k}=0.02$.

Firstly, we note that the OU2 process is closely related to the problem of
Brownian motion in an intermittent harmonic potential [21, 22], where the
potential switches stochastically between two states with finite stiffnesses
$k_{1}$ and $k_{2}$ in the manner of a telegraph process. The case where
$k_{1}=0$ and $k_{2}>0$ has recently received some attention as it represents
a realistic implementation of stochastic resetting [23, 24, 25, 26, 16];
interestingly, it was shown in this case that any degree of intermittency
leads to the establishment of a nonequilibrium stationary probability density
for the trapped particle’s position displaying a Gaussian bulk which
eventually crosses over into exponential tails.

Furthermore, the OU2 process introduced here constitutes a continuous time
extension of a class of so-called “random difference equations” (see for
instance [28]), i.e. recurrence relations involving random parameters. For
instance, the model introduced in Ref. [29] is effectively a time-discretised
version of the OU process where the confining potential exhibits a fluctuating
stiffness, whose fluctuations are uncorrelated in time. More recently, Morita
[30] focused on a version of these processes where the random parameter is
allowed to have correlations in time. Nevertheless, in this work, the author
consider the much simpler case where $k(t)$ is governed by a Poisson jump
process, arguing that studying the case where $k(t)$ is governed by an OU
process, which is precisely a discrete-time analog of OU2 process, is
particularly challenging.

In the study of transport in inhomogeneous environments, recent models of
diffusing diffusivity have been introduced, which allow for stochastic
fluctuations in diffusivity of a free Brownian particle. These models display
Fickian diffusion (characterised by a linear time dependence of the mean-
square displacement) in the presence of non-Gaussian displacement
distributions [36, 37, 38].

Moreover, establishing a connection with non-equilibrium thermodynamics, the
OU2 process can be seen as a stochastically breathing harmonic potential.
Standard breathing potentials, whereby the stiffness is modulated
deterministically in time according to a pre-defined protocol, are often
studied in the context of heat engines operating in finite time cycles [39,
40, 5]. Interestingly, generic results have been obtained in the slow driving
regime for the full distribution of the stochastic work [41].

An alternative inspiration for the OU2 process may be found in the context of
motile active matter, particularly in the canonical active particle model
known as the Active Ornstein-Uhlenbeck particle [42, 43] (AOUP). Here, out-of-
equilibrium self-propulsion is introduced via a forcing term in the Langevin
equation of motion, whose statistics are those of a zero-mean OU process.
However, one may equivalently interpret this term as a linear potential whose
amplitude is modulated stochastically in time. The OU2 process is thus a
natural, non-motile counterpart to the (much more studied) AOUP, offering a
minimal example of dissipative trapping and single-particle irreversibility
beyond active motility.

Finally, we will see shortly that the so-called random acceleration model can
be a seen as a special case of the OU2 process [31, 32, 33, 34, 35].

The paper is organised as follows: we begin in Section 2 with a preliminary
calculation of the marginal probability density function of the position for
an OU2 process with vanishing positional noise, $D_{x}=0$, highlighting some
non-trivial characteristics of the associated statistics. In Section 3, we
move away from this limit and derive the conditional and full Green’s
functions of the full OU2 process, clarifying the condition for the stability
of the dynamics. Section 4 deals with the moments of the marginal, stationary
probability density function for coordinate $x$ in Eq. (2a). In particular, we
obtain the necessary conditions for the existence of the even moments $\langle
x^{2n}\rangle$ for arbitrary $n\geq 1$ in the form of an upper bound on
$\bar{k}\mu^{2}/D_{k}$ which decreases monotonically to zero with increasing
$n$, as well as closed form expressions for the second and fourth moments. In
Section 5, we study two limits of fast stiffness dynamics by means of
homogenisation [44], solving analytically the resulting coarse-grained Fokker-
Planck equation for the slow dynamics. In Section 6, we draw on known
heuristic arguments from extreme value statistics of weakly correlated time
series to conjecture the distribution of the maximum of a finite time OU2
process, verifying our proposed classification via numerical simulations.
Owing to the algebraic nature of the tails of the marginal probability
density, we observe an unexpected transition from the Gumbel to the Frechet
universality class. Section 7 is dedicated to investigating the impact of
stiffness fluctuations on the mean first passage time (mFPT) to a stationary
target. Finally, we re-derive previous results for the steady-state entropy
production rate [16, 17] in a compact way in Section 8, thus offering a simple
thermodynamic characterisation of the model. Finally, some remarks and
potential directions for future research are discussed in the Conclusion.

##  2 Exact solution in the limit of vanishing positional noise

As a preliminary analysis, we consider the limiting case of vanishing
positional noise, $D_{x}=0$ in Eq. (2a), for which a series of exact results
can be derived. In this case the dynamics of the particle are constrained to
the half-line to which the initial condition belongs. Here, by symmetry, we
set $x(0)>0$ so that $x(t)>0$ without loss of generality.

###  2.1 Probability density function

Let us denote $x_{0}\equiv x(0)$ and define $z\equiv\ln x$ such that we can
recast the dynamics as

$\dot{z}(t)=-\bar{k}-k(t)$ (3)

or, equivalently,

$\ddot{z}(t)=\mu k(t)-\sqrt{2D_{k}}\zeta_{k}(t)~{},$ (4)

which reduces to the random acceleration process [31, 32, 33, 34, 35] in the
limit $\mu\to 0$.

Figure 2: Marginal probability density function $P(x,t)$ in the limit of
vanishing positional noise, as given by Eq. (9), at various time in linear and
double-logarithmic scales. Here, we set $\bar{k}=1$, $x_{0}=1$ and
$2D_{k}/\mu^{2}=1$.

As initial condition at $t=0$, we choose $z(0)=z_{0}$ (and correspondingly,
$x_{0}=e^{z_{0}}$) and assume that the fluctuating stiffness has been evolving
from $t\to-\infty$ such that at $t=0$ the particle experiences a value of $k$
randomly drawn from the steady-state distribution. The solution of Eq. (3) is
thus written

$z(t)=z_{0}-\bar{k}t-\sqrt{2D_{k}}\int_{0}^{t}dt^{\prime}\int_{-\infty}^{t^{\prime}}dt^{\prime\prime}e^{-\mu(t^{\prime}-t^{\prime\prime})}\zeta_{k}(t^{\prime\prime})~{},$
(5)

of which the time-dependent mean and variance are computed straightforwardly
to be

$\displaystyle\bar{z}(t)$ $\displaystyle\equiv\langle
z(t)\rangle=z_{0}-\bar{k}t,$ (6a) $\displaystyle\sigma_{z}^{2}(t)$
$\displaystyle\equiv\langle z^{2}(t)\rangle-\langle
z(t)\rangle^{2}=\frac{2D_{k}}{\mu^{2}}t-\frac{2D_{k}}{\mu^{3}}(1-e^{-\mu
t})~{},$ (6b)

where $\langle\cdots\rangle$ means an average over realisations of the noise
$\zeta_{k}$. It is interesting to note that at long times, i.e.
$t\gg\mu^{-1}$, the variance scales linearly with time as
$\sigma_{z}^{2}(t)\simeq 2D_{k}t/\mu^{2}$ confirming that the process remains
diffusive even in the absence of positional noise. Since (5) is Gaussian, the
first two moments are sufficient to determine the time-dependent moment
generating function,

$Z(q,t)\equiv\langle
e^{-iqz(t)}\rangle=\exp\left[-iq\bar{z}(t)-\frac{\sigma_{z}^{2}(t)}{2}q^{2}\right]$
(7)

from which the corresponding probability density is easily obtained,

$P(z,t)=\frac{1}{\sqrt{2\pi\sigma^{2}_{z}(t)}}\exp\left[-\frac{(z-\bar{z}(t))^{2}}{2\sigma^{2}_{z}(t)}\right].$
(8)

Via a straightforward transformation of probability, we then obtain that the
probability density for the original variable $x(t)$ is given by the following
log-normal distribution

$\displaystyle P(x,t)$
$\displaystyle=\frac{dz}{dx}P(z(x),t)=\frac{1}{x\sqrt{2\pi\sigma^{2}_{z}(t)}}\exp\left[-\frac{(\ln(x)-\bar{z}(t))^{2}}{2\sigma^{2}_{z}(t)}\right],$
(9)

The exact probability density (9) is shown for different values of $t$ in Fig.
2.

###  2.2 Growth and trapping

Equipped with the time-dependent probability density (9), we can start to
study the main qualitative features of the system. Namely, studying three key
statistics of the process — the median, mean and mode of $x(t)$ — each of
which offers a different perspective on the dynamics, we discuss the
conditions under which the process is said to be: (i) trapped in the sense
that the associated statistic reverts back to the center of the potential
(here, $x=0$) or (ii) growing in the sense that the associated statistic grows
exponentially in time.

##### Median —

The median is defined as the value $x_{M}$ in the support of the distribution
such that

$\int_{0}^{x_{M}}dx\ P(x,t)=\frac{1}{2}~{}.$ (10)

For the log-normal distribution (9), an explicit expression for the median can
be written and simply gives

$x_{M}(t)=\exp[\bar{z}(t)]=x_{0}\exp(-\bar{k}t)$ (11)

Interestingly, we note that the behavior of the median changes as $\bar{k}$
changes sign. Indeed for $\bar{k}>0$, we can easily see that the median
$x_{M}(t)$ approaches zero exponentially as $t$ increases. Conversely, for
$\bar{k}<0$, we find that $x_{M}(t)$ grows exponentially with $t$. Note that
the same result holds for every finite percentile of the distribution. In
other words, as $\bar{k}$ changes sign from positive to negative, most of the
particles go from being trapped around $x=0$ to seeing their position grow
exponentially in time.

##### Mean —

Secondly, let us consider the mean of the distribution (9), which can be
computed exactly to be

$\bar{x}(t)=\exp\left[\bar{z}(t)+\frac{\sigma_{z}^{2}}{2}\right]=x_{0}e^{\frac{D_{k}}{\mu^{3}}(1-e^{-\mu
t})}\,e^{-\left(\bar{k}-\frac{D_{k}}{\mu^{2}}\right)t}\mathrel{\mathop{\propto}\limits_{t\to\infty}}\exp\left[-\left(\bar{k}-\frac{D_{k}}{\mu^{2}}\right)t\right].$
(12)

The mean thus decays exponentially with time for $\bar{k}>D_{k}/\mu^{2}$, a
more stringent condition compared to that of median trapping.

##### Mode —

Finally, we turn our attention to the mode $x_{m}$, which is defined as the
location of the maximum of the probability density,
$\left.\partial_{x}P(x,t)\right|_{x_{m}}=0$. For the log-normal distribution
in Eq. (9), the mode is given by

$x_{m}(t)=\exp\left[\bar{z}(t)-\sigma_{z}^{2}\right]=x_{0}e^{-\frac{2D_{k}}{\mu^{3}}(1-e^{-\mu
t})}\,e^{-\left(\bar{k}+\frac{2D_{k}}{\mu^{2}}\right)t}\mathrel{\mathop{\propto}\limits_{t\to\infty}}\exp\left[-\left(\bar{k}+\frac{2D_{k}}{\mu^{2}}\right)t\right].$
(13)

Here, we notice that the behavior of the mode for the OU2 process changes from
trapped to growing as $\bar{k}+2D_{k}/\mu^{2}$ changes sign. In other words,
for $\bar{k}<-2D_{k}/\mu^{2}$, the distance of the most likely outcome grows
exponentially.

#####

Remarkably, these results imply the existence of two non-trivial regimes:

  * (i)

for $0<\bar{k}<D_{k}/\mu^{2}$, the mean position of an ensemble of particles
undergoing OU2 processes without positional noise grows indefinitely even if
most of the particles approach zero exponentially; this is due to rare
trajectories with exceptionally large displacements.

  * (ii)

for $-2D_{k}/\mu^{2}<\bar{k}<0$, most of the particles escape exponentially to
infinity, however the most likely outcome remains trapped at the proximity of
the origin.

###  2.3 Optimal trapping and condition for growth

It is interesting to note that, while the parameters $\mu$ and $D_{k}$ might
not be accessible to direct experimental control in many physical
implementations of this model, one may still be able to control the strength
of the couplings between the particle and the potential via a medium-dependent
parameters $u$ such that $\dot{x}=-u(\bar{k}+k(t))$. In this case, the mean
position becomes
$\bar{x}(t;u)\propto\exp{\left[\left(-\bar{k}u+u^{2}\frac{D_{k}}{\mu^{2}}\right)t\right]}$
for large enough $t$. For example in an experiment with optical tweezers, $u$
can be tuned by changing the viscosity or dielectric properties of the
colloid, while it may not be possible to improve the properties of the
confining potential by increasing its mean $\bar{k}$ or reducing the noise
strength $D_{k}$. The exponent is now negative for
$u<\frac{\mu^{2}\bar{k}}{D_{k}}$. In other words, provided that $\bar{k}$ is
positive we can always induce mean trapping by reducing the coupling of the
particles with the potential itself. On the other hand, an excessive reduction
of the coupling $u$ results in a weak confinement. The optimal value of $u$
minimising the exponent, hence the magnitude of the mean, and thus providing
the best confinement of the latter is $u_{\rm opt}={\mu^{2}}/{2D_{k}}$.

Another application of the OU2 process can be found in finance, where Eq. (2a)
can be used as a simplified model for the growth of the capital $x$ of a
company. In this case, exponential growth rather than trapping is the desired
outcome. In this model, the stiffness $k$ represents instead (minus) the
return on investment of its operations, which is itself subject to stochastic
market fluctuations. While the mean return $-\bar{k}$ and volatility $D_{k}$
of investment can be difficult to improve, $u$ can be tuned by simply
reinvesting more capital, while $u>1$ can be obtained by using leverage. The
expected capital will grow for $u>u_{\rm{th}}=\frac{\mu^{2}\bar{k}}{D_{k}}$,
which can happen also for companies with negative mean returns.

##  3 Conditional and full Green’s function

We now return to the full model, including noise in the displacement, and
calculate its Green’s function, first conditional on a particular value of $k$
at the time of perturbation, then averaged over the corresponding steady-state
distribution. From Eq. (2), we can write that the formal solution for $x$ is
given by

$x(t)=\sqrt{2D_{x}}\int_{-\infty}^{t}dt^{\prime}\
\zeta_{x}(t^{\prime})\exp\left[-\bar{k}(t-t^{\prime})-\int_{t^{\prime}}^{t}dt^{\prime\prime}\
k(t^{\prime\prime})\right]~{}.$ (14)

From Eq. (14), we identify the conditional Green’s function of the process,
which describes the typical evolution of a noise-generated perturbation,

$\mathcal{G}(t;k_{0})=\left\langle\exp\left[-\bar{k}t-\int_{0}^{t}dt^{\prime}\
k(t^{\prime})\right]\right\rangle_{k_{0}}\Theta(t)$ (15)

where $\langle\bullet\rangle_{k_{0}}$ denotes an average with respect to the
possible realisations of the process $k(t)$ conditioned on the initialisation
$k(0)=k_{0}$, and the Heaviside theta function $\Theta(t)$ ensures causality.
The conditional Green’s function $\mathcal{G}(t;k_{0})$ quantifies the typical
temporal evolution of a perturbation generated at $t=0$ by the noise
$\zeta_{x}$. Exploiting the relation between the moment and cumulant
generating functions and the fact that the cumulants of order 3 and above
vanish for the OU process due to it being Gaussian, we write

$\displaystyle\mathcal{G}(t;k_{0})$
$\displaystyle=e^{-\bar{k}t}\exp\sum_{n=1}^{\infty}\frac{(-1)^{n}}{n!}\left\langle\left(\int_{0}^{t}dt^{\prime}k(t^{\prime})\right)^{n}\right\rangle_{c,k_{0}}$
(16a) $\displaystyle=e^{-\bar{k}t}\exp\left[-\int_{0}^{t}dt^{\prime}\langle
k(t^{\prime})\rangle_{c,k_{0}}+\frac{1}{2}\int_{0}^{t}dt^{\prime}dt^{\prime\prime}\
\langle k(t^{\prime})k(t^{\prime\prime})\rangle_{c,k_{0}}\right]~{},$ (16b)

where $\langle\bullet\rangle_{c,k_{0}}$ denotes the conditional cumulants. The
conditional cumulants of first- and second-order can easily be calculated
independently

$\displaystyle\langle k(t^{\prime})\rangle_{c,k_{0}}$
$\displaystyle=k_{0}e^{-\mu t^{\prime}},$ (17a) $\displaystyle\langle
k(t^{\prime})k(t^{\prime\prime})\rangle_{c,k_{0}}$
$\displaystyle=\frac{D_{k}}{\mu}\left(e^{-\mu|t^{\prime}-t^{\prime\prime}|}-e^{-\mu(t^{\prime}+t^{\prime\prime})}\right).$
(17b)

With this result in hand, we can perform the integral in (16) to obtain the
conditional propagator

$\mathcal{G}(t;k_{0})=\exp\left[-\left(\bar{k}-\frac{D_{k}}{\mu^{2}}\right)t-\frac{k_{0}}{\mu}(1-e^{-\mu
t})+\frac{D_{k}}{2\mu^{3}}(4e^{-\mu t}-e^{-2\mu t}-3)\right]~{}.$ (18)

Interestingly, the propagator decays to zero at long times only if
$\bar{k}>D_{k}/\mu^{2}$, while fluctuations grow exponentially otherwise. This
highlights the importance of the competition between the two timescales in the
problem, namely $\tau_{x}=1/\bar{k}$ — the typical mean reversion time for the
particle position — and $\tau_{k}=\mu^{2}/D_{k}$ — the typical timescale for
the stiffness fluctuations.

Also note that the dependence on the initial condition for the stiffness,
$k_{0}$, is rather simple. Expanding the exponent to leading order in small
times $t\ll 1$, we find

$\mathcal{G}(t;k_{0})=\exp[-(\bar{k}+k_{0})t+\mathcal{O}(t^{2})]$ (19)

indicating that at short times the growth/decay of fluctuations is controlled
by the initial condition $k_{0}$, such that fluctuations might initially grow
exponentially (when $k_{0}<-\bar{k}$), even when they are eventually
suppressed on average at long times. In other words, the conditional
propagator is not necessarily monotonic.

We might also be interested in a situations where the initial value of the
potential stiffness $k_{0}$ is unknown. Assuming that the statistics of $k(t)$
have reached steady-state by the time we perturb our system, we can calculate
the full propagator by averaging Eq. (18) over $k_{0}$ which has a known
steady-state Gaussian probability density function, i.e.

$\mathcal{G}_{\rm
full}(t)=\exp\left[-\left(\bar{k}-\frac{D_{k}}{\mu^{2}}\right)t+\frac{D_{k}}{2\mu^{3}}(4e^{-\mu
t}-e^{-2\mu t}-3)\right]\left\langle\exp\left[-\frac{k_{0}}{\mu}(1-e^{-\mu
t})\right]\right\rangle$ (20)

where $\langle\bullet\rangle$ now denotes an expectation with respect to the
steady-state probability density function of $k_{0}$. Using the fact that the
last term in the above is simply the moment generating function of a zero-mean
normal distribution with conjugate variable $s=-(1-e^{-\mu t})/\mu$, we
eventually arrive at

$\displaystyle\mathcal{G}_{\rm full}(t)$
$\displaystyle=\exp\left[-\left(\bar{k}-\frac{D_{k}}{\mu^{2}}\right)t+\frac{D_{k}}{2\mu^{3}}(4e^{-\mu
t}-e^{-2\mu t}-3)+\frac{D_{k}}{2\mu^{3}}(1-e^{-\mu t})^{2}\right]$
$\displaystyle=\exp\left[-\left(\bar{k}-\frac{D_{k}}{\mu^{2}}\right)t+\frac{D_{k}}{\mu^{3}}(e^{-\mu
t}-1)\right]~{}.$ (21)

We note that in the limit where $D_{k}\to 0$, we recover the Green’s function
for the standard OU process, c.f. Eq. (62). The long time behaviour is the
same as for the conditional case, however expanding again at small times $t\ll
1$, we now find $\mathcal{G}_{\rm
full}(t)=\exp\left[-\bar{k}t+\mathcal{O}(t^{2})\right]$. This is to be
expected since the average of the exponential converges to the exponential of
the average when $t\to 0$.

##  4 Positional moments

We now discuss the steady-state moments of the marginal probability density
function for the coordinate $x$. In particular, we derive conditions for the
existence of finite moments. Already in the case of a Brownian particle
confined in a potential whose stiffness switches stochastically between two
finite values $k_{1}$ and $k_{2}$ following a two-state Markov jump process,
it was shown that the condition for existence of a moment of order $s$ is more
restrictive than merely ensuring that the stiffness is positive on average,
$\langle k\rangle_{t}>0$. Indeed, the existence of $\langle|x|^{s}\rangle$
requires that $k_{1}P(k_{1})+k_{2}P(k_{2})-sk_{1}k_{2}<0$ with $P$ the
stationary probability mass function of the jump process [27]. We will see in
this section that similar conditions can be derived in the OU2 case.

Starting from the formal solution, Eq. (14), we first argue by symmetry that
all odd moments are expected to vanish, $\langle x^{2n+1}(t)\rangle=0$ for all
$n\in\mathbb{N}$. The even moments are on the other hand given by

$\displaystyle\langle
x^{2n}(t)\rangle=(2D_{x})^{n}\left\langle\left(\int_{-\infty}^{t}dt^{\prime}\
\zeta_{x}(t^{\prime})e^{-\bar{k}(t-t^{\prime})}\exp\left[-\int_{t^{\prime}}^{t}dt^{\prime\prime}\
k(t^{\prime\prime})\right]\right)^{2n}\right\rangle~{}.$ (22)

Using the fact that $\zeta_{x}(t)$ and $k(t)$ are uncorrelated stochastic
processes, we write

$\displaystyle\Bigg{\langle}\prod_{i=1}^{2n}\zeta_{x}(t_{i}^{\prime})\exp\Bigg{[}-\bar{k}$
$\displaystyle(t-t_{i}^{\prime})-\int_{t_{i}^{\prime}}^{t}dt_{i}^{\prime\prime}k(t_{i}^{\prime\prime})\Bigg{]}\Bigg{\rangle}$
$\displaystyle=\left\langle\prod_{j=1}^{2n}\zeta_{x}(t_{j}^{\prime})\right\rangle\left\langle\prod_{i=1}^{2n}\exp\left[-\bar{k}(t-t_{i}^{\prime})-\int_{t_{i}^{\prime}}^{t}dt_{i}^{\prime\prime}k(t_{i}^{\prime\prime})\right]\right\rangle~{}.$
(23)

The first expectation on the right-hand side can be simplified by Wick-
Isserlis theorem to a sum of product of white noise correlators, i.e. Dirac
delta functions, with $(2n-1)!!=(2n)!/(2^{n}n!)$ summands corresponding to all
possible pairings $\mathcal{P}_{2n}$ of the random variables:

$\left\langle\prod_{j=1}^{2n}\zeta_{x}(t_{j}^{\prime})\right\rangle=\sum_{p\in\mathcal{P}_{2n}}\prod_{\\{i,j\\}\in
p}\left\langle\zeta_{x}(t^{\prime}_{i})\zeta_{x}(t^{\prime}_{j})\right\rangle~{}.$
(24)

Since the overall integral is invariant under permutation of the indices, all
summands give the same contribution. The expression for the moments thus
simplifies to

$\langle x^{2n}(t)\rangle=\mathcal{N}_{n}\int_{-\infty}^{t}dt_{1<...<n}\
\exp\left[{-2\bar{k}\sum_{i=1}^{n}(t-t_{i})}\right]\left\langle\exp\left[-2\sum_{j=1}^{n}\int_{t_{i}}^{t}dt_{i}^{\prime}k(t_{i}^{\prime})\right]\right\rangle$
(25)

with $\mathcal{N}_{n}=(2D_{x})^{n}(2n-1)!!n!$, where we have additionally
imposed the arbitrary time ordering $t_{1}<t_{2}<...<t_{n}$ in the multiple
integrals, compensated by the combinatorial prefactor $n!$, without loss of
generality. Next, we exploit the identity relating the moment generating
function and the exponential of the corresponding cumulant generating function

$\displaystyle\Bigg{\langle}\exp\Bigg{[}-2\sum_{i=1}^{n}$
$\displaystyle\int_{t_{i}}^{t}dt_{i}^{\prime}\
k(t_{i}^{\prime})\Bigg{]}\Bigg{\rangle}$
$\displaystyle=\exp\sum_{m=1}^{\infty}\frac{1}{m!}\left\langle\left(-2\sum_{i=1}^{n}\int_{t_{i}}^{t}dt_{i}^{\prime}\
k(t_{i}^{\prime})\right)^{m}\right\rangle_{c}$ $\displaystyle=\exp\left\langle
2\left(\int_{t_{1}}^{t}dt_{1}^{\prime}\
k(t_{1}^{\prime})+\int_{t_{2}}^{t}dt_{2}^{\prime}\
k(t_{2}^{\prime})+...+\int_{t_{n}}^{t}dt_{n}^{\prime}\
k(t_{n}^{\prime})\right)^{2}\right\rangle_{c}$ (26)

where we have used the fact that cumulants of order $m>2$ vanish for the
equilibrium OU process governing $k(t)$, Eq. (2b), while the first order
cumulant is zero at steady state. The right-hand side of (26) can be evaluated
as

$\displaystyle\Bigg{\langle}2\Big{(}\int_{t_{1}}^{t}$ $\displaystyle
dt_{1}^{\prime}\ k(t_{1}^{\prime})+\int_{t_{2}}^{t}dt_{2}^{\prime}\
k(t_{2}^{\prime})+...+\int_{t_{n}}^{t}dt_{n}^{\prime}\
k(t_{n}^{\prime})\Bigg{)}^{2}\Bigg{\rangle}_{c}$
$\displaystyle=\frac{2D_{k}}{\mu}\left[\left(\sum_{i=1}^{n}\int_{t_{1}}^{t}dt_{i}^{\prime}\int_{t_{1}}^{t}dt_{i}^{\prime\prime}\
e^{-\mu|t_{i}^{\prime}-t_{i}^{\prime\prime}|}\right)+2\left(\sum_{i<j}^{n}\int_{t_{i}}^{t}dt_{i}^{\prime}\int_{t_{j}}^{t}dt_{j}^{\prime}\
e^{-\mu|t_{i}^{\prime}-t_{j}^{\prime}|}\right)\right]$
$\displaystyle=\frac{4D_{k}}{\mu^{3}}\left[\left(\sum_{i=1}^{n}\mu(t-t_{i})-1+e^{-\mu(t-t_{i})}\right)\right.$
$\displaystyle\quad\qquad\qquad\left.+2\left(\sum_{i<j}^{n}e^{-\mu(t-t_{i})}+e^{-\mu(t-t_{j})}-e^{-\mu|t_{j}-t_{i}|}+2\mu(t-t_{j})-1\right)\right]~{}.$
(27)

Here we have used the result for the double integral

$\displaystyle\int_{t_{i}}^{t}dt_{i}^{\prime}\int_{t_{j}}^{t}dt_{j}^{\prime}\
e^{-\mu|t_{i}^{\prime}-t_{j}^{\prime}|}$
$\displaystyle=\frac{1}{\mu^{2}}\left(e^{-\mu(t-t_{i})}+e^{-\mu(t-t_{j})}-e^{-\mu|t_{j}-t_{i}|}+2\mu(t-{\rm
max}(t_{i},t_{j}))-1\right)~{}.$ (28)

We now compute (27) for different values of $n$.

###  4.1 Variance ($n=1$)

First consider the particular case $n=1$ for which compact expressions can be
obtained. In this case the right-hand side of (26) using (27) simplifies to

$\left.\left\langle\exp\left[-2\sum_{i=1}^{n}\int_{t_{i}}^{t}dt_{i}^{\prime}\
k(t_{i}^{\prime})\right]\right\rangle\right|_{n=1}=\exp\left[\frac{4D_{k}}{\mu^{3}}\left(e^{-\mu(t-t_{1})}+\mu(t-t_{1})-1\right)\right]~{}.$
(29)

Using this result, we can then rewrite (25) for $n=1$ as

$\langle
x^{2}\rangle=2D_{x}\int_{-\infty}^{t}dt_{1}\exp\left[-\left(2\bar{k}-\frac{4D_{k}}{\mu^{2}}\right)(t-t_{1})+\frac{4D_{k}}{\mu^{3}}\left(e^{-\mu(t-t_{1})}-1\right)\right]~{}.$
(30)

It is clear by inspection that the second moment exists if and only if
${D_{k}}/{\mu^{2}}<{\bar{k}}/{2}$. Note that this is a stricter condition
compared to that found in Sec. 3 for the exponential decay of the Green’s
function, suggesting the existence of parameter regions for which the steady
state exists but not the variance. We now write the double exponential term in
Eq. (30) as a power series,

$\exp\left[\frac{4D_{k}}{\mu^{3}}e^{-\mu(t-t_{1})}\right]=\sum_{\ell=0}^{\infty}\frac{1}{\ell!}\left(\frac{4D_{k}}{\mu^{3}}\right)^{\ell}e^{-\mu\ell(t-t_{1})}~{}.$
(31)

Substituting back into (30), swapping integral and sum, performing the simple
exponential integral and rearranging terms, we eventually arrive at the
expression

$\langle
x^{2}\rangle=\frac{2D_{x}}{\mu}e^{-\xi}\sum_{\ell=0}^{\infty}\frac{\xi^{\ell}}{\ell!(\sigma-\xi+\ell)}$
(32)

with $\xi={4D_{k}}/{\mu^{3}}$ and $\sigma={2\bar{k}}/{\mu}$, which reduces to
$\langle x^{2}\rangle=D_{x}/\bar{k}$ for $D_{k}=0$, as expected of the
standard OU process. This analytical result is plotted against numerical
simulation in Fig. 3, showing good agreement. Formally, the right hand side of
Eq. (32) can be written more compactly in terms of the lower incomplete Gamma
function $\gamma(a,b)$, which has the following series expansion [45]

$\gamma(a,b)=b^{a}\sum_{\ell=0}^{\infty}\frac{(-b)^{\ell}}{\ell!(a+\ell)}~{},$
(33)

allowing us to reduce Eq. (32) to $\langle
x^{2}\rangle=2D_{x}e^{-\xi}(-\xi)^{-\sigma+\xi}\gamma(\sigma-\xi,-\xi)$ or, in
the original notation,

$\langle
x^{2}\rangle=2D_{x}\left(-\frac{4D_{k}}{\mu^{3}}\right)^{\frac{4}{\mu}\left(\frac{D_{k}}{\mu^{2}}-\frac{\bar{k}}{2}\right)}\exp\left[-\frac{4D_{k}}{\mu^{3}}\right]\gamma\left(\frac{4}{\mu}\left(\frac{\bar{k}}{2}-\frac{D_{k}}{\mu^{2}}\right),-\frac{4D_{k}}{\mu^{3}}\right)$
(34)

Figure 3: Comparison between analytical results and numerical simulations for
the second and quartic steady-state moments of the marginal probability
density function for the coordinate $x$ of the OU2 process, showing good
agreement between the two. The analytical result for the second moment
$\langle x^{2}\rangle$, given by Eq. (32) is defined for
$D_{k}/\mu^{2}<\bar{k}/2$, while that for the quartic moment $\langle
x^{4}\rangle$, given by Eq. (37) is defined for $D_{k}/\mu^{2}<\bar{k}/4$.
Here, we set $D_{x}=1$ and $\bar{k}=1$. Grey vertical lines indicate the
predicted radius of convergence.

###  4.2 Quartic moment ($n=2$)

In this case the right-hand side of (26) becomes, for $n=2$ and using (27),

$\displaystyle\langle x^{4}\rangle=6(2D_{x})^{2}$
$\displaystyle\int_{-\infty}^{t}dt_{1}\int_{t_{1}}^{t}dt_{2}\
e^{-2[\bar{k}(t-t_{1})+\bar{k}(t-t_{2})]}$ $\displaystyle\times\exp\Bigg{[}$
$\displaystyle\frac{4D_{k}}{\mu^{3}}\left(2e^{-\mu(t-t_{1})}+2e^{-\mu(t-t_{2})}-e^{-\mu(t_{2}-t_{1})}-3+\mu(t-t_{1})+3\mu(t-t_{2})\right)\Bigg{]}~{}.$
(35)

Expanding once again the double exponential as power series we then obtain

$\displaystyle\langle x^{4}\rangle=$ $\displaystyle
6(2D_{x})^{2}e^{-\frac{12D_{k}}{\mu^{3}}}\sum_{n,m,\ell=0}^{\infty}\frac{\left(-\frac{1}{2}\right)^{k}\left(\frac{8D_{k}}{\mu^{3}}\right)^{m+\ell+k}}{m!k!\ell!}$
$\displaystyle\times\int_{-\infty}^{t}dt_{1}\int_{t_{1}}^{t}dt_{2}\
\exp\Bigg{[}-2\left[\bar{k}-\frac{2D_{k}}{\mu^{2}}(t-t_{1})-2(\bar{k}-\frac{6D_{k}}{\mu^{2}})(t-t_{2})\right]$
$\displaystyle\qquad\qquad\qquad\quad\qquad-\mu
m(t-t_{1})-\mu\ell(t-t_{2})-\mu k(t_{2}-t_{1})\Bigg{]}~{},$ (36)

indicating that the quartic moment exists if and only if
${D_{k}}/{\mu^{2}}<{\bar{k}}/{4}$. The double integral can now be performed in
closed form,

$\langle
x^{4}\rangle=24D_{x}^{2}e^{-\frac{12D_{k}}{\mu^{3}}}\sum_{k,m,\ell=0}^{\infty}\frac{\left(-\frac{1}{2}\right)^{k}\left(\frac{8D_{k}}{\mu^{3}}\right)^{m+\ell+k}}{m!k!\ell!\left(\frac{4D_{k}}{\mu^{2}}-2\bar{k}-\mu(m+k)\right)\left(\frac{16D_{k}}{\mu^{2}}-4\bar{k}-\mu(m+\ell)\right)}~{},$
(37)

giving us the most compact expression for the quartic moment. This is plotted
against numerical simulations in Fig. 3, showing good agreement in the domain
of convergence.

###  4.3 Higher values of $n$

For the general case $n>2$, we focus on determining the criteria for
convergence of the moments. As shown in B, one can generalise the arguments
developed above and obtain the following criteria of convergence for the
moment of order $2n$

$\frac{D_{k}}{\mu^{2}}<\frac{\bar{k}}{2n}$ (38)

which is in agreement with the results we just obtained for the particular
cases $n=1$ and $n=2$.

Figure 4: Marginal probability density functions for the coordinate $x$ of the
OU2 process obtained from numerical simulations. The asymptotic scaling of the
empirical probability is studied as a function of the dimensionless parameter
combination $D_{k}/\mu^{2}$ controlling the convergence of the marginal
moments, as studied in Sec. 4, showing good agreement with the exponents
stated in Eq. (39), solid lines. For $D_{k}/\mu^{2}=0$, corresponding to the
standard OU limit, the distribution is a simple Gaussian (dashed gray line).
For any finite $D_{k}/\mu^{2}>0$, we instead observe a transition to algebraic
scaling. Here, we work in units such that $D_{x}=\bar{k}=1$ without loss of
generality.

From this criterion of convergence, we can argue for the asymptotic scaling of
the marginal probability density function $P(x)$. Indeed, assuming that the
asymptotic scaling exponent is a continuous function of $D_{k}/\mu^{2}$, we
expect that the value of $D_{k}/\mu^{2}$ at which the moment of order $2n$
becomes divergent, namely $D_{k}/\mu^{2}=\bar{k}/(2n)$, corresponds to that at
which the marginal probability density scales asymptotically as $x^{-2n-1}$ at
large $x$. We conclude that, asymptotically,

$P(x)\sim|x|^{-1-\frac{\bar{k}\mu^{2}}{D_{k}}}~{}.$ (39)

This result in agreement with the closed form expression for $P(x)$ derived
analytically in Sec. 5 below by means of homogenisation in the fast-slow
limit, as well as with numerical simulations for the full model, as shown in
Fig. 4.

##  5 Fast stiffness limit

In many applications, such as the optical tweezer example discussed in the
introduction, the potential stiffness fluctuations can be reasonably assumed
to occur on a comparatively fast timescale. In this regime, analytical results
for the marginal probability density function $P(x)$ can be obtained by
enforcing a formal separation of timescales between the slow dynamics of the
particle position $x(t)$ and the fast dynamics of the stiffness fluctuations
$k(t)$. The elimination of the fast stiffness dynamics can subsequently be
carried out following a multiscale approach [44]. In the following, we
consider two fast-slow regimes: (i) a naïve adiabatic limit, where the
characteristic timescale of the $k$ dynamics is sent to zero at fixed variance
$\sigma_{k}^{2}=D_{k}/\mu$, and (ii) a nontrivial limit, where the variance
$\sigma_{k}^{2}$ diverges as the inverse of the characteristic timescale,
leading to $k(t)$ in Eq. (2a) becoming statistically equivalent to a Gaussian
white noise.

###  5.1 Naïve adiabatic limit

First, we consider a naïve adiabatic limit, in which stiffness fluctuations
are expected to be irrelevant. Formally, we introduce a real dimensionless
coefficient $\varepsilon$ and proceed to the rescaling
$\mu\to\tilde{\mu}/\varepsilon^{2}$ and
$D_{k}\to\tilde{D}_{k}/\varepsilon^{2}$ in Eq. (2b) which governs the dynamics
of $k(t)$. We then take the limit $\varepsilon\to 0$, keeping the variance of
the stiffness $\sigma_{k}^{2}=D_{k}/\mu=\tilde{D}_{k}/\tilde{\mu}$ constant.
In this limit, one finds that the effective dynamics of the slow variable
$x_{\varepsilon}$ are obtained by replacing $k(t)\to 0$ in Eq. (2b) by its
mean value. Consequently, $x_{\varepsilon}$ is shown to be governed by the OU
dynamics,

$\dot{x}_{\varepsilon}(t)=-\bar{k}x_{\varepsilon}(t)+\sqrt{2D_{x}}\zeta_{x}(t),$
(40)

A rigorous derivation of this result is presented in C.2. We conclude that in
this trivial limit, we do not retain any signature of the stiffness
fluctuations.

###  5.2 White noise limit

Intuitively, we can go beyond this first trivial limit by replacing the
Gaussian process $k(t)$ not by its mean value but by a Gaussian white noise
with appropriate mean and standard deviation. Formally, this second regime is
obtained by performing the alternative rescaling
$\mu\to\tilde{\mu}/\varepsilon^{2}$ and
$D_{k}\to\tilde{D}_{k}/\varepsilon^{4}$, before taking the limit
$\varepsilon\to 0$ which keeps the ratio $D_{k}/\mu^{2}$ constant. In
situations where $x(t)$ denotes the position of an overdamped particle, this
regime can be understood physically as a low viscosity, high temperature
limit. Indeed, given an effective friction coefficient $\gamma$ and bath
temperature $T$, we have by the Stokes-Einstein relation that
$\mu\propto\gamma^{-1}$ while $D_{x}\propto\gamma^{-1}T$. Taking
$\gamma=\tilde{\gamma}\varepsilon^{2}$ and $T=\tilde{T}\varepsilon^{-2}$
produces the desired rescaling.

Mathematically, this amounts to $k(t)$ in Eq. (2b) becoming statistically
equivalent to a Gaussian white noise with covariance $\langle
k(t)k(t^{\prime})\rangle=2\tilde{D}_{k}/\tilde{\mu}^{2}\delta(t-t^{\prime})$.
Importantly, the term $k(t)x(t)$ appearing when integrating Eq. (2a) should
now be treated as a Stratonovich product [44, 46]. Taking care of the
Stratonovich-to-Itô conversion [60], we find that the Itô form of the
resulting Langevin equation in the limit $\varepsilon\to 0$ reads

$\dot{x}_{\varepsilon}(t)=-\left(\bar{k}-\frac{D_{k}}{\mu^{2}}\right)x_{\varepsilon}(t)+\sqrt{2\left(D_{x}+\frac{D_{k}x_{\varepsilon}^{2}(t)}{\mu^{2}}\right)}\zeta_{x}(t)$
(41)

where we have used that ${\tilde{D}_{k}}/{\tilde{\mu}^{2}}={D_{k}}/{\mu^{2}}$
in this case. Interestingly, this shows a clear instability as
$D_{k}/\mu^{2}>\bar{k}$ due to a renormalisation of the confining potential
stiffness and, unlike the original dynamics, is characterised by
multiplicative noise. The Fokker-Planck representation of Eq. (41),

$\partial_{t}P(x_{\varepsilon},t)=\partial_{x_{\varepsilon}}\left\\{\partial_{x_{\varepsilon}}\left[\left(D_{x}+\frac{D_{k}x_{\varepsilon}^{2}}{\mu^{2}}\right)P(x_{\varepsilon},t)\right]+\left(\bar{k}-\frac{D_{k}}{\mu^{2}}\right)x_{\varepsilon}P(x_{\varepsilon},t)\right\\}~{},$
(42)

can equivalently be derived by multiscale methods (see C.2). Interestingly,
Eq. (42) can be mapped onto an associated Legendre differential equation, see
D. In the rest of this subsection we drop the subscript of $\varepsilon$ and
define the shorthands $h\equiv D_{k}/\mu^{2}$ and $\kappa\equiv\bar{k}-h$ for
the sake of simplicity.

We now proceed to determining the steady state probability density function
$P(x)$ associated with the Langevin equation (41) for the slow $x$ dynamics.
To do so, we introduce a variable $z(x)$ whose stochastic dynamics do not
involve multiplicative noise [47]

$z(x)=\int^{x}d\xi\left(D_{x}+\frac{D_{k}{\xi}^{2}}{\mu^{2}}\right)^{-\frac{1}{2}}=\frac{1}{\sqrt{h}}\tanh^{-1}\left(\frac{\sqrt{h}x}{\sqrt{D_{x}+hx^{2}}}\right).$
(43)

By Itô’s lemma, the dynamics for $z$ take the form

$\displaystyle\dot{z}(t)=-\frac{(\kappa+h)x}{\sqrt{(D_{x}+hx^{2})}}+\sqrt{2}\zeta_{x}(t)=-\frac{(\kappa+h)\tanh\left[\sqrt{h}z(t)\right]}{\sqrt{h}}+\sqrt{2}\zeta_{x}(t)$
(44)

where the noise is now additive as we anticipated.

The dynamics for $z$ are exactly those of a passive Brownian particle in a
static potential

$V(z)=\left(\frac{\kappa+h}{h}\right)\ln\left[\cosh\left(\sqrt{h}z\right)\right]~{},$
(45)

whence the steady state probability distribution for $z$ is given by the
Boltzmann measure

$P_{Z}(z)=e^{-V(z)}/\mathcal{Z},\quad\mathcal{Z}=\sqrt{\frac{\pi}{h}}\frac{\Gamma[(h+\kappa)/2h]}{\Gamma[(2h+\kappa)/2h]}~{}.$
(46)

Finally, we perform a transformation of probability distributions to obtain a
simple expression for the probability density of $x$,

$P_{X}(x)=\frac{dz}{dx}P_{Z}(z(x))=\sqrt{\frac{h}{\pi}}\Gamma\left(\frac{\kappa}{2h}+1\right)\left[\Gamma\left(\frac{h+\kappa}{2h}\right)\right]^{-1}D_{x}^{\frac{h+\kappa}{2h}}\left(D_{x}+hx^{2}\right)^{-\frac{\kappa}{2h}-1}~{}.$
(47)

Asymptotically, we again find that

$P_{X}(x)\sim|x|^{-\bar{k}\mu^{2}/D_{k}-1}$ (48)

It can be checked by direct substitution that $P(x)=P_{X}(x)$ solves the
Fokker-Planck equation (42) at steady state. Notice also that Eq. (47) is in
agreement with the asymptotic scaling of the marginal probability density for
the full OU2 process, Eq. (39). In particular, it is straightforward to check
that the existence of moments of order $2n$ demands
$D_{k}/\mu^{2}<\bar{k}/(2n)$.

##  6 Statistics of maxima

Despite their relevance in many domains, such as climate modelling, there is
currently no general framework to study the extreme value statistics (EVS) of
correlated random variables [48]. Amongst the few exceptions is the standard
OU process, for which the EVS can be compute exactly and can be shown to
belong to the Gumbel universality class [49]. Accordingly, the asymptotic
probability density of the maximum
$X(t)\equiv\,\mathrel{\mathop{\max}\limits_{\tau\in[0,t]}}\\{x(\tau)\\}$ is
given by

$\Phi_{G}(X;m,s)=s^{-1}e^{-(z+e^{-z})}~{},\qquad z=\frac{X-m}{s}$ (49)

where the first two moments of $X(t)$ can be expressed as

$\langle X\rangle=m+s\gamma,\quad\langle X^{2}\rangle=\frac{\pi^{2}s^{2}}{6}$
(50)

respectively, where $\gamma=0.5772...$ is the Euler-Mascheroni constant. It is
thus natural to wonder how the EVS of the OU2 process compare to those of the
standard OU process.

While a fully analytical approach is beyond the scope of this work, we can
draw on the renormalisation group heuristic introduced in Ref. [48] to argue
that, as long as correlations decay over a finite time, the EVS for a weakly
correlated stochastic process are still expected to converge to one of the
three limiting distributions for uncorrelated random variables, namely Gumbel,
for exponentially decaying parent distributions, Fréchet, for fat-tailed
parent distributions, and Weibull, for parent distributions with compact
support [48, 49]. Combining this heuristic with the finding of Sec. 4 that the
probability density of $x$ decays algebraically at large $x$, specifically as
$-\ln P(x)\sim 1+\bar{k}\mu^{2}/D_{k}\equiv 1+\alpha$, leads us to conjecture
that the EVS for the OU2 process should converge to a Fréchet distribution
with $D_{k}$-dependent characteristic exponent. In particular, we consider the
Fréchet probability density for
$X(t)=\,\mathrel{\mathop{\max}\limits_{\tau\in[0,t]}}\\{x(\tau)\\}$ given by

$\Phi_{F}(X;m,s)=\begin{cases}0&\text{for}\quad z\leq m\\\
\frac{\alpha}{s}z^{-1-\alpha}e^{-z^{-\alpha}}&\text{for}\quad
z>m\end{cases},\qquad\mbox{with}\quad z=\frac{X-m}{s}$ (51)

where the first two moments of $X(t)$ are expressed as

$\langle X\rangle=m+s\Gamma\left(1-\frac{1}{\alpha}\right),\quad\langle
X^{2}\rangle=s^{2}\left[\Gamma\left(1-\frac{2}{\alpha}\right)-\Gamma^{2}\left(1-\frac{1}{\alpha}\right)^{2}\right]$
(52)

the latter being defined only for $\alpha>2$. We find this conjecture to be in
good agreement with numerical simulations of the full model, as shown in Fig.
5. As expected, Gumbel EVS are recovered upon setting $D_{k}=0$, i.e. in the
absence of potential fluctuations.

Figure 5: Extreme value statistics for the OU2 process. In the absence of
potential fluctuations, $D_{k}=0$, we recover the standard OU process, for
which the limiting distribution of the maximum is known to be of the Gumbel
type, Eq. (49) (a). In the presence of potential fluctuations, $D_{k}>0$, we
conjecture a transition from Gumbel to Fréchet universality class, Eq. (51),
with a non-universal exponent $\alpha\equiv\bar{k}\mu^{2}/D_{k}$ (b). In both
cases, we find good agreement between numerical and analytical standardised
distributions.

##  7 First-passage time statistics

We now focus on the impact of the continuous fluctuations of the potential
stiffness present on the first-passage time statistics of the standard
Ornstein-Uhlenbeck process, as characterised in [21, 50]. The problem of
finding the mean first-passage times to an absorbing boundary can be mapped
onto that of the escape of particle over a fluctuating energy barrier. The
rate at which a Brownian particle escapes over an energy barrier is a central
problem in statistical physics dating back to Kramers [51], finding
applications across disciplines through reaction-rate theory [52].

Calculating the escape rate over a fluctuating energy barrier for a Brownian
particle has attracted some attention in the past [53, 54, 55, 56]. At low
temperatures, zero-mean fluctuations in the energy barrier height lead to so-
called resonant activation and to a reduction of the mean first-passage time,
effectively aiding the escape process [53, 54, 55]; resonant activation has
been observed for general confining potentials [52, 56]. In general, solving
for the first-passage time distribution and its moments for the coupled
dynamics of the particle position and potential stiffness constitutes a
formidable task. To the best of our knowledge, only approximate results can be
derived in the limit where the timescale associated with the potential
stiffness fluctuations is negligible compared to that of the particle position
dynamics, i.e. in the fast stiffness limit described above. The analysis of
the high temperature case is particularly challenging: even for static
harmonic potentials Kramers’ theory has been shown to breakdown in this limit
[50]; the high temperature limit in the case of a fluctuating potential is an
open problem.

Here, we instead tackle this problem numerically. Specifically, we numerically
integrate Eq. (2) using the Euler-Maruyama method with timestep $dt=10^{-4}$.
For all results presented here (see Fig. 6), we simulate $m=10^{5}$
realisations of the coupled dynamics in which the particle is initialised at
$x_{0}=0$ and we place an absorbing boundary condition in $x_{a}>0$, fixing
$D_{x}=\bar{k}=1$ with $\mu=0.1$. We have confirmed that all realisations lead
to a finite first-passage time to the absorbing boundary condition. We probe a
wide range of values for the stiffness fluctuations strength, such that for
high enough values of $D_{k}/\mu^{2}$, the steady-state distribution for the
particle position may not exist. However, for all values of the stiffness
fluctuations strength $D_{k}/\mu^{2}$ studied here, we have checked that the
standard deviation of the mean-first passage time (mFPT) is finite and
independent of the number of realisations for $m$ sufficiently large ensuring
the convergence of the mFPT.

Figure 6: First-passage time statistics as a function of stiffness
fluctuations strength. (a) Mean first-passage time as a function of
fluctuation strength $D_{k}/\mu^{2}$ for various absorbing target locations
$x_{a}>0$. (b) Mean first-passage time normalised by the mean first-passage
time of the standard Ornstein-Uhlenbeck process (i.e. the $D_{k}/\mu^{2}\to 0$
limit). We observe that stiffness fluctuations can significantly reduce the
mean first-passage time to targets located far enough from the initial
conditions (here, in $x_{0}=0$). (c) Coefficient of variation for the first-
passage statistics showing non-monotonic behavior for all target locations.
(d) Empirical histograms of the first passage time normalised by its standard
deviation, showing an exponential decay of the first passage probability at
long times. All results presented here are obtained via the integration of Eq.
(2) using an Euler-Maruyama method with $D_{x}=\bar{k}=1$ over an ensemble of
$m=10^{5}$ realisations of the coupled dynamics

A first look at Fig. 6 shows that the mFPT independently of the target
location is generically a non-trivial non-monotonic function of the stiffness
fluctuation strength, $D_{k}/\mu^{2}$. In particular, we find that the mFPT to
the absorbing target can be significantly reduced (more than one order of
magnitude) by strong enough stiffness fluctuations $D_{k}/\mu^{2}$ for targets
which are not too close to the initial conditions, which is consistent with
previous studies. Interestingly, we show that instead stiffness fluctuations
can increase the mFPT compared to the standard Ornstein-Uhlenbeck case (i.e.
the zero stiffness fluctuations limit, $D_{k}/\mu^{2}\to 0$) if $x_{a}$ is
small. Said differently, for targets very close to the initial conditions,
stiffness fluctuations can be detrimental and a system with constant stiffness
$\bar{k}$ will instead be optimal. Furthermore, we observe that in all cases
the coefficient of variation of the first-passage times,
$\sigma_{\tau}/\bar{\tau}$, first strongly increases before reaching a maximum
at the same value of $D_{k}/\mu^{2}$; at large enough fluctuations strength,
the coefficient of variation decays monotonously with fluctuation strength for
all target locations. Interestingly, we find that the first-passage times are
exponentially distributed at both low and high values of $D_{k}/\mu^{2}$ for
the most distant absorbing target location.

##  8 Entropy production

We now consider the thermodynamic implications of the stiffness fluctuations
characterising the OU2 process. Indeed, we expect a non-zero rate of entropy
production at steady-state [57, 16, 17] as the system performs work to change
the stiffness. Furthermore, the existence of steady-state divergence-free
probability currents in the $(x,k)$-plane, as shown on Fig. 7, implies a
breaking of time-reversal symmetry. It is in this sense that we have
previously referred to the OU2 process as a minimal model of dissipative
confinement. As detailed in Ref. [16] and E, such dissipation can be written
in terms of the second moment of the steady-state marginal probability density
function for the coordinate $x$, which we computed in Sec. 4. We now rederive
this result using a shortcut and call upon the results of Sec. 4 to provide a
clearer picture of the thermodynamic properties of the OU2 process.

As a preliminary step, let us introduce the Fokker-Planck formulation [58] of
the Langevin dynamics (2a) and (2b),

$\partial_{t}P(x,k)=[D_{x}\partial_{x}^{2}+D_{k}\partial_{k}^{2}]P(x,k)+(\bar{k}+k)\partial_{x}[xP(x,k)]+\mu\partial_{k}[kP(x,k)]~{}.$
(53)

Multiplying both sites of Eq. (53) by $x^{2}$ and subsequently integrating
with respect to both $x$ and $k$, using integration by parts where necessary,
gives the remarkably simple relation

$\left\langle(k(t)+\bar{k})x^{2}(t)\right\rangle=D_{x}~{}.$ (54)

Note however that, when carrying out this procedure, one encounters the
following integral

$\displaystyle\int dxdk\ x^{2}(k+\bar{k})\partial_{x}[xP(x,k)]$
$\displaystyle=\bar{k}\int dxdk\ x^{2}\partial_{x}[xP(x,k)]$
$\displaystyle=\bar{k}\left[x^{3}\int
dkP(x,k)\right]_{x=-\infty}^{+\infty}-2\bar{k}\int dxdk\ x^{2}P(x,k)~{}.$ (55)

For the boundary term on the right-hand side to vanish, it is required that
the marginal probability density $P(x)\equiv\int dkP(x,k)$ decays faster than
$P(x)\sim x^{-3}$ as $x\to\infty$. This is also the condition for the
existence of the second moment, such that the validity of Eq. (54) is
contingent upon the existence of $\langle x^{2}\rangle$.

The mean rate of entropy production, denoted $\dot{S}_{i}$, is related to the
Jarzynski stochastic work [41, 20]

$W_{\tau}=\int_{0}^{\tau}dt\
\dot{x}\circ[(k(t)+\bar{k})x(t)]=\frac{1}{2}\int_{0}^{\tau}dt\ \dot{k}(t)\circ
x^{2}(t)$ (56)

via

$\dot{S}_{i}=\lim_{\tau\to\infty}\frac{1}{D_{x}}\langle
W_{\tau}\rangle=-\frac{\mu}{2D_{x}}\langle
kx^{2}\rangle=\frac{\mu\bar{k}}{2D_{x}}\left(\langle
x^{2}\rangle-\frac{D_{x}}{\bar{k}}\right)$ (57)

where $\circ$ denotes a Stratonovich product and we have used (54) to replace
$\langle kx^{2}\rangle$ in Eq. (57). Note that $\langle
x^{2}\rangle>D_{x}/\bar{k}$ for all $D_{k}>0$ [16], such that the second law
of thermodynamics $\dot{S}_{i}\geq 0$ is always satisfied. Now, using Eq. (32)
for the variance, we can write the steady-state entropy production rate as

$\dot{S}_{i}=\frac{\mu}{2}\left(se^{-\xi}\sum_{\ell=0}^{\infty}\frac{\xi^{\ell}}{\ell!(s-\xi+\ell)}-1\right)$
(58)

where again we have defined $\xi={4D_{k}}/{\mu^{3}}$ and $s={2\bar{k}}/{\mu}$.
Expanding to leading order in weak stiffness noise,
$\dot{S}_{i}={\mu\xi}/[{2s(1+s)}]+\mathcal{O}(\xi^{2})=D_{k}/[\bar{k}(\mu+\bar{k})]+\mathcal{O}(\xi^{2})$,
which vanishes at $D_{k}=0$. Remarkably, the entropy production rate diverges
together with $\langle x^{2}\rangle$ as $D_{k}/\mu^{2}$ approaches $\bar{k}/2$
from below.

Figure 7: Steady state probability density (contours) and out-of-plane
component of the curl of the steady-state probability current. The existence
of steady state divergence-free currents implies time-reversal symmetry
breaking and thus a non-zero entropy production, which we calculate
analytically in Section 8.

##  9 Conclusion

We have examined a number of statistical, dynamic and thermodynamic properties
of the OU2 process: a generalised Ornstein-Uhlenbeck (OU) process obtained by
allowing for the associated stiffness coefficient to undergo OU-like
stochastic fluctuations around a positive mean. To the best of our knowledge,
this is the first systematic exploration of such a model, which was originally
introduced in two recent works by some of the authors on the thermodynamics of
Brownian motion in fluctuating potentials [16, 17]. This process finds
physical relevance in contexts where effective harmonic confining potentials
are generated by a non-ideal contextual processes, e.g. when colloid are
manipulated by optical tweezers with realistic laser power stability [12, 13].
We also argued that the OU2 model is closely related to models in stochastic
search with resetting [23, 16, 24, 25, 26] and active matter [42, 43, 17], and
that it constitutes a stochastic counterpart to the breathing harmonic
potentials studied in the thermodynamics literature [41, 39, 40, 5].

We started out analysis by considering the limit of vanishing positional
noise, for which the time-dependent probability density can be obtained in
closed form, Eq. (8). From this, we derived exact expressions for the median
(Eq. (11)), mode (Eq. (13)) and mean (Eq. (12)), whose time-dependence was
found to transition from exponentially decreasing (trapped regime) to
exponentially increasing (growing regime) at different non-trivial critical
values of the non-dimensional parameter $\mathcal{F}\equiv
D_{k}/(\mu^{2}\bar{k})$.

Having reintroduced a finite positional noise, we computed the conditional and
full Green’s function of the process, Eqs. (18) and (21), showing that both
grow exponentially with time when $\mathcal{F}$ exceeds unity, indicating loss
of ergodicity. Starting from the formal solution of the OU2 dynamics, Eq.
(14), we subsequently computed the second and fourth moments of the steady-
state positional probability density function, Eqs. (32) and (37), as well as
a necessary condition for the existence of moments of order $2n$ in terms of
an n-dependent upper bound on $\mathcal{F}$, Eq. (38). From this condition we
inferred an algebraic asymptotic decay of the positional probability density,
with scaling exponent $\eta=-1-\mathcal{F}^{-1}$. We subsequently considered
two limiting regimes of fast stiffness dynamics: while the naïve adiabatic
limit produced trivial OU dynamics for the slow degree of freedom, the second
limit, obtained by taking $\mu\to\tilde{\mu}/\epsilon^{2}$ and
$D_{k}\to\tilde{D}_{k}/\epsilon^{4}$ with $\epsilon\to 0$, led to non-trivial
coarse grained dynamics for the position involving a renormalised stiffness
and multiplicative noise, Eq. (41), for which the steady-state probability
density was obtained in closed form, Eq. (47).

Borrowing from extreme value theory heuristics, we then conjectured that in
the presence of finite stiffness fluctuations, $\mathcal{F}>0$, the
standardised distribution of the running maximum should converge at long times
to a Frechet form with $\mathcal{F}$-dependent exponent, in good agreement
with numerical experiments (cf. the standard OU process, whose maximum is
Gumbel distributed at long times). Further, the dependence on $\mathcal{F}$ of
the mean first passage time to a positive target was studied numerically and
we show that sufficiently strong stiffness fluctuations can aid the particle
by reducing drastically the mean first-passage time to that target. A formal
analytical treatment of this problem remains an intriguing open question and
is left for future studies.

Finally, we presented a compact derivation of the steady-state entropy
production rate which calls upon results for Jarzinsky’s stochastic work, Eq.
(58), showing that it depends solely on the second steady-state moment of the
positional probability density. Remarkably, the entropy production diverges
for $\mathcal{F}>1/2$. It would be interesting to explore higher order
statistics of the stochastic work (56), similarly to what was done in Ref.
[19] for the stochastically sliding potential and more recently in Ref. [59]
for an AOU particle confined in a harmonic potential, and to compare any such
result to the full distribution of the stochastic work obtained in Ref. [41]
for a breathing harmonic potential in the slow driving regime.

Taken together, our results point to the OU2 process as a widely applicable
minimal model of dissipative confinement. It’s rich phenomenology, emerging
from the dynamical establishment of a non-equilibrium steady-state analogous
to that of Brownian motion under stochastic resetting, renders it a valuable
non-motile counterpart to other minimal models of single-particle out-of-
equilibrium dynamics, such as the AOU particle, that have been explored
extensively in recent years.

## References

  * [1] Uhlenbeck G E and Ornstein L S 1930 Phys. Rev. 36(5) 823–841
  * [2] Einstein A et al. 1905 Annalen der physik 17 208
  * [3] Le Bellac M and Barton G 1991 Quantum and statistical field theory (Oxford: Clarendon Press)
  * [4] Hohenberg P C and Halperin B I 1977 Rev. Mod. Phys. 49(3) 435–479
  * [5] Martínez I A, Roldán É, Dinis L and Rica R A 2017 Soft matter 13 22–36
  * [6] Saha T K, Lucero J N, Ehrich J, Sivak D A and Bechhoefer J 2022 Physical Review Letters 129 130601
  * [7] Saha T K, Lucero J N, Ehrich J, Sivak D A and Bechhoefer J 2021 Proceedings of the National Academy of Sciences 118 e2023356118
  * [8] Saha T K, Ehrich J, Gavrilov M c v, Still S, Sivak D A and Bechhoefer J 2023 Phys. Rev. Lett. 131(5) 057101
  * [9] Tal-Friedman O, Pal A, Sekhon A, Reuveni S and Roichman Y 2020 The journal of physical chemistry letters 11 7350–7355
  * [10] Ashkin A, Dziedzic J M, Bjorkholm J E and Chu S 1986 Optics letters 11 288–290
  * [11] Dufresne E R, Spalding G C, Dearing M T, Sheets S A and Grier D G 2001 Review of Scientific Instruments 72 1810–1816
  * [12] Bustamante C J, Chemla Y R, Liu S and Wang M D 2021 Nature Reviews Methods Primers 1 25
  * [13] Pesce G, Jones P H, Maragò O M and Volpe G 2020 The European Physical Journal Plus 135 1–38
  * [14] Chiang R, Liu P and Okunev J 1995 Journal of Banking & Finance 19 1327–1340
  * [15] Bartoszek K, Glémin S, Kaj I and Lascoux M 2017 Journal of theoretical biology 429 35–45
  * [16] Alston H, Cocconi L and Bertrand T 2022 Journal of Physics A: Mathematical and Theoretical 55 274004
  * [17] Cocconi L, Alston H and Bertrand T 2023 Phys. Rev. Res. 5(4) 043032
  * [18] Cox J C, Ingersoll J E and Ross S A 1985 Econometrica 53 385–407 ISSN 00129682, 14680262
  * [19] Manikandan S K and Krishnamurthy S 2017 The European Physical Journal B 90 1–19
  * [20] Lee H K, Kwon Y and Kwon C 2023 Journal of the Korean Physical Society 83 331–336
  * [21] Santra I, Das S and Nath S K 2021 Journal of Physics A: Mathematical and Theoretical 54 334001
  * [22] Zhang Z and Wang W 2017 Communications in Statistics-Simulation and Computation 46 4783–4794
  * [23] Evans M R, Majumdar S N and Schehr G 2020 Journal of Physics A: Mathematical and Theoretical 53 193001
  * [24] Gupta D, Plata C A, Kundu A and Pal A 2020 Journal of Physics A: Mathematical and Theoretical 54 025003
  * [25] Gupta D, Plata C A and Pal A 2020 Physical review letters 124 110608
  * [26] Jerez M J Y, Bonachita M A and Confesor M N P 2021 Physical Review E 104 044609
  * [27] Guyon X, Iovleff S and Yao J F 2004 ESAIM: Probability and Statistics 8 25–35
  * [28] Kesten H 1973 Acta Mathematica 131 207–248
  * [29] Morita S 2016 Europhysics Letters 113 40007
  * [30] Morita S 2018 Journal of the Physical Society of Japan 87 053801
  * [31] Burkhardt T W 2000 Journal of Physics A: Mathematical and General 33 L429
  * [32] Majumdar S N and Bray A J 2001 Physical Review Letters 86 3700
  * [33] Majumdar S N, Rosso A and Zoia A 2010 Journal of Physics A: Mathematical and Theoretical 43 115001
  * [34] Boutcheng H J O, Bouetou T B, Burkhardt T W, Rosso A, Zoia A and Crepin K T 2016 Journal of Statistical Mechanics: Theory and Experiment 2016 053213
  * [35] Singh P 2020 Journal of Physics A: Mathematical and Theoretical 53 405005
  * [36] Chubynsky M V and Slater G W 2014 Phys. Rev. Lett. 113(9) 098302
  * [37] Banks D S, Tressler C, Peters R D, Höfling F and Fradin C 2016 Soft matter 12 4190–4203
  * [38] Sposini V, Chechkin A and Metzler R 2018 Journal of Physics A: Mathematical and Theoretical 52 04LT01
  * [39] Blickle V and Bechinger C 2012 Nature Physics 8 143–146
  * [40] Martínez I A, Roldán É, Dinis L, Petrov D, Parrondo J M and Rica R A 2016 Nature physics 12 67–70
  * [41] Speck T 2011 Journal of Physics A: Mathematical and Theoretical 44 305001
  * [42] Martin D, O’Byrne J, Cates M E, Fodor E, Nardini C, Tailleur J and van Wijland F 2021 Phys. Rev. E 103(3) 032607
  * [43] Bothe M and Pruessner G 2021 Phys. Rev. E 103(6) 062105
  * [44] Pavliotis G and Stuart A 2008 Multiscale methods: averaging and homogenization (Springer Science & Business Media)
  * [45] Arfken G B, Weber H J and Harris F E 2011 Mathematical methods for physicists: a comprehensive guide (Academic press)
  * [46] de Pirey T A, Cugliandolo L F, Lecomte V and van Wijland F 2023 Advances in Physics 1–85
  * [47] Rubin K J, Pruessner G and Pavliotis G A 2014 Journal of Physics A: Mathematical and Theoretical 47 195001
  * [48] Majumdar S N, Pal A and Schehr G 2020 Physics Reports 840 1–32
  * [49] Haan L and Ferreira A 2006 Extreme value theory: an introduction vol 3 (Springer)
  * [50] Abkenar M, Gray T H and Zaccone A 2017 Phys. Rev. E 95(4) 042413
  * [51] Kramers H 1940 Physica 7 284–304 ISSN 0031-8914
  * [52] Hänggi P, Talkner P and Borkovec M 1990 Rev. Mod. Phys. 62(2) 251–341
  * [53] Reimann P 1995 Phys. Rev. E 52(2) 1579–1600
  * [54] Reimann P 1995 Phys. Rev. Lett. 74(23) 4576–4579
  * [55] Stein D L, Doering C R, Plamer R G, van Hemmen J L and McLaughlin R M 1990 Journal of Physics A: Mathematical and General 23 L203
  * [56] Doering C R and Gadoua J C 1992 Phys. Rev. Lett. 69(16) 2318–2321
  * [57] Cocconi L, Garcia-Millan R, Zhen Z, Buturca B and Pruessner G 2020 Entropy 22 1252
  * [58] Risken H and Risken H 1996 Fokker-planck equation (Springer)
  * [59] Semeraro M, Gonnella G, Suma A and Zamparo M 2023 Phys. Rev. Lett. 131(15) 158302
  * [60] Gardiner C W et al. 1985 Handbook of stochastic methods vol 3 (springer Berlin)
  * [61] Abramowitz M, Stegun I A and Romer R H 1988 Handbook of mathematical functions with formulas, graphs, and mathematical tables

##  Appendix A Key results for the Ornstein-Uhlenbeck process

In this appendix, we recall some useful key results for the original Ornstein-
Uhlenbeck process as described by Eq. (1). First, the transition probability
$P(x,t|x^{\prime},t^{\prime})$ for the process takes the form

$P(x,t|x^{\prime},t^{\prime})=\sqrt{\frac{\bar{k}}{2\pi
D_{x}(1-e^{-2\bar{k}(t-t^{\prime})})}}\exp\left[-\frac{\bar{k}(x-x^{\prime}e^{-\bar{k}(t-t^{\prime})})^{2}}{2D_{x}(1-e^{-2\bar{k}(t-t^{\prime})})}\right],$
(59)

where $t^{\prime}<t$. The steady-state probability density function is thus a
Gaussian distribution with zero mean and variance, $\langle
x^{2}\rangle=D_{x}/\bar{k}$ [60]. The exponential decay at large $x$ of the
steady-state distribution ensures that all even moments $\langle
x^{2n}\rangle$ are finite for $n\in\mathbb{Z}_{\geq 0}$, with odd moments
vanishing due to the $x\rightarrow-x$ symmetry. The associated cumulant
generating function then takes the form

$K(t)=\log\langle e^{ax}\rangle=\frac{a^{2}D_{x}}{2\bar{k}}$ (60)

where cumulants of order 3 and above vanish. The solution $x(t)$ for a given
realisation of the noise can be written as

$x(t)=\sqrt{2D_{x}}\int_{-\infty}^{t}dt^{\prime}\>\zeta(t^{\prime})\exp\left[-\bar{k}(t-t^{\prime})\right],$
(61)

from which we read off the Green’s function for the process $\mathcal{G}(t)$
of the form

$\mathcal{G}(t)=\exp\left[-\bar{k}t\right]\Theta(t)$ (62)

where $\Theta(t)$ is the Heaviside function.

##  Appendix B Convergence criterion for higher moments

In this section, we discuss the derivation of the convergence criterion for
moments of order $n>2$. Looking at Eq. (27), we argue that since the function
is finite for all $(t_{1},...,t_{n})$, any divergence of the integral in the
right-hand side of (25) must be controlled by the behaviour of the integrand
in the regime $t_{i}\to-\infty$. Keeping only leading order terms in this
limit, we thus rewrite (26) using (27) as

$\displaystyle\left\langle\exp\left[-2\sum_{i=1}^{n}\int_{t_{i}}^{t}dt_{i}^{\prime}\
k(t_{i}^{\prime})\right]\right\rangle$
$\displaystyle\simeq\exp\left[\frac{4D_{k}}{\mu^{2}}\left(\sum_{i=1}^{n}(t-t_{i})+2\sum_{i<j}(t-t_{j})\right)\right]$
$\displaystyle=\exp\left[\frac{4D_{k}}{\mu^{2}}\sum_{i=1}^{n}(1+2(i-1))(t-t_{i})\right]$
(63)

where $\simeq$ denotes approximate equality at $t_{i}\to-\infty$. Combining
this with (25), we argue that the moment of order $2n$ converges when the
following integral also converges

$\mathcal{I}_{n}=\int
dt_{1<...<n}\exp\left[-\sum_{i=1}^{n}\left(2\bar{k}-\frac{4D_{k}}{\mu^{2}}(1+2(i-1))\right)(t-t_{i})\right]~{},$
(64)

where $\int
dt_{1<...<n}\equiv\int_{-\infty}^{t}dt_{1}\dots\int_{t_{n-1}}^{t}dt_{n}$. We
now define $a_{i}=2\bar{k}-{4D_{k}}[1+2(i-1)]/{\mu^{2}}$ and re-write the
integral (64) as

$\mathcal{I}_{n}=\int dt_{1<...<n}\prod_{i=1}^{n}e^{-a_{i}(t-t_{i})}~{}.$ (65)

We now call upon the following result for generic $m$ iteratively

$\int_{t_{m-1}}^{t}dt_{m}e^{-a(t-t_{m})}=\frac{1-e^{-a(t-t_{m-1})}}{a}$ (66)

to evaluate $\mathcal{I}_{n}$ up to a multiplicative constant. We integrate
over $dt_{n}$ and re-arrange to write

$\displaystyle\mathcal{I}_{n}$ $\displaystyle\propto\int
dt_{1<...<(n-1)}\prod_{i=1}^{n-1}e^{-a_{i}(t-t_{i})}\left[1-e^{-a_{n}(t-t_{n-1})}\right]$
(67) $\displaystyle=\int
dt_{1<...<(n-1)}\prod_{i=1}^{n-2}e^{-a_{i}(t-t_{i})}\left[e^{-a_{n-1}(t-t_{n-1})}-e^{-(a_{n-1}+a_{n})(t-t_{n-1})}\right].$
(68)

Integrating over $t_{n-1}$ we subsequently derive

$\displaystyle\mathcal{I}_{n}$ $\displaystyle\propto\int
dt_{1<...<(n-2)}\prod_{i=1}^{n-2}e^{-a_{i}(t-t_{i})}[$
$\displaystyle\lambda_{n-1}\left(1-e^{-a_{n-1}(t-t_{n-2})}\right)-1+e^{-(a_{n-1}+a_{n})(t-t_{n-2})}]$
$\displaystyle\propto\int
dt_{1<...<(n-2)}\prod_{i=1}^{n-3}e^{-a_{i}(t-t_{i})}[$
$\displaystyle(\lambda_{n-1}-1)e^{-a_{n-2}(t-t_{n-2})}-\lambda_{n-1}e^{-(a_{n-2}+a_{n-1})(t-t_{n-2})}$
$\displaystyle-e^{-(a_{n-2}+a_{n-1}+a_{n})(t-t_{n-2})}].$ (69)

where $\lambda_{n+1}=(a_{n-1}+a_{n})/a_{n-1}$.

Treating the remaining time variables in the same manner, we conclude that the
convergence of the integral $\mathcal{I}_{n}$ is determined by the positivity
of the partial sums $A_{m}=\sum_{i=1}^{m}a_{i}$ for $m=1\dots n.$ More
specifically, the most strict condition is always the positivity of the full
sum, $m=n$, as the sequence $\\{a_{i}\\}$ is monotonically decreasing:
$a_{i+1}<a_{i}$. The condition for the convergence of the $n$-th moment is
thus

$0<A_{n}=\sum_{i=1}^{n}\left[2\bar{k}-\frac{4D_{k}}{\mu^{2}}(1+2(i-1))\right]=2\bar{k}n-\frac{4D_{k}}{\mu^{2}}n^{2}\implies\frac{D_{k}}{\mu^{2}}<\frac{\bar{k}}{2n},$
(70)

in agreement with the results obtained in Section 4 for the particular cases
$n=1$ and $n=2$.

##  Appendix C Homogenisation procedure for coupled dynamics

In this appendix, we first review multiscale methods for the coarse graining
of the coupled dynamics of two stochastic variables with multiple timescales
and then apply it to the homogeneisation of the OU2 process.

###  C.1 General theory

In this section, we follow the treatment of the homogenisation procedure
presented in Chapter 11 of Ref. [44]. We begin from the most general form for
the coupled dynamics of two stochastic processes $x(t)$ and $k(t)$,

$\displaystyle\dot{x}$
$\displaystyle=\frac{1}{\varepsilon}f_{0}(x,k)+f_{1}(x,k)+\alpha(x,k)\zeta_{x}(t)$
(71) $\displaystyle\dot{k}$
$\displaystyle=\frac{1}{\varepsilon^{2}}g(x,k)+\frac{1}{\varepsilon}\beta(x,k)\zeta_{k}(t)$
(72)

where $\varepsilon>0$ is a dimensionless factor that will be taken to zero to
enforce a separation of timescales between the “slow” $x$ and “fast” $k$
dynamics. Note that both $\dot{x}$ and $\dot{k}$ involve fast contributions to
the dynamics, but the dynamics for $k$ are an order of $\varepsilon$ faster
than those of $x$. Combining Eqs. (71) and (72), we construct the backward
Kolmogorov equation for $v(x,k,t)$, which takes the form

$\partial_{t}v=\mathcal{L}v=\frac{1}{\varepsilon^{2}}\left(g\partial_{k}v+\frac{1}{2}B\partial_{k}^{2}v\right)+\frac{1}{\varepsilon}\left(f_{0}\partial_{x}v\right)+\left(f_{1}\partial_{x}v+\frac{1}{2}A\partial_{x}^{2}v\right),$
(73)

where we have defined the diffusivities $A(x,y)=\alpha(x,y)\bar{\alpha}(x,y)$
and $B(x,y)=\beta(x,y)\bar{\beta}(x,y)$, with $\bar{\bullet}$ denoting the
conjugate transpose. We write the order $\mathcal{O}(\varepsilon^{-2})$
contribution to the backward Fokker-Planck operator acting on $v$ as

$\mathcal{L}_{-2}v\equiv g\partial_{k}v+\frac{1}{2}B\partial_{k}^{2}v.$ (74)

Clearly, $\mathcal{L}_{-2}w(x)=0$ for any function $w(x)$ that does not depend
on $k$. Additionally, we define $\rho^{\infty}(k;x)$ as the normalised measure
obtained by solving the associated steady-state forward Kolmogorov equation
for a fixed value of the slow variable $x$,
$\mathcal{L}^{*}_{-2}\rho^{\infty}(k;x)=0$. For the limit $\varepsilon\to 0$
of this problem to be well-posed, we require that $f_{0}(x,k)$ satisfies the
so-called centering condition, namely that its average with respect to
$\rho^{\infty}$ vanishes

$\int dkf_{0}(x,k)\rho^{\infty}(k;x)=0.$ (75)

We now construct a perturbative solution to Eq. (73) of the form
$v=v_{0}+\varepsilon v_{1}+\varepsilon^{2}v_{2}+\mathcal{O}(\varepsilon^{3})$.
Matching terms of order $\mathcal{O}(\varepsilon^{-2})$, we conclude that
$v_{0}(x)$ is independent of $k$. At $\mathcal{O}(\varepsilon^{-1})$ and
relying on the centering condition on $f_{0}$ introduced above, we find an
expression for $v_{1}(x,k)$ in terms of $v_{0}(x)$ and a function $\Phi(x,k)$
which solves

$\displaystyle-\mathcal{L}_{-2}\Phi({x,k})=f_{0}(x,k)\quad\text{ with
}\quad\int dk\ \Phi(x,k)\rho^{\infty}(k;x)=0.$ (76)

Finally, at $\mathcal{O}(\varepsilon^{0})$ we derive a closed equation for
$\partial_{t}v_{0}$ from the ergodicity assumption that $\int
dk\rho^{\infty}\mathcal{L}_{-2}v_{2}=0$, namely

$\partial_{t}v_{0}=F(x)\partial_{x}v_{0}+\frac{1}{2}A(x)A(x)^{T}\partial_{x}^{2}v_{0}$
(77)

where we have defined the vector fields

$F(x)=\int dk\
\left(f_{1}(x,k)+(\partial_{x}\Phi(x,k))f_{0}(x,k)\right)\rho^{\infty}(k;x)$
(78)

and

$A(x)A^{T}(x)=A_{1}(x)+\frac{1}{2}(A_{0}(x)+A_{0}^{T}(x))~{},$ (79)

with

$\displaystyle A_{0}(x)$ $\displaystyle=2\int dk\
f_{0}(x,k)\Phi(x,k)\rho^{\infty}(k;x)$ (80) $\displaystyle A_{1}(x)$
$\displaystyle=\int dk\ A(x,k)\rho^{\infty}(k;x)~{}.$ (81)

The Itô Langevin equation corresponding to the backward Kolmogorov equation
(77) for $v_{0}$ constitutes our slow variable dynamics in the regime
$\varepsilon\ll 1$ and reads

$\dot{x}(t)=F(x(t))+A(x(t))\eta_{x}(t)~{}.$ (82)

This is the key result that we draw on in Section 5, as detailed in the rest
of this Appendix.

###  C.2 Homogenisation for OU2 Process

We now discuss the treatment of dynamics on multiple timescales in the
specific context of the OU2 process, Eqs. (2a) and (2b), which we restate for
convenience here:

$\displaystyle\dot{x}$
$\displaystyle=-\left(\bar{k}+k\right)x+\sqrt{2D_{x}}\zeta_{x}$ (83)
$\displaystyle\dot{k}$ $\displaystyle=-\mu k+\sqrt{2D_{k}}\zeta_{k}.$ (84)

####  C.2.1 Adiabatic limit —

First, following Sec. 5.1, we consider a naïve separation of timescales
between the two dynamics, akin to the adiabatic limit in thermodynamics,
obtained by introducing a small dimensionless parameter $\varepsilon$ via the
rescaling $\mu\to\tilde{\mu}/\varepsilon^{2}$ and
$D_{k}\to\tilde{D}_{k}/\varepsilon^{2}$,

$\displaystyle\dot{x}$
$\displaystyle=-\left(\bar{k}+k\right)x+\sqrt{2D_{x}}\zeta_{x}$ (85)
$\displaystyle\dot{k}$
$\displaystyle=-\frac{\tilde{\mu}}{\varepsilon^{2}}k+\sqrt{\frac{2\tilde{D}_{k}}{\varepsilon^{2}}}\zeta_{k}.$
(86)

Comparing these equations with Eqs. (71) and (72) and matching terms by their
order in $\varepsilon$, we conclude that in this limit $f_{0}(x,k)=0$,
$f_{1}(x,k)=-\left(\bar{k}+k\right)x$ and $\alpha(x,k)=\sqrt{2D_{x}}$, while
$g(x,k)=-\tilde{\mu}k$ and $\beta(x,k)=\sqrt{2\tilde{D}_{k}}$. Employing the
procedure outlined in the previous section, it is then straightforward to
verify that $F(x)=-\bar{k}x$ and $A(x)=\sqrt{2D_{x}}$ leading to

$\dot{x}(t)=-\bar{k}x(t)+\sqrt{2D_{x}}\zeta_{x}(t),$ (87)

In other words, no signature of the coupling between $x$ and $k$ survives in
the effective dynamics for the slow variable $x$ in the limit
$\varepsilon\rightarrow 0.$

####  C.2.2 White noise limit —

For a non-trivial contribution to appear in the slow dynamics, we subsequently
consider a second fast-slow regime, which we refer to as the white noise limit
in Sec. 5.2. This time, we perform the rescaling
$\mu\to\tilde{\mu}/\varepsilon^{2}$ and
$D_{k}\to\tilde{D}_{k}/\varepsilon^{4}$, whereby

$\displaystyle\dot{x}$
$\displaystyle=-\left(\bar{k}+k\right)x+\sqrt{2D_{x}}\zeta_{x}~{},$ (88)
$\displaystyle\dot{k}$
$\displaystyle=-\frac{\tilde{\mu}}{\varepsilon^{2}}k+\sqrt{\frac{2\tilde{D}_{k}}{\varepsilon^{4}}}\zeta_{k}~{}.$
(89)

Now, let $\chi\equiv\varepsilon k$, such that

$\displaystyle\dot{x}$ $\displaystyle=-\frac{1}{\varepsilon}\chi
x-\bar{k}x+\sqrt{2D_{x}}\zeta_{x}~{},$ (90) $\displaystyle\dot{\chi}$
$\displaystyle=-\frac{1}{\varepsilon^{2}}\tilde{\mu}\chi+\frac{1}{\varepsilon}\sqrt{2\tilde{D}_{k}}\zeta_{k}~{}.$
(91)

Comparing once again these equations with Eqs. (71) and (72), we identify
$f_{0}(x,\chi)=-\chi x$, $f_{1}(x,\chi)=-\bar{k}x$ and
$\alpha(x,\chi)=\sqrt{2D_{x}}$, while $g(x,\chi)=-\tilde{\mu}\chi$ and
$\beta(x,\chi)=\sqrt{2\tilde{D}_{k}}$. Following the procedure outlined above,
we find that the effective dynamics for the slow variable $x$ take the form

$\dot{x}(t)=-\left(\bar{k}-\frac{D_{k}}{\mu^{2}}\right)x(t)+\sqrt{2\left(D_{x}+\frac{D_{k}x^{2}(t)}{\mu^{2}}\right)}\zeta_{x}(t)~{},$
(92)

where we have used that ${\tilde{D}_{k}}/{\tilde{\mu}^{2}}={D_{k}}/{\mu^{2}}$
in this case. We thus recover the effective slow dynamics (41) studied in the
main text.

##  Appendix D Connection to the associated Legendre equation

Upon enforcing the separation of timescales detailed in Sec. 5.2, the Fokker-
Planck equation (42) for the marginal distribution of the slow variable $x(t)$
can be mapped on to an associated Legendre differential equation through a
change of variable, which we detail here. Indeed, the marginal steady-state
distribution $P(x)$ is solution to

$0=\frac{\partial}{\partial_{x}}\left\\{\frac{\partial}{\partial_{x}}\left[\left(D_{x}+\frac{\tilde{D}_{k}x^{2}}{\tilde{\mu}^{2}}\right)P(x,t)\right]+\left(\bar{k}-\frac{\tilde{D}_{k}}{\tilde{\mu}^{2}}\right)xP(x,t)\right\\}~{}.$
(93)

Introducing $h\equiv\tilde{D}_{k}/\tilde{\mu}^{2}$, we define

$P(x)=a_{0}\left(D_{x}+hx^{2}\right)^{-\frac{1}{4}\left(1+\frac{\bar{k}}{h}\right)}\Phi(x)~{},$
(94)

which we substitute into Eq. (93) to obtain

$0=\left(1+\frac{h}{D_{x}}x^{2}\right)\Phi^{\prime\prime}(x)+\frac{2h}{D_{x}}x\Phi^{\prime}(x)+\frac{h+\bar{k}}{4D_{x}}\left(1+\frac{D_{x}-\bar{k}x^{2}}{D_{x}+hx^{2}}\right)\Phi(x).$
(95)

Finally, we perform the change of variable $s=ix\sqrt{h/D_{x}}$ and derive the
equation for $\Phi(s)$ with imaginary argument

$0=(1-s^{2})\Phi^{\prime\prime}(s)-2s\Phi^{\prime}(x)+\left[\left(\frac{\bar{k}}{2h}-\frac{1}{2}\right)\left(\frac{\bar{k}}{2h}-\frac{1}{2}+1\right)-\left(\frac{h+\bar{k}}{2h}\right)^{2}\frac{1}{1-s^{2}}\right]\Phi(s),$
(96)

which is now in the form of an associated Legendre differential equation [61].
The formal solution of the original equation, up to a normalisation factor, is
thus

$P(x)\propto\left(D_{x}+hx^{2}\right)^{-\frac{\mu}{2}}\left[c_{1}P_{\lambda}^{\mu}\left(i\sqrt{\frac{h}{D_{x}}}x\right)+c_{2}Q_{\lambda}^{\mu}\left(i\sqrt{\frac{h}{D_{x}}}x\right)\right]$
(97)

where

$\mu=\frac{\bar{k}+h}{2}\quad\text{and}\quad\lambda=\frac{\bar{k}-h}{2},$ (98)

and $P_{\lambda}^{\mu}$ and $Q_{\lambda}^{\mu}$ denote the associated Legendre
functions of the first and second kind, respectively, which are precisely the
two linearly independent solutions of the Legendre equation.

##  Appendix E Derivation of the entropy production rate from the Gibbs-
Shannon entropy

In this appendix, we summarise the full derivation of the entropy production
rate $\dot{S}$ found in Ref. [16]. We begin from the Fokker-Planck equation
for the joint probability distribution for the particle position and stiffness
governed by Eq. (2):

$\partial_{t}P(x,k,t)=D_{x}\partial_{x}^{2}P+(\bar{k}+k)\partial_{x}(xP)+D_{k}\partial_{k}^{2}P+\mu\partial_{k}(kP)=-\partial_{x}J(x,k,t)-\mathcal{J}(x,k,t)$
(99)

where we have defined the probability currents for the positional and
stiffness variables respectively as

$J(x,k,t)=-(\bar{k}+k)xP-
D_{x}\partial_{x}P\quad\text{and}\quad\mathcal{J}(x,k,t)=-\mu kP-
D_{k}\partial_{k}P.$ (100)

Following the standard approach for the thermodynamic treatment of diffusive
systems with fluctuating potentials [16], we differentiate the Gibbs-Shannon
entropy with respect to time to write

$\dot{S}(t)=-\int dx\int
dk\partial_{t}P(x,k,t)\log\left(\frac{P(x,k,t)}{\bar{P}}\right),$ (101)

where $\bar{P}$ is an arbitrary constant for dimensional consistency and we
work in units such that $k_{B}=1.$ We identify two equal and opposite
contributions to this rate that we write as

$\dot{S}(t)=\dot{S}_{i}(t)+\dot{S}_{e}(t)$ (102)

where have defined the internal (or total) entropy production rate

$\dot{S}_{i}(t)=\int dk\int
dx\frac{1}{P(x,k,t)}\left[\frac{J^{2}(x,k,t)}{D_{x}}+\frac{\mathcal{J}^{2}(x,k,t)}{D_{k}}\right]$
(103)

and the external entropy production (or entropy flow)

$\dot{S}_{e}(t)=\int dk\int
dx\frac{(\bar{k}+k)xJ(x,k,t)}{D_{x}}+\frac{1}{D_{k}}\int dk\int dx\mu
k\mathcal{J}(x,k,t).$ (104)

At steady-state, $\partial_{t}P=0$ and hence
$\lim_{t\rightarrow\infty}\dot{S}(t)$ vanishes, which implies the steady-state
relation
$\lim_{t\rightarrow\infty}\dot{S}_{i}(t)=-\lim_{t\rightarrow\infty}\dot{S}_{e}(t)$.
In general, the internal entropy production is the quantity of interest in the
thermodynamic characterization of non-equilibrium stochastic processes. In
what follows, we evaluate the integrals for the external entropy production at
steady-state due to their simpler form, then employ the steady-state relation
to evaluate the more classical thermodynamic quantity.

The dynamics for the stiffness are given by the equilibrium Ornstein-Uhlenbeck
process, thus at steady-state the current $\mathcal{J}(k)=\int
dx\mathcal{J}(x,k)$ vanishes. Thus the only contribution to $\dot{S}_{e}(t)$
at steady-state is the first term, that we can re-write as

$\lim_{t\rightarrow\infty}\dot{S}_{e}(t)=-\frac{\langle((\bar{k}+k)x)^{2}\rangle}{D_{x}}+\bar{k},$
(105)

where the average $\langle\cdot\rangle$ is taken over the joint probability
distribution $P(x,k).$ Multiplying the Fokker-Planck equation by $x^{2}$, we
integrate over $x$ to derive an equation for the dynamics of the marginal
variance $\Xi(k,t)$:

$\partial_{t}\Xi(k,t)=2D_{x}P^{\rm
tot}(k,t)-2(\bar{k}+k)\Xi(k,t)+\partial_{k}\big{[}D_{k}\partial_{k}\Xi(k,t)+\mu
k\Xi(k,t)\big{]}.$ (106)

Integrating this last equation at steady-state with respect to $k$ leads to
$\langle kx^{2}\rangle=D_{x}$.

We then multiply (106) by $k$ before again integrating over $k$ to obtain

$\int
dk\bigg{[}\frac{(\bar{k}+k)^{2}}{D_{x}}\Xi(k)\bigg{]}=\bar{k}+\frac{1}{2D_{x}}\int
dk\bigg{[}D_{k}\partial_{k}\Xi(k)+k\partial_{k}\big{[}\mu
k\Xi(k)\big{]}\bigg{]}.$ (107)

Finally, we argue that the second term on the right-hand side of (107)
equation can be written as

$\frac{1}{2D_{x}}\int
dk\bigg{[}D_{k}\partial_{k}\Xi(k)+k\partial_{k}\big{[}\mu\bar{k}\Xi(k)\big{]}\bigg{]}=-\frac{\mu}{2D_{x}}\big{[}\langle(\bar{k}+k)x^{2}\rangle-\bar{k}\langle
x^{2}\rangle\big{]},$ (108)

noticing that the term proportional to $D_{k}$ vanishes by imposing a
sufficiently fast decay of $\partial_{x}P$ at $x\to\pm\infty$. Using
$\langle(\bar{k}+k)x^{2}\rangle=D_{x}$, we conclude that the internal entropy
production rate at steady-state can be expressed as

$\lim_{t\rightarrow\infty}\dot{S}_{i}=\frac{\mu\bar{k}}{2D_{x}}\bigg{(}\langle
x^{2}\rangle-\frac{D_{x}}{\bar{k}}\bigg{)}.$ (109)