Field theory of survival probabilities, extreme values, first passage times, and mean span of non-Markovian stochastic processes

We provide a perturbative framework to calculate extreme events of non-Markovian processes, by mapping the stochastic process to a two-species reaction diffusion process in a Doi-Peliti field theory combined with the Martin-Siggia-Rose formalism. This field theory treats interactions and the effect of external, possibly self-correlated noise in a perturbation about a Markovian process, thereby providing a systematic, diagrammatic approach to extreme events. We apply the formalism to Brownian Motion and calculate its survival probability distribution subject to self-correlated noise.


I. INTRODUCTION
Many non-equilibrium systems are studied by projecting out a single slow degree of freedom which evolves stochastically and often displays non-negligible memory effects [1][2][3]. A classic object of study is its survival probability which describes the probability of the degree of freedom not having reached a threshold yet [4,5]. The survival probability defines not only the persistence exponents [6], but is also closely linked to the distribution of first-passage times, running maxima, and spans via some simple relations. All of these extreme events aptly characterise the non-equilibrium nature of complex systems and have been studied separately over the last hundred years (for classic references see [7][8][9][10][11], for recent overviews [12,13]).
Survival probabilities of non-Markovian processes are, however, notoriously hard to compute as they depend on the entire trajectory whose distribution is usually impossible to obtain [13][14][15][16]. Perturbative schemes such as [17,18] have proven to be successful in characterising the behaviour of the survival probabilities for large times in classical non-equilibrium models. More recently, a similar perturbation theory has been applied to fractional Brownian Motion to further access the full survival probability in a perturbation theory about the Hurst parameter [15,19]. These techniques, however, heavily rely on Gaussianity and do not readily translate to more general non-Gaussian non-Markovian processes.
In this article, we compute the full survival probability of processes subject to both uncorrelated and selfcorrelated noise in a perturbation theory in the strength of the self-correlated noise. A physical realisation of such * bwalter@sissa.it processes is a particle immersed in a heat bath and subject to a random self-propelling force. Our perturbative framework is valid in the regime where the self-correlated noise is small compared to thermal fluctuations stemming from the heat bath. More generally, this type of process is central to the study of active matter [20][21][22][23][24] and nonequilibrium phenomena [25][26][27].
In our recent work [28], we presented a scheme to calculate first-passage distributions for the same class of non-Markovian processes. These results relied on a perturbative functional expansion of a renewal type equation inspired by the classical work of [10,29]. Using a field-theoretic approach that draws on both the Doi-Peliti [30,31] as well as the Martin-Siggia-Rose formalism [32], in the present work we generalise these results to a broader class of extreme events.
Using a field theory has some notable advantages. Firstly, the diagrammatics give a clear intuition of the underlying microscopic processes otherwise hidden within cumbersome expressions. Secondly, the field theory provides a systematic perturbative framework naturally drawing on renormalisation techniques. Thirdly, it is easily extended to incorporate further interactions, such as reactions, external and pair potentials.
This article is structured as follows. First, we map a Markovian process onto a field theory. Secondly, we introduce a field-theoretic mechanism which is designed to keep track of the space already visited by the process. This defines the visit probability Q(x 0 , x, t), the probability that the process started at x 0 has been at x prior to time t, and which is the complement of the survival probability. Thirdly, we add the self-correlated driving noise, thus breaking Markovianity, and compute the corrections induced in Q. Finally, we illustrate the approach by computing the correction to the survival probability of Brownian Motion driven by self-correlated noise.

II. FIELD THEORY FOR MARKOVIAN VISIT PROBABILITIES A. Markovian transition probabilities
In this article, we construct a perturbation theory around Markovian processes characterised by a Langevin Equation [33],ẋ where V (x t ) is the gradient of a potential and ξ t Gaussian white noise with correlator ξ t ξ t = 2D x δ(t − t ). Further, we introduce T (x, t) ≡ T (x 0 , x; t 0 , t) as the transition probability for the walker to travel from x 0 at time t 0 to x at time t. This probability density is also known as Green's function or propagator in related fields of mathematics. We will state x 0 and t 0 only where needed for clarity. The transition probability satisfies a Fokker-Planck equation [34] ∂ with initial condition T (x, t 0 ) = δ(x − x 0 ). As is detailed in [35], the process (1) can be mapped to a Doi-Peliti field theory [35][36][37] containing two fields, the annihilator field ϕ(x, t) and the creator field ϕ † (x, t) = 1 + ϕ(x, t), which are jointly distributed according to Here the action S ϕ [ϕ, ϕ] is constructed as Moreover, the transition probability satisfies T (x, t) = ϕ(x, t)(1 + ϕ(x 0 , t 0 )) Sϕ , where • Sϕ denotes the expectation over the measure (3). Constructing a solution to the partial differential equation in Eq. (2) via a path integral can in principle be done with the Feynman-Kac theorem [38,39]. Here, however, we use a non-equilibrium field theory following [32] (see Sec. V A 2 for further discussion).

B. Visit probability and extreme events
The key problem we address here is how to approximate the distribution of first-passage times, running maxima, and mean volume explored of the process defined in Eq. (1). These extreme events are all mutually related via the visit probability which we define as Q(x 0 , x, t 0 , t) = Q(x, t) = P[x s = x at some time t 0 ≤ s ≤ t], i.e., the complement of the survival probability x 0 x 1x Figure 1. A continuous random path xt (red line) evolves in time. The visited area is shaded in blue. In this article, we study the visit probability Q(x, t), i.e., the probability that a point (x, t) lies within the blue area. As discussed in Sec. II B, the visit probability allows one to compute the distribution of (i) first-passage times τx 0 ,x 1 (blue dashed line, Eq. (6)) and (ii) running maximaxt (green line, Eq. (7)), as well as (iii) the average volume explored (purple arrow, Eq. (8)). P surv = 1 − Q. This measures the probability that the particle has been at x at or before time t. In Fig. 1, we show a single realisation of x t , together with its visited area.
The visit probability contains various informations about the process: When taking the derivative ∂ t Q(x, t), one measures the weight of those paths which visit x at t for the first time. The latter is the first-passage time, shown in Fig. 1 as a blue dashed line, and thus its distribution satisfies Analogously, taking the derivative −∂ x Q(x, t), weighs those paths who at time t visit x for the first time, or alternatively, the distribution of the maximumx t = max s≤t x s , shown as a green line in Fig. 1, i.e., for x > x 0 . Moreover, integrating over dx Q(x, t), gives the average of the volume explored, which is defined as the difference between running maximum and minimum, Vol([x t ], t) = max s≤t x s − min s≤t x s , as illustrated in Fig. 1. Its mean is given by Higher moments of the volume explored are considered in [40].

C. Overview of main results
In the following, we build a framework to compute the visit probability Q(x, t) for a specific class of non-Markovian processes. These are given by the solution to the stochastic differential equatioṅ This equation extends the Markovian Langevin equation (1) by adding a second independent noise term y t which is assumed to be stationary, of zero mean, but not necessarily Gaussian. The driving noise y t further carries dimensions of a velocity, leaving g as a dimensionless coupling constant which we suppose to be small. The main result of this article then is a perturbative expansion of the visit probability of x t to leading order in g 2 , i.e., assuming the visit probability allows for an analytical expansion around g = 0 as Q(x, t) = , we find formulas for the correction terms. As is further detailed below, the assumption that Q(x, t) be analytic in g restricts the possible choice of driving noises depending on the choice of potential V (x t ). Together, Eq. (28) and Eq. (29) provide a general formula for the leading perturbative correction, Q (2) (x, t) , which is expressed in terms of the Markovian transition probability T (x, t) and the two-time correlation function of y t , where • denotes the average with respect to the path measure of y t . In principle, the framework also allows to compute higher-order corrections, i.e., g n Q (n) (x, t), using the n point correlations of the driving noise. In the presentation of the results, however, we restrict ourselves to the leading order perturbation only.
The processes described by Eq. (9) do not satisfy a (generalised) fluctuation-dissipation relation, and hence cannot be brought into the form of a generalised Langevin Equation. Instead, these processes are often used to model active matter in thermal environments [20][21][22][23][24][25][26][27], which typically operate away from equilibrium.
Finally, although the expression for the visit probability applies to all potentials V (x), and can be employed numerically to study these, an analytically closed expression can only be expected in cases where an analytic solution to the Markovian Fokker Planck equation (2) is known. This effectively reduces the class of potentials for which we obtain analytical results to harmonic or flat potentials, i.e., perturbations of Brownian Motion or Ornstein-Uhlenbeck processes.

D. Markovian visit probabilities
In this section, we present a field theory of visit probabilities for Markovian processes. Whilst in the case of transition probabilities it is well known that the solution to the Fokker Planck equation (2) can be expressed as a path integral, Eq. (5), this has not yet been established for the visit probability Q(x, t). Our aim is to construct a field theory whose correlation functions equal Q(x, t), in close analogy to Eq. (5). In difference to the case for the transition probability, however, this field theory cannot be straightforwardly constructed for Q(x, t), since no evolution equation for Q(x, t), comparable to Eq. (2), exists to our knowledge.
As is explained in great detail in [40][41][42], and briefly discussed in App. A, the visit probability Q(x, t) can be expressed as a field-theoretic expectation value under a Doi-Peliti field theory by introducing two additional auxiliary ("trace"-) fields ψ(x, t) and ψ(x, t) with a joint distribution such that the visit probability can be written as where • S is understood as the average with respect to the measure in Eq. (11). Here, we have introduced a normalising density n 0 which is further detailed below and in App. A. The pair of fields ψ, ψ are a stochastic auxiliary variable which tracks the volume explored by the process x t (see Fig. 1). This is in analogy to ϕ, ϕ whose correlation tracks the current position of the process x t (cf. Eq. (5)). Hence, measuring the average field density ψ(x, t), and normalising by a unit density n 0 , amounts to computing the probability that the process visited x up to time t. The average in Eq. (12) then corresponds to the probability density that x has been visited prior to t conditioned on the process having been initialised at x 0 at time t 0 ; thus matching our definition of the visit probability.
The field action S ϕ + S ψ − γS γ consists of three actions which model (i) the diffusion of the process x t (cf. Eq. (4)), (ii) the non-interacting dynamics of the auxiliary fields tracking the volume explored by the particle, and (iii) the interaction between the random process x t and the auxiliary fields tracking its explored volume. Clearly, the explored volume depends on all the previous positions of the process x t and therefore the third contribution contains all four fields ϕ, ϕ, ψ, ψ. It is multiplied with a rate γ to be taken to ∞. This rate is interpreted as the rate with which the process x t "traces", i.e., marks as explored, a given point x. Taking γ → ∞ amounts to the field ψ tracking every visited point. 1 In what follows, we outline the components of the action in Eq. (11) and refer the reader to the appendices for technical details.
The first term in the exponential of Eq. (11) is the diffusion action S ϕ introduced in (4). The second term S ψ denotes the trace action Comparing with S ϕ , Eq. (4), the trace action can be interpreted as corresponding to the process (1) with V (x) ≡ 0, D x = 0, i.e., a deterministic immobile particle. This reflects the fact that a point, once visited by the process x t , remains visited forever. To ensure convergence of the path integral over (11), we have included a positive parameter ε > 0. Latter further ensures causality [40,44], i.e., that a point is only marked visited after it has been visited by the process x t . In a field-theoretic context, the parameter ε is referred to as a mass or an infrared regulator as it suppresses the divergencies otherwise arising in Eq. (11) from the contributions of Eq. (13) at large times. The parameter ε is to be taken to zero at the end of the calculation. The third term of the exponential of Eq. (11) is the deposition action. It describes the growth of the volume explored due to fluctuations of the process x t and is derived in [41]. In the continuum limit, it reads The four dimensionfull couplings τ , σ, κ and λ, are introduced as differently for independent renormalisation. However, the rates τ and σ and the densities κ, λ are each equal and are related to each other via see App. A for details. Overall, S γ is multiplied with a dimensionless constant γ which needs to be taken to γ → ∞.
Each of the four vertices in Eq. (14) is diagrammatically represented as and enters into the action multiplied by γ, see Eq. (11).
Here, straight red lines represent the propagators of ϕ, ϕ, and green wriggly lines those of ψ, ψ. As a convention, we read vertices/diagrams from right to left. It follows that the diagrammatic expansion of the trace function (12) is where the central dot stands for the renormalised coupling τ R and the limit in γ → ∞ stems from the definition of the action in Eq. (11). This renormalisation is given by the diagrammatic expansion The only diagrams contributing to this expansions are chains of the loop-diagram . We introduce the return probability R (x, t) = T (x, x, t) and likewise its Fourier transform R (x, ω). As shown in App. B the fully renormalised vertex can be evaluated using a geometric sum as which is an exact result for all γ, such that the effective trace function, Fourier transformed, is where we made use of time-translational invariance to write the Fourier-transform in one frequency only, tacitly took the limit ε → 0, and used τ /κ = n 0 , see Eq. (15). Multiplying this result with (−iω) gives the first-passage time moment generating function. Then, the result in (20) agrees with classical results given in [10,29].
In summary, by introducing an extended field theory via the four-field action in Eq. (11), the visit probability of a Markovian process can be written as the fieldtheoretic average (12), in direct analogy to the (simpler) transition probability which can be represented with two fields, cf. Eq. (5). This effective field theory for visit probabilities comes at the cost of the additional fields ψ, ψ which need to be related to the fields ϕ, ϕ by a nonlinear interaction S γ . The effect of this interaction on the growth of the volume explored by x t can be captured by evaluating the effective ("renormalised") deposition vertex τ R . Formally, this amounts to evaluating an infinite series of correction terms, diagrammatically represented in Eq. (18). Using field-theoretic tools, the entire sum can be exactly evaluated, see Eq. (19). The result, Eq. (20), is in agreement with previous classical results. So far, we therefore have constructed a field theory for Markovian visit probabilities that reproduces known results. In the next section we consider non-Markovian processes, and use this field-theoretic formulation to compute the perturbative corrections to the exact result (20).

III. NON-MARKOVIAN VISIT PROBABILITIES: A PERTURBATIVE APPROACH
The process introduced in (1) is driven by δ-correlated noise and hence is Markovian [45]. As a perturbative generalisation towards non-Markovian processes, we thus extend the class of processes given by (1) tȯ The additional driving noise y t is assumed to be stationary with zero mean and a general autocorrelation function By • we denote in the following averages over the path distribution of y t . Importantly, y t is not required to be Gaussian such that higher non-trivial cumulants C n (t 2 − t 1 , ..., t n − t n−1 ) = y t1 ...y tn c may exist. Such higher order cumulants enter only at perturbative order g n . The correlation function C 2 (t − s) of the driving noise y t may decay exponentially, such as for run-and-tumble particles in noisy environments [20,22], or algebraically. In the following we perform a diagrammatic expansion of the visit probability Q(x, t) in the dimensionless coupling constant g. Therefore, g is assumed to be small (g 1). A condition for the validity of the expansion is that cumulants of the noise C n (t 2 − t 1 , ..., t n − t n−1 ) = y t1 ...y tn c are such that diagrams of the expansion are finite, for instance the diagram in Eq. (27). In what follows, we derive the visit probability Q(x, t) averaged over both, the Gaussian white noise ξ t and the driving noise y t to first leading perturbative order g 2 . For the case of Brownian motion driven by self-correlated noise, we find expressions that depend on double integral of the cumulant C 2 (t), Eq. (31). Existence of this double integral allows for a broad class of correlation functions, even when C 2 (t) decays algebraically slowly at large times. Such algebraic decay occurs in processes driven by fractional Gaussian noise [21].
As is outlined in the App. C, the visit probability conditioned on a fixed realisation of the driving noise can still be obtained by Eq. (12) using the average with respect to the modified path measure As is shown in the App. C, one can compute averages with respect to this y t -dependent path-average, to then integrate over the path measure of y t . As it turns out, one does not need to know the full path measure of y t , but can instead rewrite the double average • using the moment generating functional of y t which is defined as Expanding the exponential to second order, and averaging over y t , gives (cf. Eq. (22)) Therefore, we may approximate expectation values such as the one appearing in Eq. (23) using the identity which to this perturbative order only contains C 2 (t − s), the correlation function of y t (Eq. (22)). Diagrammatically speaking, this amounts to computing the same diagrams as in Eq. (18), but additionally decorating them with driving noise correlation loops represented as As in Eq. (16), red solid lines denote bare transition probabilities. The blue dashed line connecting two internal vertices, (cf. Eq. (25) and App. C), represents the correlation kernel C 2 (s 2 , s 1 ). The two vertical bars inserted to the right of each such vertex represent the gradient operator acting on the target point of the incoming transition probability. Combinatorally, there are four ways in which these loops decorate the diagrams of the (Markovian) visit probability, which are displayed in Fig. 2.
As is shown in the App. E, the new visit probability, , acquires a functionally similar form to the Markovian case (20) This is a central result of the present article. Here, as we do from now on, we denote the (Fourier transformed) Markovian transition probability density of the undriven process in Eq. (1) as . Further, we introduced the g 2 correction to the y t -averaged transition probability  (27)). The Markovian diagrams (cf. Eq. (18)) are corrected by the external driving noise g · yt, giving rise to the expansion Eq. (C11). These four types of diagrams result in four different correction terms explicitly derived in App. E, and which together give the final result for the driving noise corrected transition probability stated in Eq. (29).
and, analogously, In the case of a Brownian Motion additionally driven by self-correlated noise, the second order corrections given by Eq. (29) drastically simplifies as shall be demonstrated below.

IV. EXAMPLE: BROWNIAN MOTION DRIVEN BY SELF-CORRELATED NOISE
We briefly discuss the example process of active thermal Brownian motion (ATBM, [28]) on a real line. In our previous article [28], we found the perturbative correction to the first-passage time moment generating function for an active thermal Ornstein Uhlenbeck process, and an active Brownian Motion on a ring. Since the distribution of first-passage times and visits are closely linked via Q(x, ω) = (−iω) −1 χ FPT (x, ω), we report the visit probabilities of the latter two processes in App. G.
To illustrate the method outlined above we use the example of ATBM on a real line defined viȧ where ξ t is white noise of strength ξ s ξ t = 2D x δ(t − s), and y t is a stationary, not necessarily Gaussian, driving noise with a non-white self-correlation function C 2 (|s − t|) = y s y t . Without loss of generality, we assume x 0 = 0 and x > 0. The transition probability and a subsequent inverse Fourier transform. For g = 0, we compute the correction T (2) (x, ω) using Eq. (29). As is detailed in App. F, the integral simplifies drastically, and we state here only the results, which can be summarised most succinctly using a "time-stretch function" Υ(t) as follows. Given the correlation function C 2 (t), Eq. (22), it is defined as (31) or, alternatively, The second order correction to the transition probability can then be written as where we made use of the Brownian can be absorbed into the t-dependence of the tree-level as This result is in agreement with the expression for the full transition probability found in [46], and is exact if y t further is Gaussian. The externally driven non-Markovian process x t has therefore transition probabilities that are, to order g 2 , equal to those of a time-stretched Brown- That, however, does not mean that the return probabilities agree between the original, externally driven non-Markovian process x t and the timestretched Brownian motion, because of the latter not accounting for the now hidden variable y t (see also the discussion in [46]). The time-stretched Brownian motion ignores correlations between that y t and x t . For example, the transition probability is no longer correctly given by the first passage and repeated return, as first passage is favoured for particular values of y t that the subsequent return does not account for. Returning to Eq. (28), and using Eq. (32), we slightly rephrase the result for driven Brownian Motion by explicitly expanding in g 2 to This relation illuminates the relation between the correction to the Fourier-transformed visit probability and the correlation function of the driving noise. To obtain the visit probabilities, and derive extreme event distributions, in real time (cf. Eqs. (6), (7), (8)), numerical integration will be necessary in most cases. For exponentially correlated driving noise, such as "coloured" or telegraphic noise, we have C 2 (t) = D x βe −β|t| , Eq. (22), for some noise strength D x and timescale of relaxation β −1 . The corresponding timestretch function is From this function alone, one can read off that for large times the driving noise effectively shifts the diffusion constant by D x → D x + g 2 D y . At short time-scales (t β −1 ), however, the correction is non-trivial. We can further verify the validity of the expansion of the trace function Q(x, ω; g), Eq. (34), by numerically inverting the Fourier-transform and compare to Monte-Carlo simulations of the process Eq. (30). To this end, we estimate numerically the probability that x t has reached x 1 , at some given t and parameters such as g, D x and D y , and subtract from it the exact result Q(x, t; g = 0).
Plotting this difference over g 2 produces an estimate of the correction of Q(x, t; 0) to Q(x, t, g), described by Eq. (34) to leading order (see Eq. (G1) for explicit result). Fig. 3 shows this numerically estimated correction together with the inverted correction Eq. (34).
The persistence exponent θ is defined as the tail exponent of the survival probability P surv (x, t) = 1−Q(x, t) ∼ t −θ [4,6]. It is also contained in the small ω expansion of the trace as (−iω) −1 − Q(x, ω) ∼ ω θ−1 . Evaluating Eq. (34) using the time-stretch function (35) shows that the Markovian result of θ = 1 2 [4] does not acquire corrections as is expected for short-range correlated driving noise y t .

V. DISCUSSION AND SUMMARY
In this section, we discuss our findings in the context of the literature of (quantum) field theory and stochastic dynamics and summarise our results.
A. Discussion

Relation to Markovian techniques
Our study is motivated by the study of complex systems comprising many interacting degrees of freedom of which we single out the slowest one as as a stochastically evolving coordinate x t (see Sec. I). The remaining degrees of freedom are subsumed into a bath, exerting a stochastic force onto the particle, here modelled by Eq. (21).
Often, the fast degrees of freedom in a complex system are assumed to evolve infinitely fast, thus rendering the stochastic evolution Markovian, since all correlations disappear within an infinitesimal time [47]. This assumption is reflected in the mathematical structure of the usual stochastic representations, such as Langevin equations [48,49] or Fokker-Planck equations [34] which evolve locally in time, i.e., with no recourse to the past evolution. In this article, we set up a field-theoretic framework for Markovian processes in Sec. II. The "Markovianness" of the field-theory can be seen from the field action (14) which is local in time. At this stage, the field-theory is fully equivalent to any other Markovian description, and in fact reproduces the known Markovian result for the visit probability in Eq. (20).
The assumption of an infinitely fast evolving bath, however, is unphysical [50] and hence in each timeincrement the future stochastic evolution depends on the, potentially entire, past of the trajectory. This implies that the time-local formalisms mentioned above need to be extended to include non-local time interactions (e.g., to generalised Langevin equations [51] or fractional diffusion equations [52]). In our work, this self-interaction of the process with its own past is encapsulated by the non-local contribution (25). The selfinteraction is most clearly seen diagrammatically by the loop-diagram given in Eq. (27).
In this article, we begin by writing the solution of the Markovian (g = 0) forward equation (2), the transition probability T (x, t), as a path-integral Eq. (5) over fields whose distribution is given by the action Eq. (4). In setting up the field theory, Eq. (11), we constructed a nonequilibrium field theory using the Doi-Peliti formalism [30,31]. In principle, the transition probability can be obtained using alternative routes, following for instance the Feynman-Kac theorem [38,39] which expresses the solution to parabolic PDEs such as Eq. (2) as path integrals over trajectories rather than fields.
These alternative techniques, however, do not extend to either (i) the transition probability of driven, non-Markovian, processes, such as Eq. (21), which render the forward equation nonlocal in time, nor (ii) to the stochastic description of the visit probability Q(x, t) which to the best of our knowledge cannot be characterised as a solution to a parabolic PDE. This then requires two respective additional technical points with respect to alternative methods suitable to study transition probabilities of Markovian processes.
In order to address the average over the driving noise y t , we introduced Eq. (26) where we replaced the path integral over the driving noise by its moment generating function (or partition function) Z y . Following fieldtheoretic standard procedures [58,59], we expand the latter in a power series in g. This is analogous to the perturbative treatment of self-interaction in classical field theories.
Secondly, in order to cast the visit probability Q(x, t) (which does not readily follow from a Fokker-Planck equation) into a field theory, we used the tracing mechanism and its Doi-Peliti formulation to track the volume explored via the auxiliary fields ψ, ψ. This formalism (introduced in [41]) does not correspond to a classical field theoretic technique, but is rooted in the study of reactiondiffusion methods with field-theoretic methods.

Stochastic Dynamics
The problem of finding the transition probability, let alone the extreme event distribution, of non-Markovian processes has been a long-standing problem in stochastic dynamics. In terms of computing the transition probability, our field-theoretic approach recovers results known in the literature [46] in which non-δ-correlated driving noises are also treated using functional methods (see also [60,61] for a similar discussion of first-passage times). However, for non-Markovian processes the visit probability does not straight-forwardly follow from the transition probability. The visit probability Q(x, t) therefore represents the central result of the present work.

B. Summary
In the present work we have established a method to compute the visit probability (the complement of the survival probability) for random motion in a onedimensional potential, as defined in Eq. (21).
By mapping the problem to a field-theory we systematically compute corrections to the Markovian result (20) to any order in the coupling g. The leading order correction is of the form (28), which involves the g 2 -corrected transition probability given in (29). Generally, to compute contribution to order g n , it is sufficient to know the Markovian transition probability and the n-point function of the driving noise y t . In the absence of an external potential, the expressions reduce to (34) which depends on the twice-integrated driving noise correlation Υ(t), Eq. (31).
By casting the problem in a field-theoretic language, we replaced the single degree of freedom by a field ϕ representing the full density. This allows us to extend the model to the case of many, potentially interacting, random walkers. Also, the Doi-Peliti framework allows for the inclusion of potentials. To the best of our knowledge, this is the first time that both Doi-Peliti and Martin-Siggia-Rose have been used simultaneously to construct an action.
Overall, we have established a method which further lays bare the interplay between non-Markovianness and extreme events in stochastic processes.
The visit probability for Markovian processes (cf. (1)), given in (12), is a result that can be derived in various ways (indirectly, via the first-passage time distribution, ∂ t Q(x, t), this result has been found in, e.g., , [10] using a renewal type approach). Here, we derived the result by taking the continuum limit of a discrete reaction-diffusion process designed to track visits, which we refer to as "tracing mechanism" and which has been introduced in [40,41].
To describe the tracing mechanism, we consider a coarse-grained version x t of the stochastic process x t of Eq. (1), which takes values only on a lattice δ a Z where δ a is the lattice-spacing, so that formally x t = δ a [δ −1 a x t ] where [x] rounds to the nearest integer. At any time, the random walker attempts to deposit a trace at x t with Poissonian rate proportional to γ. Each lattice site, however, has a "carrying capacity" n 0 which limits the number of trace particles that can be deposited at this site. If at a site n 0 trace particles have already been deposited, any further deposition is suppressed. Hence, the number of trace particles deposited at any lattice site is an integer bound above by n 0 . Taking γ → ∞, the particle deterministically deposits trace particles at any site visited for the first time. In the limit of δ a → 0, the process x t tends to x t , and the expected number of trace particles at a site, divided by n 0 , converges to Q(x, t).
Cast into the language of reaction diffusion processes, we have at each lattice site i where particles of species W ("walkers") deposit particles of species T ("traces") at rate γ, provided their number does not surpass the carrying capacity n 0 . Meanwhile W diffuses according to (a discretised form of) (1), T remains at a given site and "evaporates" with rate ε, later sent to zero. The Doi-Peliti formalism [30,31] describes the continuum limit of particle densities in a reaction-diffusion system, such as defined in Eqs. (A1), (A2), by mapping the problem onto a non-equilibrium field theory, where each particle species corresponds to a pair of fields (see [41] for details). The local density of walkers δ −1 a W i is mapped to ϕ(x, t), ϕ(x, t) (referred to as annihilation and creation fields, respectively), and the trace particle density lim δa→0 δ −1 a T i to ψ(x, t), ψ(x, t). As shown in [40,41], the joint distribution of the four fields then follows from the Doi Peliti framework to be distributed according to the action given in Eqs. (4), (11), and (14). In order to turn the density of tracers into a probability of visit, it needs to be divided by the normalising density n 0 corresponding to the continuum limit of n 0 . Hence, the field-theoretic formula for the visit probability, Eq. (12), contains a prefactor of n −1 0 .

Appendix B: Field-theoretic Calculation of Markovian visit probability
We derive Eq. (20), the expression for the visit probability in the Markovian case. First, we consider the field theory in the case of γ = 0, when the probability measure (11) is Gaussian.
We introduce the forward and backward operators L, L † , associated to (1) and generating the corresponding Fokker Planck Equation with T (x, t) the transition probability of the process. The forward and backward operator have a set of eigenfunctions with eigenvalues 0 ≤ λ 0 < λ 1 .... The eigenfunctions are L 2 -normalised to satisfy the orthonormal relation and since they form a complete set in L 2 they further satisfy [34] n Thus, every field ϕ, ϕ, ψ, ψ has a unique decomposition into the u n (x), v n (x) in space. Together with a Fourier transform in time, we then write and analogously for ψ, ψ with coefficients ψ k (ω), ψ k (ω ), respectively. This mode transform diagonalises the nonperturbative parts of the action, S ϕ and S ψ (cf. Eq. (11)), which read For γ = 0, the measure in Eq. (11) is Gaussian and the (Fourier transformed) bare propagators of both fields, ϕ ϕ and ψ ψ , therefore immediately follow from Eqs. (B10) and (B11) using standard path integral techniques [62], where we introduced a diagrammatic representation for both bare propagators. These propagators are commonly interpreted as the (linear) response functions [62]. Using Eq. (5), and transforming back into real space and time using Eqs. (B8) and (B9), we obtain the transition probability where Θ(t) is the Heaviside Θ-function. Crucially, we made use of the property [44,62] ϕ S = 0 (B15) such that ϕ(1 + ϕ) = ϕ ϕ . We next consider the expectation (5) in the case of γ = 0 when the non-linear contributions of S γ , Eq. (14), enter. Each of the four vertices is diagrammatically represented as and enters into the action multiplied by γ. Following Refs. [40,41], we introduce the carrying capacity density n 0 which is the continuum limit of δ −1 a n 0 , where n 0 is the maximal number of trace particles which can simultaneously be deposited at a single site. At bare level, the carrying capacity enters in the couplings via the relation It follows that the diagrammatic expansion of the trace function (12), which counts the average number of tracer particles at a site, needs to be normalised with n −1 0 in order to indicate the visit probability, and hence where the central dot in the diagram stands for the renormalised coupling τ R , and T (x 0 , x, ω) is the Fourier transform of the transition probability. This renormalisation of τ is given by the diagrammatic expansion of the amputated vertex x t x 0 x Time Figure 4. The random walker xt (red solid path) travels from x0 < x1 to xt > x1, thereby passing x1 (green solid line) infinitely often (black dots). If the random walker is additionally driven by self-correlated noise (cf. Eq. (21)), this induces correlations between increments at any different times t1, t2 (blue dashed line). Meanwhile the perturbative expansion in γ (cf. Eq. (18)) tracks the probability of all possible transitions T (x1, t) and subsequent returns from x0 to x1, the second perturbative expansion in g includes the effect of correlated increments. If g = 0, the Markovian case, the process undergoes renewal at every return at x1, thus rendering the results such as Eq. (B23) exact.
The only diagrams contributing to this expansions are chains of the loop-diagram , refered to in the following as a "bubble". Considering the expansion in Fourier and mode transform, Eqs. (B8) and (B9) (App. D for details), each diagram factorises into a product over the bubbles and can hence be evaluated using a geometric or Dyson sum (App. D for derivation) resulting in such that the effective trace function, Fourier transformed, is where we made use of time-translational invariance to write the Fourier-transform in one frequency only. The couplings τ /κ in Eq. (B22) cancel with n −1 0 following their bare values, Eq. (B17). This then leads to the central Markovian result for the trace function where we have tacitly taken the limit ε → 0.

Appendix C: Visit probability for driven process
In this appendix, we provide some technical details to the derivation of the key result 28 and 29 which together provide the visit probability for non-Markovian processes of the form (21).
There are three technical steps: First, we consider the perturbative correction of the visit probability in the presence of a fixed, but random, realisation of the driving noise y t . Secondly, we average over all such realisations of y t . This then leads to a large set of correction terms which we interpret diagrammatically and which, thirdly, we evaluate to leading perturbative order.

Averaging general observables over driving noise
The presence of the autocorrelated driving noise y t affects the transition and return probabilities of the random walker x t (see Fig. 4). Conditioned on a fixed realistion of y t , the forward operator (cf. Eqs. (1), (B3), (B1) and (21)) is shifted by Hence, the y t -conditioned random walker action Eq. (4), which is essentially of the form ϕ(∂ t − L)ϕ, is shifted also, To obtain the y-averaged joint distribution of the four fields, one needs to evaluate (Eqs. (11), (23)) where gy t ϕ∂ x ϕ may be interpreted as a new vertex, (see Eq. (27). Whilst we were able to compute the result exactly in γ, we need to resort to perturbative methods to approximate in small couplings g. Although the path-integral in (23) appears to require complete knowledge of the full path measure P[y t ] of the driving noise y t , we can relax this requirement by introducing the normalised partition function (or moment generating function) of y t , where we ignore terms of higher perturbative order. We may thus replace the y t -integration in Eq. (23) by inserting Z y g dx dt ϕ∂ x ϕ into expectations over the fields and consequently evaluate double averages via For example, the transition probability, averaged over all driving noises y t , acquires a correction term which up to g 2 reads The remaining field average in (C7) can now be evaluated normally, using standard Wick product rules (see App. B), Eq. (5), and the non-driven path measure (11). A straightforward calculation shows that the only non-vanishing correction to order g 2 in Eq. (C7) is using ϕ(x 1 , t) ϕ(y 1 , s 1 ) = T (y 1 , x 1 , t − s 1 ), ∂ y1 ϕ(y 1 , t) ϕ(y 2 , s 1 ) = ∂ y1 T (y 2 , y 1 , s 1 − s 2 ) etc. Here, and in particular for more cumbersome expressions, it is advantageous to use diagrammatics to keep track of perturbative correction terms. The correction of the transition probability in Eq. (C8) is represented as As in Eq. (B12), red solid lines denote bare transition probabilities. The blue dashed line connecting two internal vertices, Eq. (23), represents the correlation kernel C 2 (s 2 , s 1 ). The two vertical bars inserted to the right of each such vertex represent the gradient operator acting on the target point of the incoming transition probability. As usual for Feynman diagrams, external fields depend on fixed parameters (x 0 , x, t), whilst internal fields depend on variables to be integrated over (e.g., y 1 , s 1 , y 2 , s 2 ). In addition to providing a better overview of terms arising in the perturbative expansion, diagrams also act as graphical cues illustrating how the driving noise induces memory into the evolution of x t .

The driving noise averaged visit probability
Following Eq. (26), we are able to evaluate driving noise averaged observables. In order to derive the central results of Eqs. (28) and (29), it remains to evaluate the driving noise averaged visit probability where we use ϕ i = ϕ(y i , s i ) for brevity. Each term appearing in the γ-perturbative expansion of the trace function (cf. (18)) is additionally corrected to order g 2 by replacing two internal ϕ ϕ propagators by two "y t -driven propagators" ϕ∂ x ϕ and connecting them with the two-point correlator C 2 of y t . As can be seen easiest diagrammatically, all possible corrections fall into one of the four categories shown in Fig. 2. They are classified according to whether both g-vertices couple to the same or different propagator of a transition or a return. It is simplest to compute the four contributions in frequency rather than direct time as the loops factorise. The calculation itself is given in App. E.

Appendix D: Renormalisation of transmutation rate
In this appendix, we derive the result for the renormalisation of the coupling τ which is stated in Eq. (19). To simplify the computation, we perform the calculation in x, ω variables, i.e., in real space and Fourier transformed time. The renormalisation is given by (cf. Eq. (18) which is a diagrammatic representation of the terms arising in the path-integrated average of the visit probability, Eq. (12), when expanding in γ, Crucially, the averages • S;γ=0 are taken over Gaussian random variables since for γ = 0 the path action (Eq. (11)) is bilinear. Thus, Wick's theorem (e.g., [58]) applies and all averages in Eq. (D2) decompose into products of two-point functions. The only such Gaussian two-point functions which do not vanish are The second correlator intuitively characterises the behaviour of the trace which, once deposited at z 0 at time t 0 , remains there for an infinitely long time, as ε → 0. Equipped with these correlators, the non-vanishing contributions to the averages appearing in Eq. (D2) are, following Wick's theorem, and to second order, where in the first equality we used the definition of the return probability to abbreviate and in the second equality used Cauchy's residue formula to solve the integral by evaluating the residue of the simple pole at ω 2 = −iε.
Since both correlators in Eq. (D3) are proportional to δ¯(ω 0 +ω 1 ), all higher order expansion terms factorise, after integrating over the internal frequencies, into a product over amputated one-loop bubble diagrams (i.e., interpreted here as a function of external parameters z 1 , ω 1 and z 2 , ω 2 , respectively) which by analogous reasoning to the calculation above evaluate as The bubble diagram graphically encodes the probability of a particle depositing a trace and then returning to it, in other words a "time-ordered" return probability: The green wriggly line may be understood as the trace which once placed remains immobile, meanwhile the red solid line represents the diffusing, and returning, walker. Likewise, the higher order diagrams in Eq. (18) may be interpreted as repeated returns to x. Returning to Eq. (18), the renormalised τ coupling, τ R , is the effective factor satisfying and collecting the factors generated by the terms in Eq. (D2), and evaluated using Eq. (D10), one obtains This series can be resummed using the geometric series, in field theory often referred to as Dyson summation [58].
Rearranging the sum gives where we made use of the bare values given in Eq. (14), i.e., replaced σ with τ and used λ = κ at bare level. This vertex interpolates the physical pictures for γ = 0, where no deposition takes place (τ R = 0), and γ → ∞, where every newly visited site gets marked immediately by a deposited trace. For γ → ∞, the coupling tends to using n 0 = τ /κ.
As usual, we employ Wick's theorem to evaluate this average over Gaussian random variables: When expanding the average n −1 0 ψ(x, ω) ϕ(x 0 , ω ) in powers of γ and g, as shown in Eq. (E2), the resulting coefficients are averages over finite products of ϕ 1 , ...ϕ j1 , ϕ 1 ...ϕ j2 , ψ 1 ...ψ j3 and ψ 1 .. ψ j4 fields (where ϕ i = ϕ(z i , ω i ) etc.) with respect to the Gaussian measure (defined by the action in Eq. (11) for γ = 0). Each of those coefficients is evaluated using Wick's theorem, i.e., where the sum runs over all possible pairwise pairings of the j 1 + j 2 + j 3 + j 4 indices, and we use φ in lieu of ϕ, ϕ, ψ, ψ to alleviate notation. Although the number of possible combinations of such pairings is very large, the right hand side of Eq. (E3) drastically simplifies because most of the pairwise averages vanish under the Gaussian average. Again, as in the case of g = 0 (cf. App. B), the only non-vanishing Gaussian correlators are those of the form ϕ ϕ , ψ ψ as given by Eqs. (D3), (D4). Thus, the sum over averages in Eq. (E2) simplifies into a sum over integrals over products of these two Gaussian propagators. The integrals run over n internal spatial variables z 1 , ..., z n which stem from the expansion in γ n , and four internal spatial variables y 1 , ..., y 4 which stem from the expansion in g 2 (and thus are related to the non-Markovian correction). Analagously, the integration also runs over n frequencies ω 1 , ..., ω n stemming from the γ n -expansion and four frequenciesω 1 , ...,ω 4 from the g 2 -expansion. The integral over the internal space variables z 1 , ..., z n , simplifies significantly, since the corresponding ψ ψ propagators (cf. Eq. (D4)) are proportional to spatial δ-functions. This results in all internal space coordinates to identify as z 1 = ... = z n = x.
The integral over the internal spatial coordinates y 1 , ..., y 4 appearing in the four ϕ and ϕ fields in Eq. (E2), however, is less trivial as it involves the correlators ("propagators") of ϕ ϕ . Likewise, the integration overω 1 , ...,ω 4 , further involves the correlator of the driving noise, C 2 (ω 3 +ω 4 ). Hence, the integrals stemming from the g 2 -expansion need to be dealt with more carefully: The four fields ϕ(y 1 ,ω 1 )∂ y2 ϕ(y 2 ,ω 2 ) ϕ(y 3 ,ω 4 )∂ y4 ϕ(y 4 ,ω 4 ) appearing in every term of Eq. (E2) have to be paired up (according to Eq. (E3)) with another creation field ϕ(z i ) or annihilation field ϕ(z i ) appearing in Eq. (E2), respectively, in order to have a non-vanishing contribution (cf. Eq. (D3)). Up to permutation of indices, each of these possible combinations (which we refer to as "Wick pairings", Eq. (E3)) can fundamentally be grouped into four different ways that are best understood diagrammatically: In any case, each of the two ϕ(y i )∂ yj ϕ(y j ) attaches via Wick-pairing to a ϕ(z k ) ϕ(z ) field appearing in the terms of Eq. (E2). Concurrently, ϕ(z k ), ϕ(z ), Wickpair with two other corresponding fields giving rise to connected ϕ ϕ -"propagators" and so on. Thereby, the Wick-pairing of ϕ(y i )∂ yj ϕ(y j ), diagrammatically speaking, splits an existing propagator in the γ n -expansion into two or, as shown below, in three.