First passage time distribution of active thermal particles in potentials

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.

Understanding the statistical properties of firstpassage times (FPTs), the time a stochastic process takes to first reach a prescribed target, has fuelled research in stochastic dynamics for over a century. Thanks to its wide-ranging applications it has enjoyed increased attention over the last decade [1][2][3]. First-passage times are often used as 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].
Historically, the first settings in which FPT-problems were studied were Markovian processes in one spatial dimension. Schrödinger approached this problem first by integrating over the probability density with absorbing boundary conditions [12]. Pontryagin et al. introduced a method which casts the mean first-passage time (MFPT) into an ODE derived from the Kolmogorov backward equation [13,14]. Further, Siegert and Darling introduced a method to obtain the moment-generating function of the FPT starting from a renewal equation [15,16]. These three key advances of the first "classical" period suggest that, despite the innocent looking simplicity of the problem, even Markovian processes do generally not allow for a closed form expression of the full distribution of FPT. Expectations alone, the MFPT, are difficult enough to compute in very simple systems.
In a parallel development, and sparked by the work of Kramers [17], much effort was put into investigating the rate with which fluctuating particles escape from metastable potentials. The forces the surrounding heat bath exerts on the particle, usually inter-molecular collisions, were modelled as white noise ξ t with correlator Active driving noise correlation β −1 x V(x) 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 τx 0 ,x 1 ; the time such a random walker (blue rough path) takes to first reach x1 starting from x0 (dashed lines). ξ s ξ t = 2D x δ(s − t). This correlator implies that the timescale on which the particle evolves is infinitely separated from the heat bath's correlations [1]. In this setup, the fluctuation-dissipation theorem (FDT) [18,19] identifies the diffusivity D x with a fixed temperature T of the surrounding heat bath. Particles subject to white noise follow Markovian (memoryless) trajectories and are immersed in a heat bath in thermodynamical equilibrium (e.g. [20]). In many systems, the paradigm of white noise is a drastic over-simplification. Coloured noise, i.e. noise with non-uniform power spectrum, was introduced to account for correlations due to the heat bath [21,22]. Coloured noise is usually assumed to have a correlator that decays exponentially with some characteristic inverse rate arXiv:2006.00116v1 [cond-mat.stat-mech] 29 May 2020 β −1 , sometimes referred to as the "colour" of the noise. White noise corresponds to the limit of β → ∞. Escape rate problems with coloured noise were a highly active area of research in the late 1980s and early 1990s, when various approximation methods to calculate the MFPT were developed [23][24][25][26][27][28]. There contradictory predictions initially led to a considerable degree of confusion [29]. Many of these developments are reviewed in [21].
More recently, coloured noise was suggested as a model for "active" swimmers which are self-propelled and whose energy-consumption is fuelled by the environment [30]. In these systems, the fluctuation-dissipation theorem does not hold (e.g. [31]) such that no thermodynamic equilibrium can be assumed. Ensembles of non-interacting "active swimmers" have been intensively studied over the last five years in the light of their nonequilibrium features [32][33][34][35]. Notably, their first-passage time distributions have been recognised as one of their characterising properties [36].
Over the last twenty years or so, first-passage times have regained interest independently of escape problems, and new methods have been developed to extend older techniques into the realm of non-Markovian stochastic processes. The classic approaches of Schrödinger (cf. [37][38][39]), Pontryagin (cf. [40,41]) and Siegert (cf. [5]) have been successfully employed to a plethora of non-Markovian systems such as generalised Langevin and Fokker-Planck equations and coupled oscillator chains. Further approaches have been found in [42][43][44][45][46]. The vast majority of that work is concerned with MFPT, also because in many cases higher moments or the full distribution contain prohibitively complicated expressions. And yet, the MFPT has recently been criticised as insufficient or misleading in characterising the timescale of a dynamics. Therefore a more precise understanding of the full distribution poses a pressing current challenge [47,48].
In this work, 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. To our knowledge, this is the first time 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 simulation.

B. Outline
In the following, we study the first-passage time distribution of non-interacting active active active active active active active active active active active active active active active active active particles subject to thermal thermal thermal thermal thermal thermal thermal thermal thermal thermal thermal thermal thermal thermal thermal thermal thermal noise in an external potential. Particles which are driven by coloured noise alone, such as those studied in [31,49], are used as a model for self-propelled particles, but do not have a coupling to a thermal heat bath; Such systems are referred to as active and athermal [50]. More recent models contain an active force term to capture self-propulsion in addition to thermal white noise which represents the heat bath, e.g. [35,51,52]. Additionally, a conservative force may be considered stemming from an external potential. Together, for a particle moving in a 1D space, this class of thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active thermal and active models is characterised by a single degree of freedom x t satisfying a Langevin-type equation [1,52] (see Fig. 1 for illustration) Here, V (x) denotes a potential, ξ t a white noise with ξ s ξ t = 2D x δ(t − s) and D x a diffusivity fixed by the FDT. The second, stochastic term y t in Eq. (1)  Since y t is self-correlated, x t by itself is no longer Markovian; Our framework therefore provides a tool to study FPT-distributions of non-Markovian processes. On the other hand, this work stands in the context of recent efforts to approach active thermal systems with field-theoretic methods. As will be reported elsewhere, the process Eq. (1) can be mapped to a non-equilibrium field-theory from which a variety of observables can be deduced. In the present work, we focus on first-passage times, for which a simpler setup is sufficient in which we use functional perturbation theory instead of a fullyfledged field theory.
The perturbation takes place in the regime where εy t is small compared to ξ t , that means where thermal fluctuations dominate. In our perturbation theory, we will emphasise this point by controlling the amplitude of the driving noise term via ε which is chosen small. A priori no assumption is made about the noise-colour β −1 and we recover all moments all moments all moments all moments all moments all moments all moments all moments all moments all moments all moments all moments all moments all moments all moments all moments all moments at once to equal perturbative order. Our framework is based on the established renewal approach [15,16], but uses functional expansions to circumvent the problem of non-Markovianity of Eq. (1).

C. Main results
The central result of the present work concerns the moment-generating function of first-passage times of x t , Theory (a) First-passage time distribution of ATOU (cf. Eq. (75)) compared to numerical Laplace inversion of analytically obtained moment generating function (97) for various values of ν. The plot marks indicate the distribution as sampled through Monte Carlo simulations with 5 × 10 6 runs. Simulation parameters are x 0 = 0, x 1 = 1, α = 1, Dx = 1, Dy = 1, β = 1 2 while ε is tuned to fix ν = Dyβε 2 /(Dxα) to values as indicated in legend. The inset shows rescaled deviations to the undriven first-passage time distribution (cf. Eq. (105)) as compared to the first-order correction (ν = 0.1 omitted). Higher-order corrections appear for growing values of ν (cf. Fig. 4). See Sec. III A 2 for discussion.
x to values as indicated in legend. The inset shows rescaled deviations to the undriven first-passage time distribution (cf. Eq. (105)) as compared to the first-order correction (ν = 0.1 omitted). Higher-order corrections appear for growing values of ν (cf. Fig. 7). See Sec. III B 2 for discussion. where τ x0,x1 is the first-passage time of x t defined as The following argument is based on two assumptions. First, we assume the moment-generating function F (s) is known for ε = 0. The case of ε = 0 corresponds to the case of a purely "passive" particle with no self propulsion. The particle behaves Markovian in this case, and already established techniques can be applied (see Introduction). We refer to this state as in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium in equilibrium. Secondly, we assume that around this state F (s) is analytic in ε or, to be more precise, in a dimensionless parameter ν which is of order O(ε 2 ). This means that the moment generating function of first-passage times has an expansion in ν of the form where M 0 and M 1 are the coefficients of expansion of F (s) in ν around ν = 0. The equilibrium component, M 0 , can be found by classical methods such as the Darling-Siegert method which is discussed in the next subsection. 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 ν for arbitrarily coloured noise (cf. Eq. (67)) as long as the autocorrelations can be integrated suitably. Further the potentials V (x) are arbitrary, as long as an associated differential operator can be diagonalised (Sec. III and [53]), 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 y t1 y t2 = D y β exp(−β|t 1 − t 2 |)) with some diffusivity D y and correlation time β −1 . In the first case, the particle is placed in a harmonic potential, V (x) = α 2 x 2 . This particular model has been studied in, e.g. , [52]. We refer to this model as active thermal Ornstein Uhlenbeck process (ATOU). While M 0 (see Eq. (98)) has been long known, M 1 (see Eq. (99)) 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 ν. 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 The paper is structured as follows. In Sec. II we give detailed account of how to calculate F (s) for small ν. First, we reproduce the Darling-Siegert argument in the equilibrium case (ν = ε = 0). Next, we introduce a perturbative version of the Darling-Siegert equation. Then, we obtain, as an intermediate result, a formula for F (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. (74) 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. (99). 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 momentgenerating function varies smoothly as ε, the coupling to the self-propelling force, is switched on. The momentgenerating function of the equilibrium version of the process (ε = 0) is assumed to be known in closed form, as is for instance justified for the Ornstein-Uhlenbeck process or Brownian Motion ( [53]). 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 ( [15,16]). Next, we outline our perturbative approach to the active case.
Before going further, we introduce some notations. The transition probability density of progressing from x 0 at t 0 to x 1 at t 1 is denoted by where the subscripts x 0 and x 1 are 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 Further, the first-passage time density to first reach x 1 starting from x 0 is denoted by In the following, we denote the Fourier transform of a function f (t) byf where F denotes the Fourier transform operation, with inverse where dω = dω 2π . 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 diagonal diagonal diagonal diagonal diagonal diagonal diagonal diagonal diagonal diagonal diagonal diagonal diagonal diagonal diagonal diagonal diagonal. The Fourier transform in two variables of these functions read

Equilibrium case: The Darling-Siegert solution
We here consider the equilibrium case of Eq. (1) defined by setting ε = 0. As x t is Markovian, the functionŝ F andT satisfy the following renewal equation: for all x 2 ∈ (x 0 , x 1 ]. Applying a Fourier transform to Eq. (11), the timehomogeneity of both T (t 0 , t 1 ) and F (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 Since the Fourier transform of a probability density equals its characteristic function, Eq. (12) recovers all moments of the first passage time providedT is known. Further, setting ω = −is for some s ∈ 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.

Out-of-equilibrium: A perturbative approach
The argument made by Darling and Siegert breaks down when ε = 0: indeed when averaging over the driving noise y(t), the renewal equation Eq. (11) no longer is true. The approach we take in this paper, consists of three steps: 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 2. As long as y(t) is fixed, the process, when un-derstood as conditioned on this particular driving, satisfies a renewal equation of the type Eq. (11). Inserting the perturbative transition and return probabilities from the previous step, gives a perturbative series for the first-passage time densitŷ F (ω 0 , ω 1 ; [y]) of x t conditioned on a particular y t .
3. Averaging over the ensemble of driving noises. For simplicity, we here assume that the correlation function of the driving noise is given by where β is the inverse correlation time. Generally, when computing the term of order ε 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 ε reads In the next sections, we derive this relation in more details.

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 F (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. (11) still applies and gives rise to: (16) 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 F (t 0 , t ; [y]) vanishes for t < t 0 and R(t , t 1 ) vanishes for t > 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 F , as done in Eq. (12), since neither terms in the integral are diagonal, i.e. they depend explicitly on both t 0 , t and t , t 1 . However, introducing the inverse functional R −1 which is defined by the implicit equation: the renewal equation can be formally solved by the relation: Our approach is then to perform a functional expansion of the quantities involved in Eq. (18) in the function y, around the Markovian case y = 0.

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 ): 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ỹ(t), according to: and we will use the fact that: for any time-shift τ . 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; such that Eq. (21) gives, choosing τ = −t 0 : Using Eq. (23), the functional derivatives of T can be related to the functional derivatives of T 0 through (24) and the functional expansion of T can now be written As a result, calculating the Fourier transform of Eq. (25) with respect to the variables t 0 , t 1 , we obtain the functional expansion: The perturbative expansion of the return probability R is simply obtained from by letting x 0 → x 1 in the expansion of T and has therefore the same structure, using R 0 (τ ; [ỹ]) = R(0, τ ; [ỹ]). To ease the notation, we now introduce a shorter notation for the set of functions: where the factor 1/ε n has been introduced for convenience, to make the expansion in small ε explicit in Eq. (27). We then arrive at the following form: where we have made a choice of signs ofω 1 , ..ω n for convenience.

Functional expansion of the inverse of the return probability density
We now turn to the expansion of R −1 in Eq. (18). As T and R, R −1 satisfy the time-shift invariance relation (21). This can be seen by time shifting Eq. (17) by the time τ : Eq. (17) defines R −1 . The analysis of the previous sub-section still applies and one obtains the expansion We then obtain by applying a Fourier transform to Eq. (17): This relation allows to relate the functions (R −1 ) (n) to the functions R (n) , by plugging Eqs. (30) and (32) and identifying order by order in ε. Limiting ourselves up to order 2, we obtain:

Second-order expansion of the first passage time distribution
Equipped with these expansions, one now can expand the first-passage time density expression (18) to obtain a functional expansion of the first-passage density in ε, involving the functionsT (n) andR (n) which are simpler to calculate. We first need to write down the Fourier-transformed version of Eq. (18): where the dependency on x 0 , x 1 is here implicit. Performing an expansion of this relation in ε, the result reads to second order: At this stage, we have obtained a perturbative expansion of the first-passage time density for a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a particular realisation of a 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.
B. 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 ε = 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. (22) and (28). The transition density T (0) solves the Kolmogorov forward equation where we introduce the forward evolution operator L x as where f is a twice differentiable test function. In Eq. (39), we denote the forward operator as L x1 to indicate that its gradient terms are acting on the x 1 dependency of T (0) . Correspondingly, the L 2 -adjoint operator L † x , also referred to as backward operator, is The forward operator L x has a countable set of eigenfunctions {u n (x)} and a non-positive spectrum 0 ≥ −λ 0 > −λ 1 > . . . [54], [55] however that the operator is self-adjoint and that therefore the family of as eigenfunctions of L x , form an orthogonal set of eigenfunctions spanning L 2 (R). We further impose a choice of normalization of the functions {u n }, such that the set as left eigenfunctions, satisfying we obtain a biorthogonal system for the pair {u n (x)}, {v n (x)}:ˆd 1 Instead, it is self-adjoint in the weighted space L 2 (u 0 (x)), [54].
This is useful to solve the forward equation; taking the Fourier transform in time of Eq. (40), one obtains Inserting the ansatẑ into Eq. (48) and using Eq. (42) leads to where we made use of the decomposition of unity, Since the u n (x 1 ) are linearly independent, their prefactors in Eq. (50) need to agree. Therefore, implying, together with Eq. (49), Turning to the case of ε = 0, the transition probability of the driven Langevin equation (1) solves the forward equation where time-homogeneity can no longer be assumed. Under Fourier transform, this forward equation becomes where the y-dependent term turns from a product into a convolution under the Fourier transform. We develop a perturbative solution ofT (x 0 , x 1 , ω 0 , ω 1 , [ŷ]) in powers of y.
We first discuss the first-order correction in ε. Following the functional expansion (30), and using the zeroth order result (53), we obtain: and we wish to determine the functionT (1) . Since the u n (x) span the L 2 -space, and in analogy to the ansatz (49), the first-order correction too can be written as a sum n (x 0 ; ω 1 ,ω 1 ) u n (x 1 ). (57) Re-inserting Eqs. (56) and (57) into Eq. (55) causes all terms to zeroth order in y to cancel, and one obtains an equation relating the contributions proportional to ε, The right hand side, which is the convolution ofT (0) (ω) and y(ω), no longer sums over u n (x 1 ) but their derivative ∂ x1 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 where we refer to the ∆ nk as derivative coupling matrix whose entries, as follows from bi-orthogonality, are Further, Eq. (58) holds true for arbitrary y(ω 1 ). In order to compare both integrations over dω 1 , we relabel ω 1 → −ω 1 in the left-hand side of the equation. Using this notation, inserting the sum (59) into Eq. (58), and resolving the ansatz (57), one obtains: In a similar way, the second order correction can be found; using the functional expansion (30) to second or-der: with the results from Eqs. (53) and (61) to zeroth and first order, inserting this ansatz into the forward equation (55) gives, following in complete analogy to the previous steps, Following this method, it is straightforward to generate the perturbative terms ofT (n) to arbitrary order in n, Finally, choosing x 0 = x 1 in any of the expressions (53), (61), (63) and (64) gives the corresponding terms for the return probability coefficientsR (0) (ω) ,R (1) (ω,ω 1 ) , .... Equipped with these expressions, we are able to compute the relevant integrals in the formula for the y-averaged first-passage time density Eq. (66).

C. 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 underlying stochastic process x t , in that sense akin to "quenched disorder averages" which replace a stochastic background field by an effective deterministic correction to observables. The expansion in Eq. (38) is a power series in orders of ε, where contributions of order ε 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 nth order response functionals of the first-passage time distribution (in s) to perturbations in the driving noise y. To calculate the y-average ofF (ω 0 , ω 2 ; [y]), each term in Eq. (38) is integrated over the path-measure of y, P[y]. The order of internal integration and y-averaging can be swapped. Consequently, since y = 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 where the last equality arises from the stationary in time of the random process y(t). Stationary in time. also implies thatĈ 2 (ω 2 ) is symmetric inω → −ω. Performing an averaging over P[ŷ] of Eq. (38), we obtain: The first term, of zeroth order, represents the Darling-Siegert solution (12). 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 δ(ω 0 + ω 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 ε are labelled (I) to (IV), and need to be calculated explicitly. The corresponding expressions are derived in the following section, for the case of Gaussian coloured noise.

D. Second-order correction with Gaussian coloured noise
So far, in our derivation of the second-order correction to the first-passage time density, Eq. (66), 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 [21]. 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 nth 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 O(ε 2 ) hold for any any any any any any any any any any any any any any any any any noise with such autocorrelation, in particular telegraphic noise as used in run-and-tumble processes [56].
Gaussian coloured noise is defined by its exponential correlator, which in Fourier space reads (Eq. (65)) With the explicit expressions (53),(61), (63), we perform the integration in (I) (see Eq. (66) for notation) in eigenfunction-representation, where in the last equality we employed Cauchy's residue theorem closing the contour in the lower half-plane containing the simple pole atω = −iβ. Likewise, we find where again the integral is evaluated by closing the contour in the lower half-plane enclosing the pole atω = −iβ.
The integrals (II) and (III), featuringω-dependent de-nominators, require some more careful analysis. We have 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 ω = iβ stemming from the correlator, as can be confirmed by inspection of Eq. (61). Besides, the denomina-torR (0) (ω 1 −ω) does not have any roots for (ω) > 0. Indeed, using the relations Eq. (53) and (45), i(ω1−ω)+λ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, λ n + (ω) > 0, assuming ω 1 ∈ 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 By analogous reasoning one obtains All four terms together give a general formula for the moment generating function of first-passage times for arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary underlying processes and driving noises arbitrary 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 This general result concludes this section. In the next section, we consider two concrete examples to demon-strate how this perturbation theory can be turned into analytical results.

A. 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 Equatioṅ with driving noise correlator (see Eq. (67)) This process reduces to the simple Ornstein Uhlenbeck process when ε = 0 which models a particle in a harmonic potential (V (x) = α 2 x 2 ) within a thermal bath. We consider, however, the process driven by an additional "active" term εy t . We therefore refer to this process as active thermal Ornstein Uhlenbeck process (ATOU). In the undriven case (ε = 0), the dynamics are characterised by the time and length-scales α −1 and = D x α −1 .
The coupling matrix is not diagonal: incoming momentum n is upgraded to outgoing momentum n + 1 by the noise coupling. The amplitude of ∆ m,n grows like |∆ m,n | ∼ √ λ m . For later use, we also note that By L 2 -adjointness, it follows that the adjoint creation and annihilation operators are In order to compute the transition and return probabilities, the following identity [57] proves to be useful. We introduce all quantities in dimensionless form, and take ω to be imaginary, such as the reduced frequency σ and the reduced autocorrelation time β as By considering the Fourier transform at σ = iα −1 ω, 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 rescale by as In the discussion that follows, we analyse all densities as densities in these dimensionless quantities to simplify notation and calculations. Following Eq. (53), one obtains for the transition probabilitŷ He n (x 0 ) He n (x 1 )e − x 1 2 2 n!(σα + αn) where we used identity (88) setting z = e −t . This integral is the Laplace transform of the Ornstein-Uhlenbeck propagator (in t ↔ σ), and its value is known in the literature to be ( [15]) where we introduced the parabolic cylinder functions D −σ (x) ( [58]). By continuity, for x 0 → x 1 it follows thatR In order to computeT (1) (−iσ, −iβ ), we use Eq. (61) and the derivative coupling matrix computed in (83) to find where we made use of relation (86) in the second equality. Letting x 0 → x 1 , one obtainsR (1) (−iσ, −iβ ). The counterpart,T (1) (−iσ − iβ , iβ ), is similarly found to bê These terms can be explicitly calculated and simplified. The rather lengthy but explicit expressions are given in appendix C. For the second order derivative term, using for-mula (63), one findŝ Again, the evaluated terms, including for x 0 → x 1 are given in appendix C. Equipped with the return and transition probabilities and its first two derivatives with respect to driving noise y (cf. (92)-(96)), we obtain the four contributions (69)-(73) which constitute the second-order correction formula (74). Whilst all the explicit expressions are given in appendix C, we here give the moment generating function in full as a undriven part and a perturbative correction, using s = ασ = iω, where we introduced the dimensionless parameter of expansion ν = ε 2 Dyβ Dxαβ . As is already known from literature (e.g. [16]), By symmetry (x 0 , x 1 ) ↔ (−x 1 , −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 x 0 < x 1 . The central result of this section then is Using for instance a computer algebra system like [60], all moments can be obtained by differentiation and eval-uating the limit of σ → 0 at which all derivatives have a removable singularity.

Numerical Validation
In order to corroborate the closed form result of the first-order correction to the moment generating function of first-passage times, Eq. For larger values of ν, 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 ν, is numerically estimated by taking the second numerical derivative, and is shown in the inset of Fig. 4. For 0.2 ≤ ν ≤ 0.8, the second-order corrections collapse, indicating that the deviations in the main figure are well accounted for by second-order corrections. For ν = 0.1,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 [59,60] of M 0 and M 1 to obtain the first-passage time distribution (cf. Eq. (7)) neglecting, as usual, terms of higher order in ν. In Fig. 2a, we show the first-order corrected distribution compared to a numerical sampling of the probability function,F (τ x0,x1 ), for values of 0 ≤ ν ≤ 3.2. The agreement for ν 1 is excellent, and plotting the rescaled deviationF against the theoretical result of F 1 (τ x0,x1 ) (see inset of Fig. 2a) shows systematic higher-order corrections consistent with our framework. Since for small ν the deviation to the undriven case is small,F (ν) suffers from statistical noise, therefore we omit showingF (ν = 0.1) in the inset. These numerical results therefore confirm the analytically obtained first-order correction to the momentgenerating 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 i , which we assume to expand in ν as T 1 (ν) = T 1 0 + νT 1 1 + ν 2 T 1 2 + . . . and T 2 (ν) = T 2 0 + νT 2 1 + ν 2 T 2 2 + . . .. The first-order corrections introduced are obtained by differentiation wrt σ using the result of Eq. (99) 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 and compare it to the result obtained from Eq. (104). 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 ∈ [0.05, 2]. The figures show a clear agreement with the theoretical result and systematic deviations for larger ν. Further, we observe that in this setting the coloured noise increases the mean-first passage time for smaller distances (x 1 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 momentsT 1 ,T 2 .

B. Active Thermal Brownian Motion on a ring (ATBM)
In this subsection, we consider the case of Brownian Motion x t driven by coloured noise with periodic boundary conditions (x ≡ x+2πr). The position of the particle satisfies the Langevin equatioṅ with    x 0 x 1 θ r y t , ξ t Figure 5. A Particle on a circle of radius r is driven by both white (thermal) and coloured (active) noise (cf. Eq. (107)). We study the first passage time distribution from x0 to x1 as a function of the angle θ ∈ [0, 2π).
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 firstpassage time distribution.

From eigenfunctions to the moment generating function of first-passage times
The Fokker Planck associated to the forward equation corresponding to Langevin Equation (107) is the Laplace equation with periodic boundary conditions x ≡ x + 2πr. Since its eigenfunctions are complex-valued square-integrable functions, we follow the identical steps of the framework introduced in Sec. II B but replace the scalar-product in Eq. (47) by its complex equivalent In what follows, we therefore consider the complex conjugate of the left eigenfunctions and obtain with corresponding eigenvalues where by allowing n ∈ Z, we enumerate the eigenfunctions efficiently (cf. Eq. (42)). The eigenfunctions are conjugate to each other since the forward operator of simple diffusion, L = D x ∂ 2 x , is self-adjoint. From this follows that the noise-coupling matrix defined in Eq. (60), is diagonal and purely imaginary. Because of scaleinvariance and rotational symmetry, we simplify the following discussion by introducing the dimensionless angle where we restrict ourselves to θ ∈ [0, 2π), and the diffusive timescale with which we rescale the Fourier-frequency ω to again effectively evaluating the Laplace transform (cf. Eq. (8)) of the respective probability densities. With this simplified notation, the transition density to zeroth order readŝ and the return probability, setting θ = 0, iŝ Assuming the y-averaged moment generating function has an expansion of where we use as in the Ornstein-Uhlenbeck case the dimensionless perturbative parameter ν = Dyε 2 β Dxα (here α has however a different definition than used in the OU case, see Eq. (116)). The zeroth order contribution is, using the classic result Eq. (12), and the results in (118), (119), which expands around σ = 0 as Being a moment generating function, the prefactors in front of σ n correspond to the (rescaled) moments of the first-passage time (−1) n n! τ n x0→x1 Dx r 2 n .
To second order, the transition density is where ∂ θ = r∂ x1 . Setting θ = 0, one obtains the second order response of the return probability, Inserting these quantities into the FPT correction (74), giveŝ This expression is the full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε full moment generating function up to second order in ε in elementary functions. In line with our previous notation, we identify the dimensionless scaling functions which together form the first-order correction of the moment-generating function An expansion around σ = 0 gives corrections to all moments. The correction to first order in ν of the mean first-passage time over an angle of θ, τ 0→θ , for instance reads where we indicated that the result can be written as the classical contribution (T 1 0 ) times the dimensionless per-turbative coefficient ν and a dimensionless scaling function T 1 1 (θ, β ). By successive derivation, any higher order moment may be obtained from Eq. (127).

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 A 2. Using Monte-Carlo simulations, we sample the first-passage times τ i of the integrated stochastic equation (107). To validate the moment-generating function, we average over N 10 6 to 10 7 iterations to samplẽ for various values of ν and across x 1 − x 0 ∈ (0, π] (with r = 1). Again, symbols with a tilde denote quantities which are measured numerically. To validate the theoretically predicted first-order correction Eq. (129), we subtract the ν = 0 contribution, M 0 (cf. Eq. (128)), and rescale by ν,M where we introduced, in analogy to Eq. (101), as shorthand for numerically measured higher-order corrections. In Fig. 7, we show the analytic result M 1 (cf. Eq. (129)) and numerically obtainedM 1 (ν) (cf. Eq. (133)) for various values of 0 ≤ ν ≤ 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 ν. The rescaled discrepancy,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. (129). By numerically computing the inverse Laplace transform, we obtain the perturbative correction to the first-passage time distribution (cf. Eq. (102)) which in Fig. 2b is compared to the numerically sampled distribution. The inset shows the numerical estimate of the first order correction (cf. Eq. (103)) which agrees well for small ν. For larger values of ν, second-order corrections play an increasing role as indicated in Fig. 7. In Figs. 6a and 6b, we show the first and second moment of first-passage times and how its deviation to the ν = 0 case is captured by the first-order correction obtained using Eq. (105) and Eq. (106) with the result of Eq. (129). 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 ν, indicating the systematic decrease for increased values of ν.

C. Limit Cases
The framework we introduce here allows to study coloured driving noises at any noise-colour β. In particular, this includes two limit cases of β → 0 and β → ∞. 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.

The white noise limit
For very small autocorrelation times β −1 the driving noise y t appears more and more as white noise. The correlator of y t tends towards In this limit, the driving noise features in the Langevin equation (1) as additional white noise and is absorbed aṡ such that effectively the diffusion constant is shifted by D x → D x + ε 2 D y . In the white noise limit, the theory is Markovian and Eq. (12) 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: One can verify that this relation is indeed satisfied for the Brownian motion case (Eqs. (128) and (129)) as well as the Ornstein-Uhlenbeck case (Eqs. (98) and (99)).

Quenched Disorder limit
In the opposite limit of β → 0, provided D y β = w 2 remains fixed, the driving noise "freezes" to a random constant since Effectively, the Langevin equation (1) therefore turns intȯ   where v is a constant driving velocity which is normal distributed according to v ∼ N (0, ε 2 w 2 ). The driving noise average e −sτ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 Our framework therefore predicts the first-order correction in 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 B, one can compute the moment generating function of first-passage times for a particular fixed drift v exactly (see Appendix A and in particular Eq. (A16) 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 N (0, ε 2 w 2 ) we obtain a resulting quenched average approximation in ... When employing our framework and letting β → 0 in our general result Eq. (127) we recover precisely . The necessary calculations are given in App. B 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 ν as (compare to Eq. (131)) such that quenched disorder lowers the mean firstpassage time for any choice of parameters. Further, for the mean squared first-passage time we obtain where ϕ = 1+ √ 5 2 is the golden ratio, and ψ ± = 1± √ 10 √ 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 self-propelled particles in a thermal environment. The self-propulsion is modelled by a noise with exponential autocorrelation and characteristic timescale β −1 , while the thermal bath is modelled by white noise with diffusion constant D x . The expansion parameter in which the perturbation takes place is a dimensionless quantity, ν, which indicates how strong the fluctuations of the self-propulsion are in comparison to the strength of thermal fluctuations.
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 (18) stands at the centre of this work. In order to solve it perturbatively, one needs to calculate the expansion terms (cf. Eqs. (30) and (32)) which involve the eigenfunctions of the Fokker Planck operator associated to the non driven process (cf. Eq. (42)). To first order in ν, we obtain an analytic result of the moment-generating function in terms of the associated eigenfunctions (cf. Eq. (74)). Higher order contributions can be obtained by further iterating the steps outlined in Sec. II B.
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 A 1, we calculate all necessary response functions to find the first order correction to the moment-generating function of firstpassage times (cf. Eq. (99)). 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 ν, 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 B). 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 β −1 , we are able to recover the limiting cases for β → 0 and ∞, respectively. The case of β → 0 is of particular interest since it recovers quenched disorder averages over processes with additional fixed and normal distributed drift (see Sec. III C).
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. (40)) can be diagonalised. This significantly extends the range of systems our framework can be applied to. In this work, we focused on Fouriermodes 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. (40) 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 + sgn(x)∆D).
Moreover, the functional expansion inŷ(ω) (cf. Eq. (27)) 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 ν, the corrections, as given in Eqs. (69) to (73), 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 Feynmandiagrams 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.
The results obtained in this work can be derived alternatively using field-theoretic methods. In fact, the field theoretic treatment allows for the study of a broader class of extreme events and will be reported in [61]. and initial condition T (x 0 , x, t = 0) = δ(x − x 0 ) (A5) The net flux of the surviving probability over the boundary is related to the first-passage time distribution via F (x 0 , x 1 , t) = −∂ tˆx x1−2πr dx T (x 0 , x, t). (A6) Using the backward equation (A2), and applying a further derivative in time one finds At the boundary, the first-passage time distribution is given by since the particle is immediately absorbed, while at time t = 0 for x 0 away from the boundary F (x 0 , x 1 , t = 0) = 0 (A10) since the particle has not yet begun to diffuse. Under Laplace transform in t, the differential equation (A8), together with these boundary and initial conditions returns The ordinary differential equation (away from the boundary) is solved by the exponential ansatz F (x 0 , x 1 , s) = Ae ω1x0 + Be ω2x0 (A12) Inserting this ansatz into (A11) enforces The boundary conditions (A11) fix the normalising constants A, B to A = −e πrω1−ω1x1 sinh(πrω 2 ) sinh(πr(ω 1 − ω 2 )) (A14) B = sinh(πrω 1 ) e πrω2−ω2x1 sinh(πr(ω 1 − ω 2 )) (A15) After some further simplifications one arrives at the v-dependent moment generating function with θ = (x 1 − x 0 )/r following the notation from the main text. We note that for v → 0, and one recovers the undriven moment generating function F (x 0 , x 1 , s; v = 0) = F (x 0 , x 1 , s) = sinh (θ √ σ) − sinh ((θ − 2π) √ σ) in agreement with the independently found expression (122). In App. B, 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 β → 0.