0.8pt

# First passage time distribution of active thermal particles in potentials

Benjamin Walter  Department of Mathematics, Imperial College London, 180
Queen’s Gate, SW7 2AZ London, United Kingdom  Centre for Complexity &
Networks, Imperial College London, SW7 2AZ London, United Kingdom  SISSA -
International School for Advanced Studies, via Bonomea 265, 34135 Trieste,
Italy  Gunnar Pruessner  Department of Mathematics, Imperial College London,
180 Queen’s Gate, SW7 2AZ London, United Kingdom  Centre for Complexity &
Networks, Imperial College London, SW7 2AZ London, United Kingdom  Guillaume
Salbreux  The Francis Crick Institute, 1 Midland Road, NW1 1AT London, United
Kingdom  Department of Genetics and Evolution, University of Geneva, Quai
Ernest-Ansermet 30, 1205 Geneva, Switzerland

(August 27, 2024; August 27, 2024)

###### Abstract

We introduce a perturbative method to calculate all moments of the first-
passage time distribution in stochastic one-dimensional processes which are
subject to both white and coloured noise. This class of non-Markovian
processes is at the centre of the study of thermal active matter, that is
self-propelled particles subject to diffusion. The perturbation theory about
the Markov process considers the effect of self-propulsion to be small
compared to that of thermal fluctuations. To illustrate our method, we apply
it to the case of active thermal particles (i) in a harmonic trap (ii) on a
ring. For both we calculate the first-order correction of the moment-
generating function of first-passage times, and thus to all its moments. Our
analytical results are compared to numerics.

##  I  Introduction and Main results

###  I.1 Introduction

Understanding the statistical properties of first-passage times (FPT), the
time a stochastic process takes to first reach a prescribed target, has
enjoyed increased attention over the last two decades [1, 2, 3], since it is a
key characteristic of complex systems, such as chemical reactions [4],
polymer-synthesis [5], intra-cellular events [6], neuronal activity [7] or
financial systems [8]. Besides their dynamical information, FPTs are helpful
to understand spatial properties of complex networks [9], extreme values of
stochastic processes [10] and characteristic observables in out-of-equilibrium
statistical physics [11].

Figure 1: A particle in a potential (orange parabola) subject to both white
and coloured noise (see Eq. (1)). While the white noise models a thermal
environment whose timescale of correlation is negligibly small, the driving
term models hidden degrees of freedom which are correlated over timescales
comparable to those of the particle’s stochastic dynamics. Those driving
forces induce correlations (pink correlation kernel) in the particle’s
increments and therefore break its Markovianity. In this work, we study first-
passage times $\tau_{x_{0},x_{1}}$; the time such a random walker (blue rough
path) takes to first reach $x_{1}$ starting from $x_{0}$ (dashed lines).

For Markovian stochastic processes, the FPT has been studied for over a
century, and often the full distribution of FPT, or their moment generating
function, can be found in closed form [12, 13, 14, 15]. For non-Markovian
processes, however, the problem of finding the FPT distribution is much more
difficult and often the focus has been on the mean first-passage time (MFPT)
alone [16, 17, 18, 19, 20, 21, 22, 23, 24]. The MFPT, however, can be
insufficient to characterise complex dynamics, emphasizing the need for a more
precise understanding of the full FPT distribution [25, 26].

In this article, we address this challenge and compute the full moment-
generating function of a class of non-Markovian stochastic processes
perturbatively. In doing so, we obtain all moments of the distribution to the
same order in the perturbative expansion. By numerical integration, we further
recover the full FPT distribution With this method the full distribution is
obtained systematically in the presence of correlated driving noise and white
thermal noise for a wide range of settings, including a potential. The
formulas we obtain order by order are exact, and the results we obtain for two
systems, as an illustration, are in excellent agreement with numerical
simulations. This allows for a detailed analysis of the qualitative changes of
the FPT distribution induced by correlations in stochastic forces.

In this article, we consider the FPT problem in one dimension. We study the
FPT of a particle at position $x_{t}$ placed in a potential $V(x)$, and
subject to white thermalwhitethermal noise $\xi_{t}$,
$\left\langle{\xi_{t_{0}}\xi_{t_{1}}}\right\rangle=2D_{x}\delta(t_{1}-t_{0})$,
modelling a surrounding heat bath at a temperature proportional to $D_{x}$
(via the fluctuation-dissipation theorem [27, 28]). Additionally, we assume
that the particle is subject to a second stochastic force which may model
either self-propulsion or hidden internal degrees of freedom. We refer to this
force as activewhiteactive, as we equate “activeness” with the presence of
coloured noise (e.g. telegraphic noise or an Ornstein-Uhlenbeck noise),
following e.g. [29, 30, 31]. The active noise $y_{t}$ is a stochastically
independent stationary stochastic force with zero mean and a non-vanishing
autocorrelation function $C_{2}(t_{1}-t_{0})=\overline{y_{t_{0}}y_{t_{1}}}$.
Since the process includes two independent stochastic forces, we introduce
$\left\langle{\bullet}\right\rangle$ to denote averages over $\xi_{t}$, and
$\overline{\bullet}$ to denote averages over $y_{t}$. The particle’s position
satisfies a Langevin equation [1] of the type

$\displaystyle\dot{x}_{t}=-V^{\prime}(x_{t})+\xi_{t}+\varepsilon y_{t};$ (1)

where $\varepsilon$ denotes a dimensionless coupling to the active term which
is chosen small enough such that $\varepsilon y_{t}\ll\xi_{t}$ in probability.
We refer to the process described in Eq. (1) as active thermalwhiteactive
thermal process, as opposed to $(i)$ purely active (“athermal”, [32])
processes and $(ii)$ purely thermal processes since $y_{t}$ breaks thermal
equilibrium ([33]). Active thermal processes have recently been considered in
e.g. [34, 35, 36, 37]. While the FPT problem in purely active matter such as
“run-and-tumble” processes ([38]) has been studied in [39]. Understanding the
FPT in active thermal matter is relevant to e.g. neural activity [40],
transport in living cells [41], and molecular motors [42, 43] where the
influence of the thermal environment cannot necessarily be neglected.

###  I.2 Main results

(a)  First-passage time distribution of active thermal Ornstein-Uhlenbeck
Process (ATOU) (cf. Eq. (76)) compared to numerical Laplace inversion of
analytically obtained moment generating function (cf. theoretical result in
Eq. (98)) for various values of $\nu$ (solid lines). The plot marks indicate
the distribution as sampled through Monte Carlo simulations with $5\times
10^{6}$ runs. Simulation parameters are
$x_{0}=0,x_{1}=1,\alpha=1,D_{x}=1,D_{y}=1,\beta=\frac{1}{2}$ while
$\varepsilon$ is tuned to fix $\nu=D_{y}\beta\varepsilon^{2}/(D_{x}\alpha)$ to
values as indicated in legend. The inset shows rescaled deviations to the
undriven first-passage time distribution (plot marks, cf. Eq. (104)) as
compared to the first-order correction (solid lines, $\nu=0.1$ omitted).
Higher-order corrections appear for growing values of $\nu$ (cf. Fig. 4). See
Sec. III.1.2 for discussion.

(b)  First-passage time distribution of active thermal Brownian Motion (ATBM)
(cf. Eq. (108)) compared to numerical Laplace inversion of analytically
obtained moment generating function (cf. theoretical result in Eq. (127)) for
various values of $\nu$ (solid lines). The plot marks indicate the
distribution as sampled through Monte Carlo simulations with $10^{6}$ runs.
Simulation parameters are
$x_{0}=0,x_{1}=\pi,r=1,D_{x}=1,D_{y}=1,\beta=\frac{1}{2}$ while $\varepsilon$
is tuned to fix $\nu=D_{y}r^{2}\beta\varepsilon^{2}/D_{x}^{2}$ to values as
indicated in legend. The inset shows rescaled deviations to the undriven
first-passage time distribution (cf. Eq. (104)) as compared to the first-order
correction (solid lines, $\nu=0.1$ omitted). Higher-order corrections appear
for growing values of $\nu$ (cf. Fig. 7). See Sec. III.2.2 for discussion.

Figure 2: First order corrections to the probability distribution of first-
passage times as found by the framework presented in this work for two example
processes.

In the following we denote $\tau_{x_{0},x_{1}}$ the first-passage time of
$x_{t}$ defined as

$\displaystyle\tau_{x_{0},x_{1}}\coloneqq\inf_{t_{1}>t_{0}}\left\\{t:x_{t_{1}}=x_{1}|x_{t_{0}}=x_{0}\right\\},$
(2)

and $p_{x_{0},x_{1}}(t_{0},t_{1})=p_{x_{0},x_{1}}(t_{1}-t_{0})$ the
probability distribution of $\tau_{x_{0},x_{1}}$, where time-invariance of the
$y_{t}$-averaged process ensures that $p_{x_{0},x_{1}}$ depends only on the
time difference $t_{1}-t_{0}$. The central result of the present work concerns
the moment-generating function of first-passage times of $x_{t}$,

$M_{x_{0},x_{1}}(s)=\left\langle{\overline{e^{-s\tau_{x_{0},x_{1}}}}}\right\rangle\\\
=1-s\left\langle{\overline{\tau_{x_{0},x_{1}}}}\right\rangle+\frac{s^{2}}{2}\left\langle{\overline{\left(\tau_{x_{0},x_{1}}\right)^{2}}}\right\rangle+...$
(3)

for $s\geq 0$. We note that $M_{x_{0},x_{1}}(s)$ is the Laplace transform of
the probability distribution $p_{x_{0},x_{1}}(t)$, and the restriction to the
negative imaginary axis of the complex function
$\hat{p}_{x_{0},x_{1}}(z)=\int_{0}^{\infty}\mathrm{d}{t}\,p_{x_{0},x_{1}}(t)e^{-izt}$.
Because $p_{x_{0},x_{1}}(t)$ is non-zero only for $t\geq 0$, and
$\int_{0}^{\infty}\mathrm{d}{t}\,p_{x_{0},x_{1}}(t)=1$, the function
$\hat{p}_{x_{0},x_{1}}(z)$ is bounded, and by the Paley-Wiener theorem
holomorphic, on the lower half complex plane [44]. Besides,
$\hat{p}_{x_{0},x_{1}}(z)$ converges to the Fourier transform of
$p_{x_{0},x_{1}}(t)$ on the real axis (defined in Eq. 7). Therefore in the
following we will use alternatively the Fourier transform and Laplace
transform for our calculations whenever convenient, and assume that conditions
are met such that a simple rotation $\omega=-is$ allows to go from one
transform of $p_{x_{0},x_{1}}(t)$ to the other.

We now assume that the moment generating function $M_{x_{0},x_{1}}(s)$ has an
expansion of the form

$\displaystyle M_{x_{0},x_{1}}(s)=M^{0}_{x_{0},x_{1}}(s)+\nu
M^{1}_{x_{0},x_{1}}(s)+\mathcal{O}(\nu^{2}).$ (4)

where $M^{0}$ and $M^{1}$ are the coefficients of expansion of
$M_{x_{0},x_{1}}(s)$ in a dimensionless parameter $\nu$ of order
$\mathcal{O}(\varepsilon^{2})$, around $\nu=0$. The case $\varepsilon=0$
corresponds to a purely thermal process, and the corresponding moment-
generating function $M^{0}$, can be found by classical methods such as the
Darling-Siegert method which is discussed in the next subsection. Here we
assume that around this state $M_{x_{0},x_{1}}(s)$ is analytic in $\nu$.

The first-order contribution, $M^{1}$, requires some deeper analysis. Much of
what follows is dedicated to the calculation of $M^{1}$, which to our
knowledge is new in the literature. In principle, the method we present here
is capable of calculating coefficients $M^{2},M^{3},...$ of arbitrarily high
order of $\nu$ for arbitrarily coloured noise (cf. Eq. (68)) as long as the
autocorrelations can be integrated suitably. Furthermore the potentials $V(x)$
are arbitrary, as long as an associated differential operator can be
diagonalised (Sec. III and [45]), otherwise it needs to be treated
perturbatively as well.

We further illustrate our framework by explicitly computing $M^{0}$ and
$M^{1}$ for two cases, each of which are additionally driven by coloured
Gaussian noise, i.e. $y_{t}$ is Gaussian and has correlator
$\overline{y_{t_{0}}y_{t_{1}}}=D_{y}\beta e^{-\beta|t_{1}-t_{0}|}$ with some
diffusivity $D_{y}$ and correlation time $\beta^{-1}$. In the first case, the
particle is placed in a harmonic potential, $V(x)=\frac{\alpha}{2}x^{2}$. This
particular model has been studied in, e.g. , [37]. We refer to this model as
active thermal Ornstein Uhlenbeck process (ATOU). While $M^{0}$ (see Eq. (99))
has been long known, $M^{1}$ (see Eq. (100)) is a new result. By numerically
computing its inverse Laplace transform, we obtain the full probability
distribution of first-passage times to first order in $\nu$. In Fig. 2a we
compare this to the numerically obtained probability density function.
Secondly, we calculate $M^{0}$ and $M^{1}$ for the case of a Brownian Motion
on a ring of radius $r$ driven by coloured Gaussian noise which we refer to as
active thermal Brownian Motion (ATBM). Again, through numerically computing
its inverse Laplace transform we obtain the corrected probability distribution
of first-passage times which is compared to numerics in Fig. 2b.

Our method is systematicwhitesystematic since it allows its user to calculate
in principle corrections to arbitrary order, and it is
controlledwhitecontrolled in the sense that the error can be made arbitrarily
small. Further all momentswhiteall moments are available at once and to equal
perturbative order. It is also valid for arbitrary noise colours $\beta^{-1}$.

The paper is structured as follows. In Sec. II we give detailed account of how
to calculate $M_{x_{0},x_{1}}(s)$ for small $\nu$. First, we reproduce the
Darling-Siegert argument in the equilibrium case ($\nu=\varepsilon=0$). Next,
we introduce a perturbative version of the Darling-Siegert equation. Then, we
obtain, as an intermediate result, a formula for $M_{x_{0},x_{1}}(s)$ which is
still a functional of the colored noise $y_{t}$. In the last step, we need to
average over the stationary distribution of $y_{t}$ to arrive at the explicit
formula Eq. (II.4) which is the main result of our work. In the subsequent
section III, we calculate all quantities required for the case of a harmonic
potential and a Brownian motion with periodic boundary conditions and arrive
at the first-order correction to the moment-generating function of first-
passage times Eq. (100). Section IV concludes with a discussion of our
findings.

##  II  Perturbation Theory

As outlined above, in this work we present a way to calculate the moment-
generating function of first-passage times of stochastic processes which are
close to an equilibrium state. The underlying assumption is that the moment-
generating function varies smoothly as $\varepsilon$, the coupling to the
self-propelling force, is switched on. The moment-generating function of the
equilibrium version of the process ($\varepsilon=0$) is assumed to be known in
closed form, as is for instance justified for the Ornstein-Uhlenbeck process
or Brownian Motion ([45]). This exact form is then corrected by terms in the
spirit of a perturbative expansion which is controlled by powers of a
dimensionless parameter describing the distance to equilibrium. First, we
revise the arguments given by Darling and Siegert for the equilibrium case
([14, 15]). Next, we outline our perturbative approach to the active case.

Before going further, we introduce some notations (see appendix A for a list).
The transition probability density of progressing from $x_{0}$ at $t_{0}$ to
$x_{1}$ at $t_{1}$ is denoted by

$\displaystyle T_{x_{0},x_{1}}(t_{0},t_{1})=T(t_{0},t_{1})$ (5)

where the arguments $x_{0}$ and $x_{1}$ will later be dropped wherever
confusion can be avoided, for the sake of easier notation. Analogously the
return probability at $x_{1}$, $T_{x_{1},x_{1}}(t_{0},t_{1})$ is denoted by

$\displaystyle
R_{x_{1}}(t_{0},t_{1})=R(t_{0},t_{1})=T_{x_{1},x_{1}}(t_{0},t_{1})$ (6)

In the following, we denote the Fourier transform of a function $f(t)$ by

$\displaystyle\hat{f}(\omega)=\int_{-\infty}^{\infty}\mathrm{d}{t}\,e^{-i\omega
t}f(t)\leavevmode\nobreak\ ,$ (7)

with inverse

$\displaystyle f(t)=\int_{-\infty}^{\infty}\mathrm{d}\bar{}\hskip
1.00006pt\omega\,e^{i\omega t}f(\omega),$ (8)

where $\mathrm{d}\bar{}\hskip
1.00006pt{\omega}=\frac{\mathrm{d}{\omega}\,}{2\pi}$.

In the following, we consider functions $f(t_{0},t_{1})$ which depend on the
difference $t_{1}-t_{0}$ only, i.e. $f(t_{0},t_{1})=f(t_{1}-t_{0})$. We refer
to functions of this type as diagonalwhitediagonal. The Fourier transform in
two variables of these functions read

$\displaystyle\hat{f}(\omega_{0},\omega_{1})$
$\displaystyle=\int_{-\infty}^{\infty}\mathrm{d}{t_{0}}\,\int_{-\infty}^{\infty}\mathrm{d}{t_{1}}\,e^{-i\omega_{0}t_{0}-i\omega_{1}t_{1}}f(t_{0},t_{1})$
$\displaystyle=\hat{f}(\omega_{1})\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{1})\leavevmode\nobreak\ ,$ (9)

where $\delta\bar{}\hskip 1.00006pt\,(\omega)=2\pi\delta(\omega)$.

####  II.0.1 Equilibrium case: The Darling-Siegert solution

We here consider the equilibrium case of Eq. (1) defined by setting
$\varepsilon=0$. As $x_{t}$ is Markovian, the functions $p$ and $T$ satisfy
the following renewal equation:

$\displaystyle
T_{x_{0},x_{1}}(t_{0},t_{1})=\int_{t_{0}}^{t_{1}}\mathrm{d}{t^{\prime}}\,p_{x_{0},x_{2}}(t_{0},t^{\prime})T_{x_{2},x_{1}}(t^{\prime},t_{1})$
(10)

for all $x_{2}\in(x_{0},x_{1}]$. Applying a Fourier transform to Eq. (10), the
time-homogeneity of both $T_{x_{0},x_{1}}(t_{0},t_{1})$ and $p(t_{0},t_{1})$
translates into diagonality in frequencies and turns the convolution into a
product, such that the result can be stated at the level of the amplitudes
alone. Rearranging the terms and choosing $x_{2}=x_{1}$ results in

$\displaystyle\hat{p}(\omega)=\frac{\hat{T}(\omega)}{\hat{R}(\omega)}\leavevmode\nobreak\
.$ (11)

Since the Fourier transform of a probability density equals its characteristic
function, Eq. (11) recovers all moments of the first passage time provided
$\hat{T}$ is known. Further, setting $\omega=-is$ for some
$s\in\mathbb{R}_{+}$ turns the Fourier transforms into Laplace transforms and
the characteristic function into the moment generating function. This recovers
the Darling-Siegert equation in its original form in which the Laplace
transform of transition and return probabilities is linked to the moment-
generating function of first-passage times.

####  II.0.2 Out-of-equilibrium: A perturbative approach

The argument made by Darling and Siegert breaks down when $\varepsilon\neq 0$:
indeed when averaging over the driving noise $y(t)$, the renewal equation Eq.
(10) no longer is true. The approach we take in this paper, consists of three
steps:

  1. 1.

Fix a particular realisation $y(t)$, and expand transition and return
probabilities of $x_{t}$ as functional expansion around $y=0$ of the form

$\hat{T}(\omega_{0},\omega_{1},[\hat{y}])=\sum_{n=0}^{\infty}\frac{\varepsilon^{n}}{n!}\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{1}}...\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{n}}\delta\bar{}\hskip
1.00006pt\,\left(\omega_{0}+\omega_{1}+\sum_{i=1}^{n}\tilde{\omega}_{i}\right)\\\
\times\hat{T}^{(n)}\left(\omega_{1},\tilde{\omega}_{1},...,\tilde{\omega}_{n}\right)\hat{y}(-\tilde{\omega}_{1})...\hat{y}(-\tilde{\omega}_{n})\leavevmode\nobreak\
.$ (12)

  2. 2.

As long as $y(t)$ is fixed, the process, when understood as conditioned on
this particular driving, satisfies a renewal equation of the type Eq. (10).
Inserting the perturbative transition and return probabilities from the
previous step, gives a perturbative series for the first-passage time density
$\hat{p}\left(\omega_{0},\omega_{1};[y]\right)$ of $x_{t}$ conditioned on a
particular $y_{t}$.

  3. 3.

Averaging over the ensemble of driving noises. For simplicity, we here assume
that the correlation function of the driving noise is given by

$\displaystyle\overline{y_{t_{0}}y_{t_{1}}}=D_{y}\beta
e^{-\beta|t_{1}-t_{0}|}\leavevmode\nobreak\ ,$ (13)

where $\beta$ is the inverse correlation time. Generally, when computing the
term of order $\varepsilon^{n}$, the first $n$ moments of $y_{t}$ need to be
known.

This procedure leads to the central result of this work: the moment-generating
function of first-passage times to second order in $\varepsilon$ reads

$\displaystyle\hat{p}(\omega)=$
$\displaystyle\frac{\hat{T}^{(0)}\left(\omega\right)}{\hat{R}^{(0)}\left(\omega\right)}+\frac{\varepsilon^{2}D_{y}\beta}{2}\left[\frac{\hat{T}^{(2)}\left(\omega;i\beta,-i\beta\right)}{\hat{R}^{(0)}\left(\omega\right)}-2\frac{\hat{T}^{(1)}\left(\omega-i\beta,i\beta\right)\hat{R}^{(1)}\left(\omega,i\beta\right)}{\hat{R}^{(0)}\left(\omega-i\beta\right)\hat{R}^{(0)}\left(\omega\right)}\right.$
(14)
$\displaystyle\qquad\qquad\qquad\left.+2\frac{\hat{T}^{(0)}\left(\omega\right)\hat{R}^{(1)}\left(\omega-i\beta,i\beta\right)\hat{R}^{(1)}\left(\omega,i\beta\right)}{\left(\hat{R}^{(0)}\left(\omega\right)\right)^{2}\hat{R}^{(0)}\left(\omega-i\beta\right)}-\frac{\hat{T}^{(0)}\left(\omega\right)\hat{R}^{(2)}\left(\omega;i\beta,-i\beta\right)}{\left(\hat{R}^{(0)}\left(\omega\right)\right)^{2}}\right]+\mathcal{O}(\varepsilon^{4}).$

In the next sections, we derive this relation in more details.

###  II.1 Perturbative Darling-Siegert equation

####  II.1.1 Expression for the first-passage time distribution

By imposing an additional driving noise $y_{t}$, the transition probability
and FPT probability density of $x_{t}$ depend on a particular realisation of
$y_{t}$. Accordingly, we introduce the transition probability density
$T(t_{0},t_{1},[y])$ and FPT probability density $p(t_{0},t_{1},[y])$ as the
densities of the process $x_{t}$ conditioned on $y$ given. The conditional
densities are explicitly dependent on $t_{0}$ and $t_{1}$ rather then their
difference $t_{1}-t_{0}$ because $y_{t}$ is an explicit function of time.

For $y$ fixed, the process remains Markovian and therefore Eq. (10) still
applies and gives rise to:

$T_{x_{0},x_{1}}(t_{0},t_{1};[y])=\int_{-\infty}^{\infty}\mathrm{d}{t^{\prime}}\,p_{x_{0},x_{1}}(t_{0},t^{\prime};[y])\\\
\times R_{x_{1}}(t^{\prime},t_{1};[y])\leavevmode\nobreak\ ,$ (15)

where the dependency of the functions on the spatial values $x_{0}$, $x_{1}$
has been made explicit for clarity, and we have used the fact that
$p(t_{0},t^{\prime};[y])$ vanishes for $t^{\prime}<t_{0}$ and
$R(t^{\prime},t_{1})$ vanishes for $t^{\prime}>t_{1}$ to integrate over the
full real axis.

It is no longer possible to directly invert this equation in Fourier space to
solve for $p$, as done in Eq. (11), since neither terms in the integral are
diagonal, i.e. they depend explicitly on both $t_{0},t^{\prime}$ and
$t^{\prime},t_{1}$. However, introducing the inverse functional $R^{-1}$ which
is defined by the implicit equation:

$\displaystyle\int_{-\infty}^{\infty}\mathrm{d}{t}\,R^{-1}(t_{0},t;[y])R(t,t_{1};[y])=\delta(t_{1}-t_{0})\leavevmode\nobreak\
,$ (16)

the renewal equation can be formally solved by the relation:

$p_{x_{0},x_{1}}(t_{0},t_{1};[y])\\\
=\int_{-\infty}^{\infty}\mathrm{d}{t}\,T_{x_{0},x_{1}}(t_{0},t;[y])R^{-1}_{x_{1}}(t,t_{1};[y])\leavevmode\nobreak\
.$ (17)

Our approach is then to perform a functional expansion of the quantities
involved in Eq. (17) in the function $y$, around the Markovian case $y=0$.

####  II.1.2 Functional expansion of the transition and return probability
densities

We start by performing a Taylor expansion of the transition probability
$T(t_{0},t_{1})$:

$T(t_{0},t_{1};[y])=\sum_{n=0}^{\infty}\frac{1}{n!}\int\mathrm{d}{\tilde{t}_{1}}\,...\mathrm{d}{\tilde{t}_{n}}\,\\\
\frac{\delta^{n}T}{\delta y(\tilde{t}_{1})...\delta
y(\tilde{t}_{n})}(t_{0},t_{1};[0])y(\tilde{t}_{1})...y(\tilde{t}_{n})\leavevmode\nobreak\
.$ (18)

We now argue that this functional expansion can be rewritten in a simpler
form, using the time-invariance properties of the dynamic equation (1). We
introduce the time-shifting operator which transforms a function $y_{t}$ into
a time-shifted function $\tilde{y}_{t}$, according to:

$\displaystyle S_{t_{0}}[y](t)=y_{t+t_{0}}\equiv\tilde{y}_{t}$ (19)

and we will use the fact that:

$\displaystyle T(t_{0},t_{1};[y])=T(t_{0}+\tau,t_{1}+\tau;[S_{-\tau}[y]])$
(20)

for any time-shift $\tau$. This relation results from the invariance by time
shifting of Eq. (1), when the function $y_{t}$ is shifted in time accordingly.
As a result we introduce the shifted function $T_{0}$, defined by;

$\displaystyle T_{0}(\tau;[\tilde{y}])=T(0,\tau;[\tilde{y}])$ (21)

such that Eq. (20) gives, choosing $\tau=-t_{0}$:

$\displaystyle
T(t_{0},t_{1};[y])=T_{0}(t_{1}-t_{0};[S_{t_{0}}[y]])\leavevmode\nobreak\ .$
(22)

Using Eq. (22), the functional derivatives of $T$ can be related to the
functional derivatives of $T_{0}$ through

$\frac{\delta^{n}T}{\delta y(\tilde{t}_{1})...\delta
y(\tilde{t}_{n})}(t_{0},t_{1};[y])=\\\
\frac{\delta^{n}T_{0}}{\delta\tilde{y}(\tilde{t}_{1}-t_{0})...\delta\tilde{y}(\tilde{t}_{n}-t_{0})}(t_{1}-t_{0};[S_{t_{0}}[y]])\leavevmode\nobreak\
,$ (23)

and the functional expansion of $T$ can now be written

$T(t_{0},t_{1};[y])=\sum_{n=0}^{\infty}\frac{1}{n!}\int\mathrm{d}{\tilde{t}_{1}}\,...\mathrm{d}{\tilde{t}_{n}}\,\\\
\times\frac{\delta^{n}T_{0}}{\delta\tilde{y}(\tilde{t}_{1}-t_{0})...\delta\tilde{y}(\tilde{t}_{n}-t_{0})}(t_{1}-t_{0};[0])y(\tilde{t}_{1})...y(\tilde{t}_{n})\leavevmode\nobreak\
,$ (24)

where we have used that $S_{t_{0}}[y=0]=[0]$. So far, the transitional
probability $T(t_{0},t_{1};[y])$ has been expressed as a functional of $y$ in
the time-domain. In order to consider the transition probability as a
functional of the Fourier-transformed driving noise $\hat{y}$, we introduce
$\hat{T}(\omega_{0},\omega_{1};[\hat{y}])$, defined by taking the Fourier
transform of $T(t_{0},t_{1},[y])$ with respect to $t_{0}$, $t_{1}$, and taking
$y$ equal to the inverse Fourier transform of $\hat{y}$ in the definition of
$T$:

$\displaystyle\hat{T}(\omega_{0},\omega_{1};[\hat{y}])$
$\displaystyle=\int\mathrm{d}{t_{0}}\,\mathrm{d}{t_{1}}\,T(t_{0},t_{1},[y])e^{-i(\omega_{0}t_{0}+\omega_{1}t_{1})}\leavevmode\nobreak\
,$ (25)

and similarly:

$\displaystyle\hat{T}_{0}(\omega;[\hat{y}])$
$\displaystyle=\int\mathrm{d}{\tau}\,T_{0}(\tau,[y])e^{-i\omega\tau}\leavevmode\nobreak\
.$ (26)

The functional derivatives of $T_{0}(\tau;[y])$ and
$\hat{T}_{0}(\omega;[\hat{y}])$ are related by

$\displaystyle\frac{\delta^{n}T_{0}(\tau;[y])}{\delta
y(\tilde{t}_{1})...\delta y(\tilde{t}_{n})}$
$\displaystyle=\int\mathrm{d}\bar{}\hskip
1.00006pt{\omega}\mathrm{d}{\tilde{\omega}_{1}}\,...\mathrm{d}{\tilde{\omega}_{n}}\,\frac{\delta^{n}\hat{T}_{0}(\omega;[\hat{y}])}{\delta\hat{y}(\tilde{\omega}_{1})...\delta\hat{y}(\tilde{\omega}_{n})}$
$\displaystyle\qquad\times
e^{-i\tilde{\omega}_{1}\tilde{t}_{1}}...e^{-i\tilde{\omega}_{n}\tilde{t}_{n}}e^{i\omega\tau}\leavevmode\nobreak\
.$ (27)

As a result, combining Eq. (24) and Eq. (25), we obtain the functional
expansion:

$\hat{T}(\omega_{0},\omega_{1},[\hat{y}])\\\
=\sum_{n=0}^{\infty}\frac{1}{n!}\int\mathrm{d}{\tilde{\omega}_{1}}\,...\mathrm{d}{\tilde{\omega}_{n}}\,\delta\bar{}\hskip
1.00006pt\,\left(\omega_{0}+\omega_{1}-\sum_{i=1}^{n}\tilde{\omega}_{i}\right)\\\
\times\frac{\delta^{n}\hat{T}_{0}}{\delta\hat{y}(\tilde{\omega}_{1})...\delta\hat{y}(\tilde{\omega}_{n})}(\omega_{1};[0])\hat{y}(\tilde{\omega}_{1})...\hat{y}(\tilde{\omega}_{n})\leavevmode\nobreak\
.$ (28)

The perturbative expansion of the return probability $R$ is simply obtained
from by letting $x_{0}\to x_{1}$ in the expansion of $T$ and has therefore the
same structure, using $R_{0}(\tau;[\tilde{y}])=R(0,\tau;[\tilde{y}])$. To ease
the notation, we now introduce a shorter notation for the set of functions:

$\displaystyle
T^{(n)}(\omega_{1},\tilde{\omega}_{1},...,\tilde{\omega}_{n})=\frac{(2\pi)^{n}}{\varepsilon^{n}}\frac{\delta^{n}\hat{T}_{0}}{\delta\hat{y}(-\tilde{\omega}_{1})...\delta\hat{y}(-\tilde{\omega}_{n})}(\omega_{1};[0])$
(29) $\displaystyle
R^{(n)}(\omega_{1},\tilde{\omega}_{1},...,\tilde{\omega}_{n})=\frac{(2\pi)^{n}}{\varepsilon^{n}}\frac{\delta^{n}\hat{R}_{0}}{\delta\hat{y}(-\tilde{\omega}_{1})...\delta\hat{y}(-\tilde{\omega}_{n})}(\omega_{1};[0])$
(30)

where the factor $1/\varepsilon^{n}$ has been introduced for convenience, to
make the expansion in small $\varepsilon$ explicit in Eq. (28). We then arrive
at the following form:

$\hat{T}(\omega_{0},\omega_{1},[\hat{y}])=\sum_{n=0}^{\infty}\frac{\varepsilon^{n}}{n!}\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{1}}...\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{n}}\delta\bar{}\hskip
1.00006pt\,\left(\omega_{0}+\omega_{1}+\sum_{i=1}^{n}\tilde{\omega}_{i}\right)\\\
\times\hat{T}^{(n)}\left(\omega_{1},\tilde{\omega}_{1},...,\tilde{\omega}_{n}\right)\hat{y}(-\tilde{\omega}_{1})...\hat{y}(-\tilde{\omega}_{n})\\\
\hat{R}(\omega_{0},\omega_{1},[\hat{y}])=\sum_{n=0}^{\infty}\frac{\varepsilon^{n}}{n!}\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{1}}...\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{n}}\delta\bar{}\hskip
1.00006pt\,\left(\omega_{0}+\omega_{1}+\sum_{i=1}^{n}\tilde{\omega}_{i}\right)\\\
\times\hat{R}^{(n)}\left(\omega_{1},\tilde{\omega}_{1},...,\tilde{\omega}_{n}\right)\hat{y}(-\tilde{\omega}_{1})...\hat{y}(-\tilde{\omega}_{n})\leavevmode\nobreak\
,$ (31)

where we have made a choice of signs of
$\tilde{\omega}_{1},..\tilde{\omega}_{n}$ for convenience.

####  II.1.3 Functional expansion of the inverse of the return probability
density

We now turn to the expansion of $R^{-1}$ in Eq. (17). As $T$ and $R$, $R^{-1}$
satisfies the time-shift invariance relation (20). This can be seen by time
shifting Eq. (16) by the time $\tau$:

$\int_{-\infty}^{\infty}\mathrm{d}{t}\,R^{-1}(t_{0}+\tau,t+\tau;S_{-\tau}[y])\\\
\times R(t+\tau,t_{1}+\tau;S_{-\tau}[y])=\delta(t_{1}-t_{0})$ (32)

and, since $R(t+\tau,t_{1}+\tau;S_{-\tau}[y])=R(t,t_{1},[y])$, one also
obtains $R^{-1}(t_{0}+\tau,t+\tau;S_{-\tau}[y])=R^{-1}(t_{0},t;[y])$ since Eq.
(16) defines $R^{-1}$. The analysis of the previous subsection still applies
and one obtains the expansion

$\displaystyle\hat{R}^{-1}(\omega_{0},\omega_{1},[\hat{y}])=\sum_{n=0}^{\infty}\frac{\varepsilon^{n}}{n!}\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{1}}...\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{n}}\delta\bar{}\hskip
1.00006pt\,\left(\omega_{0}+\omega_{1}+\sum_{i=1}^{n}\tilde{\omega}_{i}\right)$
(33)
$\displaystyle\times(\hat{R}^{-1})^{(n)}(\omega_{1},\tilde{\omega}_{1},...,\tilde{\omega}_{n})\hat{y}(-\tilde{\omega}_{1})...\hat{y}(-\tilde{\omega}_{n})\leavevmode\nobreak\
.$

We then obtain by applying a Fourier transform to Eq. (16):

$\displaystyle\int\mathrm{d}\bar{}\hskip
1.00006pt\omega\hat{R}^{-1}(\omega_{0},\omega;[\hat{y}])\hat{R}(-\omega,\omega_{1};[\hat{y}])=\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{1})\leavevmode\nobreak\ .$ (34)

This relation allows to relate the functions $(R^{-1})^{(n)}$ to the functions
$R^{(n)}$, by plugging Eqs. (31) and (33) and identifying order by order in
$\varepsilon$. Limiting ourselves up to order $2$, we obtain:

$\displaystyle(\hat{R}^{-1})^{(0)}(\omega)$
$\displaystyle=\frac{1}{\hat{R}^{(0)}(\omega)}$ (35)
$\displaystyle(\hat{R}^{-1})^{(1)}(\omega,\tilde{\omega}_{1})$
$\displaystyle=-\frac{\hat{R}^{(1)}\left(\omega,\tilde{\omega}_{1}\right)}{\hat{R}^{(0)}\left(\omega\right)\hat{R}^{(0)}\left(\omega+\tilde{\omega}_{1}\right)}$
(36)
$\displaystyle(\hat{R}^{-1})^{(2)}\left(\omega,\tilde{\omega}_{1},\tilde{\omega}_{2}\right)$
$\displaystyle=\frac{1}{\hat{R}^{(0)}\left(\omega+\tilde{\omega}_{1}+\tilde{\omega}_{2}\right)\hat{R}^{(0)}\left(\omega\right)}\Bigg{\\{}\frac{2\hat{R}^{(1)}\left(\omega+\tilde{\omega}_{2},\tilde{\omega}_{1}\right)\hat{R}^{(1)}\left(\omega,\tilde{\omega}_{2}\right)}{\hat{R}^{(0)}\left(\omega+\tilde{\omega}_{2}\right)}-\hat{R}^{(2)}\left(\omega,\tilde{\omega}_{1},\tilde{\omega}_{2}\right)\Bigg{\\}}.$
(37)

####  II.1.4 Second-order expansion of the first passage time distribution

Equipped with these expansions, one now can expand the first-passage time
density expression (17) to obtain a functional expansion of the first-passage
density in $\varepsilon$, involving the functions $\hat{T}^{(n)}$ and
$\hat{R}^{(n)}$ which are simpler to calculate. We first need to write down
the Fourier-transformed version of Eq. (17):

$\displaystyle\hat{p}(\omega_{0},\omega_{1};[\hat{y}])=\int_{-\infty}^{\infty}\mathrm{d}\bar{}\hskip
1.00006pt\omega\,\hat{T}(\omega_{0},\omega;[\hat{y}])\,\hat{R}^{-1}(-\omega,\omega_{1};[\hat{y}])\leavevmode\nobreak\
.$ (38)

where the dependency on $x_{0}$, $x_{1}$ is here implicit. Performing an
expansion of this relation in $\varepsilon$, the result reads to second order:

$\displaystyle\hat{p}(\omega_{0},\omega_{1};[\hat{y}])=\frac{\hat{T}^{(0)}(\omega_{1})}{\hat{R}^{(0)}\left(\omega_{1}\right)}\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{1})+\varepsilon\left[-\frac{\hat{T}^{(0)}\left(-\omega_{0}\right)\hat{R}^{(1)}\left(\omega_{1},-\omega_{0}-\omega_{1}\right)}{\hat{R}^{(0)}\left(-\omega_{0}\right)\hat{R}^{(0)}\left(\omega_{1}\right)}+\hat{T}^{(1)}\left(\omega_{1},-\omega_{0}-\omega_{1}\right)\frac{1}{\hat{R}^{(0)}\left(\omega_{1}\right)}\right]\hat{y}(\omega_{0}+\omega_{1})$
$\displaystyle+\frac{1}{2}\varepsilon^{2}\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\bigg{[}\frac{\hat{T}^{(0)}\left(-\omega_{0}\right)}{\hat{R}^{(0)}\left(-\omega_{0}\right)\hat{R}^{(0)}\left(\omega_{1}\right)}\bigg{\\{}\frac{2\hat{R}^{(1)}\left(-\tilde{\omega}-\omega_{0};\tilde{\omega}\right)\hat{R}^{(1)}\left(\omega_{1},-\omega_{0}-\omega_{1}-\tilde{\omega}\right)}{\hat{R}^{(0)}\left(-\tilde{\omega}-\omega_{0}\right)}-\hat{R}^{(2)}\left(\omega_{1},\tilde{\omega},-\omega_{0}-\omega_{1}-\tilde{\omega}\right)\bigg{\\}}$
$\displaystyle+\frac{\hat{T}^{(2)}\left(\omega_{1},\tilde{\omega},-\omega_{0}-\omega_{1}-\tilde{\omega}\right)}{\hat{R}^{(0)}\left(\omega_{1}\right)}-\frac{2\hat{T}^{(1)}\left(-\tilde{\omega}-\omega_{0};\tilde{\omega}\right)\hat{R}^{(1)}\left(\omega_{1},-\omega_{0}-\omega_{1}-\tilde{\omega}\right)}{\hat{R}^{(0)}\left(\omega_{1}\right)\hat{R}^{(0)}\left(-\tilde{\omega}-\omega_{0}\right)}\bigg{]}\hat{y}(-\tilde{\omega})\hat{y}(\omega_{0}+\omega_{1}+\tilde{\omega})+...$
(39)

At this stage, we have obtained a perturbative expansion of the first-passage
time density for a particular realisation ofwhitea particular realisation of
$y$, and only in terms of the expansion coefficients of transition and return
probability $T$ and $R$. We now give a more explicit expression of these
expansion coefficients.

###  II.2 Finding the coefficient terms for probability densities in the
functional expansion

In this section we show how the functional expansion of transition and return
probability are obtained perturbatively in terms of some suitable
eigenfunctions.

The transition probability $T_{x_{0},x_{1}}(t_{0},t_{1})$ of the undriven
process, characterised by Langevin equation (1) for $\varepsilon=0$, depends
on the time-difference only and can therefore be written
$T_{x_{0},x_{1}}(t_{0},t_{1})=T^{(0)}_{x_{0},x_{1}}(t_{1}-t_{0})$; following
definitions given in Eqs. (21) and (29). The transition density $T^{(0)}$
solves the Kolmogorov forward equation

$\displaystyle\begin{cases}\partial_{t}T^{(0)}_{x_{0},x_{1}}(t)&=\mathcal{L}_{x_{1}}T^{(0)}_{x_{0},x_{1}}(t)\qquad
t>0\\\ T^{(0)}_{x_{0},x_{1}}(t=0)&=\delta(x_{1}-x_{0})\\\
T^{(0)}_{x_{0},x_{1}}(t)&=0\qquad t<0\end{cases}$ (40)

where we introduce the forward evolution operator $\mathcal{L}_{x}$ as

$\displaystyle\mathcal{L}_{x}f$
$\displaystyle=\partial_{x}(V^{\prime}(x)f)+D_{x}\partial_{x}^{2}f,$ (41)

where $f$ is a twice differentiable test function. In Eq. (40), we denote the
forward operator as $\mathcal{L}_{x_{1}}$ to indicate that its gradient terms
are acting on the $x_{1}$ dependency of $T^{(0)}$. Correspondingly, the
$L^{2}$-adjoint operator $\mathcal{L}_{x}^{\dagger}$, also referred to as
backward operator, is

$\displaystyle\mathcal{L}_{x}^{\dagger}f$
$\displaystyle=-V^{\prime}(x)\partial_{x}f+D_{x}\partial_{x}^{2}f.$ (42)

The forward operator $\mathcal{L}_{x}$ has a countable set of eigenfunctions
$\left\\{u_{n}(x)\right\\}$ and a non-positive spectrum
$0\geq-\lambda_{0}>-\lambda_{1}>\ldots$ [46],

$\displaystyle\mathcal{L}_{x}u_{n}(x)=-\lambda_{n}u_{n}(x)\qquad n\geq 0$ (43)

but is a priori not self-adjoint in $L^{2}(\mathbb{R})$ 111Instead,
$\mathcal{L}_{x}^{\dagger}$ is self-adjoint and non-positive in the weighted
space $L^{2}(u_{0}(x))$, [46].. It is straightforward to show [47] however
that the operator

$\displaystyle\mathfrak{L}_{x}=e^{\frac{V(x)}{2D_{x}}}\mathcal{L}_{x}e^{-\frac{V(x)}{2D_{x}}}$
(44)

is self-adjoint and that therefore the family of

$\displaystyle\left\\{e^{\frac{V(x)}{2D_{x}}}u_{n}(x)\right\\},n\geq 0$ (45)

as eigenfunctions of $\mathfrak{L}_{x}$, form an orthogonal set of
eigenfunctions spanning $L^{2}(\mathbb{R})$. We further impose a choice of
normalization of the functions $\\{u_{n}\\}$, such that the set
$\left\\{e^{\frac{V(x)}{2D_{x}}}u_{n}(x)\right\\}$ is orthonormal. Defining
$u_{n}(x)$ as right eigenfunctions, and

$\displaystyle v_{n}(x)=e^{\frac{V(x)}{D_{x}}}u_{n}(x)$ (46)

as left eigenfunctions, satisfying

$\displaystyle\mathcal{L}_{x}^{\dagger}v_{n}(x)=-\lambda_{n}v_{n}(x),$ (47)

we obtain a biorthogonal system for the pair $\\{u_{n}(x)\\}$,
$\\{v_{n}(x)\\}$:

$\displaystyle\int\mathrm{d}{x}\,v_{m}(x)u_{n}(x)=\delta_{m,n}.$ (48)

This is useful to solve the forward equation; taking the Fourier transform in
time of Eq. (40), one obtains

$\displaystyle
i\omega\hat{T}^{(0)}_{x_{0},x_{1}}(\omega)=\mathcal{L}_{x}\hat{T}^{(0)}_{x_{0},x_{1}}(\omega)+\delta(x_{1}-x_{0}).$
(49)

Inserting the ansatz

$\displaystyle\hat{T}^{(0)}_{x_{0},x_{1}}(\omega)=\sum_{n\geq
0}\hat{T}^{(0)}_{n,x_{0}}(\omega)u_{n}(x_{1})$ (50)

into Eq. (49) and using Eq. (43) leads to

$\displaystyle\sum_{n\geq
0}(i\omega+\lambda_{n})\hat{T}^{(0)}_{n,x_{0}}(\omega)u_{n}(x_{1})=\sum_{n\geq
0}v_{n}(x_{0})u_{n}(x_{1})\leavevmode\nobreak\ ,$ (51)

where we made use of the decomposition of unity,

$\displaystyle\delta(x_{1}-x_{0})=\sum_{n\geq
0}v_{n}(x_{0})u_{n}(x_{1})\leavevmode\nobreak\ .$ (52)

Since the $u_{n}(x_{1})$ are linearly independent, their prefactors in Eq.
(51) need to agree. Therefore,

$\displaystyle\hat{T}^{(0)}_{n,x_{0}}(\omega)=\frac{v_{n}(x_{0})}{i\omega+\lambda_{n}}$
(53)

implying, together with Eq. (50),

$\displaystyle\hat{T}^{(0)}_{x_{0},x_{1}}(\omega)=\sum_{n\geq
0}\frac{v_{n}(x_{0})u_{n}(x_{1})}{i\omega+\lambda_{n}}\leavevmode\nobreak\ .$
(54)

Turning to the case of $\varepsilon\neq 0$, the transition probability of the
driven Langevin equation (1) solves the forward equation

$\displaystyle\partial_{t_{1}}T_{x_{0},x_{1}}(t_{0},t_{1};[y])$
$\displaystyle=$ $\displaystyle\left(\mathcal{L}_{x_{1}}+\varepsilon
y(t_{1})\partial_{x_{1}}\right)$ $\displaystyle
T_{x_{0},x_{1}}(t_{0},t_{1};[y]),\quad t_{1}>t_{0}\ $ $\displaystyle
T_{x_{0},x_{1}}(t_{0},t_{0};[y])$ $\displaystyle=\delta(x_{1}-x_{0})$
$\displaystyle\ T_{x_{0},x_{1}}(t_{0},t_{1};[y])$ $\displaystyle=0\quad
t_{1}<t_{0}$ (55)

where time-homogeneity can no longer be assumed. Under Fourier transform, this
forward equation becomes

$\left(i\omega_{1}-\mathcal{L}_{x_{1}}\right)\hat{T}_{x_{0},x_{1}}(\omega_{0},\omega_{1},[\hat{y}])\\\
=\delta(x_{1}-x_{0})\delta\bar{}\hskip 1.00006pt\,(\omega_{1}+\omega_{0})\\\
+\varepsilon\int\mathrm{d}\bar{}\hskip
1.00006pt\tilde{\omega}_{1}\,\hat{y}(\tilde{\omega}_{1})\partial_{x_{1}}\hat{T}_{x_{0},x_{1}}(\omega_{0},\omega_{1}-\tilde{\omega}_{1},[\hat{y}])$
(56)

where the $y$-dependent term turns from a product into a convolution under the
Fourier transform. We develop a perturbative solution of
$\hat{T}_{x_{0},x_{1}}(\omega_{0},\omega_{1},[\hat{y}])$ in powers of
$\hat{y}$.

We first discuss the first-order correction in $\varepsilon$. Following the
functional expansion (31), and using the zeroth order result (54), we obtain:

$\hat{T}_{x_{0},x_{1}}(\omega_{0},\omega_{1},[\hat{y}])=\sum_{n\geq
0}\frac{v_{n}(x_{0})u_{n}(x_{1})\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{1})}{i\omega_{1}+\lambda_{n}}\\\
+\varepsilon\int\mathrm{d}\bar{}\hskip
1.00006pt\tilde{\omega}_{1}\hat{T}^{(1)}_{x_{0},x_{1}}(\omega_{1},\tilde{\omega}_{1})\hat{y}(-\tilde{\omega}_{1})\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{1}+\tilde{\omega}_{1})+...$ (57)

and we wish to determine the function $\hat{T}^{(1)}$. Since the $u_{n}(x)$
span the $L^{2}$-space, and in analogy to the ansatz (50), the first-order
correction too can be written as a sum

$\displaystyle\hat{T}^{(1)}_{x_{0},x_{1}}(\omega_{1},\tilde{\omega}_{1})=\sum_{n\geq
0}\hat{T}^{(1)}_{n,x_{0}}(\omega_{1},\tilde{\omega}_{1})u_{n}(x_{1}).$ (58)

Re-inserting Eqs. (57) and (58) into Eq. (56) causes all terms to zeroth order
in $y$ to cancel, and one obtains an equation relating the contributions
proportional to $\varepsilon$,

$\sum_{n\geq 0}(i\omega_{1}+\lambda_{n})\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{1}}\hat{T}^{(1)}_{n,x_{0}}(\omega_{1},\tilde{\omega}_{1})\hat{y}(-\tilde{\omega}_{1})u_{n}(x_{1})\\\
=\int\mathrm{d}\bar{}\hskip 1.00006pt{\tilde{\omega}_{1}}\sum_{n\geq
0}\frac{v_{n}(x_{0})\partial_{x_{1}}u_{n}(x_{1})}{i(\omega_{1}-\tilde{\omega}_{1})+\lambda_{n}}\hat{y}(\tilde{\omega}_{1}).$
(59)

The right hand side, which is the convolution of
$\hat{T}^{(0)}\left(\omega\right)$ and $y(\omega)$, no longer sums over
$u_{n}(x_{1})$ but their derivative $\partial_{x_{1}}u_{n}(x_{1})$. In order
to compare both left and right terms, we need to express this sum as a sum
over the linearly independent $u_{n}(x_{1})$ again. The decomposition of the
derivative in terms of $u_{n}(x_{1})$ is given by

$\displaystyle\partial_{x_{1}}u_{n}(x_{1})=\sum_{k}\Delta_{nk}u_{k}(x_{1})$
(60)

where we refer to the $\Delta_{nk}$ as derivative coupling matrix whose
entries, as follows from bi-orthogonality, are

$\displaystyle\Delta_{nk}=\int\mathrm{d}{x}\,v_{k}(x)\partial_{x}u_{n}(x).$
(61)

Further, Eq. (59) holds true for arbitrary $y(\tilde{\omega}_{1})$. In order
to compare both integrations over $\mathrm{d}\bar{}\hskip
1.00006pt\tilde{\omega}_{1}$, we relabel
$\tilde{\omega}_{1}\mapsto-\tilde{\omega}_{1}$ in the left-hand side of the
equation. Using this notation, inserting the sum (60) into Eq. (59), and
resolving the ansatz (58), one obtains:

$\displaystyle\hat{T}^{(1)}_{x_{0},x_{1}}(\omega_{1},\tilde{\omega}_{1})=\sum_{n,k\geq
0}\frac{v_{k}(x_{0})\Delta_{kn}u_{n}(x_{1})}{(i(\omega_{1}+\tilde{\omega}_{1})+\lambda_{k})(i\omega_{1}+\lambda_{n})}\leavevmode\nobreak\
.$ (62)

In a similar way, the second order correction can be found; using the
functional expansion (31) to second order:

$\hat{T}_{x_{0},x_{1}}(\omega_{0},\omega_{1},[\hat{y}])=\sum_{n\geq
0}\frac{v_{n}(x_{0})u_{n}(x_{1})}{\omega_{1}+\lambda_{n}}\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{1})\\\ +\varepsilon\sum_{n,k\geq
0}\frac{v_{k}(x_{0})\Delta_{kn}u_{n}(x_{1})}{(-i\omega_{0}+\lambda_{k})(i\omega_{1}+\lambda_{n})}\hat{y}(\omega_{0}+\omega_{1})\\\
+\frac{\varepsilon^{2}}{2}\iint\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{1}}\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}_{2}}\hat{T}^{(2)}_{x_{0},x_{1}}(\omega_{1},\tilde{\omega}_{1},\tilde{\omega}_{2})\hat{y}(-\tilde{\omega}_{1})\hat{y}(-\tilde{\omega}_{2})\\\
\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{1}+\tilde{\omega}_{1}+\tilde{\omega}_{2})+\cdots\leavevmode\nobreak\
,$ (63)

with the results from Eqs. (54) and (62) to zeroth and first order, inserting
this ansatz into the forward equation (56) gives, following in complete
analogy to the previous steps,

$\hat{T}^{(2)}_{x_{0},x_{1}}(\omega_{1},\tilde{\omega}_{1},\tilde{\omega}_{2})=\\\
\sum_{n,k,m\geq
0}\frac{2v_{n}(x_{0})\Delta_{nk}\Delta_{km}u_{m}(x_{1})}{(i(\omega_{1}\\!+\\!\tilde{\omega}_{1}\\!+\\!\tilde{\omega}_{2})+\lambda_{n})(i(\omega_{1}\\!+\\!\tilde{\omega}_{1})+\\!\lambda_{k})(i\omega_{1}\\!+\\!\lambda_{m})}.$
(64)

Following this method, it is straightforward to generate the perturbative
terms of $\hat{T}^{(n)}$ to arbitrary order in $n$,

$\hat{T}^{(n)}_{x_{0},x_{1}}(\omega_{1},\tilde{\omega}_{1},...,\tilde{\omega}_{n})\\\
=n!\sum_{k_{0},...,k_{n}\geq
0}\frac{v_{k_{0}}(x_{0})\Delta_{k_{0}k_{1}}\cdot...\cdot\Delta_{k_{n-1}k_{n}}u_{k_{n}}(x_{1})}{\prod_{j=0}^{n}\left(i(\omega_{1}+\sum_{1\leq\ell\leq
n-j}\tilde{\omega}_{\ell})+\lambda_{k_{j}}\right)}.$ (65)

Finally, choosing $x_{0}=x_{1}$ in any of the expressions (54), (62), (64) and
(65) gives the corresponding terms for the return probability coefficients
$\hat{R}^{(0)}\left(\omega\right),\hat{R}^{(1)}\left(\omega,\tilde{\omega}_{1}\right),...$.
Equipped with these expressions, we are able to compute the relevant integrals
in the formula for the $y$-averaged first-passage time density Eq. (67).

###  II.3 Driving noise averaging

As was set out initially, the quantity of interest is the first-passage time
density when averaged over all driving noises (even if the quantity given
above might be of interest in itself). The average over driving noise
realisations $y_{t}$ is an average different to the average over the
stochastic process $x_{t}$. This double average is reminiscent of “quenched
disorder averages” where a system subject to thermal fluctuations is embedded
into a random disorder potential. The expansion in Eq. (II.1.4) is a power
series in orders of $\varepsilon$, where contributions of order
$\varepsilon^{n}$ contain an internal integration over $n-1$ free frequencies.
The expansion terms which stand in front of the $y$ terms, those denoted
within square brackets, are independent of $y$. They may be interpreted as the
$n$th order response functionals of the first-passage time distribution to
perturbations in the driving noise $y$. To calculate the $y$-average of
$\hat{p}\left(\omega_{0},\omega_{1};[y]\right)$, each term in Eq. (II.1.4) is
integrated over the path-measure of $y$, $\mathcal{P}[y]$. The order of
internal integration and $y$-averaging can be swapped. Consequently, since
$\overline{y_{t}}=0$ by assumption, all terms in first order in $y$ vanish. To
second order, correlations of $y$ come into play. We introduce the correlation
function

$\displaystyle\overline{\hat{y}(\omega_{0})\hat{y}(\omega_{1})}=\int\mathcal{D}[\hat{y}]\mathcal{P}[\hat{y}]\hat{y}(\omega_{0})\hat{y}(\omega_{1})=\hat{C}_{2}(\omega_{1})\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{1}),$ (66)

where the last equality arises from the stationary in time of the random
process $y(t)$. Stationarity in time also implies that $\hat{C}_{2}(\omega)$
is symmetric in $\omega\mapsto-\omega$. Performing an averaging over
$\mathcal{P}[\hat{y}]$ of Eq. (II.1.4), we obtain:

$\overline{\hat{p}(\omega_{0},\omega_{1};[\hat{y}])}=\frac{\hat{T}^{(0)}\left(\omega_{1}\right)}{\hat{R}^{(0)}\left(\omega_{1}\right)}\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{1})\\\
+\frac{1}{2}\varepsilon^{2}\left[-\underbrace{\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{\hat{T}^{(0)}\left(\omega_{1}\right)\hat{R}^{(2)}\left(\omega_{1},\tilde{\omega},-\tilde{\omega}\right)}{\left(\hat{R}^{(0)}\left(\omega_{1}\right)\right)^{2}}\hat{C}_{2}(\tilde{\omega})}_{=:({\rm
I})}+\underbrace{2\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{\hat{T}^{(0)}\left(\omega_{1}\right)\hat{R}^{(1)}\left(\omega_{1}-\tilde{\omega},\tilde{\omega}\right)\hat{R}^{(1)}\left(\omega_{1},-\tilde{\omega}\right)}{\left(\hat{R}^{(0)}\left(\omega_{1}\right)\right)^{2}\hat{R}^{(0)}\left(\omega_{1}-\tilde{\omega}\right)}\hat{C}_{2}(\tilde{\omega})}_{=:({\rm
II})}\right.\\\ \left.-\underbrace{2\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{\hat{T}^{(1)}\left(\omega_{1}-\tilde{\omega},\tilde{\omega}\right)\hat{R}^{(1)}\left(\omega_{1},-\tilde{\omega}\right)}{\hat{R}^{(0)}\left(\omega_{1}\right)\hat{R}^{(0)}\left(\omega_{1}-\tilde{\omega}\right)}\hat{C}_{2}(\tilde{\omega})}_{=:(\rm{III)}}+\underbrace{\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{\hat{T}^{(2)}\left(\omega_{1},\tilde{\omega},-\tilde{\omega}\right)}{\hat{R}^{(0)}\left(\omega_{1}\right)}\hat{C}_{2}(\tilde{\omega})}_{=:({\rm
IV})}\right]\delta\bar{}\hskip 1.00006pt\,(\omega_{0}+\omega_{1})+...$ (67)

The first term, of zeroth order, represents the Darling-Siegert solution (11).
This is consistent with our expansion around the base-point of no driving
noise (fully Markovian process). Once averaged, the second-order contribution
is again diagonal (i.e. proportional to $\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{2})$) indicating that the $y$-averaged first-
passage distribution is again invariant under time-shifts. The four correction
terms featuring in the second order expansion in $\varepsilon$ are labelled
$\rm(I)$ to $\rm(IV)$, and need to be calculated explicitly.

The corresponding expressions are derived in the following section, for the
case of Gaussian coloured noise. Generally, however, the entire expression Eq.
(67) remains valid for driving noise correlations $\hat{C}_{2}(\omega)$ which
keep the integrals finite.

###  II.4 Second-order correction with Gaussian coloured noise

So far, in our derivation of the second-order correction to the first-passage
time density, Eq. (67), we only demanded the active driving noise to be
stationary, with finite correlations and vanishing mean. In what follows, we
specify $y_{t}$ to be Gaussian coloured noise. This choice is almost canonical
in the study of coloured noise [23]. In our case it greatly simplifies the
necessary integrals. It is, however, possible to use any other correlation
functions as long as the integrals remain manageable. Generally, to compute
the perturbative contribution of $n$th order in $y$, the $n$-point correlation
function of $y_{t}$ needs to be known; For Gaussian processes all higher
moments follow from the two-point correlation function which simplifies the
calculation of potential higher order corrections. Since the correlation
function of coloured noise is an exponential, the results obtained in this
section to order $\varepsilon^{2}$ hold for anywhiteany noise with such
autocorrelation, in particular telegraphic noise as used in run-and-tumble
processes [38]. Run-and-tumble and other purely active processes correspond to
$D_{x}=0$, or alternatively $\nu\to\infty$, and are therefore not captured by
our perturbative approach valid for $\nu\lesssim 1$.

Gaussian coloured noise is defined by its exponential correlator,

$\displaystyle\overline{y_{t_{0}}y_{t_{1}}}=D_{y}\beta
e^{-\beta|t_{1}-t_{0}|},$ (68)

which in Fourier space reads (Eq. (66))

$\displaystyle\hat{C}_{2}(\omega)=\frac{2D_{y}\beta^{2}}{\omega^{2}+\beta^{2}}\leavevmode\nobreak\
.$ (69)

With the explicit expressions (54),(62), (64), we perform the integration in
${\rm(I)}$ (see Eq. (67) for notation) in eigenfunction-representation,

$\displaystyle{\rm(I)}$ $\displaystyle=\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{\hat{T}^{(0)}\left(\omega_{1}\right)\hat{R}^{(2)}\left(\omega_{1},\tilde{\omega},-\tilde{\omega}\right)}{\left(\hat{R}^{(0)}\left(\omega_{1}\right)\right)^{2}}\hat{C}_{2}(\tilde{\omega})$
$\displaystyle=2D_{y}\beta^{2}\,\frac{\hat{T}^{(0)}\left(\omega_{1}\right)}{\left(\hat{R}^{(0)}\left(\omega_{1}\right)\right)^{2}}\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{\hat{R}^{(2)}\left(\omega_{1},\tilde{\omega},-\tilde{\omega}\right)}{\beta^{2}+\tilde{\omega}^{2}}$
$\displaystyle=2D_{y}\beta^{2}\,\frac{\hat{T}^{(0)}\left(\omega_{1}\right)}{\left(\hat{R}^{(0)}\left(\omega_{1}\right)\right)^{2}}\sum_{n,k,m\geq
0}$ $\displaystyle\cdot 2\int_{-\infty}^{\infty}\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{v_{n}(x_{1})\Delta_{nk}\Delta_{km}u_{m}(x_{1})}{(i\omega_{1}+\lambda_{n})(i(\omega_{1}+\tilde{\omega})+\lambda_{k})(i\omega_{1}+\lambda_{m})}\frac{1}{\tilde{\omega}^{2}+\beta^{2}}$
$\displaystyle=D_{y}\beta\frac{\hat{T}^{(0)}\left(\omega_{1}\right)\hat{R}^{(2)}\left(\omega_{1},-i\beta,i\beta\right)}{\left(\hat{R}^{(0)}\left(\omega_{1}\right)\right)^{2}}.$
(70)

where in the last equality we employed Cauchy’s residue theorem closing the
contour in the lower half-plane containing the simple pole at
$\tilde{\omega}=-i\beta$. Likewise, we find

$\displaystyle{\rm(IV)}$ $\displaystyle=\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{\hat{T}^{(2)}\left(\omega_{1},\tilde{\omega},-\tilde{\omega}\right)}{\hat{R}^{(0)}\left(\omega_{1}\right)}\hat{C}_{2}(\tilde{\omega})$
$\displaystyle=2D_{y}\beta^{2}\frac{1}{\hat{R}^{(0)}\left(\omega_{1}\right)}\sum_{n,k,m\geq
0}$ $\displaystyle\cdot 2\int_{\infty}^{\infty}\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{v_{n}(x_{0})\Delta_{nk}\Delta_{km}u_{m}(x_{1})}{(i\omega_{1}+\lambda_{n})(i(\omega_{1}+\tilde{\omega})+\lambda_{k})(i\omega_{1}+\lambda_{m})}\frac{1}{\tilde{\omega}^{2}+\beta^{2}}$
$\displaystyle=D_{y}\beta\frac{\hat{T}^{(2)}\left(\omega_{1},-i\beta,i\beta\right)}{\hat{R}^{(0)}\left(\omega_{1}\right)}$
(71)

where again the integral is evaluated by closing the contour in the lower
half-plane enclosing the pole at $\tilde{\omega}=-i\beta$.

The integrals ${\rm(II)}$ and ${\rm(III)}$, featuring
$\tilde{\omega}$-dependent denominators, require some more careful analysis.
We have

$\displaystyle{\rm(II)}$ $\displaystyle=2\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{\hat{T}^{(0)}\left(\omega_{1}\right)\hat{R}^{(1)}\left(\omega_{1}-\tilde{\omega},\tilde{\omega}\right)\hat{R}^{(1)}\left(\omega_{1},-\tilde{\omega}\right)\hat{C}_{2}(\tilde{\omega})}{\left(\hat{R}^{(0)}\left(\omega_{1}\right)\right)^{2}\hat{R}^{(0)}\left(\omega_{1}-\tilde{\omega}\right)}\leavevmode\nobreak\
.$ (72)

As before we calculate this integral by application of Cauchy’s residue
theorem, but now closing the contour in the upper half-plane. The numerator’s
poles all lie in the lower half-plane with the exception of the pole at
$\tilde{\omega}=i\beta$ stemming from the correlator, as can be confirmed by
inspection of Eq. (62). Besides, the denominator
$\hat{R}^{(0)}\left(\omega_{1}-\tilde{\omega}\right)$ does not have any roots
for $\Im(\tilde{\omega})>0$. Indeed, using the relations Eq. (54) and (46),
$\hat{R}^{(0)}(\omega_{1}-\tilde{\omega})=\sum_{n\geq
0}\frac{v_{n}(x_{1})u_{n}(x_{1})}{i(\omega_{1}-\tilde{\omega})+\lambda_{n}}=e^{-V(x_{1})/D}\sum_{n\geq
0}\frac{\left(u_{n}(x_{1})\right)^{2}}{i(\omega_{1}-\tilde{\omega})+\lambda_{n}}$.
Since the $\\{u_{n}(x)\\}$ span $L^{2}$, there cannot be a $x_{1}$ for which
all $u_{n}(x_{1})=0$. Besides the real part of the denominator in the sum is
strictly positive, $\lambda_{n}+\Im(\tilde{\omega})>0$, assuming
$\omega_{1}\in\mathbb{R}$. In the upper half plane the sum therefore only
contains terms whose real part is non-negative and at least once strictly
positive; hence the sum is free of roots in the upper half plane. Invoking
Cauchy’s residue formula the integral is then given by

$\displaystyle{\rm(II)}=2D_{y}\beta\frac{\hat{T}^{(0)}\left(\omega_{1}\right)\hat{R}^{(1)}\left(\omega_{1}-i\beta,i\beta\right)\hat{R}^{(1)}\left(\omega_{1},-i\beta\right)}{\left(\hat{R}^{(0)}\left(\omega_{1}\right)\right)^{2}\hat{R}^{(0)}\left(\omega_{1}-i\beta\right)}\
.$ (73)

By analogous reasoning one obtains

$\displaystyle{\rm(III)}$ $\displaystyle=2\int\mathrm{d}\bar{}\hskip
1.00006pt{\tilde{\omega}}\frac{\hat{T}^{(1)}\left(\omega_{1}-\tilde{\omega},\tilde{\omega}\right)\hat{R}^{(1)}\left(\omega_{1},-\tilde{\omega}\right)}{\hat{R}^{(0)}\left(\omega_{1}\right)\hat{R}^{(0)}\left(\omega_{1}-\tilde{\omega}\right)}\hat{C}_{2}(\tilde{\omega})$
$\displaystyle=2D_{y}\beta\frac{\hat{T}^{(1)}\left(\omega_{1}-i\beta,i\beta\right)\hat{R}^{(1)}\left(\omega_{1},-i\beta\right)}{\hat{R}^{(0)}\left(\omega_{1}\right)\hat{R}^{(0)}\left(\omega_{1}-i\beta\right)}\leavevmode\nobreak\
.$ (74)

All four terms together give a general formula for the moment generating
function of first-passage times for arbitrary underlying processes and driving
noiseswhitearbitrary underlying processes and driving noises provided the
eigenfunctions and correlators are known. In the case of driving noise with
exponentially decaying auto-correlation, the full formula for
$\overline{\hat{p}\left(\omega_{0},\omega_{1}\right)}=\hat{F}({\omega_{1}})\delta\bar{}\hskip
1.00006pt\,(\omega_{0}+\omega_{1})$, reads

$\displaystyle\hat{p}({\omega})$
$\displaystyle=\frac{\hat{T}^{(0)}\left(\omega\right)}{\hat{R}^{(0)}\left(\omega\right)}+\frac{\varepsilon^{2}D_{y}\beta}{2\hat{R}^{(0)}\left(\omega\right)}\left[\frac{\hat{T}^{(0)}\left(\omega\right)}{\hat{R}^{(0)}\left(\omega\right)}\left(2\frac{\hat{R}^{(1)}\left(\omega-i\beta,i\beta\right)\hat{R}^{(1)}\left(\omega,-i\beta\right)}{\hat{R}^{(0)}\left(\omega-i\beta\right)}-\hat{R}^{(2)}\left(\omega,-i\beta,i\beta\right)\right)\right.$
$\displaystyle\left.\qquad\qquad\qquad+\hat{T}^{(2)}\left(\omega,-i\beta,i\beta\right)-2\frac{\hat{T}^{(1)}\left(\omega-i\beta,i\beta\right)\hat{R}^{(1)}\left(\omega,-i\beta\right)}{\hat{R}^{(0)}\left(\omega-i\beta\right)}\right]+\mathcal{O}(\varepsilon^{4})\leavevmode\nobreak\
.$ (75)

This general result concludes this section. In the next section, we consider
two concrete examples to demonstrate how this perturbation theory can be
turned into analytical results.

##  III  Results for simple potentials

###  III.1 Active Thermal Ornstein Uhlenbeck Process (ATOU)

In this example we study the case of a particle in a harmonic potential driven
by white and coloured noise described by the Langevin Equation

$\displaystyle\dot{x}_{t}=-\alpha x_{t}+\xi_{t}+\varepsilon y_{t}$ (76)

with driving noise correlator (see Eq. (68))

$\displaystyle\overline{y_{t_{0}}y_{t_{1}}}=D_{y}\beta
e^{-\beta|t_{1}-t_{0}|}.$ (77)

This process reduces to the simple Ornstein Uhlenbeck process when
$\varepsilon=0$ which models a particle in a harmonic potential
$(V(x)=\frac{\alpha}{2}x^{2}$) within a thermal bath. We consider, however,
the process driven by an additional “active” term $\varepsilon y_{t}$. We
therefore refer to this process as active thermal Ornstein Uhlenbeck process
(ATOU). In the undriven case ($\varepsilon=0$), the dynamics are characterised
by the time and length-scales $\alpha^{-1}$ and

$\ell=\sqrt{D_{x}\alpha^{-1}}.$

####  III.1.1 From eigenfunctions to the moment generating function of first-
passage times

The Fokker-Planck equation associated (cf. Eq. (41)) to the Langevin Equation
(76) is

$\displaystyle\partial_{t}T^{(0)}_{x_{0},x_{1}}=\mathcal{L}_{x_{1}}T^{(0)}_{x_{0},x_{1}}$
$\displaystyle=\left(D_{x}\partial_{x_{1}}^{2}+\alpha
x_{1}\partial_{x_{1}}+\alpha\right)T^{(0)}_{x_{0},x_{1}}\leavevmode\nobreak\
.$ (78)

with $\mathcal{L}_{x}=D_{x}\partial_{x}^{2}+\alpha x\partial_{x}+\alpha$ the
forward evolution operator. The operator has eigenvalues [46]

$\displaystyle\lambda_{n}=\alpha n$ (79)

and is diagonalised by the (normalised) eigenfunctions

$\displaystyle v_{n}(x)$ $\displaystyle=\frac{1}{\sqrt{\sqrt{2\pi}\ell\cdot
n!}}\operatorname{He}_{n}\left(\frac{x}{\ell}\right)$ (80) $\displaystyle
u_{n}(x)$ $\displaystyle=\frac{1}{\sqrt{\sqrt{2\pi}\ell\cdot
n!}}\operatorname{He}_{n}\left(\frac{x}{\ell}\right)\cdot
e^{-\frac{x^{2}}{2\ell^{2}}}$ (81)

where we introduced Hermite polynomials using the convention

$\displaystyle\operatorname{He}_{n}(x)=(-1)^{n}e^{\frac{{x}^{2}}{2}}\frac{\mathrm{d}^{n}}{\mathrm{d}x^{n}}e^{-\frac{{x}^{2}}{2}}$
(82)

which satisfy the relation

$\displaystyle\partial_{x}\operatorname{He}_{n}(x)=x\operatorname{He}_{n}(x)-\operatorname{He}_{n+1}(x)$
(83)

such that the coupling matrix (cf. Eq. (61)) resolves to

$\displaystyle\Delta_{mn}=-\frac{\sqrt{n}}{\ell}\delta_{m+1,n}\leavevmode\nobreak\
.$ (84)

The coupling matrix is not diagonal: incoming momentum $n$ is upgraded to
outgoing momentum $n+1$ by the noise coupling. The amplitude of $\Delta_{m,n}$
grows like $|\Delta_{m,n}|\sim\sqrt{\lambda_{m}}$. For later use, we also note
that

$\displaystyle\left(\frac{x}{\ell}-\ell\partial_{x}\right)v_{n}(x)$
$\displaystyle=\sqrt{n+1}\,v_{n+1}(x)$ (85)
$\displaystyle\ell\partial_{x}v_{n}(x)$
$\displaystyle=\sqrt{n}\,v_{n-1}(x)\leavevmode\nobreak\ .$ (86)

By $L^{2}$-adjointness, it follows that the adjoint creation and annihilation
operators are

$\displaystyle-\ell\partial_{x}u_{n}(x)$
$\displaystyle=\sqrt{n+1}\,u_{n+1}(x)$ (87)
$\displaystyle\left(\frac{x}{\ell}+\ell\partial_{x}\right)u_{n}(x)$
$\displaystyle=\sqrt{n}u_{n-1}(x)\leavevmode\nobreak\ .$ (88)

In order to compute the transition and return probabilities, the following
identity [48]

$\sum_{k=0}^{\infty}\frac{\operatorname{He}_{k}(x)\operatorname{He}_{k}(y)e^{-\frac{y^{2}}{2}}}{k!}z^{k}=\frac{1}{\sqrt{1-z^{2}}}e^{-\frac{(y-zx)^{2}}{2(1-z^{2})}}$
(89)

proves to be useful.

We introduce all quantities in dimensionless form, and take $\omega$ to be
imaginary, such as the reduced frequency $\bar{s}$ and the reduced
autocorrelation time $\bar{\beta}$ as

$\displaystyle\bar{s}=i\alpha^{-1}\omega\qquad\bar{\beta}=\alpha^{-1}\beta$
(90)

By considering the Fourier transform at $\bar{s}=i\alpha^{-1}\omega$, we are
effectively studying the Laplace transform. This is intended since our final
observable is the moment generating function, the Laplace transform of the
first-passage time distribution. Further, we denoted the lengths rescaled by
$\ell$ as

$\displaystyle\bar{x}_{1}=\ell^{-1}x_{1}\qquad\bar{x}_{0}=\ell^{-1}x_{0}.$
(91)

In the discussion that follows, we analyse all densities as densities in these
dimensionless quantities to simplify notation and calculations. Following Eq.
(54), one obtains for the transition probability

$\displaystyle\hat{T}^{(0)}\left(\omega=-i\bar{s}\alpha\right)$ (92)
$\displaystyle=\frac{1}{\sqrt{2\pi}\ell}\sum_{n=0}\frac{\operatorname{He}_{n}(\bar{x}_{0})\operatorname{He}_{n}(\bar{x}_{1})e^{-\frac{{\bar{x}_{1}}^{2}}{2}}}{n!(\bar{s}\alpha+\alpha
n)}$
$\displaystyle=\frac{1}{\sqrt{2\pi}\ell\alpha}\int_{0}^{\infty}\mathrm{d}{t}\,\,e^{-\bar{s}t}\sum_{n=0}^{\infty}\frac{\operatorname{He}_{n}(\bar{x}_{0})\operatorname{He}_{n}(\bar{x}_{1})e^{-\frac{{\bar{x}_{1}}^{2}}{2}}}{n!}\left(e^{-t}\right)^{n}$
$\displaystyle=\frac{1}{\sqrt{2\pi}\ell\alpha}\int_{0}^{\infty}\mathrm{d}{t}\,\frac{e^{-\bar{s}t}}{\sqrt{1-e^{-2t}}}e^{-\frac{(\bar{x}_{1}-\bar{x}_{0}e^{-t})^{2}}{2(1-e^{-2t})}}$

where we used identity (89) setting $z=e^{-t}$. This integral is the Laplace
transform of the Ornstein-Uhlenbeck propagator (in $t\leftrightarrow s$), and
its value is known in the literature to be ([14])

$\displaystyle\hat{T}^{(0)}\left(-i\bar{s}\alpha\right)$ (93)
$\displaystyle=\begin{cases}\frac{\Gamma(\bar{s})}{\sqrt{2\pi}\ell\alpha}e^{\frac{{\bar{x}_{0}}^{2}-{\bar{x}_{1}}^{2}}{4}}D_{-\bar{s}}\left(-\bar{x}_{0}\right)D_{-\bar{s}}\left(\bar{x}_{1}\right)&\bar{x}_{0}<\bar{x}_{1}\\\
\frac{\Gamma(\bar{s})}{\sqrt{2\pi}\ell\alpha}e^{\frac{{\bar{x}_{0}}^{2}-{\bar{x}_{1}}^{2}}{4}}D_{-\bar{s}}\left(\bar{x}_{0}\right)D_{-\bar{s}}\left(-\bar{x}_{1}\right)&\bar{x}_{0}>\bar{x}_{1}\end{cases}$

where we introduced the parabolic cylinder functions $D_{-\bar{s}}(x)$ ([49]).
By continuity, for $\bar{x}_{0}\to\bar{x}_{1}$ it follows that

$\displaystyle\hat{R}^{(0)}\left(-i\bar{s}\alpha\right)=\frac{\Gamma(\bar{s})}{\sqrt{2\pi}\ell\alpha}D_{-\bar{s}}(\bar{x}_{1})D_{-\bar{s}}(-\bar{x}_{1}).$
(94)

In order to compute
$\hat{T}^{(1)}\left(-i\bar{s}\alpha,-i\bar{\beta}\alpha\right)$, we use Eq.
(62) and the derivative coupling matrix computed in (84) to find

$\displaystyle\hat{T}^{(1)}\left(-i\bar{s}\alpha,-i\bar{\beta}\alpha\right)$
(95)
$\displaystyle=\frac{1}{\ell\alpha^{2}}\sum_{n}\frac{v_{n}(\bar{x}_{0})(-\sqrt{n+1})u_{n+1}(\bar{x}_{1})}{(\bar{s}+\bar{\beta}+n)(\bar{s}+n+1)}$
$\displaystyle=\frac{1}{\ell\alpha^{2}(1-\bar{\beta})}$
$\displaystyle\qquad\times\partial_{\bar{x}_{1}}\sum_{n}v_{n}(\bar{x}_{0})u_{n}(\bar{x}_{1})\left[\frac{1}{\bar{s}+\bar{\beta}+n}-\frac{1}{\bar{s}+n+1}\right]$
$\displaystyle=\frac{\partial_{\bar{x}_{1}}}{\ell\alpha(1-\bar{\beta})}\left[\hat{T}^{(0)}\left(-i(\bar{s}+\bar{\beta})\alpha\right)-\hat{T}^{(0)}\left(-i(\bar{s}+1)\alpha\right)\right]$

where we made use of relation (87) in the second equality. Letting
$\bar{x}_{0}\to\bar{x}_{1}$, one obtains
$\hat{R}^{(1)}\left(-i\bar{s}\alpha,-i\bar{\beta}\alpha\right)$. The
counterpart,
$\hat{T}^{(1)}\left(-i\bar{s}\alpha-i\bar{\beta}\alpha,i\bar{\beta}\alpha\right)$,
is similarly found to be

$\hat{T}^{(1)}\left(-i\bar{s}\alpha-i\bar{\beta}\alpha,i\bar{\beta}\alpha\right)\\\
=\frac{\partial_{\bar{x}_{1}}}{\ell\alpha(1+\bar{\beta})}\left[\hat{T}^{(0)}(-i\bar{s}\alpha)-\hat{T}^{(0)}(-i(\bar{s}+\bar{\beta}+1)\alpha)\right]$
(96)

These terms can be explicitly calculated and simplified. The rather lengthy
but explicit expressions are given in appendix D.

For the second order derivative term, using formula (64), one finds

$\displaystyle\hat{T}^{(2)}\left(-i\bar{s}\alpha,-i\bar{\beta}\alpha,i\bar{\beta}\alpha\right)$
(97)
$\displaystyle=\frac{2}{\ell^{2}\alpha^{3}}\sum_{n=0}^{\infty}\frac{v_{n}(\bar{x}_{0})\sqrt{n+1}\sqrt{n+2}u_{n+2}(\bar{x}_{1})}{(\bar{s}+n)(\bar{s}+\bar{\beta}+n+1)(\bar{s}+n+2)}$
$\displaystyle=\frac{2\partial_{\bar{x}_{1}}^{2}}{\ell^{2}\alpha^{3}}\sum_{n=0}^{\infty}\frac{v_{n}(\bar{x}_{0})u_{n}(\bar{x}_{1})}{(\bar{s}+n)(\bar{s}+\bar{\beta}+n+1)(\bar{s}+n+2)}$
$\displaystyle=\frac{2\partial_{\bar{x}_{1}}^{2}}{\ell^{2}\alpha^{2}}\left[\frac{1}{2(\bar{\beta}+1)}\hat{T}^{(0)}\left(-i\bar{s}\alpha\right)-\frac{1}{2(\bar{\beta}-1)}\hat{T}^{(0)}\left(-i(\bar{s}+2)\alpha\right)\right.$
$\displaystyle\qquad\left.+\frac{1}{({\bar{\beta}}^{2}-1)}\hat{T}^{(0)}\left(-i(\bar{s}+\bar{\beta}+1)\alpha\right)\right]$

Again, the evaluated terms, including for $\bar{x}_{0}\to\bar{x}_{1}$ are
given in appendix D.

Equipped with the return and transition probabilities and its first two
derivatives with respect to driving noise $y$ (cf. (93)-(97)), we obtain the
four contributions (II.4)-(74) which constitute the second-order correction
formula (II.4). Whilst all the explicit expressions are given in appendix D,
we here give the moment generating function in full as a undriven part and a
perturbative correction, using $s=\alpha\bar{s}=i\omega$,

$M_{x_{0},x_{1}}(s)=M^{0,{\rm
OU}}_{\frac{x_{0}}{\ell},\frac{x_{1}}{\ell}}\left(\frac{s}{\alpha}\right)+\underbrace{\frac{D_{y}\varepsilon^{2}\beta}{D_{x}\alpha}}_{=:\nu}M^{1,{\rm
OU}}_{\frac{x_{0}}{\ell},\frac{x_{1}}{\ell}}\left(\frac{s}{\alpha}\right)+\mathcal{O}(\nu^{2})$
(98)

where we introduced the dimensionless parameter of expansion
$\nu=\frac{\varepsilon^{2}D_{y}\beta}{D_{x}\alpha}.$

As is already known from literature (e.g. [15]),

$\displaystyle M^{0,{\rm
OU}}_{\bar{x}_{0},\bar{x}_{1}}(\bar{s})=\begin{cases}e^{\frac{{\bar{x}_{0}}^{2}-{\bar{x}_{1}}^{2}}{4}}\frac{D_{-\bar{s}}(-\bar{x}_{0})}{D_{-\bar{s}}(-\bar{x}_{1})}&\bar{x}_{0}<\bar{x}_{1}\\\
e^{\frac{{\bar{x}_{0}}^{2}-{\bar{x}_{1}}^{2}}{4}}\frac{D_{-\bar{s}}(\bar{x}_{0})}{D_{-\bar{s}}(\bar{x}_{1})}&\bar{x}_{0}>\bar{x}_{1}.\end{cases}$
(99)

By symmetry
$(\bar{x}_{0},\bar{x}_{1})\leftrightarrow(-\bar{x}_{1},-\bar{x}_{0})$ of the
problem and symmetry of driving noise, it suffices to regard one case only,
such that without loss of generality we assume $\bar{x}_{0}<\bar{x}_{1}$. The
central result of this section then is

$\displaystyle M^{1,{\rm OU}}_{\bar{x}_{0},\bar{x}_{1}}(\bar{s})$
$\displaystyle=\frac{\bar{s}e^{\frac{{\bar{x}_{0}}^{2}-{\bar{x}_{1}}^{2}}{4}}}{2\left({\bar{\beta}}^{2}-1\right)D_{-\bar{s}}(-\bar{x}_{1}){}^{2}D_{-{\bar{\beta}}-\bar{s}}(-\bar{x}_{1})}$
$\displaystyle\times\left[\left({\bar{\beta}}+1\right)(\bar{s}+1)D_{-{\bar{\beta}}-\bar{s}}(-\bar{x}_{1})\left(D_{-\bar{s}}(-\bar{x}_{0})D_{-\bar{s}-2}(-\bar{x}_{1})-D_{-\bar{s}-2}(-\bar{x}_{0})D_{-\bar{s}}(-\bar{x}_{1})\right)\right.$
$\displaystyle\qquad\left.-2({\bar{\beta}}+\bar{s})D_{-\bar{s}-1}(-\bar{x}_{1})\left(D_{-\bar{s}}(-\bar{x}_{0})D_{-{\bar{\beta}}-\bar{s}-1}(-\bar{x}_{1})-D_{-{\bar{\beta}}-\bar{s}-1}(-\bar{x}_{0})D_{-\bar{s}}(-\bar{x}_{1})\right)\right]$
(100)

Using for instance a computer algebra system like [50], all moments can be
obtained by differentiation and evaluating the limit of $\bar{s}\to 0$ at
which all derivatives have a removable singularity.

####  III.1.2 Numerical Validation

(a) Correction to mean first passage time of ATOU (cf. Eq. (76)) as obtained
from Eq. (106) versus target positions $x_{1}$, $x_{0}=0$ fixed, and various
values of $\nu$ (plot marks) compared to theoretical result to first order in
$\nu$ (black line) using Eq. (105) and the result obtained in (100). The inset
shows the mean first-passage time $\tau_{x_{0},x_{1}}$ as measured vs $x_{1}$
for values of $\nu=0$ to $0.8$. Correction due to active driving noise
increases MFPT for $x_{1}\lesssim 1.6$ and decreases MFPT for $x_{1}\gtrsim
1.6$. This behaviour is fully captured by the analytic result.

(b) Correction to mean squared first passage time of ATOU (cf. Eq. (76)) as
obtained from Eq. (107) versus target positions $x_{1}$, $x_{0}=0$ fixed, and
various values of $\nu$ (plot marks) compared to theoretical result to first
order in $\nu$ (black line) using Eq. (105) and the result obtained in (100).
The inset shows the mean squared first-passage time $\tau_{x_{0},x_{1}}$ as
measured vs $x_{1}$ for values of $\nu=0$ to $0.8$.

Figure 3: First order correction to first and second moment of ATOU.
Simulation parameters are $x_{0}=0,D_{x}=1,D_{y}=1,\alpha=1,\beta=\frac{1}{2}$
and $\varepsilon$ suitably chosen to fix
$\nu=D_{y}\beta\varepsilon^{2}/(D_{x}\alpha)$. Averages were taken over
$10^{6}$ samples.

In order to corroborate the closed form result of the first-order correction
to the moment generating function of first-passage times, Eq. (100), we employ
Monte Carlo simulations integrating the driven Langevin equation (76) $N\simeq
10^{6}$ times, numerically find the first-passage time $\tilde{\tau}_{i}$, and
average the moment generating function
$\tilde{M}=\frac{1}{N}\sum_{i=1}^{N}e^{-s\tilde{\tau}_{i}}$ for values of $s$
within the range $s\in[0,5]$. In the following, $\tilde{M}$ and
$\tilde{\mathcal{T}}$ denote quantities that have been numerically obtained.
Since we assume an expansion of the form $M=M^{0}+\nu M^{1}+\nu^{2}M^{2}+...$,
we take the numerical first derivative

$\displaystyle\tilde{M}_{1}|_{\nu}=\frac{\tilde{M}|_{\nu}-\tilde{M}|_{\nu=0}}{\nu}$
(101)

to verify our analytic prediction of $M^{1}$. In Fig. 4, the numerical
estimate $\tilde{M}_{1}$ is shown for various values of $\nu$, together with
the analytic expression Eq. (100) of the scaling function $M^{1}$. For small
$\nu$, the agreement is excellent. For larger values of $\nu$, higher-order
corrections become more visible. The next-higher contribution, which we did
not calculate analytically but which can be found by following the framework
to second order in $\nu$, is numerically estimated by taking the second
numerical derivative,

$\displaystyle\tilde{M}^{2}|_{\nu}=\frac{\tilde{M}^{1}|_{\nu}-M^{1}}{\nu},$
(102)

and is shown in the inset of Fig. 4. For $0.2\leq\nu\leq 0.8$, the second-
order corrections collapse, indicating that the deviations in the main figure
are well accounted for by second-order corrections. For $\nu=0.1$,
$\tilde{M}^{2}$ deviates slightly due to the statistical noise, since the
second order correction is very small.

For better physical intuition of the process, we numerically compute the
inverse Laplace transform [51, 50] of $M^{0}$ and $M^{1}$ to obtain the first-
passage time distribution:

$\displaystyle p(t)=p^{0}(t)+\nu p^{1}(t)$ (103)

neglecting, as usual, terms of higher order in $\nu$. In Fig. 2a, we show the
first-order corrected distribution compared to a numerical sampling of the
probability function, $\tilde{p}(t)$, for values of $0\leq\nu\leq 3.2$. The
agreement for $\nu\lesssim 1$ is excellent, and plotting the rescaled
deviation

$\displaystyle\tilde{p}^{1}(t)|_{\nu}=\frac{\tilde{p}(t)|_{\nu}-\tilde{p}(t)|_{\nu=0}}{\nu}$
(104)

against the theoretical result of $p^{1}(y)$ (see inset of Fig. 2a) shows
systematic higher-order corrections consistent with our framework. Since for
small $\nu$ the deviation to the undriven case is small, $\tilde{p}$ suffers
from statistical noise, therefore we omit showing $\tilde{p}(\nu=0.1)$ in the
inset.

These numerical results therefore confirm the analytically obtained first-
order correction to the moment-generating function; consequently, the
correction to all moments has been gained. As an illustration, we further show
the first and second moment of the Ornstein-Uhlenbeck process driven by
coloured noise in Fig. 3a and Fig. 3b. In analogy to the moment-generating
function, we measure the mean and mean square first passage times
$\tilde{\mathcal{T}}_{1}=\frac{1}{N}\sum_{i=1}^{N}\tilde{\tau}_{i}$,
$\tilde{\mathcal{T}}_{2}=\frac{1}{N}\sum_{i=1}^{N}\tilde{\tau}_{i}^{2}$, which
we assume to expand in $\nu$ as
$\mathcal{T}_{1}(\nu)=\mathcal{T}^{0}_{1}+\nu\mathcal{T}_{1}^{1}+\nu^{2}\mathcal{T}^{2}_{1}+\ldots$
and
$\mathcal{T}_{2}(\nu)=\mathcal{T}^{0}_{2}+\nu\mathcal{T}^{1}_{2}+\nu^{2}\mathcal{T}_{2}^{2}+\ldots$.
The first-order corrections introduced are obtained by differentiation wrt
$\bar{s}$

$\displaystyle\mathcal{T}^{1}_{n}=\left.\frac{\mathrm{d}\,{}^{n}}{\mathrm{d}{(-\bar{s})^{n}}\,}\right|_{\bar{s}=0}M^{1}$
(105)

using the result of Eq. (100) which is performed by a computer algebra system
and evaluated exactly. Due to their lengthiness, we do not give their full
expression here. In order to numerically confirm these predictions, we measure
the first order derivatives

$\displaystyle\tilde{\mathcal{T}}^{1}_{1}|_{\nu}$
$\displaystyle=\frac{\tilde{\mathcal{T}}_{1}|_{\nu}-\tilde{\mathcal{T}}_{1}|_{\nu=0}}{\nu}$
(106) $\displaystyle\tilde{\mathcal{T}}^{1}_{2}|_{\nu}$
$\displaystyle=\frac{\tilde{\mathcal{T}}_{2}|_{\nu}-\tilde{\mathcal{T}}_{2}|_{\nu=0}}{\nu}$
(107)

and compare it to the result obtained from Eq. (105). In Fig. 3a and Fig. 3b,
we show the resulting moments of first-passage times obtained for fixed start
position $x_{0}=0$ (at the minimum of the potential) but varied
$x_{1}\in[0.05,2]$. The figures show a clear agreement with the theoretical
result and systematic deviations for larger $\nu$. Further, we observe that in
this setting the coloured noise increases the mean-first passage time for
smaller distances ($x_{1}\lesssim 1.6$) and decreases it for larger distances.
This also holds true for the mean squared first-passage time. This example
therefore further illustrates that the effect of coloured driving (or memory)
on the Langevin dynamics is highly non-trivial, yet our framework is able to
capture this effect. The insets in both figures show the measured moments
$\tilde{\mathcal{T}}^{1},\tilde{\mathcal{T}}^{2}$.

Figure 4: Numerical validation of first order correction to FPT moment
generating function $M^{1}$ (cf. Eq. (4)) of ATOU (see Eq. (76) and Sec. III.1
for discussion). The result is calculated in Eq. (100). Numerical Simulations
are shown for various values of $0.1\leq\nu\leq 0.8$ (plot marks). The moment
generating functions were sampled for
$x_{0}=0,x_{1}=1,D_{x}=1,D_{y}=1,\alpha=1,\beta=\frac{1}{2}$ and $\varepsilon$
suitably chosen to fix $\nu=D_{y}\varepsilon^{2}\beta/(D_{x}\alpha)$. For
small values of $\nu$ agreement with theoretical first-order correction (black
line) is excellent. For larger values of $\nu$ the deviation increases. The
rescaled deviation, $\tilde{M}_{2}$ (see Eq. (102)), (inset) collapse and thus
confirm that these deviations are systematic higher-order corrections. See
(III.1.2) for further results and discussion.

###  III.2 Active Thermal Brownian Motion on a ring (ATBM)

In this subsection, we consider the case of a Brownian particle at $x_{t}$
driven by coloured noise which is placed on a ring of radius $r$, thus
satisfying periodic boundary conditions ($x\equiv x+2\pi r$). The position of
the particle satisfies the Langevin equation

$\displaystyle\dot{x}_{t}=\xi_{t}+\varepsilon y_{t}$ (108)

with

$\displaystyle\left\langle{\xi_{t_{0}}\xi_{t_{1}}}\right\rangle$
$\displaystyle=2D_{x}\delta(t_{1}-t_{0})$ (109)
$\displaystyle\overline{y_{t_{0}}y_{t_{1}}}$ $\displaystyle=D_{y}\beta
e^{-\beta|t_{1}-t_{0}|}.$ (110)

We refer to this system as Active Thermal Brownian Motion (ATBM) on a ring
(see Fig. 5 for illustration). In analogy to the previous subsection, we
derive the correction to the moment generating function of the first-passage
time distribution.

####  III.2.1 From eigenfunctions to the moment generating function of first-
passage times

Figure 5: A Particle on a circle of radius $r$ is driven by both white
(thermal) and coloured (active) noise (cf. Eq. (108)). We study the first
passage time distribution from $x_{0}$ to $x_{1}$ as a function of the angle
$\theta\in[0,2\pi)$.

The Fokker Planck associated to the forward equation corresponding to Langevin
Equation (108) is the Laplace equation with periodic boundary conditions
$x\equiv x+2\pi r$. Since its eigenfunctions are complex-valued square-
integrable functions, we follow the identical steps of the framework
introduced in Sec. II.2 but replace the scalar-product in Eq. (48) by its
complex equivalent

$\displaystyle\int\mathrm{d}{x}\,v_{m}^{*}(x)u_{n}(x)=\delta_{mn}.$ (111)

In what follows, we therefore consider the complex conjugate of the left
eigenfunctions and obtain

$\displaystyle u_{n}(x)=\frac{1}{\sqrt{2\pi r}}e^{i\frac{nx}{r}}$ (112)
$\displaystyle v_{n}^{*}(x)=\frac{1}{\sqrt{2\pi r}}e^{-i\frac{nx}{r}}$ (113)

with corresponding eigenvalues

$\displaystyle\lambda_{n}=D_{x}r^{-2}n^{2}\qquad(n\in\mathbb{Z}).$ (114)

where by allowing $n\in\mathbb{Z}$, we enumerate the eigenfunctions
efficiently (cf. Eq. (43)). The eigenfunctions are conjugate to each other
since the forward operator of simple diffusion,
$\mathcal{L}=D_{x}\partial_{x}^{2}$, is self-adjoint. From this follows that
the noise-coupling matrix defined in Eq. (61),

$\displaystyle\Delta_{mn}$
$\displaystyle=\int\mathrm{d}x\,v^{*}_{m}\partial_{x}u_{n}(x)$
$\displaystyle=\frac{1}{2\pi r}\int_{0}^{2\pi
r}\mathrm{d}{x}\,\,i\frac{n}{r}e^{i(n-m)\frac{x}{r}}=i\frac{n}{r}\delta_{mn}$
(115)

is diagonal and purely imaginary. Because of scale-invariance and rotational
symmetry, we simplify the following discussion by introducing the
dimensionless angle

$\displaystyle\theta\coloneqq\frac{x_{1}-x_{0}}{r},$ (116)

where we restrict ourselves to $\theta\in[0,2\pi)$, and the diffusive
timescale

$\displaystyle\alpha^{-1}=\frac{r^{2}}{D_{x}}$ (117)

with which we rescale the Fourier-frequency $\omega$ to

$\displaystyle\bar{s}=i\alpha^{-1}\omega,$ (118)

again effectively evaluating the Laplace transform (cf. Eq. (7)) of the
respective probability densities. With this simplified notation, the
transition density to zeroth order reads

$\displaystyle\hat{T}^{(0)}\left(-i\bar{s}\alpha\right)$
$\displaystyle=\frac{1}{2\pi\alpha
r}\sum_{n=-\infty}^{\infty}\frac{e^{-in\theta}}{\bar{s}+n^{2}}$
$\displaystyle=\frac{1}{2\alpha
r}\frac{\cosh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\sqrt{\bar{s}}\sinh\left(\pi\sqrt{\bar{s}}\right)},$
(119)

and the return probability, setting $\theta=0$, is

$\displaystyle\hat{R}^{(0)}\left(-i\bar{s}\alpha\right)$
$\displaystyle=\frac{1}{2\alpha
r}\frac{\cosh\left(\pi\sqrt{\bar{s}}\right)}{\sqrt{\bar{s}}\sinh\left(\pi\sqrt{\bar{s}}\right)}.$
(120)

Assuming the $y$-averaged moment generating function has an expansion of

$\displaystyle
M\left(s\right)=M^{0,\rm{BM}}_{\frac{x_{1}-x_{0}}{r}}\left(\frac{r^{2}s}{D_{x}}\right)+\underbrace{\frac{D_{y}\varepsilon^{2}\beta}{D_{x}\alpha}}_{\nu}M^{1,\rm{BM}}_{\frac{x_{1}-x_{0}}{r}}\left(\frac{r^{2}s}{D_{x}}\right)+\mathcal{O}(\nu^{2}),$
(121)

where we use as in the Ornstein-Uhlenbeck case the dimensionless perturbative
parameter $\nu=\frac{D_{y}\varepsilon^{2}\beta}{D_{x}\alpha}$ (here $\alpha$
has however a different definition than used in the OU case, see Eq. (117)).
The zeroth order contribution is, using the classic result Eq. (11), and the
results in (119), (120),

$\displaystyle
M^{0,\rm{BM}}_{\theta}(\bar{s})=\frac{\cosh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\cosh\left(\pi\sqrt{\bar{s}}\right)}$
(122)

which expands around $\bar{s}=0$ as

$\displaystyle M^{0,\rm{BM}}_{\theta}(\bar{s})=$ $\displaystyle
1+\left(\frac{\theta^{2}}{2}-\pi\theta\right)\bar{s}$ (123)
$\displaystyle+\frac{1}{24}\left(\theta^{4}-4\pi\theta^{3}+8\pi^{3}\theta\right)\bar{s}^{2}$
$\displaystyle+\frac{1}{720}\left(\theta^{6}-6\pi\theta^{5}+40\pi^{3}\theta^{3}-96\pi^{5}\theta\right)\bar{s}^{3}+\ldots$

Being a moment generating function, the prefactors in front of $\bar{s}^{n}$
correspond to the (rescaled) moments of the first-passage time
$\frac{(-1)^{n}}{n!}\left\langle{\tau_{x_{0}\to
x_{1}}^{n}}\right\rangle\left(\frac{D_{x}}{r^{2}}\right)^{n}$.

We turn to higher orders in $\varepsilon$. To first order, the transition
density

$\displaystyle\hat{T}^{(1)}\left(-i\bar{s}\alpha,-i\bar{\beta}\alpha\right)$
$\displaystyle=\frac{1}{2\pi\alpha^{2}r^{2}}\sum_{n=-\infty}^{\infty}\frac{ine^{-in\theta}}{(\bar{s}+\bar{\beta}+n^{2})(\bar{s}+n^{2})}$
$\displaystyle=\frac{1}{2\pi\alpha^{2}r^{2}}\sum_{n=0}^{\infty}\frac{-2n{\sin{n\theta}}}{(\bar{s}+\bar{\beta}+n^{2})(\bar{s}+n^{2})}$
(124)

Here, we introduced as before $\bar{\beta}=\alpha^{-1}\beta$ as the
dimensionless correlation time-scale of the coloured noise. The function
$\hat{R}^{1}(-i\bar{s}\alpha,-i\bar{\beta}\alpha)$, which corresponds to
$\hat{T}^{1}(-i\bar{s}\alpha,-i\bar{\beta}\alpha)$ evaluated at $\theta=0$,
vanishes. This implies that the contribution of integrals $({\rm II})$ and
$({\rm III})$, as given in Eq. (73) and (74), vanish, leaving only $({\rm I})$
and $({\rm IV})$ as correction terms.

To second order, the transition density is

$\displaystyle\hat{T}^{(2)}\left(-i\bar{s}\alpha,-i\bar{\beta}\alpha,i\bar{\beta}\alpha\right)$
$\displaystyle=-\frac{1}{\pi\alpha^{3}r^{3}}\sum_{n=-\infty}^{\infty}\frac{n^{2}e^{-in\theta}}{(\bar{s}+n^{2})^{2}(\bar{s}+\bar{\beta}+n^{2})}$
$\displaystyle=-\frac{1}{\pi\alpha^{3}r^{3}}\sum_{n=-\infty}^{\infty}\left[\frac{1}{{\bar{\beta}}^{2}}\frac{n^{2}e^{-in\theta}}{\bar{s}+\bar{\beta}+n^{2}}-\frac{1}{{\bar{\beta}}^{2}}\frac{n^{2}e^{-in\theta}}{\bar{s}+n^{2}}+\frac{1}{{\bar{\beta}}}\frac{n^{2}e^{-in\theta}}{(\bar{s}+n^{2})^{2}}\right]$
$\displaystyle=\frac{2}{\alpha^{2}r^{2}}\cdot\frac{1}{{\bar{\beta}}^{2}}\partial_{\theta}^{2}\left[\hat{T}^{(0)}\left(\bar{s}+\bar{\beta}\right)-\hat{T}^{(0)}\left(\bar{s}\right)-\bar{\beta}\partial_{\bar{s}}\hat{T}^{(0)}\left(\bar{s}\right)\right]$
$\displaystyle=\frac{2}{\alpha^{3}r^{3}}\cdot\frac{1}{4{\bar{\beta}}^{2}\sqrt{\bar{s}}}\left[\frac{\cosh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\sinh\left(\pi\sqrt{\bar{s}}\right)}\left(\pi\bar{\beta}\sqrt{\bar{s}}\coth\left(\pi\sqrt{\bar{s}}\right)-2\bar{s}-\bar{\beta}\right)\right.$
$\displaystyle\qquad\qquad\qquad+\sqrt{\bar{s}}\left.\left\\{2\sqrt{\bar{\beta}+\bar{s}}\frac{\cosh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\sinh(\pi\sqrt{\bar{\beta}+\bar{s})}}+\bar{\beta}(\pi-\theta)\frac{\sinh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\sinh(\pi\sqrt{\bar{s}}}\right\\}\right]$
(125)

where $\partial_{\theta}=r\partial_{x_{1}}$. Setting $\theta=0$, one obtains
the second order response of the return probability,

$\displaystyle\hat{R}^{(2)}\left(-i\bar{s}\alpha,-i\bar{\beta}\alpha,i\bar{\beta}\alpha\right)=\frac{1}{2\alpha^{3}{\bar{\beta}}^{2}r^{3}}\left[\pi\bar{\beta}\coth^{2}\left(\pi\sqrt{\bar{s}}\right)-\frac{(\bar{\beta}+2\bar{s})\coth\left(\pi\sqrt{\bar{s}}\right)}{\sqrt{\bar{s}}}+2\sqrt{\bar{\beta}+\bar{s}}\coth\left(\pi\sqrt{\bar{\beta}+\bar{s}}\right)-\pi\bar{\beta}\right],$
(126)

Inserting these quantities into the FPT correction (II.4), gives

$\displaystyle
M(s=\alpha\bar{s})=\frac{\cosh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\cosh\left(\pi\sqrt{\bar{s}}\right)}+\frac{D_{y}\varepsilon^{2}}{2D_{x}}\cdot\frac{\sqrt{\bar{s}}\tanh\left(\pi\sqrt{\bar{s}}\right)}{{\bar{\beta}}^{2}}\left[\frac{\cosh\left((\theta-\pi)\sqrt{\bar{\beta}+\bar{s}}\right)}{\sinh\left(\pi\sqrt{\bar{\beta}+\bar{s}}\right)}2\sqrt{\bar{\beta}+\bar{s}}\right.$
$\displaystyle\qquad+\left.\frac{\cosh\left((\theta-\pi)\sqrt{\bar{\beta}+\bar{s}}\right)}{\cosh\left(\pi\sqrt{\bar{s}}\right)}\left(\pi\bar{\beta}-2\sqrt{\bar{\beta}+\bar{s}}\coth\left(\pi\sqrt{\bar{\beta}+\bar{s}}\right)\right)+\frac{\sinh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\sinh\left(\pi\sqrt{\bar{s}}\right)}\bar{\beta}\left(\pi-\theta\right)\right]+\mathcal{O}(\varepsilon^{4})$
(127)

This expression is the full moment generating function up to second order in
$\varepsilon$whitefull moment generating function up to second order in
$\varepsilon$ (first order in $\nu$) in elementary functions. In line with our
previous notation, we identify the dimensionless scaling functions

$\displaystyle M^{0,\rm{BM}}_{\theta}(\bar{s})$
$\displaystyle=\frac{\cosh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\cosh\left(\pi\sqrt{\bar{s}}\right)}$
(128) $\displaystyle M^{1,\rm{BM}}_{\theta}(\bar{s})$
$\displaystyle=\frac{\sqrt{\bar{s}}\tanh\left(\pi\sqrt{\bar{s}}\right)}{2{\bar{\beta}}^{2}}\left[\frac{\cosh\left((\theta-\pi)\sqrt{\bar{\beta}+\bar{s}}\right)}{\sinh\left(\pi\sqrt{\bar{\beta}+\bar{s}}\right)}2\sqrt{\bar{\beta}+\bar{s}}\right.$
$\displaystyle\left.+\bar{\beta}(\pi-\theta)\frac{\sinh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\sinh\left(\pi\sqrt{\bar{s}}\right)}+\left(\pi\bar{\beta}-2\sqrt{\bar{\beta}+\bar{s}}\coth\left(\pi\sqrt{\bar{\beta}+\bar{s}}\right)\right)\frac{\cosh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\cosh\left(\pi\sqrt{\bar{s}}\right)}\right]$
(129)

which together form the first-order correction of the moment-generating
function

$\displaystyle
M(s)=M^{0,\rm{BM}}_{\frac{x_{1}-x_{0}}{r}}\left(\frac{r^{2}s}{D_{x}}\right)+\frac{D_{y}\varepsilon^{2}\beta}{D_{x}\alpha}M^{1,\rm{BM}}_{\frac{x_{1}-x_{0}}{r}}\left(\frac{r^{2}s}{D_{x}}\right)+\mathcal{O}(\varepsilon^{4})\
.$ (130)

An expansion around $\bar{s}=0$ gives corrections to all moments. The
correction to first order in $\nu$ of the mean first-passage time over an
angle of $\theta$, $\left\langle{\tau_{0\to\theta}}\right\rangle$, for
instance reads

$\displaystyle\alpha\left\langle{\tau_{0,\theta}}\right\rangle=\underbrace{\pi\theta-\frac{\theta^{2}}{2}}_{\mathcal{T}^{1}_{0}(\theta)}-\nu\underbrace{\frac{1}{2{\bar{\beta}}^{3/2}}\left(\sqrt{\bar{\beta}}(2\pi-\theta)\theta-2\pi\coth\left(\pi\sqrt{\bar{\beta}}\right)+2\pi\frac{\cosh\left(\sqrt{\bar{\beta}}(\theta-\pi)\right)}{\sinh(\pi\sqrt{\bar{\beta}})}\right)}_{=:\mathcal{T}^{1}_{1}(\theta,\bar{\beta})}+\mathcal{O}(\nu^{2})$
(131)

where we indicated that the result can be written as the classical
contribution ($\mathcal{T}^{1}_{0}$), plus the dimensionless perturbative
coefficient $\nu$ times a dimensionless scaling function
$\mathcal{T}^{1}_{1}(\theta,\bar{\beta})$. By successive derivation, any
higher order moment may be obtained from Eq. (127).

In the limit of an infinitely large radius $r\to\infty$ one obtains, under
suitable rescaling, the moment generating function for an active thermal
Brownian motion on the real line. We have confirmed that the perturbative
approximation Eq. (127) agrees, up to second order in $\varepsilon$, with the
analytic solution for a run-and-tumble process subject to additional white
noise in the limit of small tumble noise which has been calculated in [34].

####  III.2.2 Numerical Validation

In order to validate the analytic result of the moment generating function of
first-passage times Eq. (127) and the mean first-passage time Eq. (131), we
follow the same steps as in the previous section III.1.2. Using Monte-Carlo
simulations, we sample the first-passage times $\tilde{\tau}_{i}$ of the
integrated stochastic equation (108). To validate the moment-generating
function, we average over $N\simeq 10^{6}$ to $10^{7}$ iterations to sample

$\displaystyle\tilde{M}=\frac{1}{N}\sum_{i=1}^{N}e^{-s\tilde{\tau}_{i}}$ (132)

for various values of $\nu$ and across $\bar{x}_{1}-\bar{x}_{0}\in(0,\pi]$
(with $r=1$). Again, symbols with a tilde denote quantities which are measured
numerically. To validate the theoretically predicted first-order correction
Eq. (III.2.1), we introduce the numerically measured corrections:

$\displaystyle\tilde{M}^{1}|_{\nu}$
$\displaystyle=\frac{\tilde{M}|_{\nu}-\tilde{M}_{\nu=0}}{\nu}$ (133)
$\displaystyle\tilde{M}^{2}|_{\nu}$
$\displaystyle=\frac{\tilde{M}^{1}|_{\nu}-M^{1}}{\nu}\leavevmode\nobreak\ .$
(134)

In Fig. 7, we show the analytic result $M^{1}$ (cf. Eq. (III.2.1)) and
numerically obtained $\tilde{M}^{1}|_{\nu}$ (cf. Eq. (133)) for various values
of $0\leq\nu\leq 0.8$. The agreement is again excellent, and the discrepancy
between simulated result and theoretical first order correction grows, as
expected, with larger values of $\nu$. The rescaled discrepancy,
$\tilde{M}^{2}$, to leading order the second-order correction $M^{2}$ (cf. Eq.
(134)), is plotted in the inset and collapses, indicating that the discrepancy
is systematic and confirming the validity of the result in Eq. (III.2.1). By
numerically computing the inverse Laplace transform, we obtain the
perturbative correction to the first-passage time distribution (cf. Eq. (103))
which in Fig. 2b is compared to the numerically sampled distribution. The
inset shows the numerical estimate of the first order correction (cf. Eq.
(104)) which agrees well for small $\nu$. For larger values of $\nu$, second-
order corrections play an increasing role as indicated in Fig. 7.

(a) Correction to mean first passage time of ATBM on a ring (cf. Eq. (108)) as
obtained from Eq. (106) versus target position $x_{1}$ and various values of
$\nu$ (plot marks). This is compared to theoretical result of Eq. (131) (solid
black line). The inset shows the measured mean first passage time versus a
varying target position $x_{1}$ and different values of $\nu$.

(b) Correction to mean squared first passage time of ATBM on a ring (cf. Eq.
(108)) as obtained from Eq. (107) versus target position $x_{1}$ and various
values of $\nu$ (plot marks) . This is compared to theoretical result of twice
differentiating Eq. (III.2.1) (solid black line). The inset shows the measured
mean squared first passage time versus a varying target position $x_{1}$ and
different values of $\nu$.

Figure 6: First order correction to first and second moment of first-passage
times of ATBM with periodic boundary conditions. Simulation parameters are
$x_{0}=0,D_{x}=1,r=1,D_{y}=1,\beta=\frac{1}{2}$ and $\varepsilon$ suitably
chosen to fix $\nu=D_{y}r^{2}\beta\varepsilon^{2}/(D_{x}^{2})$.

In Figs. 6a and 6b, we show the first and second moment of first-passage times
and how their deviations to the $\nu=0$ case is captured by the first-order
correction obtained using Eq. (106) and Eq. (107) with the result of Eq.
(III.2.1). For the first and second moment, the agreement is again excellent,
showing that the correction induced by the active driving noise is accurately
captured to leading order. The insets of Figs. 6a and 6b show the respective
moments of the FPT for various $\nu$, indicating the systematic decrease for
increased values of $\nu$.

Figure 7: Numerical validation of first order correction to FPT moment
generating function $M^{1}$ of ATBM (cf. Eq. (108) and Sec. III.2 for
discussion) for various values of $0.1\leq\nu\leq 0.8$. The moment generating
functions were sampled for $x_{0}=0,x_{1}=\pi$, $D_{x}=1$, $\alpha=1$,
$D_{y}=1$, $\beta=\frac{1}{2}$ and $\varepsilon$ suitably chosen to fix
$\nu=D_{y}r^{2}\beta\varepsilon^{2}/(D_{x}^{2})$ (plot marks). For small
values of $\nu$ agreement with theoretical first-order correction (black line)
is very good. For larger values of $\nu$ the deviation increases. The rescaled
deviations, $\tilde{M}_{2}$ (see Eq. (134)), (inset) collapse and thus confirm
that these deviations are systematic higher-order corrections. See (III.2.2)
for further results and discussion.

###  III.3 Limit Cases

The framework we introduce here allows to study coloured driving noises at any
correlation time $\beta^{-1}$. In particular, this includes two limit cases of
$\beta\to 0$ and $\beta\to\infty$. For appropriate re-scaling of $D_{y}$, the
former limit corresponds to a particular quenched disorder model, and the
latter to additional white noise. In what follows we discuss these limit cases
in more detail.

####  III.3.1 The white noise limit

For very small autocorrelation times $\beta^{-1}$ the driving noise $y_{t}$
appears more and more as white noise. The correlator of $y_{t}$ tends towards

$\displaystyle\overline{y_{t_{0}}y_{t_{1}}}=D_{y}\beta
e^{-\beta|t_{1}-t_{0}|}\to 2D_{y}\delta(t_{1}-t_{0}),\qquad\beta\to\infty.$
(135)

In this limit, the driving noise features in the Langevin equation (1) as
additional white noise and is absorbed as

$\displaystyle\dot{x}=-V^{\prime}(x_{t})+\sqrt{2\left(D_{x}+\varepsilon^{2}D_{y}\right)}\xi,$
(136)

such that effectively the diffusion constant is shifted by $D_{x}\mapsto
D_{x}+\varepsilon^{2}D_{y}$. In the white noise limit, the theory is Markovian
and Eq. (11) may be applied using the shifted diffusivity to obtain exact
results.The perturbation theory presented here then correspond to an expansion
of the Markovian result in a perturbation to the diffusion constant $D_{x}$;
and as a result the following identity must hold:

$\displaystyle\lim_{\beta\to\infty}\beta M^{1}_{x_{0},x_{1}}(s)=\alpha
D_{x}\frac{\partial}{\partial D_{x}}M^{0}_{x_{0},x_{1}}(s).$ (137)

One can verify that this relation is indeed satisfied for the Brownian motion
case (Eqs. (128) and (III.2.1)) as well as the Ornstein-Uhlenbeck case (Eqs.
(99) and (100)).

####  III.3.2 Quenched Disorder limit

In the opposite limit of $\beta\to 0$, provided $D_{y}\beta=w^{2}$ remains
fixed, the driving noise “freezes” to a random constant since

$\displaystyle\overline{y_{t_{0}}y_{t_{1}}}=D_{y}\beta
e^{-\beta|t_{1}-t_{0}|}\to w^{2},\qquad\beta\to 0.$ (138)

Effectively, the Langevin equation (1) therefore turns into

$\displaystyle\dot{x}=-V^{\prime}(x_{t})+\xi_{t}+v$ (139)

where $v$ is a constant driving velocity which is normal distributed according
to $v\sim\mathcal{N}(0,\varepsilon^{2}w^{2})$. The driving noise average
$\overline{e^{-s\tau_{x_{0},x_{1}}}}$ then corresponds to a quenched average
over the ensemble of normal distributed velocities $v$. If we treat $M^{1}$ as
a functional of the potential $V(x)$ in which the particle is embedded, then
formally

$\displaystyle\frac{w^{2}\varepsilon^{2}}{D_{x}\alpha}\lim_{\beta\to
0}M^{1}_{x_{0},x_{1}}(s;[V(x)])$ (140)
$\displaystyle=\frac{1}{2}\underbrace{\varepsilon^{2}w^{2}}_{=\overline{v^{2}}}\left.\partial_{v}^{2}\right|_{v=0}M_{x_{0},x_{1}}(s;[V(x)+vx]).$
(141)

Our framework therefore predicts the first-order correction in
$\overline{v^{2}}$.

##### Example: Brownian Motion with periodic boundary conditions and a random
drift

For Brownian Motion with periodic boundary conditions, as studied in Sec.
III.2, one can compute the moment generating function of first-passage times
for a particular fixed drift $v$ exactly (see Appendix B and in particular Eq.
(159) for the result. We could not find this result elsewhere in the
literature.) On expanding this result in orders of the drift $v$ and averaging
$v^{2}$ over its distribution $\mathcal{N}(0,\varepsilon^{2}w^{2})$ we obtain
a resulting quenched average approximation in orders of $\overline{v^{2}}$,
${\overline{M\left(s;v\right)}}={{M\left(s;v=0\right)}}+\frac{1}{2}{{\left.\partial_{v}^{2}\right|_{v=0}M\left(s;v\right)}}\cdot\overline{v^{2}}+...$.
When employing our framework and letting $\beta\to 0$ in our general result
Eq. (127) we recover precisely
$\frac{1}{2}{{\partial_{v}^{2}|_{v=0}\hat{p}\left(s;v\right)}}$. The necessary
calculations are given in App. C and show that this is indeed the case. By way
of this relation, our framework for instance returns the correction to the
mean-first passage time of a Brownian motion with quenched disordered drift to
first order in $\nu$ as (compare to Eq. (131))

$\displaystyle\alpha\left\langle{\overline{\tau_{0,\theta}}}\right\rangle=\frac{\theta\left(2\pi-\theta\right)}{2}-\frac{\varepsilon^{2}w^{2}}{\alpha^{2}r^{2}}\cdot\frac{\theta^{2}(\theta-2\pi)^{2}}{24}+...$
(142)

such that quenched disorder lowers the mean first-passage time for any choice
of parameters. Further, for the mean squared first-passage time we obtain

$\displaystyle\alpha^{2}\left\langle{\overline{\tau_{0,\theta}^{2}}}\right\rangle=\frac{1}{12}\theta\left(\theta-2\pi\right)\left(\theta-2\varphi\pi\right)\left(\theta-2\pi+2\varphi\pi\right)$
$\displaystyle+\frac{\varepsilon^{2}w^{2}}{\alpha^{2}r^{2}}\frac{\theta(\theta-2\pi)(\theta-\pi(1+\psi^{+}))(\theta-\pi(1-\psi^{+}))}{120}$
$\displaystyle\qquad\cdot(\theta-\pi(1+\psi^{-}))(\theta-\pi(1-\psi^{-})),$
(143)

where $\varphi=\frac{1+\sqrt{5}}{2}$ is the golden ratio, and
$\psi^{\pm}=\frac{1\pm\sqrt{10}}{\sqrt{3}}$.

##  IV  Conclusion

In the present work, we introduce a perturbative approach to study the first-
passage time distribution of stochastic processes which are driven both by
white and coloured noise. This class of stochastic processes lies at the heart
of the study of active particles in a thermal environment. The activeness,
which can either represent self-propulsion or hidden degrees of freedom, is
modelled by coloured noise $y_{t}$ with characteristic timescale $\beta^{-1}$
while the thermal bath is modelled by white noise $\xi_{t}$ with diffusion
constant $D_{x}$. The expansion parameter in which the perturbation takes
place is a dimensionless quantity, $\nu$, which indicates how strong the
fluctuations of the active noise are in comparison to the strength of thermal
fluctuations. We study the regime in which $\nu$ is small.

Setting out from a renewal equation which gives the moment generating function
of first-passage times, we employ a functional expansion to obtain its
perturbative corrections. This key equation (17) stands at the centre of this
work. In order to solve it perturbatively, one needs to calculate the
expansion terms (cf. Eqs. (31) and (33)) which involve the eigenfunctions of
the Fokker Planck operator associated to the non driven process (cf. Eq.
(43)). To first order in $\nu$, we obtain an analytic result of the moment-
generating function in terms of the associated eigenfunctions (cf. Eq.
(II.4)). Higher order contributions can be obtained by further iterating the
steps outlined in Sec. II.2.

To illustrate the capabilities of our framework, we study two systems. First,
we consider an active thermal particle in a harmonic potential, the Active
Thermal Ornstein-Uhlenbeck Process. In Sec. III.1.1, we calculate all
necessary response functions to find the first order correction to the moment-
generating function of first-passage times (cf. Eq. (100)). By taking
derivatives, we could in principle obtain closed form expressions for the
first order correction to any moment of the first-passage time distribution.
We confirm these analytical results by numerical simulations. Sampling the
experimental moment generating function, we obtain an excellent agreement with
the first-order correction (see Fig. 4). For larger values of $\nu$, the
perturbative parameter, the deviations systematically indicate higher-order
corrections. Further, we compare the theoretically predicted correction to the
first two moments of the first-passage times to numerical results (see Figs.
3a and 3b) which are in excellent agreement. Secondly, we study Active Thermal
Brownian Motion on a ring (see Sec. III.2). Again, we illustrate our framework
by finding the first-order correction to the moment generating function (cf.
Eq. (127)). Numerical simulations show excellent agreement and systematic
higher-order corrections (see Fig. 7). Both first and second moment of the
first-passage time are obtained from Eq. (127) and show good agreement with
numerical simulations.

Further, since the perturbation theory we present makes no assumption on
$\beta^{-1}$, we are able to recover the limiting cases for $\beta\to 0$ and
$\infty$, respectively. The case of $\beta\to 0$ is of particular interest
since it recovers quenched disorder averages over processes with additional
fixed and normal distributed drift (see Sec. III.3).

The framework requires to find the eigenfunctions of a differential operator,
and to express all transition and return densities as sums over these
eigenfunctions. This often requires certain calculations that for more unusual
eigenfunctions may be difficult to perform.

Our approach further allows for the presence of an external potential provided
the associated differential operator (Eq. (41)) can be diagonalised. This
significantly extends the range of systems our framework can be applied to. In
this work, we focused on Fourier-modes and Hermite-polynomials which are
suitable for flat and harmonic potentials. It is, however, also possible to
study piece-wise combinations of the potentials using these eigenfunctions.
This may be relevant when studying bi-stable processes for instance. Further,
as long as Eq. (41) can be diagonalised, this framework also allows for a
space-dependent thermal diffusivity by letting $D_{x}=D(x)$. For future work,
for instance, it would be interesting to study first-passage time behaviour of
particles at the boundary between two heat baths at different temperature
(e.g. $D(x)=D_{0}+\operatorname{sgn}(x)\Delta D$).

Moreover, the functional expansion in $\hat{y}(\omega)$ (cf. Eq. (28))
drastically simplifies in the case of $y_{t}$ being a periodic driving force.
This framework therefore would not only be able to capture stochastic $y$, but
also oscillating deterministic driving forces. This will possibly be addressed
in future work.

To first order in $\nu$, the corrections, as given in Eqs. (II.4) to (74),
involve simple complex integrals which can be solved by the residue theorem.
For higher-order corrections, however, the integration runs over more than one
free internal variable and will require more work. This corresponds to the
problems of typical Feynman-diagrams of higher order in statistical field
theories which often involve non-trivial integrals. To study higher orders,
field theory therefore would provide the necessary toolbox to solve the
required correction terms. Already, the general result obtained in Eq. (II.4)
can be analogously derived using field-theoretic methods (see [52]).

In conclusion, by casting a broad class of non-Markovian processes into a
perturbative language, our framework will prove itself useful in tackling a
diverse range of future challenges.

###### Acknowledgements.

BW and GP would like to thank the Francis Crick Institute for hospitality and
Andy Thomas for invaluable computing support. BW thanks Kay Wiese, Grigorios
Pavliotis, Gregory Schehr and Raphael Voituriez for insightful discussions. GS
was supported by the Francis Crick Institute which receives its core funding
from Cancer Research UK (FC001317), the UK Medical Research Council
(FC001317), and the Wellcome Trust (FC001317).

## References

  * Kampen [2007] N. G. v. Kampen, Stochastic processes in physics and chemistry, 3rd ed. (Elsevier, 2007).
  * Redner [2001] S. Redner, A guide to first-passage processes (Cambridge University Press, 2001).
  * Metzler et al. [2014] R. Metzler, G. Oshanin, and S. Redner, First-Passage Phenomena and Their Applications (World Scientific, 2014).
  * Szabo et al. [1980] A. Szabo, K. Schulten, and Z. Schulten, First passage time approach to diffusion controlled reactions, J. Chem. Phys. 72, 4350 (1980).
  * Sokolov [2003] I. M. Sokolov, Cyclization of a Polymer: First-Passage Problem for a Non-Markovian Process, Phys. Rev. Lett. 90, 080601 (2003).
  * Ghusinga et al. [2017] K. R. Ghusinga, J. J. Dennehy, and A. Singh, First-passage time approach to controlling noise in the timing of intracellular events, Proc. Natl. Acad. Sci. 114, 693 (2017).
  * Taillefumier and Magnasco [2013] T. Taillefumier and M. O. Magnasco, A phase transition in the first passage of a Brownian process through a fluctuating boundary with implications for neural coding, Proc. Natl. Acad. Sci. 110, E1438 (2013).
  * Chicheportiche and Bouchaud [2014] R. Chicheportiche and J.-P. Bouchaud, Some Applications of First-Passage Ideas to Finance, in First-Passage Phenomena and Their Applications, edited by R. Metzler, G. Oshanin, and S. Redner (World Scientific, 2014) pp. 447–476.
  * Tejedor et al. [2009] V. Tejedor, O. Bénichou, and R. Voituriez, Global mean first-passage times of random walks on complex networks, Phys. Rev. E 80, 065104 (2009).
  * Hartich and Godec [2019] D. Hartich and A. Godec, Extreme value statistics of ergodic Markov processes from first passage times in the large deviation limit, J. Phys. A 52, 244001 (2019).
  * Bray et al. [2013] A. J. Bray, S. N. Majumdar, and G. Schehr, Persistence and first-passage properties in nonequilibrium systems, Adv. Phys. 62, 225 (2013).
  * Schrödinger [1915] E. Schrödinger, Zur Theorie der Fall- und Steigversuche an Teilchen mit Brownscher Bewegung, Z. Phys. 16, 289 (1915).
  * Pontryagin et al. [1989] L. Pontryagin, A. Andronov, and A. Vitt, Appendix: On the statistical treatment of dynamical systems, in Noise in Nonlinear Dynamical Systems, Vol. 1, edited by F. Moss and P. V. E. McClintock (Cambridge University Press, 1989) pp. 329–348.
  * Siegert [1951] A. J. F. Siegert, On the First Passage Time Probability Problem, Phys. Rev. 81, 617 (1951).
  * Darling and Siegert [1953] D. A. Darling and A. J. F. Siegert, The First Passage Problem for a Continuous Markov Process, Ann. Math. Stat. 24, 624 (1953).
  * Masoliver et al. [1986] J. Masoliver, K. Lindenberg, and B. J. West, First-passage times for non-Markovian processes: Correlated impacts on bound processes, Phys. Rev. A 34, 2351 (1986).
  * Kuś et al. [1990] M. Kuś, E. Wajnryb, and K. Wódkiewicz, Mean first-passage time in the presence of colored noise: An exact theory, Phys. Rev. A 42, 7500 (1990).
  * Kuś et al. [1991] M. Kuś, E. Wajnryb, and K. Wódkiewicz, Mean first-passage time in the presence of colored noise: A random-telegraph-signal approach, Phys. Rev. A 43, 4167 (1991).
  * Ramrez-Piscina and Sancho [1991] L. Ramrez-Piscina and J. M. Sancho, First-passage times for a marginal state driven by colored noise, Phys. Rev. A 43, 663 (1991).
  * Dykman et al. [1993] M. I. Dykman, R. Mannella, P. V. E. McClintock, N. D. Stein, and N. G. Stocks, Probability distributions and escape rates for systems driven by quasimonochromatic noise, Phys. Rev. E 47, 3996 (1993).
  * Jung [1997] P. Jung, Colored noise in dynamical systems: Some exact solutions, in Stochastic Dynamics, Lecture Notes in Physics, Vol. 484, edited by L. Schimansky-Geier and T. Pöschel (Springer Berlin Heidelberg, 1997) pp. 23–31.
  * Fox [1988] R. F. Fox, Mean first-passage times and colored noise, Phys. Rev. A 37, 911 (1988).
  * Hanggi and Jung [1995] P. Hanggi and P. Jung, Colored Noise in Dynamical Systems, Adv. Chem. Phys. 89, 239 (1995).
  * Guérin et al. [2016] T. Guérin, N. Levernier, O. Bénichou, and R. Voituriez, Mean first-passage times of non-Markovian random walkers in confinement, Nature 534, 356 (2016).
  * Godec and Metzler [2016a] A. Godec and R. Metzler, First passage time distribution in heterogeneity controlled kinetics: going beyond the mean first passage time, Sci. Rep. 6, 20349 (2016a).
  * Godec and Metzler [2016b] A. Godec and R. Metzler, Universal Proximity Effect in Target Search Kinetics in the Few-Encounter Limit, Phys. Rev. X 6, 041037 (2016b).
  * Callen and Welton [1951] H. B. Callen and T. A. Welton, Irreversibility and generalized noise, Phys. Rev. 83, 34 (1951).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Rep. Prog. Phys. 29, 255 (1966).
  * Maggi et al. [2015] C. Maggi, U. M. B. Marconi, N. Gnan, and R. Di Leonardo, Multidimensional stationary probability distribution for interacting active particles, Sci. Rep. 5, 10742 (2015).
  * Fodor et al. [2016] É. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco, and F. van Wijland, How far from equilibrium is active matter?, Phys. Rev. Lett. 117, 038103 (2016).
  * Mandal et al. [2017] D. Mandal, K. Klymko, and M. R. DeWeese, Entropy production and fluctuation theorems for active matter, Phys. Rev. Lett. 119, 258001 (2017).
  * Szamel [2019] G. Szamel, Stochastic thermodynamics for self-propelled particles, Phys. Rev. E 100, 050603 (2019).
  * Łuczka [2005] J. Łuczka, Non-Markovian stochastic processes: Colored noise, Chaos 15, 026107 (2005).
  * Malakar et al. [2018] K. Malakar, V. Jemseena, A. Kundu, K. Vijay Kumar, S. Sabhapandit, S. N. Majumdar, S. Redner, and A. Dhar, Steady state, relaxation and first-passage properties of a run-and-tumble particle in one-dimension, J. Stat. Mech. Theory Exp. 2018, 043215 (2018).
  * Caprini et al. [2019] L. Caprini, U. M. B. Marconi, A. Puglisi, and A. Vulpiani, The entropy production of Ornstein–Uhlenbeck active particles: a path integral method for correlations, J. Stat. Mech. Theory Exp. 2019, 053203 (2019).
  * Sevilla et al. [2019] F. J. Sevilla, R. F. Rodríguez, and J. R. Gomez-Solano, Generalized ornstein-uhlenbeck model for active motion, Phys. Rev. E 100, 032123 (2019).
  * Dabelow et al. [2019] L. Dabelow, S. Bo, and R. Eichhorn, Irreversibility in active matter systems: Fluctuation theorem and mutual information, Phys. Rev. X 9, 021009 (2019).
  * Cates [2012] M. E. Cates, Diffusive transport without detailed balance in motile bacteria: does microbiology need statistical physics?, Rep. Prog. Phys. 75, 042601 (2012).
  * Dhar et al. [2019] A. Dhar, A. Kundu, S. N. Majumdar, S. Sabhapandit, and G. Schehr, Run-and-tumble particle in one-dimensional confining potentials: Steady-state, relaxation, and first-passage properties, Phys. Rev. E 99, 032132 (2019).
  * Lindner [2009] B. Lindner, Some unsolved problems relating to noise in biological systems, J. Stat. Mech. Theory Exp. , P01008 (2009).
  * Fodor et al. [2015] É. Fodor, M. Guo, N. S. Gov, P. Visco, D. A. Weitz, and F. van Wijland, Activity-driven fluctuations in living cells, EPL 110, 48005 (2015).
  * Barbara J Gabryś and Bartkiewicz [2004] K. P. Barbara J Gabryś and S. J. Bartkiewicz, Brownian motion, molecular motors and ratchets, Physica A 336, 112 (2004).
  * Pavliotis [2005] G. A. Pavliotis, A multiscale approach to Brownian motors, Phys. Lett. A 344, 331 (2005).
  * Hörmander [1963] L. Hörmander, Linear Partial Differential Operators, Grundlehren der mathematischen Wissenschaften, Vol. 116 (Springer, 1963).
  * Bakry et al. [2013] D. Bakry, I. Gentil, and M. Ledoux, Analysis and geometry of Markov diffusion operators, Grundlehren der mathematischen Wissenschaften, Vol. 348 (Springer, 2013).
  * Pavliotis [2014] G. A. Pavliotis, Stochastic Processes and Applications, Texts in Applied Mathematics, Vol. 60 (Springer, New York, 2014).
  * Risken [1984] H. Risken, The Fokker-Planck Equation (Springer, 1984).
  * Magnus et al. [1966] W. Magnus, F. Oberhettinger, and R. P. Soni, Formulas and Theorems for the Special Functions of Mathematical Physics, Grundlehren der mathematischen Wissenschaften, Vol. 52 (Springer, 1966).
  * Gradshteyn and Ryzhik [2007] I. S. Gradshteyn and I. Ryzhik, Table of Integrals, Series, and Products, 7th ed. (Academic Press, 2007).
  * [50] Wolfram Research, Mathematica, Version 12.1.
  * Valkó and Abate [2004] P. Valkó and J. Abate, Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion, Comput. Math. with Appl. 48, 626 (2004).
  * Walter [2020] B. Walter, Field theories for stochastic processes, Ph.D. thesis, Imperial College London (2020).
  * Cox and Miller [1965] D. R. Cox and H. D. Miller, The theory of stochastic processes, 1st ed. (Chapman and Hall, 1965).

##  Appendix A Notations

For reference, we list the most commonly used notations in Table 1.

Symbol | Description
---|---
$p_{x_{0},x_{1}}(t_{0},t_{1})$ | First-passage time probability density at $(x_{1},t_{1})$ starting from $(x_{0},t_{0})$
$\hat{p}_{x_{0},x_{1}}(\omega)$ | Fourier transform (characteristic function) of FPT density
$M_{x_{0},x_{1}}(s)$ | Laplace transform (moment generating function) of FPT density
$T_{x_{0},x_{1}}(t_{0},t_{1})$ | Probability density of progress from $(x_{0},t_{0})$ to $(x_{1},t_{1})$
$\hat{T}_{x_{0},x_{1}}(\omega)$ | Fourier transform of transition probability (cf. Eq. (9))
$R_{x_{1}}(t_{0},t_{1})$ | Probability density of return at $x_{1}$ at time $t_{1}$ starting at $t_{0}$
$\hat{R}_{x_{1}}(\omega)$ | Fourier transform of return probability (cf. Eq. (9))
$\hat{T}^{(n)}_{x_{0},x_{1}}(\omega_{1},\tilde{\omega}_{1},...,\tilde{\omega}_{n})$ | Expansion functionals of transition probability in $\varepsilon y_{t}$ (cf. Eq. (29))
$\hat{R}^{(n)}_{x_{1}}(\omega_{1},\tilde{\omega}_{1},...,\tilde{\omega}_{n})$ | Expansion functionals of return probability in $\varepsilon y_{t}$ (cf. Eq. (30))
$M_{x_{0},x_{1}}^{0},M_{x_{0},x_{1}}^{1},...$ |  $y_{t}$-Averaged coefficients to moment generating function $M_{x_{0},x_{1}}(s)$ in $\nu$-expansion
Table 1: Overview of symbols used for various probability densities and their
functional transforms. When possible, the lower subscripts $x_{0}$, $x_{1}$
are ommitted.

##  Appendix B First passage times of Brownian Motion on a ring with drift

We calculate the first-passage time distribution of a Brownian Motion with
drift $v$ on a ring of radius $r$ departing from the methods outlined in [53,
Chp. V] although the explicit formula is not given there. Instead of
considering the first passage event of transition $x_{0}\to x_{1}$, we
calculate the exit probability of a Brownian Motion on the real line over the
absorbing boundaries at $x_{1}$ and $x_{1}-2\pi r$ where without loss of
generality we chose $x_{1}-2\pi r<x_{0}<x_{1}$. The transition density
$T\left(t\right)\equiv T_{x_{0},x_{1}}(t)$ satisfies the Fokker-Planck
Equation

$\displaystyle\partial_{t}T\left(t\right)=D_{x}\partial_{x_{1}}^{2}T\left(t\right)-v\partial_{x_{1}}T\left(t\right)$
(144)

and the Kolmogorov backward equation

$\displaystyle\partial_{t}T\left(t\right)=D_{x}\partial_{x_{0}}^{2}T\left(t\right)+v\partial_{x_{0}}T\left(t\right)$
(145)

with absorbing boundary conditions for $t\geq 0$

$\displaystyle T_{x_{0},x_{1}}(t)$ $\displaystyle=0$ (146) $\displaystyle
T_{x_{0},x_{1}-2\pi r}(t)$ $\displaystyle=0$ (147)

and initial condition

$\displaystyle T_{x_{0},x}(t=0)=\delta(x-x_{0})$ (148)

The net flux of the surviving probability over the boundary is related to the
first-passage time distribution via

$\displaystyle p_{x_{0},x_{1}}(t)=-\partial_{t}\int_{x_{1}-2\pi
r}^{x}\mathrm{d}{x}\,T_{x_{0},x}(t).$ (149)

Using the backward equation (145), and applying a further derivative in time
one finds

$\displaystyle\partial_{t}p_{x_{0},x_{1}}(t)$
$\displaystyle=-(D_{x}\partial_{x_{0}}^{2}+v\partial_{x_{0}})\partial_{t}\int_{x_{1}-2\pi
r}^{x_{1}}\mathrm{d}{x}\,T_{x_{0},x}(t)$ (150)
$\displaystyle=(D_{x}\partial_{x_{0}}^{2}+v\partial_{x_{0}})p_{x_{0},x_{1}}(t).$
(151)

At the boundary, the first-passage time distribution is given by

$\displaystyle p_{x_{1},x_{1}}(t)=p_{x_{1}-2\pi r,x_{1}}(t)=\delta(0^{+})$
(152)

since the particle is immediately absorbed, while at time $t=0$ for $x_{0}$
away from the boundary

$\displaystyle p_{x_{0},x_{1}}(t=0)=0$ (153)

since the particle has not yet begun to diffuse. Under Laplace transform in
$t$, the differential equation (151), together with these boundary and initial
conditions returns

$\displaystyle\begin{cases}sp_{x_{0},x_{1}}(s)=(D_{x}\partial_{x_{0}}^{2}+v\partial_{x_{0}})p_{x_{0},x_{1}}(s)&x_{1}-2\pi
r<x_{0}<x_{1}\\\ p_{x_{0},x_{1}}(s)=p_{x_{0},x_{1}-2\pi
r}(s)=1&x_{0}=x_{1}-2\pi r\text{ or }x_{0}=x_{1}.\end{cases}$ (154)

The ordinary differential equation (away from the boundary) is solved by the
exponential ansatz

$\displaystyle M_{x_{0},x_{1}}(s)=Ae^{\omega_{1}x_{0}}+Be^{\omega_{2}x_{0}}$
(155)

Inserting this ansatz into (154) enforces

$\displaystyle\omega_{1,2}=-\frac{v}{2D_{x}}\pm\frac{\sqrt{v^{2}+4D_{x}s}}{2D_{x}}$
(156)

The boundary conditions (154) fix the normalising constants $A,B$ to

$\displaystyle A$ $\displaystyle=-e^{\pi
r\omega_{1}-\omega_{1}x_{1}}\frac{\sinh(\pi r\omega_{2})}{\sinh(\pi
r(\omega_{1}-\omega_{2}))}$ (157) $\displaystyle B$ $\displaystyle=\sinh(\pi
r\omega_{1})\frac{e^{\pi r\omega_{2}-\omega_{2}x_{1}}}{\sinh(\pi
r(\omega_{1}-\omega_{2}))}$ (158)

After some further simplifications one arrives at the $v$-dependent moment
generating function

$\displaystyle
M_{x_{0},x_{1}}(s;v)=\frac{e^{\frac{rv}{2D_{x}}\left(\theta-2\pi\right)}\sinh\left(\frac{r\sqrt{v^{2}+4D_{x}s}}{2D_{x}}\theta\right)-e^{\frac{rv}{2D_{x}}\theta}\sinh\left(\frac{r\sqrt{v^{2}+4D_{x}s}}{2D_{x}}(\theta-2\pi)\right)}{\sinh\left(\pi\frac{r\sqrt{v^{2}+4D_{x}s}}{D_{x}}\right)}$
(159)

with $\theta=(x_{1}-x_{0})/r$ following the notation from the main text. We
note that for $v\to 0$,

$\displaystyle\frac{r\sqrt{v^{2}+4D_{x}s}}{2D_{x}}\to\sqrt{\bar{s}}$ (160)

and one recovers the undriven moment generating function

$\displaystyle
M_{x_{0},x_{1}}(s;v=0)=M_{x_{0},x_{1}}(s)=\frac{\sinh\left(\theta\sqrt{\bar{s}}\right)-\sinh\left((\theta-2\pi)\sqrt{\bar{s}}\right)}{\sinh\left(2\pi\sqrt{\bar{s}}\right)}=\frac{\cosh\left((\theta-\pi)\sqrt{\bar{s}}\right)}{\cosh\left(\pi\sqrt{\bar{s}}\right)}$
(161)

in agreement with the independently found expression (122). In App. C, we show
that to second order in $v$ this result is identical to the first-order
correction $M^{1}$ from Eq. (127) in the limit of $\beta\to 0$.

##  Appendix C Equivalence of quenched averages

In this section, we provide a more detailed proof showing that Eq. (141)
indeed holds for the case of Brownian Motion driven by coloured noise, i.e.
that our framework perturbatively gives the correct moment generating function
of first-passage times when taking the quenched average over Eq. (139) with
normal distributed drift $v\sim\mathcal{N}(0,\varepsilon^{2}w^{2})$. To that
end, we take the $\beta\to 0$ limit of the analytically found $M^{1}$ (cf. Eq.
(III.2.1)),

$\displaystyle\lim_{\beta\to 0}M^{1}=$
$\displaystyle\frac{1}{8\sqrt{\bar{s}}}\Big{\\{}\theta\cosh\left(\theta\sqrt{\bar{s}}\right)\left(\theta\sqrt{\bar{s}}-\tanh\left(\pi\sqrt{\bar{s}}\right)\right)$
(162)
$\displaystyle+\sinh\left(\theta\sqrt{\bar{s}}\right)\left[-\left(\theta^{2}-2\pi\theta+2\pi^{2}\right)\sqrt{\bar{s}}\tanh\left(\pi\sqrt{\bar{s}}\right)+2\pi(\pi-\theta)\sqrt{\bar{s}}\coth\left(\pi\sqrt{\bar{s}}\right)+\pi\tanh^{2}\left(\pi\sqrt{\bar{s}}\right)+\theta-\pi\right]\Big{\\}}.$

In Eq. (141), it is claimed that this equals

$\displaystyle\frac{D_{x}\alpha}{2\varepsilon^{2}w^{2}}\int_{-\infty}^{\infty}\mathrm{d}{v}\,\frac{e^{-\frac{v^{2}}{2\varepsilon^{2}w^{2}}}}{\sqrt{2\pi\varepsilon^{2}w^{2}}}v^{2}\left.\partial_{v}^{2}\right|_{v=0}M_{x_{0},x_{1}}(s;v)\
.$ (163)

Evaluating this expression using the result from Eq. (159) and setting $v=0$
results in

$\displaystyle\frac{D_{x}\alpha}{2}\left.\partial_{v}^{2}\right|_{v=0}M_{x_{0},x_{1}}(s;v)=$
$\displaystyle\frac{1}{8\sqrt{\bar{s}}}\Big{\\{}\frac{1}{\text{sinh}\left(2\pi\sqrt{\bar{s}}\right)}\Big{[}\sinh\left((\theta-2\pi)\sqrt{\bar{s}}\right)\left(2\pi\coth\left(2\pi\sqrt{\bar{s}}\right)-\theta^{2}\sqrt{\bar{s}}\right)+\theta\cosh\left(\theta\sqrt{\bar{s}}\right)$
$\displaystyle+(2\pi-\theta)\cosh\left((\theta-2\pi)\sqrt{\bar{s}}\right)+\sinh\left(\theta\sqrt{\bar{s}}\right)\left((\theta-2\pi)^{2}\sqrt{\bar{s}}-2\pi\coth\left(2\pi\sqrt{\bar{s}}\right)\right)\Big{]}\Big{\\}}\
.$

The expressions in (162) and (C) are indeed equal as can be shown with the
help of Mathematica [50].

##  Appendix D Explicit expressions for functional derivatives of transition
probability densities of Ornstein-Uhlenbeck processes

Following the notation from Sec. III.1.1, we here give the explicit
expressions of $\hat{T}^{(1)},\hat{T}^{(2)}$ as implicitly given in Eqs.
(95)-(97). We confine ourselves to the case of $\bar{x}_{0}<\bar{x}_{1}$.
Starting from the Fourier-transformed transition probability density (cf. Eq.
(93)), all other functional derivatives are given as partial derivatives of
this density. From formula (95), one obtains

$\displaystyle\hat{T}^{(1)}\left(-i\bar{s}\alpha,-i\beta\right)=-\frac{e^{\frac{{\bar{x}_{0}}^{2}-{\bar{x}_{1}}^{2}}{4}}}{\sqrt{2\pi}\ell^{2}\alpha^{2}(1-\bar{\beta})}\left[\Gamma(\bar{\beta}+\bar{s})D_{-\bar{\beta}-\bar{s}}(-\bar{x}_{0})D_{-\bar{\beta}-\bar{s}+1}(\bar{x}_{1})-\Gamma(\bar{s}+1)D_{-\bar{s}-1}(-\bar{x}_{0})D_{-\bar{s}}(\bar{x}_{1})\right]$
(164)

where $\Gamma(\bar{s})$ is the usual Gamma-function. Letting
$\bar{x}_{0}\to\bar{x}_{1}$, gives the first functional derivative of the
return probability at $\bar{x}_{1}$,

$\displaystyle\hat{R}^{(1)}\left(-i\bar{s}\alpha,-i\beta\right)=-\frac{\Gamma(\bar{\beta}+\bar{s})D_{-\bar{\beta}-\bar{s}}\left(-\bar{x}_{1}\right)D_{-\bar{\beta}-\bar{s}+1}\left(\bar{x}_{1}\right)-\Gamma(\bar{s}+1)D_{-\bar{s}-1}\left(-\bar{x}_{1}\right)D_{-\bar{s}}\left(\bar{x}_{1}\right)}{\sqrt{2\pi}\ell^{2}\alpha^{2}(1-\bar{\beta})}$
(165)

These results further imply

$\displaystyle\hat{T}^{(1)}\left(-i\bar{s}\alpha-i\beta,i\beta\right)=-\frac{e^{\frac{\bar{x}_{0}^{2}-\bar{x}_{1}^{2}}{4}}(\Gamma(\bar{s})D_{-\bar{s}}(-\bar{x}_{0})D_{1-\bar{s}}(\bar{x}_{1})-\Gamma(\bar{\beta}+\bar{s}+1)D_{-\bar{\beta}-\bar{s}-1}(-\bar{x}_{0})D_{-\bar{\beta}-\bar{s}}(\bar{x}_{1}))}{\sqrt{2\pi}(1+\bar{\beta})\ell^{2}\alpha^{2}}$
(166)

and

$\displaystyle\hat{R}^{(1)}\left(-i\bar{s}\alpha-i\beta,i\beta\right)=-\frac{\Gamma(\bar{s})D_{-\bar{s}}\left(-\bar{x}_{1}\right)D_{1-\bar{s}}\left(\bar{x}_{1}\right)-\Gamma(\bar{\beta}+\bar{s}+1)D_{-\bar{\beta}-\bar{s}-1}\left(-\bar{x}_{1}\right)D_{-\bar{\beta}-\bar{s}}\left(\bar{x}_{1}\right)}{\sqrt{2\pi}(1+\bar{\beta})\ell^{2}\alpha^{2}}$
(167)

The second order derivative of the transition probability is

$\displaystyle\hat{T}^{(2)}\left(-i\bar{s}\alpha,-i\beta,i\beta\right)=\frac{e^{\frac{\bar{x}_{0}^{2}-\bar{x}_{1}^{2}}{4}}}{\sqrt{2\pi}\left({\bar{\beta}}^{2}-1\right)\ell^{3}\alpha^{3}}$
$\displaystyle\left[2\Gamma(\bar{\beta}+\bar{s}+1)D_{-\bar{\beta}-\bar{s}-1}(-\bar{x}_{0})D_{-\bar{\beta}-\bar{s}+1}(\bar{x}_{1})\right.$
$\displaystyle\left.+(\bar{\beta}-1)\Gamma(\bar{s})D_{-\bar{s}}(-\bar{x}_{0})D_{2-\bar{s}}(\bar{x}_{1})-(\bar{\beta}+1)\Gamma(\bar{s}+2)D_{-\bar{s}-2}(-\bar{x}_{0})D_{-\bar{s}}(\bar{x}_{1}))\right]$
(168)

and of the return probability

$\displaystyle\hat{R}^{(2)}\left(-i\bar{s}\alpha,-i\beta,i\beta\right)=\frac{1}{\sqrt{2\pi}\left({\bar{\beta}}^{2}-1\right)\ell^{3}\alpha^{3}}$
$\displaystyle\left[2\Gamma(\bar{\beta}+\bar{s}+1)D_{-\bar{\beta}-\bar{s}-1}(-\bar{x}_{1})D_{-\bar{\beta}-\bar{s}+1}(\bar{x}_{1})\right.$
$\displaystyle\left.+(\bar{\beta}-1)\Gamma(\bar{s})D_{-\bar{s}}(-\bar{x}_{1})D_{2-\bar{s}}(\bar{x}_{1})-(\bar{\beta}+1)\Gamma(\bar{s}+2)D_{-\bar{s}-2}(-\bar{x}_{1})D_{-\bar{s}}(\bar{x}_{1})\right]$
(169)