{'fraction_non_alphanumeric': 0.05313923900822534, 'fraction_numerical': 0.016608477819200346, 'mean_word_length': 4.162203441604724, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interlacing Relaxation and First-Passage Phenomena in Reversible Discrete and Continuous Space Markovian Dynamics

We uncover a duality between relaxation and first passage processes in ergodic reversible Markovian dynamics in both discrete and continuous state-space. The duality exists in the form of a spectral interlacing -- the respective time scales of relaxation and first passage are shown to interlace. Our canonical theory allows for the first time to determine the full first passage time distribution analytically from the simpler relaxation eigenspectrum. The duality is derived and proven rigorously for both discrete state Markov processes in arbitrary dimension and effectively one-dimensional diffusion processes, whereas we also discuss extensions to more complex scenarios. We apply our theory to a simple discrete-state protein folding model and to the Ornstein-Uhlenbeck process, for which we obtain the exact first passage time distribution analytically in terms of a Newton series of determinants of almost triangular matrices.


Introduction
In his seminal work [1] Kramers analyzed the kinetics of chemical reactions in terms of diffusive barrier crossing, assuming that the kinetic rate of a chemical reaction corresponds to the inverse of the mean first crossing time. Ever since, first passage theory is at the heart of theoretical descriptions of kinetics of chemical reactions [2][3][4][5][6][7]; see e.g. [8][9][10][11] for comprehensive reviews.
First passage times play an important role in quantifying persistence properties in non-equilibrium interacting many-body systems [32][33][34]. More recent applications of first passage concepts also include stochastic thermodynamics [35][36][37], in particular, fluctuation relations for stopping time statistics and stochastic entropy production in driven molecular systems [38] and in stochastic resetting processes [39,40], as well as uncertainty relations for first passage time statistics of fluctuating currents [41,42] (see also [43]).
Moreover, our current understanding of the speed and precision of transcription regulation in biological cells, and in particular of the role of the so-called proximity effect in the co-regulation of genes, [44,45] builds on first passage time ideas. The corresponding physical principles underlying these proximity effects were explained in [46][47][48]. Notably, these works revealed the inherent insufficiency of the mean first passage time and traditional rate-based concepts for a quantitative description of biophysical dynamics in the so-called few encounter limit [48]. As a result, a quantitative understanding of phenomena such as gene regulation [44][45][46][47][48][49][50] and the misfolding-triggered pathological aggregation of proteins [51][52][53][54][55][56], which are discussed in more detail in a related study [57], requires the consideration of the full statistics of first passage time.
Existing studies of the full first passage statistics in physical systems typically focus on systems with continuous state-space dynamics, whereas much less emphasis is put on discrete-space dynamics [58]. Recent investigations of such discrete-state dynamics include, for example, simple models of enzyme kinetics [59][60][61] and novel numerical approximation schemes for studying first-passage statistics based on Bayesian inference [62] (see also [63] for a recent review).
Complementary to first passage processes are relaxation dynamics, which by contrast do not terminate upon reaching a given threshold for the first time. Relaxation phenomena in reversible diffusive dynamical systems are nowadays well understood in terms of the eigenmodes and eigenvalues of the underlying Fokker-Planck operators, which provide a generic and very intuitive understanding of the dynamics of complex stochastic systems [64][65][66]. Conversely, despite for allowing an analogous spectral representation, a similar intuitive understanding of the full first passage statistics and its physical implications remains elusive. Notwithstanding, an important approximate link between the mean first passage time for escaping the deepest potential basin and the corresponding slowest relaxation mode in the potential was established in the seminal works of Matkovsky and Schuss [67,68], which has ever since been used routinely in explaining relaxation phenomena in condensed matter systems. Nevertheless, a deeper and more generic connection between the two paradigms to date was not yet established.
Here, we present the complete duality between relaxation and first passage phenomena, which holds for all ergodic Markov processes obeying detailed balance in both, continuous and discrete state-space, in which the absorbing target is effectively one-dimensional. The duality emerges in the form of a spectral interlacing, which we prove rigorously by combining spectral-theoretic, matrix-algebraic and Greens functiontheoretic concepts. On the one hand the duality allows for an intuitive generic understanding of first passage phenomena in terms of relaxation eigenmodes. On the other hand, it enables us to determine the full first passage time statistics exactly from the corresponding relaxation eigensystem. The formalism is exact and holds for all reversible Markovian systems governed by a master equation in arbitrary dimensions or by a Fokker-Planck equation, and therefore unifies the theoretical treatment of discrete and continuous space phenomena. We note the spectral interlacing in the case of a discrete state dynamics has also recently be deduced from a 'lumping' of the state dynamics [69].
To illustrate the predictive power of the formalism in practice, we here predominantly focus on systems with discrete state-space dynamics, whereas continuous space dynamics are treated in more detail in a related study [57]. In particular, we here apply our theory to a simple discrete-state protein folding model and to diffusion in a harmonic potential, also know as the Ornstein-Uhlenbeck process. Notably, we obtain, to the best of our knowledge, for the first time an exact analytical solution for the full first passage time distribution of the Ornstein-Uhlenbeck process in the time domain.
The paper is organized as follows. In Sec. 2 we present a canonical formulation the first passage problem applicable to both discrete states-pace and continuous Fokker-Planck dynamics. Sections 3 and 4 provide a step-by-step explanation of how one can exactly determine the first-passage distribution from the corresponding relaxation process, and also contain rigorous proofs of the duality in discrete and continuous statespace dynamics, respectively. We apply the duality framework in Sec. 5 to determine the first passage statistics for a simple protein folding model and for the Ornstein-Uhlenbeck process. A concluding perspective is provided in Sec. 6. In Appendix A we derive a compact representation of the long-time asymptotics of the first passage time distribution, which inter alia extends our results for the long time asymptotics from equilibrium systems to irreversibly driven systems.

Relaxation and first passage
We assume that the probability density to find the system in state x at time t upon evolving from an initial state x 0 according to microscopically reversible dynamics, P (x, t|x 0 ), is governed by where L is a linear reversible operator, which will be specified below. We consider two classes of operators: (DS) discrete state Markov jump process, where x and x 0 assume only a finite number of states, and (FP) continuous Markovian diffusion governed by a Fokker-Planck equation.
For discrete Markov state models of class (DS) the dynamics is governed by is the total rate of leaving state x guaranteeing conservation of probability ( x ∂ t P (x, t|x 0 ) = 0). In order to have reversible dynamics we need to additionally impose detailed balance, i.e. the constraint L(x, x )/L(x , x) = exp[βU (x ) − βU (x)] (see., e.g., [70]), which assures that the system will relax to a Boltzmann distribution in a potential U (x) on ergodic timescales P eq (x) ∝ e −βU (x) , where β = 1/k B T is the inverse thermal energy. We call such a reversible ergodic process that conserves probability a relaxation process. If we add an absorbing point at x = a we call the resulting process a first passage process or in short absorption, which we introduce in the following way. First, we modify the generator (L → L a ) such that all transitions corresponding to jumps out of the absorbing state a are removed, i.e., the elements of the first passage generator read Using a bra-ket matrix notation [71] we rewrite this equation as where |a = (0, . . . , 0, 1, 0, . . . , 0) = (|a ) is a vector with all entries except the ath one; consequently, we identify L a (x, x ) = x|L a |x . The first passage time density to reach state a at time t starting from x 0 is then formally defined by which is nothing but the normalized probability flux into state a with ∞ 0 ℘ a (t|x 0 )dt = 1. Note that with Eq. (4) we use the convention that |a is the unique stationary solution with L a |a = 0.
For a continuous space Markovian diffusion the transition probability density function (the 'propagator') instead obeys the Fokker-Planck equation (2) LP where is a force field generated by the potential U (x) at position x, and β is the inverse temperature, which we set to β ≡ 1 to express energies in units of k B T from now on. The scenario with reflecting barriers at x = b ± with j(b ± , t|x 0 ) = 0 [9], we term a relaxation process, where b ± = ±∞ correspond to so-called natural boundary conditions ‡.
Conversely, an absorbing boundary at x = a enters the Fokker-Planck equation via the Dirichlet boundary condition P (a, t|x 0 ) = 0, without altering the partial differential equation (6), i.e., the first passage operator still reads L a = ∂ x D[U (x) + ∂ x ]. However, here the first passage time density becomes the probability flux into state a. For ‡ For natural boundary conditions the current and the probability density both vanish, i.e., convenience we use the operator L a as shorthand for Eq. (6) under the boundary condition P (a, t|x 0 ) = 0.
We note that without an absorbing point both, dynamics governed by the master equation (2) and the Fokker-Planck equation (6) relax to the Boltzmann distribution P eq (x) ∝ exp[−U (x)], whereas with the absorbing boundary condition the particle will eventually reach the target with probability 1.

Eigendecomposition
Since L is assumed to generate a reversible Markov process, we can expand the generator L in a bi-orthogonal eigenbasis [72]. Denoting the eigenvalues of the relaxation process by λ k and the corresponding left (right) eigenvectors by ψ L k | (|ψ R k ), respectively, the generators from Eqs. (2) and (6) become in the respective eigenbases where λ 0 = 0 ≤ λ 1 ≤ . . ., and ψ L k |ψ R l = δ kl . We assume the eigenvalues to be ordered such that λ k ≤ λ k+1 , and the generator to be irreducible λ 0 = 0 < λ 1 , which means that there is a unique equilibrium state [70]. Note that for a Fokker-Planck equation with reflecting barriers at . The generator with the absorbing point at state a, can similarly be expanded in a bi-orthogonal set of eigenfunctions where µ k is the k-th eigenvalue and φ L k | (|φ R k ) denote the corresponding left (right) eigenfunctions of the first passage process. Without loss generality we use an ordered labeling such that µ k ≤ µ k+1 , where 0 < µ 1 .
The left and right eigenvectors of the absorption (at position x = a) as well as of the relaxation process are related via x|φ L k ∝ e U (x) x|φ R k and x|ψ L k ∝ e U (x) x|ψ R k , respectively. In the case of a discrete number of states, the lowest eigenvalue of the generator (4) will be µ 0 = 0 with the right eigenfunction |φ R 0 = |a , whereas for Fokker-Planck dynamics one imposes the boundary condition a|φ R k = 0. In a previous work an explicit Newton series expression for µ 1 in terms of a series of almost triangular matrices was derived [48], which corresponds to a large deviation limit t → ∞. One of our main goals here is to obtain the full first passage statistics ℘ a (t|x 0 ) explicitly in terms of relaxation eigenmodes. Our theory builds on the renewal theorem, which we briefly review in the following subsection.

Renewal theorem
The classical renewal theorem provides a well known implicit connection between first passage and relaxation processes. It relates the probability density of the freely propagating system to be in state x at time t upon starting from a state x 0 , to the first passage distribution ℘ a (t|x 0 ) from x 0 to a: where both P (x, t|x 0 ) and ℘ a (t|x 0 ) admit a spectral representation and respectively. In other words, a system starting from state x 0 must pass through state a before reaching the final state x, which for an effectively 1-dimensional Fokker-Planck necessarily means x 0 < a ≤ x or x 0 > a ≥ x. In Eq. (11) we introduced in the first passage weights for discrete state (DS) and Fokker-Planck (FP) dynamics, respectively, which must satisfy k w k = 1 with the first nonzero weight being strictly positive w 1 (x 0 ) > 0, and where we introduced σ ± ≡ sign(a − x 0 ). Note that the first line of Eq. (12), i.e. the DS case, is equivalent to w k (x 0 ) = x =a x|φ R φ L |x 0 with L a from Eq. (4). In the case of FP dynamics the second line of Eq. (12) is equivalent to w k (x 0 ) ≡ x|φ R φ L |x 0 dx, which follows from a partial integration using both Eq. (6) and Eq. (8).
In the case of x = a the renewal theorem (9) has the simple interpretation: a system being in state a at time t must have arrived at that point at some earlier time τ for the first time (τ ≤ t), and then returned to the same position again at time t, where τ = t corresponds to the time of first arrival.
Laplace transforming the renewal theorem (9), where a generic function f is transformed according tof (s) ≡ e −st f (t)dt, yields [73] Based on this well known renewal theorem we construct in the following section a method that allows to determine explicitly the first passage time statistics ℘ a (t|x 0 ) exactly in terms of the relaxation process, i.e., we render Eq. (13) explicit in the time domain.

Principal result for discrete state systems
Starting from the renewal theorem (13), we now derive an expression for the first passage time density for discrete state Markov processes in terms of relaxation modes in the following three steps. The first step involves a crucial relation between the eigenvalues of the relaxation process λ k and absorption process µ k , which are here shown to interlace for k = 1, . . . , M . For effectively one dimensional finite lattice models with the target at an outer edge these inequalities become strict which will also apply identically to Fokker-Planck dynamics discussed in Sec. 4 in which case we formally assume M = ∞. In the second step we exactly express the first passage eigenvalues µ k in the form of a Newton series of determinants of almost triangular matrices, which generalizes the result for the slowest mode µ 1 from [48] to all first passage modes. The third and final step corresponds to a straightforward application of the residue theorem, which is used to determine the first passage weights w k (x 0 ).

Interlacing of eigenmodes (step 1)
For a discrete system with M + 1 states the eigenvalues λ k and µ k correspond to the roots of the respective characteristic polynomials i.e., χ(−λ k ) = 0 and χ a (−µ k ) = 0. Inserting Eq. (4), which is L a = L − L|a a|, into the second characteristic polynomial (16) and using the matrix determinant lemma establishes a link between the two characteristic polynomials where adj(A) is called the adjugate of a matrix A satisfying Cramer's rule A adj(A) = det(A)1. We note that the same mathematical concepts have been used recently to determine the stalling distribution of irreversibly driven systems (cf. deletioncontraction formula in [74,75]). The adjugate of a diagonal matrix D with elements which inserted into Eq. (17) gives . . . model: where we used the eigenvalue equation constitutes an essential step in our calculations, which allows us to express the diagonal of the relaxation propagatorP (a, s|a) solely in terms of eigenvalues µ k and λ k (see the following subsection for more details).
Moreover, the characteristic polynomials of the first passage process χ a (s) and relaxation process χ(s) change sign one after the other, since detailed balance imposes a|ψ R i ψ L i |a ≥ 0 for all i = 0, . . . , M , which proves that the eigenvalues of the first passage process µ k and eigenvalues of the relaxation process λ k interlace according to Eq. (14). We note that this result is directly related to the interlacing of eigenvalues generated from a "lumping" of states which is proven in [69]. In the following paragraph we briefly discuss the scenario, in which the interlacing of eigenvalues becomes strict (15), which will be the case for systems with Fokker-Planck dynamics discussed in Sec. 4.2.
The stronger condition (15) holds if all eigenfunctions are nonzero at the target | a|ψ R k | > 0 and all relaxation eigenvalues are non-degenerate, that is, λ i−1 < λ i for all i = 1, . . . , M . One can show that this condition is always trivially satisfied for 1dimensional models (L(x, x ) = 0 if |x − x | > 1), in which the target a is placed at the border (e.g., a = M or a = 0); see inset of Fig. 1a for such an exemplary 3-state system.
Inserting the relaxation eigenvalues s = −λ k into the characteristic polynomial of the first passage process (19) yields where we used the relations λ i − λ k < 0 for all i < k and λ i − λ k > 0 for all i > k, as well as χ(−λ k ) = 0. Since for k ≥ 1 each eigenvalue is positive (λ k > 0), the characteristic polynomial of the first passage process χ a (−λ k ) is equal to (−1) k−1 multiplied by a positive constant. Consequently, χ a changes sign exactly once between any two consecutive relaxation modes −λ k < s < −λ k−1 . The fact that χ a and χ are polynomials of the same degree M + 1 forbids more than a single root, and hence implies the strict interlacing of eigenvalues from Eq. (15), which completes the proof. The aforementioned reasoning is illustrated in Fig. 1a for a simple three state model in which the vertical arrows represent Eq. (20). For fine-tuned systems in which the target is not located at the very outer position (see e.g., Fig. 1b) or systems that are not effectively one dimensional (see e.g., Fig. 1c) the strict interlacing theorem (15) can be violated, whereas the "slightly weaker" interlacing condition (14) still holds.

Diagonal of the relaxation propagator in terms of bare eigenvalues
Using the results from the previous subsection we are now in the position to represent P (a, s|a) (i.e. Eq. (10) with x 0 = x = a), using only the eigenvalues of both the first passage process and the relaxation process, µ k and λ k , respectively. Laplace transforming the eigenmode expansion in Eq. (10) assuming where we identified the equilibrium probability density a|ψ R 0 ψ L 0 |a = P eq (a) in the first term. Comparing nowP (a, s|a) in (21) with χ a from Eq. (19) and χ from Eq. (16) yields after some algebrã The second equality in Eq. (22) follows from Eq. (16). Hence,P (a, s|a) encodes the eigenvalues of both, the relaxation and the first passage processes. Due to Eq. (21) P (a, s|a) contains only simple poles and decreases monotonically in s between any two consecutive poles since a|ψ R , 1d models with the target at at the border) each root λ k ofP (a, s|a) represents a first passage eigenvalue s = −µ k (k ≥ 1), which is located in between two relaxation modes λ k−1 < µ k < λ k , thus providing an alternative proof of relation (15) [76]. In the following section we determine the roots of the diagonal of the propagator explicitly, which due to Eq. (22) correspond to first passage eigenvalues µ k .
Let us briefly reformulateP (a, s|a) in a way that can also be applied to continuous systems with an infinite number of states. Isolating the equilibrium probability, which is the first term in Eq. (21), from the product formula (22) yields Since µ k , λ k increase monotonically with k we will later be able to adopt these results to systems governed by Fokker-Planck dynamics, which formally corresponds to the limit M → ∞ for which the product in Eq. (23) still converges.
3.3. From the relaxation spectrum to the first passage time spectrum (steps 2 and 3) Based on the interlacing theorem presented in Eq. (14), which is also given in Eq. (4) in a related work [57], we can determine the full first passage time spectrum {µ k , w k (x 0 )} from the corresponding relaxation spectrum, {λ k , |ψ R k , ψ L k |}. For simplicity we first consider the eigenvalues to be both ordered λ k < λ k+1 and non-degenerate, and also assume that a|ψ R k ψ L k |a > 0 holds for all values of k. The extension to situations with a|ψ R k ψ L k |a = 0, which also includes degenerate eigenvalues, for some k is straightforward and will be dealt with at the end of this subsection.
Before determining the weights w k , we first determine the first passage eigenvalues µ k , which were shown to be encoded in the roots ofP (a, s|a) in Eq. (22). We introduce the k * th "modified diagonal of the propagator" which still encodes all of the first passage eigenvalues {µ k } according to Eq. (22), i.e., it has exactly the same roots asP (a, s|a). However, in contrast toP (a, s|a) the modified function F k * (s) is strictly concave within the interval −λ k * +1 < s < −λ k * −1 , which can easily be confirmed by taking the second derivative and realizing thaẗ F k * (s) ≡ ∂ 2 s F k * (s) < 0 holds within the region of interest −λ k * +1 < s < −λ k * −1 . For k * = k and k * = k−1 the modified functions F k (s) and F k−1 (s) both are strictly concave within the interval −λ k < s < −λ k−1 and, consequently, also locally concave around the kth first passage eigenvalue The method we present in the following is an analytical technique based on the principles of Newton iteration, which is a simple root finding algorithm that is guaranteed to work for functions that are both negative and concave between the starting point and the first root. Hence, to determine the kth eigenvalue we accordingly choose the modified function such that which guarantees both negativity f (s, k) ≤ 0 and concavity ∂ 2 s f (s, k) ≤ 0 between s = −μ k and s = −µ k .
According to the interlacing theorem (15) s = −µ k is the only zero f (−µ k , k) = 0 within the interval −λ k < s < −λ k−1 . With the midpoint starting conditionμ k = (λ k + λ k−1 )/2 the kth first passage eigenvalue can be represented exactly in a series of determinants of almost triangular matrices where f n (k) is the nth derivative of f (s, k) as defined in (25) with respect to s at s = −μ k , and A n (k) stands for an almost triangular matrix with elements [48] with Θ(l) denoting the Heaviside step function (Θ(l) = 1 if l ≥ 0) and i, j = 1, 2, . . . , n − 1. Moreover, we adopt the convention det A 1 (k) = 1. We note that this method generalizes the method recently derived to determine the slowest first passage mode µ 1 [48] to all first passage eigenmodes µ k . Let us briefly repeat the two crucial steps towards Eq. (27). First, the interlacing theorem (15) guarantees that the Taylor series f (s, k) = i f i (k)(s +μ k ) i around the midpointμ k = (λ k + λ k−1 )/2 converges in the entire spectral interval −λ k < s < −λ k−1 , which also contains the first passage eigenvalue s = −µ k . Second, due to F k * (s) in Eqs. (24)-(26) the function f (s, k) is strictly concave and negative between s = −μ k and s = −µ k , which in turn guarantees the convergence of the explicit Newton series (27).
Eqs. (24)-(26) provide a universal method for determining explicitly first passage eigenvalues from the corresponding relaxation spectrum and constitute the central result of this work. We show in the Appendix A a simpler derivation of µ 1 as well as a compact approximation of the principal first passage eigenvalue µ 1 , which is particularly useful in the case of time scale separation µ 1 λ 1 (or λ 1 λ 2 ). Furthermore, Appendix A provides a generalization of the long time asymptotics from systems with reversible dynamics to irreversibly driven systems.
In the following we briefly comment on the practical implementation of the exact result for µ k to render Eqs. (24)-(27) fully explicit. The weights w k will be determined afterwards in this subsection. The nth derivative of F k * (s) with respect to s at s = −μ k , f n (k) ≡ ∂ n s f (s, k)| s=−μ k , can be written explicitly as where k * = k or k * = k − 1 is chosen according to Eq. (26). Note that condition (26) is equivalent to the condition f 0 (k) ≤ 0, implying the first line of Eq. (29) to be either negative for k * = k or for k * = k − 1, i.e. one has to evaluate the first line of Eq. (29) for k * = k: if f 0 (k) > 0 one must to change k * to k * = k − 1 and reevaluate f 0 (k). Once one has determined k * (k) and f 0 (k) one can proceed with the second line of Eq. (29) to determine f n (k)/n! and insert the result in the almost triangular matrix (28). The determinant of almost triangular matrices can be calculated elegantly using the simple recursion relation from [77], see also [78] for an efficient numerical implementation.
Having obtained the first passage eigenvalues, the weights of the first passage time distribution can be calculated using the standard residue theorem. The Laplace transform of the spectral expansion of the first passage time density (11) reads Using the residue theorem to invert the Laplace transformed renewal theorem (13) yields whereṖ (a, s|a) = ∂ sP (a, s|a) is taken at s = −µ k . The explicit Newton series (27) along with the first passage weights (31) fully characterize the first passage time distribution ℘ a (t|x 0 ) = k w k (x 0 )µ k e −µ k t in terms of relaxation eigenmodes {λ k , ψ R k }. This completes our third and final step, which allows, for the first time, to analytically deduce first passage time statistics directly from relaxation eigenmodes. We call this relation the explicit forward duality between first passage and relaxation. This completes the central result of this paper.
The spectral representation is very useful for determining the moments of the first passage time, t n ≡ t n ℘ a (t|x 0 )dt = n! k w k (x 0 )µ −n k . Moreover, as explained in more detail in a related work [57], the full spectral expansion is required for a correct explanation of kinetics in the so-called few encounter limit, where N molecules starting from position x 0 are searching for the target at a. The probability density that the first molecule out of N arrives at time t at a for the first time for this case becomes ℘ (N ) , which can be understood as follows. The probability that the first N − 1 molecule have not yet reached the target will be given by [ ∞ t ℘ a (τ |x 0 )dτ ] N −1 , while the N th particle arrives at a with a rate ℘ a (t|x 0 ); hence the probability density that any particle out of N molecules arrives at the target for the first time according to ℘ (N ) a (t|x 0 ). Further details of the N -particle problem and in particular the physical implications of the few-encounter limit are discussed in a related study [57].
Let us now briefly generalize the method to systems with degenerate eigenvalues or vanishing relaxation modes. An eigenfunction that vanishes at the target a|ψ R k = 0 will have a vanishing spectral weight as a result of Eq. (31). Hence, 'manually' removing such modes will not affect the first passage time distribution ℘ a . Moreover, if a relaxation eigenvalue λ k is degenerate we define and replace a|ψ R , a) and take the sums in Eq. (13) over all different values of λ k . After renumbering all distinct contributing eigenvalues we obtain a strict interlacing (15). Therefore, we can apply our standard forward duality also to degenerate eigensystems. In the next subsection we will briefly derive a formal backward duality after which we reformulate the results from this subsection to continuous Fokker-Planck dynamics.

Backward duality
In contrast to the explicit forward duality, which was presented in the previous subsection, an explicit reverse relation in the time-domain could not be established. In Laplace space, however, the forward duality can be inverted to give a backward duality as follows. Inserting the first passage generator L a = L − L|a a| from (4) into the Laplace transform of the propagatorP (a, s|x 0 ) = a|(1s − L) −1 |x 0 and using the Sherman-Morrison-Woodbury formula yields Let us now insert the expression for the first passage time distribution from Eq. (5), which can be written as ℘ a (s|x 0 ) = s a|(1s − L a ) −1 |x 0 , into Eq. (33) to obtaiñ where L(x, a) = x|L|a is the generator of the relaxation process. Notably, this is expression corresponds to the backward duality and is the formal inverse of the renewal theorem, where℘ a (s|x 0 )/s is the Laplace transform of the cumulative first passage time distribution t 0 ℘ a (τ |x 0 )dτ .

Greens function with natural boundaries
We restrict our discussion to effectively 1-dimensional dynamics, which include diffusion in d dimensions in an isotropic potential as discussed in [48], where d may also be fractal.
Introducing an absorbing target at position a splits the first passage problem into two cases (I) x 0 < a and (II) x 0 > a. Case (I) corresponds to an absorption from the left, and case (II) to an absorption from the right. In the following paragraph we demonstrate that all first passage modes µ k of both distinct cases (I) and (II) are entirely encoded inP (a, s|a), which allows to formulate the results from Sec. 3.3 also for systems with Fokker-Planck dynamics.
Laplace transforming the Fokker-Planck equation (6) yields where L = −∂ x D[βU (x) + ∂ x ] and x 0 is the initial position of the relaxation process. Eq. (35) is a inhomogeneous linear differential equation which can be solved using the standard Green's function approach. First, we find the two independent solutions v ± (x, s) of the homogeneous problem (L−s)v ± (x, s) = 0, where we use the label "−" and "+" for the solution satisfying the left and right boundary condition, respectively. That is, a diffusion process within an interval b − < x < b + imposes the probability current s) to vanish at the boundaries, i.e. j ± (b ± , s) = 0. The special case of so-called natural boundary conditions correspond to the limit lim x→±∞ v ± (x, s) = 0 or analogously lim x→±∞ j ± (x, s) = 0, that is, b ± = ±∞. The full solutionP (x, s|x 0 ) of (35) is a continuous function in x with a discontinuity of its first derivative at x = x 0 . Using the scaled Wronskian § the propagator, which satisfies the proper jump condition of the first derivative (current function) at x = x 0 , becomes We note that the Wronskian (36) is proportional to the Boltzmann factor (see, e.g., Ref. [76]), i.e., W s (x) = W s (x 0 ) exp[βU (x 0 ) − βU (x)]. Hence using the renewal theorem (13) andP (a, s|a) as well asP (a, s|x 0 ) from Eq. (37) yields the Laplace transform of the first passage time distribution The two independent functions v ± (x, s) are entire functions without any poles in s [79], and in turn encode in their roots all first passage eigenvalues s = −µ k . In particular v − encodes all first passage modes from case (I) x 0 < a, and v + encodes all first passage modes from case (II), in which the particle is absorbed from the right x 0 > a. Due to Eq. (37) the zeros ofP (a, s|a) at s = −µ k determine the first passage spectrum. Hence, § For convenience we defined with the scaled Wronskian with the current function j ± (x, s) instead of the first derviative ∂ x v ± (x, s), i.e., W s (x)/D would represent the standard definition of the Wronskian.
all results from Sec. 3.3 hold identically for continuous systems as well. However, the sums are here not finite, i.e., M = ∞. For example,P (a, s|a) becomes where ψ R l is the lth right eigenfunction of the Fokker-Planck operator satisfying Lψ R l (x) = −λ l ψ R l (x) with the corresponding left eigenfunction ψ L l (a) ∝ e βU (a) ψ R l (a) and normalization The first passage modes µ k can then be determined with Eqs. (27)-(29), where k * (k * = k or k * = k − 1) must be chosen such that f 0 (k) < 0 holds in Eq. (29) with a|ψ R l ψ L l |a ≡ ψ L l (a)ψ R l (a). Concurrently, the first passage weights follow from Eq. (31).
The formal backward duality from Sec. 3.4, however, must be adopted as follows. After some tedious algebra we obtain formally the exact inverse duality in the form of where sign σ ± = −1 if x 0 < x and σ ± = +1 if x 0 > x; Eq. (40) can easily be verified by inserting the Wronskian (36) and the first passage time distribution (38) into the right hand side of Eq. (40), and comparing the result with the propagator from Eq. (37). Notably, this inverse duality is the continuous version of Eq. (34).

Relaxation under reflecting boundary conditions and strict spectral interlacing
In the previous subsection the target a divided the phase space into two regions, which implies that the first passage modes for the cases (I) and (II) separate into "left" and "right" modes as well. For example, if x 1 < a and x 2 > a one of the first passage weights w k (x 1 ) or w k (x 2 ) must typically be zero for all values of k. If one uses just the first M modes to approximate the propagator [cf. Eq. (39)] the zeros of the right hand side of become approximations of the first passage modes and, hence, the weights w M k (x 1 ) and ) for finite M , i.e., modes from case (I) and (II) mix. Such a mixing can be avoided entirely if the relaxation process is analyzed with a reflecting boundary at the target position a (b + = a or b − = a).
The result for reflecting boundary conditions j + (a) = 0 (case (I)) and j − (a) = 0 (case (II)) is automatically obtained by the following replacement: respectively, which inserted into the scaled Wronskian (36) Utilizing the Wronskian for the reflecting boundary condition in Eq. (37) yields the diagonal of the propagator in the form For two linearly independent functions v ± with nonzero Wronskian (36) the zeros of v ± and j ± are different. Hence, the zeros of v − in s contain only the first passage modes for the case (I) x 0 < a, whereas the zeros of v − do not contain zeros of first passage modes corresponding to the case (II) x 0 > a. Let us from now on just focus on a case (I), in which x 0 < a, since case (II) follows by analogy. For case (I) we consider the Fokker-Planck operator L from Eq. (6) with zero current condition at x = a and natural boundary condition for x → −∞.
To that end we first determine the relaxation eigenvalues λ k and eigenmodes ψ R k , ψ L k . Note that we consider the eigensystem in the presence of a reflecting wall. Without loss of generality we here explicitly treat only the "left" problem −∞ < x ≤ a, since the opposite "right" problem (denoted later on with †) follows by analogy. As before we have the normalization } we now determine the first passage eigenvalues µ k as explained in the previous subsection. The resulting first passage eigenvalues µ k will automatically contain only first passage modes corresponding to the "left" problem. This procedure remarkably simplifies the numerical determination of the first passage distribution, especially of those modes that are faster than the slowest mode of the "right" problem (i.e. absorption from the right), µ k > µ † 1 , since a small number of modes M in Eq. (41) might otherwise be confused with 'fantom' modes from the opposite case †. In a related work [57] we investigated the "left" first passage problem (case (I)) for a triple well potential in the presence of a reflecting boundary, and found an excellent agreement between the analytical first passage time distribution and computer simulations extending over many orders of magnitude in time using merely M = 40 relaxation modes. Finally, we have to point out that solving an eigenvalue problem {λ k , ψ R k } with reflecting boundary condition is numerically easier than without reflecting boundary, i.e., natural boundaries albeit theoretically easier are numerically harder.
In the following section we apply these theoretical results to a discrete-state protein folding model and for the Ornstein-Uhlenbeck process.

Discrete protein folding model
We consider a simple continuous-time Markov state model for a protein with three structural elements as shown in Fig. 2. The protein starts from an initially unfolded state x 0 = (0, 0, 0) ≡ I, from which it is searching for the native state a = (1, 1, 1) through intermediate states II-VII (see e.g., [80,81]). Each arrow in Fig. 2 indicates a possible transition x → x (x, x = I, . . . , VIII, x = x ) that occurs with a Arrhenius type rate To test the power of the method from Sec. 3.3 we set up the 8 × 8 transition matrix L with elements x |L|x = L(x , x) for a given set of energy barriers B xx and free energies F x . Then we carry out the eigendecomposition of L, for which we first determine the eigenvalues 0 < λ 1 , . . . , λ 7 (with λ 0 = 0) corresponding to the zeros of the characteristic function (16), χ(−λ k ) = 0. We then determine the right eigenvectors |ψ R k by solving L|ψ R k = −λ k |ψ R k for k = 0, . . . , 7. The corresponding left eigenvectors, which solve ψ L k |L = − ψ L k |λ k , have components ψ L k |x = N −1 e Fx x|ψ R k , where N = VIII x=I e Fx | ψ R k |x | 2 is a normalization factor. We take the function F k * as defined in Eq. (24), where Ψ k = a|ψ R k ψ L k |a andμ k = (λ k + λ k−1 )/2 with k = 0, . . . , 7, and choose k * (k) = k, k − 1 according to Eq. (26), which guarantees f (s, k) = F k * (k) (s) to be negative at s = −μ k . The truncated Newton series (27) involving the first N terms is then given by where det A n (k) is the determinant of the almost triangular matrix from Eq. (28) and f i (k) is the ith derivative of f (s, k) at s = −µ, with explicit formulas given in Eq. (29). The weights w k (x 0 ) are determined using Eq. (31), i.e. by inserting µ k → µ N k . The calculations are performed for 100 randomly generated folding landscapes chosen as described in the caption to Fig. 2. In Fig. 3 we present the results for one particular realization of the folding landscape. Fig. 3a displays the first passage time distribution for N = 6 and N = 12 on a doubly-logarithmic scale. The solid line represents the first passage time distribution obtained via a numerical diagonalization of L a . The corresponding duality solutions nicely overlap with the numerical result even on relatively short time scales (see inset for a plot with linear scales).
Having obtained the full distribution of first passage times is important for  understanding kinetics in the so-called few encounter limit [48], in which for example 100 molecules are simultaneously searching for a state a. This scenario is indeed biologically relevant, for example, in the misfolding-triggered protein aggregation, which in turn leads to numerous diseases (see [57] for a more detailed discussion). Namely, as soon as the first protein molecule spontaneously misfolds it creates a nucleation site for further downhill misfolding and aggregation events ultimately leading to a macroscopic insoluble toxic aggregate. In such a scenario the typical timescale of first arrivals will naturally be shifted towards shorter timescales, thus requiring an accurate determination of the full first passage statistics. Standard approaches focusing on the mean first passage time alone, would therefore fail in the few encounter limit, whereas our new framework provides an accurate and consistent result (see also [57] and Fig. 2 therein for more details).
In Fig. 4 we systematically analyze the deviation of the truncated Newton series (45) with respect to corresponding numerically obtained first passage eigenvalues µ k for 100 randomly generated folding landscapes. For a given landscape the relative error is quantified in terms of the dimensionless quantity ≡ |µ N k − µ k |/µ k , and Fig. 4a depicts the typical error characterized by the median of the individual errors for all seven modes, respectively. Note that for N = 20 the relative error of the finite Newton series is typically below 10 −6 . Fig. 4b displays the maximal error out of 100 randomly picked landscapes. We observe that larger errors can occur if a first passage eigenvalue is located immediately after a gap in the relaxation spectrum. The smaller error of the slowest first passage mode µ 1 is due to the fact that µ 1 cannot be located after such a gap due to the interlacing theorem (14), which implies µ 1 ≤ λ 1 . In this specific example the maximum relative error out of 100 models randomly generated models is found for the fifth mode (µ 5 ); the inset of Fig. 4b shows the correspondingP (a, s|a) as well as the result form the finite Newton series with N = 20 (see diamonds in the inset of Fig. 4b). In this extreme scenario the weight of the fifth relaxation mode Ψ 5 Ψ l =5 is almost negligible compared to other weights, leading to an almost vanishing weight w 5 , which would in turn require an increased number of terms N entering the Newton series. Increasing the number of terms in the truncated Newton series from N = 20 to N = 200 reduces the deviation from 10 −1 to 10 −2 , the result for N = 200 is marked by the star in the inset of Fig. 4b. Fig. 4 readily demonstrates that the our duality can be robustly and reliably applied to all Markov state models.

Ornstein-Uhlenbeck process
Let us now consider a linear Ornstein-Uhlenbeck, which corresponds to a diffusion process in a harmonic potential βU (x) = ωx 2 /2. The corresponding Fokker-Planck operator reads L = D∂ x ωx + D∂ 2 x . The eigendecomposition of the relaxation process in the absence of reflecting boundaries is well known. The respective eigenvalues are given by λ k = Dωk with the corresponding eigenfunctions [72] where H k is the kth Hermite polynomial. Although this process is extremely well studied, a closed-form analytical result for the first passage time distribution ℘ a (t|x 0 ) for any non-centered target position a = 0 remained elusive [82][83][84]. We note that the well known analytical solution of the Laplace transform of the probability density℘ a (s|x 0 ) in terms of Hermite polynomials [73,82,85,86] until now could only be analytically inverted to ℘ a (t|x 0 ) for the special case a = 0 [82]. Furtheremore, the exact large deviation limit ℘ a (t|x 0 ) w 1 (x 0 )e −µ 1 t was just recently derived in [48]. To obtain the full first passage time distribution we here use Eqs.
The kth first passage eigenvalue µ k is determined by using the finite Newton series (45), where the almost triangular matrix is taken from Eq. (28), and the corresponding first passage weights w k (x 0 ) are determined using the residue theorem (31). Note that our theory allows for the first time to determine analytically all first passage eigenvalues {µ k } as well as the weights {w k } and, therefore, also provides a complete solution to the first passage time density ℘ a (t|x 0 ).  by the vertical arrows in Fig. 5a. These first passage modes correspond to an absorption from the right, where the corresponding weights vanish w k (x 0 ) = 0 (here k = 4, 7, 9, . . .) completely for all x 0 ≤ a, which, however, is only obtained in the limit M → ∞.
This numerical truncation problem can be avoided completely if the relaxation process is considered with a reflecting boundary condition as explained in Sec. 4.2, which automatically removes beforehand all zeros marked by the arrows in Fig. 5a (s = −µ 4 , −µ 7 , −µ 9 , . . .). Nevertheless, to illustrate the power and robustness of our duality approach we proceed here without a reflecting wall and use M = 2000. In Fig. 5b we show the first passage time distribution on a log-log scale (see inset for a linear scale) for three different starting positions x 0 = 0, 0.5, 1 (absorbing point a = 2). The lines represent the first passage time distribution which is determined using our new method (with M = 2000 relaxation modes) and the symbols represent the results ℘ a of a numerical Laplace inversion of the renewal theorem (see figure caption for more details). We find a perfect agreement between our new analytical method (lines) and the numerical solution. For comparison, we imposed a reflecting wall at the target in a related article [57] and obtained a similarly excellent agreement between the duality solution and the simulated first passage time density using a total of M = 40 relaxation modes to quantify the first passage time statistics for diffusion in a multi-well potential. In either case, our new duality framework is exact for infinite M , and hence the desired precision can be tuned at will.

Concluding perspectives
We rigorously established a duality between the relaxation and the corresponding first passage processes in terms of an interlacing of eigenvalues. In other words, the timescales at which a particle is absorbed into the target are proven to interlace with the corresponding relaxation timescales. This duality allows us to understand first passage processes, both qualitatively and quantitatively, in terms of relaxation eigenmodes. For example, spectral gaps in the relaxation spectrum translate directly into spectral gaps in the first passage spectrum. More explicitly, in effectively one dimensional systems N gaps in the relaxation spectrum, arising from N local (free) energy basins, translate into N − 1 gaps in the first passage time spectrum corresponding to the N − 1 barriers separating the minima. Most importantly, we established a duality that allows, for the first time, to determine exactly the first passage time distribution from the corresponding relaxation spectrum.
Our theory is developed end tested on both, continuous reversible Fokker-Planck dynamics and Markov state jump processes in arbitrary dimensions. For convenience and without loss of generality, we restricted the applications of the duality for systems obeying Fokker-Planck dynamics to effectively one dimensional problems. An extension to more general models, for example, to diffusion on graphs would be straightforward, albeit rendering the calculations more cumbersome.
We tested and applied our theory to a discrete Markov state model of a simple protein folding landscape and the Ornstein-Uhlenbeck process, while a continuous analogue of a folding landscape are discussed elsewhere [57]. Notably, we have derived, to the best of our knowledge, for the first time an exact and explicit analytical expression for the first passage time distribution of the Ornstein-Uhlenbeck process.
Looking forward it will be interesting and relevant to apply the duality to the analysis of first passage processes on graphs. Applications of the duality to narrow escape problems in arbitrary dimensions [23][24][25][26][27][28] will also be carried out in future studies.
Finally, an extension of the framework to periodically or constantly driven systems (i.e., irreversible Markovian dynamics), which goes beyond the long time limit that is presented in Appendix A, will be particularly challenging. Namely, there the interlacing theorem cannot be expected to hold anymore, since both eigenvalue spectra {λ k } and {µ k } can become complex valued.

Appendix A. Explicit formula for principal eigenvalue
In this appendix we simplify Eqs. (24)-(26) in the limit of a time-scale separation and for rare-event asymptotics for the principal first passage eigenvalue µ 1 . We obtain a compact asymptotic expression of the principal first passage eigenvalueμ 1 µ 1 , which is particular accurate if the time-scale of the slowest first passage eigenvalue is well separated from the time-scale of the slowest relaxation mode (µ 1 λ 1 ), which inter alia refines a previously proposed approximate link between the mean first passage time and the slowest relaxation mode [67,68].
First, we redefine Eq. (25)  If we now set f 3 = f 4 = . . . = 0 in the almost triangular matrices (A.3), that isÃ n ≡ A n | f 3 =f 4 =...=0 , the resulting matrixÃ n becomes triangular, implying that its determinant is simply given by the product of the diagonal elements Eq. (A.6) is nothing but the root of the second order Taylor expansion of f (s) around s = 0 (i.e., the parabolic equation). This approximation is quite accurate whenever µ 1 λ 1 . If the target is located at a high energy barrier, such that slowest first passage eigenvalue is exponentially suppressed by the (free) energy at the target (i.e., µ 1 ∝ e −U (a) ), Eq. (A.6) will lead to a quite accurate approximationμ 1 , which can be seen in Fig. A1. More precisely, in Fig. A1a we depict the relative error |μ 1 − µ 1 |/µ 1 as function of the target-site energy ∆E = a 2 /2 for the Ornstein-Uhlenbeck process from Sec. 5.2 with U (x) = x 2 /2 (ω = D = β = 1). Conversely, Fig. A1b displays the results the first passage time distribution to four different target positions a = 1.5, 2, 2.5, 3 for a particle starting from x 0 = 1. The symbols represent histograms for ℘ a (t|x 0 ) deduced from 10 6 Brownian dynamics trajectories, and the lines correspond to the large deviation asymptotic ℘ a (t|x 0 ) w 1 (x 0 )μ 1 e −μ 1 t deduced from (A.6). We conclude that the limit µ 1 λ 1 lead to both, a quite accurate approximationμ 1 µ 1 and to an effectively single exponential decay µ 1 µ 2 of the first passage statistics, which extends previous results [87] (see also [88]).
Moreover, a spectral gap such as λ 1 λ 2 will also render f n≥2 from Eq. (A.2) to be negligibly small if the target is not located at the global minimum of the potential. For example, a multi-barrier crossing, as the one studied in Ref. [57] (see Fig. 5 therein), the principal first passage eigenvalue from Eq. (A.6) deviates less than two percent from the exact value µ 1 , i.e., |µ 1 −μ 1 |/µ 1 < 0.02. Notably, the approximation (A.6) refines previous conjectures that the mean first passage time to escape from the deepest potential basin corresponds to the first nonzero relaxation mode [67,68].
We note that Eq. (A.6) can be reformulated to givẽ where we inserted f 0 = P eq (a) and defined σ n = l≥1 a|ψ R l ψ L l |a /λ n l . This relation is equivalent to Eq. (17) from a related article [57].
Finally, we emphasize that Eqs. (A.4) and (A.6) apply also to homogeneous irreversible Markov processes, i.e., the relations from this appendix are not restricted to hold just for reversible Markov chains.