# First-passage on disordered intervals

James Holehouse<EMAIL_ADDRESS>The Santa Fe Institute, 1399 Hyde
Park Road, Santa Fe, NM, 87510, USA  S. Redner<EMAIL_ADDRESS>The Santa
Fe Institute, 1399 Hyde Park Road, Santa Fe, NM, 87510, USA

###### Abstract

We investigate the first-passage properties of nearest-neighbor hopping on a
finite interval with disordered hopping rates. We develop an approach that
relies on the backward equation, in conjunction with probability generating
functions, to obtain all moments, as well as the distribution of first-passage
times. Our approach is simpler than previous approaches that are based on
either the forward equation or recursive method, in which the $m^{\rm th}$
moment requires all preceding moments. For the interval with two absorbing
boundaries, we elucidate the disparity in the first-passage times between
different realizations of the hopping rates and also unexpectedly find that
the distribution of first-passage times can be _bimodal_ for certain
realizations of the hopping rates.

First-passage problems concern the distribution of times for a stochastic
process to first reach a defined state [1, 2, 3, 4, 5]. A paradigmatic setting
for first-passage phenomena is the nearest-neighbor random walk on the finite
interval. Here, it is of fundamental importance to understand the statistics
of the time required for the random walk to _first_ reach the boundaries of
the interval. The applications of this type of first-passage problem are vast,
including biological processes, such as the Moran model [6, 7], migration
phenomena [8, 9], the behavioral dynamics of ant recruitment [10, 11], as well
as the dynamics of many types of financial instruments [12, 13].

First-passage phenomena are much richer and less well understood when the
hopping rates of the random walk are _spatially disordered_. While much
progress has been made in determining the average first-passage time and its
low-order moments [2, 4, 14, 15, 16, 17, 18, 19, 20], understanding the
properties of the full _distribution_ of first-passage times is still
incomplete. Here we tackle this problem by focusing on the generating function
for the first-passage probability. The finite interval provides a particularly
instructive platform to investigate many basic disorder-controlled physical
phenomena, such as the diffusion of a particle in a random potential [21, 19,
20, 22, 23, 24, 25], DNA translocation through a nanopore [26], channel
transport [27], and percolation [28, 29]. Our formalism allows us to readily
evaluate first-passage probabilities and their moments for individual
realizations of disorder and to uncover some unexpected disorder-controlled
phenomena.

The key to our solution lies in first writing the backward equation for the
moment generating function of the first-passage probability, and then
reformulating the resulting recursion as a linear algebra problem. This
approach leads to analytic expressions for the generating function in terms of
the elements of the inverse of a tridiagonal matrix [30]. Related methods have
previously provided a semi-analytic solution to the probability distribution
from one-step master equation in time [31], and to solve the three-term
recurrence for the moment generating function [32]. Other pertinent studies
derive semi-analytic solutions for the first-passage time distribution, again
for arbitrary hopping rates [33, 34, 35]. These studies start with the Laplace
transform of the formal solution to the master equation, and obtain their
result in terms of undetermined eigenvalues of the master operator.

However, these past investigations all focused on the forward master equation,
whereas we use the backward equation due to its utility in the study of first-
passage problems [4, 27]. The benefits of our approach are its versatility in
elucidating first-passage properties for absorbing boundaries, reflecting
boundaries, and conditional waiting times, and its conciseness through the use
of linear algebraic results that simplify the generating function.
Additionally, our approach does not rely on finding the eigenvalues of the
master operator numerically, leading to analytic, as opposed to semi-analytic
results.

_Formalism._ The backward equation is the adjoint of the commonly used master
equation for the time evolution of the probability distribution. The backward
equation is especially useful in first-passage problems for which the final
state is prescribed and the initial state becomes the fundamental dependent
variable. Let $P_{i,j}(t)$ be the probability for a random walk to reach a
final state $j$ for the _first_ time at time $t$ when starting from state $i$
at $t=0$. Here, $i$ denotes a one-dimensional coordinate. The backward
equation for this first-passage probability is [2, 4, 14],

$\displaystyle P_{i,j}(t+\Delta t)$ $\displaystyle=b_{i}\Delta
tP_{i+1,j}(t)+d_{i}\Delta tP_{i-1,j}(t)$ $\displaystyle\hskip
36.98866pt+[1-(b_{i}+d_{i})\Delta t]P_{i,j}(t),$ (1)

where $b_{i}$ and $d_{i}$ are arbitrary rates of hopping to the right and left
from site $i$, respectively. This equation states that the probability to
first arrive at $j$ starting from $i$ at time $t+\Delta t$ equals the sum of
the hopping probabilities in a time $\Delta t$ from $i$ to either $i-1$, $i$,
or $i+1$ times the first-passage probabilities to $j$ at time $t$ from $i-1$,
$i$, or $i+1$. We set $\Delta t=1$ so that the hopping probabilities in a time
$\Delta t$ are $b_{i}\Delta t$ and $d_{i}\Delta t$, and $t$ now counts the
number of hopping events. For $\Delta t=1$, the hopping rates must satisfy the
constraints $0<b_{i},d_{i}<1$ and $b_{i}+d_{i}<1$.

Figure 1: Interval with (a) two absorbing boundaries, and (b) one reflecting
and one absorbing boundary.

We now calculate the unconditional first-passage probability for the interval
with two absorbing boundaries (Fig. 1(a)). Define $f_{i}(t)=P_{0,i}+P_{N,i}$
as the probability density function to first reach either $0$ or $N$ at time
$t$ when starting from site $i$. In close analogy with Eq. (First-passage on
disordered intervals), its governing equation is

$\displaystyle
f_{i}(t\\!+\\!1)=b_{i}f_{i+1}(t)+d_{i}f_{i-1}(t)+[1\\!-\\!(b_{i}\\!+\\!d_{i})]f_{i}(t).$
(2)

We solve this equation by the generating function technique. Multiplying (2)
by $z^{t}$, summing over $t$, and introducing the generating function
$F_{i}(z)=\sum_{t=0}^{\infty}z^{t}f_{i}(t)$, gives the following three-term
recurrence

$\displaystyle\begin{split}&b_{i}F_{i+1}(z)+d_{i}F_{i-1}(z)+\beta_{i}(z)F_{i}(z)=0,\end{split}$
(3)

for $1\leq i\leq N-1$, where for notational simplicity we introduce
$\beta_{i}(z)\equiv 1-z^{-1}-b_{i}-d_{i}$. This recursion obeys the boundary
conditions $F_{0}(z)=F_{N}(z)=1$, since $f_{0}(t)=f_{N}(t)=\delta_{t,0}$.

To solve (3), it is helpful to write it as a matrix equation [32]. We first
define the column vector $\mathbf{F}(z)=(F_{1}(z),F_{2}(z),\dots,F_{N-1}(z))$
and using $F_{0}(z)=F_{N}(z)=1$, this recursion becomes

$\displaystyle\mathbb{A}(z)\cdot\mathbf{F}(z)=-\mathbf{w},$ (4)

where $\mathbf{w}=(d_{1},0,0,\dots,0,b_{N-1})$ and $\mathbb{A}(z)$ is the
tridiagonal matrix of dimension $(N-1)\times(N-1)$

$\displaystyle\mathbb{A}(z)=\begin{pmatrix}\beta_{1}(z)&b_{1}&0&\dots&0&0\\\
d_{2}&\beta_{2}(z)&0&\dots&0&0\\\ \vdots&&&&\ddots&&\\\
0&0&0&\dots&d_{N-1}&\beta_{N-1}(z)\\\ \end{pmatrix}.$

The formal solution to Eq. (4) is,

$\displaystyle\mathbf{F}(z)=-\mathbb{A}(z)^{-1}\cdot\mathbf{w}.$ (5)

In what follows we write the $ij^{\rm th}$ elements of $\mathbb{A}(z)^{-1}$ as
$\alpha_{i,j}(z)$. Performing the matrix multiplication gives

$\displaystyle
F_{i}(z)=-[d_{1}\alpha_{i,1}(z)+b_{\bar{N}}\alpha_{i,\bar{N}}(z)]\,,$ (6)

with $\bar{N}\equiv N-1$. Our task now is to find closed-from expressions for
the elements of the inverse $\alpha_{i,j}(z)$. To this end we exploit the
tridiagonal form of $\mathbb{A}(z)$, as was done in [31] for the forward
master equation, to give the elements $\alpha_{i,j}(z)$ in terms of
computationally simple products of polynomials via Cramer’s rule [30]. For
later use, we define the following products

$\displaystyle B_{i}\equiv\prod_{k=i}^{\bar{N}}b_{k},\qquad
D_{i}\equiv\prod_{k=1}^{i}d_{k}.$ (7)

To find $F_{i}(z)$, we only require the elements $\alpha_{i,1}(z)$ and
$\alpha_{i,\bar{N}}(z)$, which are given by

$\displaystyle\alpha_{i,1}=(-1)^{i+1}\frac{D_{i}p_{i}(z)}{d_{1}p_{0}(z)},\quad\alpha_{i,\bar{N}}=(-1)^{i+\bar{N}}\frac{B_{i}q_{i-1}(z)}{b_{\bar{N}}p_{0}(z)},$
(8)

where we recursively define the polynomials

$\displaystyle\begin{split}p_{N}(z)&=1,\\\
p_{\bar{N}}(z)&=\beta_{\bar{N}}(z),\\\
p_{i}(z)&=\beta_{i+1}(z)p_{i+1}(z)-b_{i+1}d_{i+2}p_{i+2}(z),\end{split}$ (9a)
for $0\leq i\leq\bar{N}-2$, and $\displaystyle\begin{split}q_{0}(z)&=1,\\\
q_{1}(z)&=\beta_{1}(z),\\\
q_{i}(z)&=\beta_{i}(z)q_{i-1}(z)-d_{i}b_{i-1}q_{i-2}(z)\,,\end{split}$ (9b)

for $2\leq i\leq\bar{N}$. Although this calculation seems complicated, it is
much faster than matrix inversion via Cramer’s rule, since it avoids
calculating the multiple zeros encountered in the evaluation of the minors of
$\mathbb{A}(z)$ [36]. More details are given in Appendix I.

We can now write the generating functions for the first-passage probability
for each starting position $i$ as

$\displaystyle\begin{split}F_{1}(z)&=\frac{1}{p_{0}(z)}\left[(-1)^{\bar{N}}B_{1}-d_{1}p_{1}(z)\right],\\\
F_{\bar{N}}(z)&=\frac{1}{p_{0}(z)}\left[(-1)^{\bar{N}}D_{\bar{N}}-b_{\bar{N}}q_{\bar{N}-1}(z)\right],\\\
F_{i}(z)&=\frac{(-1)^{i}}{p_{0}(z)}\left[p_{i}(z)D_{i}+(-1)^{\bar{N}-1}q_{i-1}(z)B_{i}\right]\,,\end{split}$
(10)

where the last equality holds for $2\leq i\leq\bar{N}-1$. Note that
$p_{\bar{N}-i}$ and $q_{i}$ are polynomials of order $i$ in $z^{-1}$, and each
$F_{i}(z)$ is a rational function that is made up of polynomials in $z$. In
practice, for an interval of length $N$ and given set of $\\{b_{i},d_{i}\\}$,
we find the polynomials $p_{i}$ and $q_{i}$ via Eq. (9), after which we can
use Eqs. 10 to compute the generating function and its series expansion.

Our formalism can be readily extended to the _conditional_ first-passage
probability, namely, the probability to first reach a specified boundary
without ever touching the other boundary (Appendix II). We can also treat the
case of a reflecting boundary. Here, the dynamics at the reflecting end of the
interval must be treated as a special case in which some of the elements of
$\mathbb{A}(z)$ are altered, although the functional forms in Eqs. 10 remain
the same. The details of this calculation will be given elsewhere.

Application to the disordered interval: A salient feature of first-passage in
the interval is the huge disparity in the mean first-passage time (FPT)
between different realizations of the hopping rates $\\{b_{j},d_{j}\\}$. We
use Eq. (10) to give the exact FPT that represents an average over all random-
walk trajectories. We choose each $b_{j}$ from a uniform distribution on
$[0,2/3]$, so that each $d_{j}=2/3-b_{j}$. Since $b_{j}+d_{j}=2/3$, there is a
1/3 probability for the walk to remain at the same site in a single event.
This choice eliminates the even-odd oscillations that arise for the nearest-
neighbor random walk with $b_{j}+d_{j}=1$.

Figure 2: (a) Mean first-passage times on a disordered interval of length
$N=100$ for four realizations of the hopping rates as a function of starting
position $i$ (solid). The dashed curve shows the first-passage time for the
homogeneous interval with $b_{i}=d_{i}=1/3$ for all $i$. (b) The smoothed
(10-point range) local bias $\ln(b_{i}/d_{i})$ for the realization in black in
(a).

Figure 2 illustrates this heterogeneity in the FPT for several realizations of
the $\\{b_{j},d_{j}\\}$; even larger realization-specific variation occurs for
higher moments of the FPT (not shown). Moreover, the dependence of the FPT on
starting position has no resemblance to the parabolic profile that arises in
the absence of disorder. Instead, the FPT is nearly independent of starting
location in certain subintervals and changes rapidly within intervening
boundary layers. This behavior stems from the existence of local potential
wells that are induced by the disordered hopping rates. For the realization
shown in black in the figure, there is a negative bias over the first quarter
of the interval so that the FPT for starting points in this range is small.
The bias suddenly becomes positive for $i/N\gtrsim 0.3$ which causes the rapid
increase in the FPT. The local bias again changes sign as $i/N$ increases
through roughly 0.65. The effective potential well for $i/N$ between 0.3 and
0.65 keeps the FPT roughly constant in this range. The change of sign in the
bias around $i/N\approx 0.7$ is the source of another sudden increase in the
FPT at this point. Generally, wherever the local bias leads to an effective
potential well, the FPT is nearly independent of starting location within this
well and then suddenly jumps when the local bias changes sign.

It is also revealing to determine how the FPT converges to its true value upon
averaging over progressively larger numbers of realizations of the hopping
rates. For this calculation, we use the dichotomous distribution in which each
$b_{j}$ takes the values 0.3 or 0.6 equiprobably, while $d_{j}=0.9-b_{j}$.
With this choice, there is a countably finite number of hopping rate
realizations so that we can average over all random walk trajectories _and_
over all realizations of the hopping rate disorder for a given (albeit short)
interval length.

Figure 3: The average deviation from Eq. (11) when a finite fraction $f$ of
all hopping rate realizations are sampled. Interval length $N=12$.

For a given realization of the hopping rates indexed by $\alpha$, we define
the mean FPT starting from site $i$ as $\langle t_{i}^{(\alpha)}\rangle$. The
true average FPT, averaged over all $\mathcal{M}=2^{N-1}$ realizations of the
hopping rates is then

$\displaystyle\overline{\langle
t_{i}\rangle}\equiv\frac{1}{\mathcal{M}}\sum_{\alpha=1}^{\mathcal{M}}\langle
t_{i}^{(\alpha)}\rangle\,.$

We now ask: how close is the average FPT over a subset of realizations of the
hopping rates to the average over _all_ realizations of the hopping rates? To
address this question, we define the partial average in which we select a
random fraction $f=M/\mathcal{M}$ of all hopping rate realizations and compute
the exact FPT for this subset, again for random walks that start at site $i$:

$\displaystyle\overline{\langle
t_{i}\rangle}_{f}\equiv\frac{1}{M}\sum_{\alpha=1}^{M}\langle
t_{i}^{(\alpha)}\rangle\,.$

An appropriate deviation measure is the relative difference between the true
average and the partial average as a function of $f$, which we define as
$E_{i}(f)$:

$\displaystyle E_{i}(f)\equiv\frac{\overline{\langle
t_{i}\rangle}-\overline{\langle t_{i}\rangle}_{f}}{\overline{\langle
t_{i}\rangle}}\,.$ (11)

Figure 4: (a) The first-passage probability $f_{2}(t)$ on the interval
$[0,10]$. The schematic shows the local bias $\ln(b_{i}/d_{i})$ at each site
(darker color = stronger bias). The inset shows the explicit values of
$\ln(b_{i}/d_{i})$. (b) The first-passage probability $f_{6}(t)$ for an
interval of length $N=20$ in which the synthetically generated bias changes
sign at $i=5$ with $|v|=|b_{i}-d_{i}|=0.3$.

In Fig. 3(a) shows this deviation $E_{i}(f)$ versus $i$ for various ensemble
fractions $f$. We see that with 1% of all realizations there is a deviation of
roughly 25%. Since a real simulation can only sample an infinitesimal fraction
of all realizations, this plot indicates that any reasonable simulation will
provide a poor approximation to the true average FPT. In particular, having
half of all possible realizations still leads to a deviation of around 10%. We
conclude that studies of first-passage in disordered systems are effectively
pointless from a simulation perspective. Again, the deviations found for the
higher-order moments are more substantial than that of the mean behavior (not
shown).

Another unexpected feature from our exact approach is that _bimodal_ first-
passage distributions arise for certain realizations of the hopping rates.
Such effects can be difficult to observe from a simulation-based approach as
multiple modes can seem apparent due to finite-sample noise. We show one such
example in Fig. 4(a) with $b_{i}$ chosen uniformly in the range $(0,1/3)$ and
$d_{i}=1/3-b_{i}$, with the walk starting at $i=2$ on an interval of length
10. Also shown in this figure are the local biases at each site $i$,
$\ln(b_{i}/d_{i})$. The essential feature of this bias profile is that it is
negative at $i=1$, almost neutral at $i=2$, and generally positive for $i>2$.
Thus a particle starting at $i=2$ exits via the left edge with appreciable
probability and does so quickly. However, if the particle initially hops to
the right, it then experiences a rightward bias, which leads to the second
peak in the first-passage probability. To verify this crude picture, we also
construct a synthetic interval of length 20 in which the segment $[0,5]$ has a
bias to the left and the segment $[5,20]$ has a bias to the right both of
magnitude $|v|=0.3$ (Fig. 4(b)). When the random walk starts at $i\leq 5$, it
is likely to remain in the region $i\leq 5$ and exit the interval on the left
side, corresponding to the early-time peak in Fig. 4(b). However, once the
random walker traverses to the region $i>5$, it is likely to remain on this
side of the interval, ultimately leading to the second, longer-time peak in
Fig. 4(b).

Summary: We derived the first-passage probability and all its moments on
disordered intervals with arbitrary nearest-neighbor hopping rates. Our
approach gives a simple way to understand various anomalies, such as the step-
like features of the first-passage time for individual hopping-rate
realizations. We also quantified how the first-passage time averaged over
finite fraction of all realizations of the hopping rates approaches the true
average when averaging over all realizations of the hopping rates. In any
realistic simulation of a disordered chain, it is practical to sample only a
tiny fraction of all realizations of the disorder, and our results show that
this limitation severely compromises the accuracy of such simulations. We
additionally found the unexpected feature of bimodal first-passage
distributions. These arise in hopping-rate realizations that have a general
outward bias from the interval. For such a system, a random walk exits via one
side of the interval quickly and more slowly from the other side, but with
comparable probabilities for each mode of exit. It would interesting to see
whether this bimodality has potentially useful applications in describing
unexplained occurrences of multimodality in other stochastic processes. In
particular, the relation between transient bimodality in the solution to the
forward equation, where bimodality can arise and dissipate on the way to the
steady state, and bimodalities in the first-passage time distribution have not
been widely explored [37, 38, 39, 40].

_Acknowledgements:_ This publication is based upon work that is supported by
the National Science Foundation under Grant No. DMR-1910736.

## References

  * [1] Feller, W.  _An introduction to probability theory and its applications, Volume 2_ , vol. 81 (John Wiley & Sons, 1991).
  * [2] Van Kampen, N. G.  _Stochastic processes in physics and chemistry_ , vol. 1 (Elsevier, 1992).
  * [3] Karlin, S.  _A first course in stochastic processes_ (Academic press, 2014).
  * [4] Redner, S.  _A guide to first-passage processes_ (Cambridge University Press, 2001).
  * [5] Bray, A. J., Majumdar, S. N. & Schehr, G.  Persistence and first-passage properties in nonequilibrium systems.  _Advances in Physics_ 62, 225–361 (2013).
  * [6] Moran, P. A. P.  Random processes in genetics.  In _Mathematical Proceedings of the Cambridge Philosophical Society_ , vol. 54, 60–71 (Cambridge University Press, 1958).
  * [7] Antal, T. & Scheuring, I.  Fixation of strategies for an evolutionary game in finite populations.  _Bulletin of Mathematical Biology_ 68, 1923–1944 (2006).
  * [8] McKane, A., Alonso, D. & Solé, R. V.  Mean-field stochastic theory for species-rich assembled communities.  _Physical Review E_ 62, 8466 (2000).
  * [9] Piñero, J., Redner, S. & Solé, R.  Fixation and fluctuations in two-species cooperation.  _Journal of Physics: Complexity_ 3, 015011 (2022).
  * [10] Kirman, A.  Ants, rationality, and recruitment.  _The Quarterly Journal of Economics_ 108, 137–156 (1993).
  * [11] Holehouse, J. & Moran, J.  Exact time-dependent dynamics of discrete binary choice models.  _Journal of Physics: Complexity_ 3, 035005 (2022).
  * [12] Perelló, J., Gutiérrez-Roig, M. & Masoliver, J.  Scaling properties and universality of first-passage-time probabilities in financial markets.  _Physical Review E_ 84, 066110 (2011).
  * [13] Chicheportiche, R. & Bouchaud, J.-P.  Some applications of first-passage ideas to finance.  In _First-Passage Phenomena and Their Applications_ , 447–476 (World Scientific, 2014).
  * [14] Ashcroft, P.  Metastable states in a model of cancer initiation.  In _The Statistical Physics of Fixation and Equilibration in Individual-Based Models_ , 91–126 (Springer, 2016).
  * [15] Murthy, K. P. N. & Kehr, K. W.  Mean first-passage time of random walks on a random lattice.  _Phys. Rev. A_ 40, 2082–2087 (1989).  URL https://link.aps.org/doi/10.1103/PhysRevA.40.2082.
  * [16] Hernández-García, E. & Cáceres, M. O.  First-passage-time statistics in disordered media.  _Phys. Rev. A_ 42, 4503–4518 (1990).  URL https://link.aps.org/doi/10.1103/PhysRevA.42.4503.
  * [17] Hernández-García, E., Cáceres, M. O. & San Miguel, M.  Characterizing strong disorder by the divergence of a diffusion time.  _Phys. Rev. A_ 41, 4562–4565 (1990).  URL https://link.aps.org/doi/10.1103/PhysRevA.41.4562.
  * [18] Pury, P. A., Cáceres, M. O. & Hernández-García, E.  First-passage time and the fluctuation of the quenched disorder in biased media.  _Phys. Rev. E_ 49, R967–R970 (1994).  URL https://link.aps.org/doi/10.1103/PhysRevE.49.R967.
  * [19] Goldhirsch, I. & Gefen, Y.  Analytic method for calculating properties of random walks on networks.  _Physical Review A_ 33, 2583 (1986).
  * [20] Noskowicz, S. & Goldhirsch, I.  First-passage-time distribution in a random random walk.  _Physical Review A_ 42, 2047 (1990).
  * [21] Sinai, Y. G.  The limiting behavior of a one-dimensional random walk in a random medium.  _Theory of Probability & Its Applications_ 27, 256–268 (1983).
  * [22] Bouchaud, J., Comtet, A., Georges, A. & Le Doussal, P.  The relaxation-time spectrum of diffusion in a one-dimensional random medium: an exactly solvable case.  _Europhysics Letters_ 3, 653 (1987).
  * [23] Fisher, D. S., Le Doussal, P. & Monthus, C.  Random walks, reaction-diffusion, and nonequilibrium dynamics of spin chains in one-dimensional random environments.  _Physical review letters_ 80, 3539 (1998).
  * [24] Le Doussal, P., Monthus, C. & Fisher, D. S.  Random walkers in one-dimensional random environments: exact renormalization group analysis.  _Physical Review E_ 59, 4795 (1999).
  * [25] Krapivsky, P. L., Redner, S. & Ben-Naim, E.  _A kinetic view of statistical physics_ (Cambridge University Press, 2010).
  * [26] Slutsky, M., Kardar, M. & Mirny, L. A.  Diffusion in correlated random potentials, with applications to dna.  _Physical Review E_ 69, 061903 (2004).
  * [27] Iyer-Biswas, S. & Zilman, A.  First-passage processes in cellular biology.  _Advances in chemical physics_ 160, 261–306 (2016).
  * [28] Lyons, R. & Pemantle, R.  Random walk in a random environment and first-passage percolation on trees.  _The Annals of Probability_ 125–136 (1992).
  * [29] Mardoukhi, Y., Jeon, J.-H., Chechkin, A. V. & Metzler, R.  Fluctuations of random walks in critical random environments.  _Physical Chemistry Chemical Physics_ 20, 20427–20438 (2018).
  * [30] Usmani, R. A.  Inversion of a tridiagonal Jacobi matrix.  _Linear Algebra and its Applications_ 212, 413–414 (1994).
  * [31] Smith, S. & Shahrezaei, V.  General transient solution of the one-step master equation in one dimension.  _Physical Review E_ 91, 062119 (2015).
  * [32] Holehouse, J.  Recurrence and Eigenfunction Methods for Non-Trivial Models of Discrete Binary Choice.  _Entropy_ 25, 996 (2023).
  * [33] Barrio, M., Leier, A. & Marquez-Lago, T. T.  Reduction of chemical reaction networks through delay distributions.  _The Journal of Chemical Physics_ 138, 104114 (2013).
  * [34] Leier, A., Barrio, M. & Marquez-Lago, T. T.  Exact model reduction with delays: closed-form distributions and extensions to fully bi-directional monomolecular reactions.  _Journal of The Royal Society Interface_ 11, 20140108 (2014).
  * [35] Ashcroft, P., Traulsen, A. & Galla, T.  When the mean is not enough: Calculating fixation time distributions in birth-death processes.  _Physical Review E_ 92, 042154 (2015).
  * [36] Higham, N. J.  _Functions of matrices: theory and computation_ (SIAM, 2008).
  * [37] Broggi, G., Lugiato, L. & Colombo, A.  Transient bimodality in optically bistable systems.  _Physical Review A_ 32, 2803 (1985).
  * [38] Lange, W., Mitschke, F., Deserno, R. & Mlynek, J.  Study of fluctuations in transient optical bistability.  _Physical Review A_ 32, 1271 (1985).
  * [39] Jia, C. & Grima, R.  Dynamical phase diagram of an auto-regulating gene in fast switching conditions.  _The Journal of Chemical Physics_ 152, 174110 (2020).
  * [40] Holehouse, J., Sukys, A. & Grima, R.  Stochastic time-dependent enzyme kinetics: Closed-form solution and transient bimodality.  _The Journal of Chemical Physics_ 153, 164113 (2020).
  * [41] Zhang, F.  _The Schur complement and its applications_ , vol. 4 (Springer Science & Business Media, 2006).
  * [42] Usmani, R.  Inversion of jacobi’s tridiagonal matrix.  _Computers & Mathematics with Applications_ 27, 59–66 (1994).

Supplementary Information

##  I Inverse elements of a tridiagonal matrix

In this appendix we succinctly reproduce the inverse elements of a general
tridiagonal matrix first given in [30]. That is, we will derive the elements
of the inverse of the following matrix,

$\displaystyle\mathbb{B}=\begin{pmatrix}b_{1}&c_{1}&0&0&\dots&0&0\\\
a_{2}&b_{2}&c_{2}&0&\dots&0&0\\\ 0&a_{3}&b_{3}&c_{3}&\dots&0&0\\\
\vdots&&&&\ddots&&\\\ 0&0&0&0&\dots&a_{N}&b_{N}\\\ \end{pmatrix},$ (S1)

First, we know that the inverse of a general $N\times N$ matrix is given by
Cramer’s rule as [36],

$\displaystyle\alpha_{ij}\equiv[\mathbb{B}^{-1}]_{i,j}=\frac{(-1)^{i+j}M_{ji}}{\textrm{det}(\mathbb{B})},$
(S2)

where $M_{ji}$ is the minor, i.e., the determinant of the matrix that results
from the removal of row $j$ and column $i$. Our aim is to exploit the
triangular nature of the matrix to give analytic expressions for the inverse
that can be computed much faster than general matrix inversion (which scales
as $N^{3}$ for Gauss–Jordan elimination). For brevity we introduce the
following determinants,

$\displaystyle U_{i}\equiv\begin{vmatrix}b_{1}&c_{1}&0&0&\dots&0&0\\\
a_{2}&b_{2}&c_{2}&0&\dots&0&0\\\ 0&a_{3}&b_{3}&c_{3}&\dots&0&0\\\
\vdots&&&&\ddots&&\\\ 0&0&0&0&\dots&a_{i}&b_{i}\\\ \end{vmatrix},\qquad\qquad
L_{i}\equiv\begin{vmatrix}b_{i}&c_{i}&0&0&\dots&0&0\\\
a_{i+1}&b_{i+1}&c_{i+1}&0&\dots&0&0\\\ 0&a_{i+2}&b_{i+2}&c_{i+2}&\dots&0&0\\\
\vdots&&&&\ddots&&\\\ 0&0&0&0&\dots&a_{N}&b_{N}\\\ \end{vmatrix},$ (S3)

with $U_{0}=1$, $U_{1}=b_{1}$ and $L_{N+1}=1$, $L_{N}=b_{N}$. The key step in
calculating the $\alpha_{ij}$ is in finding the matrix minors. Let’s calculate
the minors for a few examples,

$\displaystyle M_{11}$ $\displaystyle=L_{2},$ $\displaystyle M_{22}$
$\displaystyle=U_{1}L_{3},$ $\displaystyle M_{33}$ $\displaystyle=U_{2}L_{4},$
$\displaystyle M_{12}$ $\displaystyle=a_{2}L_{3},$ $\displaystyle M_{13}$
$\displaystyle=a_{2}a_{3}L_{4},$ $\displaystyle M_{14}$
$\displaystyle=a_{2}a_{3}a_{4}L_{5},$ $\displaystyle M_{N,N}$
$\displaystyle=U_{N-1},$ $\displaystyle M_{N,N-1}$
$\displaystyle=c_{N-1}U_{N-2},$ $\displaystyle M_{N,N-2}$
$\displaystyle=c_{N-1}c_{N-2}U_{N-3}.$

Here we have made judicious use of Schur’s formula [41] for the determinants
of block matrices, which states that for the block matrix

$\displaystyle\mathbb{E}=\begin{pmatrix}\mathbf{A}&\mathbf{B}\\\
\mathbf{C}&\mathbf{D}\end{pmatrix},$ (S4)

with invertible $\mathbf{A}$ and $\mathbf{D}$, we have

$\displaystyle\mathrm{det}(\mathbb{E})=\mathrm{det}(\mathbf{A})\cdot\mathrm{det}(\mathbf{D}-\mathbf{C}\cdot\mathbf{A}^{-1}\cdot\mathbf{B}).$
(S5)

When either $\mathbf{B}$ or $\mathbf{C}$ consist entirely of zeros, the above
formula reduces to,

$\displaystyle\mathrm{det}(\mathbb{E})=\mathrm{det}(\mathbf{A})\cdot\mathrm{det}(\mathbf{D}).$
(S6)

We observe that the minors follow the pattern,

$\displaystyle
M_{ji}=\begin{cases}\left(\prod_{k=i}^{j-1}c_{k}\right)U_{i-1}L_{j+1}&j>i,\\\\[5.69054pt]
U_{j-1}L_{j+1}&j=i,\\\\[5.69054pt]
\left(\prod_{k=j+1}^{i}a_{k}\right)U_{j-1}L_{j+1}&j<i,\\\ \end{cases}$ (S7)

which is proved by induction in [42]. Substituting $M_{ji}$ into Eq. (S2) then
gives our analytic expression for the $\alpha_{ij}$. Note that
$\det(\mathbb{B})=L_{1}=U_{N}$.

The introduction of orthogonal polynomials comes from the recursive
relationship defining the determinants of the sequences
$U_{0},U_{1},\ldots,U_{N}$ and $L_{1},L_{2},\ldots,L_{N+1}$. This recurrence
is trivially given by Leibniz’s rule as,

$\displaystyle\begin{split}L_{N+1}&=1,\;L_{N}=b_{N},\\\
L_{i}&=b_{i}L_{i+1}-c_{i}a_{i+1}L_{i+2},\\\ U_{0}&=1,\;U_{1}=b_{1},\\\
U_{i}&=b_{i}U_{i-1}-c_{i-1}a_{i}U_{i-2}.\end{split}$ (S8)

Upon re-labeling $b_{i}\Rightarrow\beta_{i}(z)$, $a_{i}\Rightarrow d_{i}$,
$c_{i}\Rightarrow b_{i}$, $N\Rightarrow\bar{N}$, $U_{i}\Rightarrow q_{i}(z)$
and $L_{i}\Rightarrow p_{i-1}(z)$ one then recovers the recursive polynomials
defined in Eq. (9) in the main text. Finally, since we only need the terms
$\alpha_{i,1}(z)$ and $\alpha_{i,\bar{N}}(z)$, we can use the properties of
$U_{0}=1$ and $L_{N+1}=1$ to arrive at Eqs. (8) in the main text.

##  II Conditional first-passage probability

We are interested in the _conditional_ FPT to reach either $0$ or $N$ without
ever reaching the opposite boundary starting from some initial site $i$. We
focus on the conditional FPT to reach $N$, and later return to the case of
reaching the origin. The probability of reaching $N$ at time $t$ starting from
$i$ at $t=0$ is given by $f_{i}(t)=P_{i,N}(t)$. Note that $f_{i}(t)$ has the
boundary conditions $f_{N}(t)=\delta_{t,0}$ and $f_{0}(t)=0$. The recurrence
relation satisfied by $f_{i}(t)$ is the same as Eq. (2), but with different
boundary conditions. Note that $f_{i}(t)$ does not strictly correspond to a
normalized probability distribution over $t$ since there is a finite
probability of being absorbed at the origin and therefore never reaching $N$.
Therefore, we further define,

$\displaystyle\widetilde{f}_{i}(t)=\frac{f_{i}(t)}{\sum_{t=0}^{\infty}f_{i}(t)}\equiv\frac{f_{i}(t)}{\phi_{i}},$
(S9)

where $\widetilde{f}_{i}(t)$ is the properly normalized probability density
function and we define $\phi_{i}$ as the probability for the walk to
ultimately reach the boundary at $N$. This eventual hitting probability is
given by

$\displaystyle\phi_{i}=\frac{1+\sum_{k=1}^{i-1}\prod_{i=1}^{k}\frac{d_{j}}{b_{j}}}{1+\sum_{k=1}^{N-1}\prod_{i=1}^{k}\frac{d_{j}}{b_{j}}}.$
(S10)

We now define the generating function of $\widetilde{f}_{i}(t)$ as,

$\displaystyle\widetilde{F}_{i}(z)=\sum_{t=0}^{\infty}z^{t}\widetilde{f}_{i}(t),$
(S11)

and taking Eq. (2) multiplying by $z^{t}$ and summing over all $t$ we get the
three-term recurrence relation that $\widetilde{F}_{i}(z)$ satisfies,

$\displaystyle
b_{i}\phi_{i+1}\widetilde{F}_{i+1}(z)+d_{i}\phi_{i-1}\widetilde{F}_{i-1}(z)+\beta_{i}(z)\phi_{i}\widetilde{F}_{i}(z)=0,\qquad
i\in{1,2,\dots,N-1},$ (S12)

with the boundary conditions now $\widetilde{F}_{0}(z)=0$ and
$\widetilde{F}_{N}(z)=1$ and we have used $\beta_{i}(z)$ as defined in the
main text. For brevity, we henceforth denote
$C_{i}(z)=\phi_{i}\widetilde{F}_{i}(z)$. As in the unconditional case, we
rewrite Eq. (S12) as a matrix equation,

$\displaystyle\mathbb{A}(z)\cdot\mathbf{C}(z)=-\mathbf{m},$ (S13)

with $\mathbf{C}(z)=(C_{1}(z),C_{2}(z),\dots,C_{\tau}(z))$,
$\mathbf{m}=(0,0,\dots,0,b_{\tau})$ and $\mathbb{A}(z)$ the $(N-1)\times(N-1)$
matrix defined by Eq. (S1). Since we already have explicit expressions for the
elements of the inverse of $\mathbb{A}$, we can readily obtain
$\widetilde{F}_{i}(z)=-b_{\bar{N}}\alpha_{i,\tau}/\phi_{i}$:

$\displaystyle\begin{split}\widetilde{F}_{i}(z)&=\frac{-(-1)^{i+\tau}q_{i-1}(z)B_{i}}{p_{0}(z)\phi_{i}},\;\qquad
i\in\\{1,2,\dots,\bar{N}-1\\},\\\\[5.69054pt]
\widetilde{F}_{\bar{N}}(z)&=\frac{-b_{\bar{N}}q_{-1}(z)}{p_{0}(z)\phi_{\bar{N}}}.\end{split}$
(S14)

By symmetry considerations one can then find the conditional generating
functions of the first-passage time to reach $n=0$ as,

$\displaystyle\begin{split}\widetilde{F}_{1}(z)&=\frac{-d_{1}p_{1}(z)}{p_{0}(z)(1-\phi_{1})},\\\\[5.69054pt]
\widetilde{F}_{i}(z)&=\frac{(-1)^{i}p_{i}(z)D_{i}}{p_{0}(z)(1-\phi_{i})},\;\qquad
i\in\\{2,3,\dots,\bar{N}\\}.\end{split}$ (S15)

In the application of the disordered interval in the main text one can use the
conditional first-passage time to observe the origin of bimodalities observed
for some realizations of disorder. This is possible since the origin of the
bimodality seen in Figure 4 comes from the competition between leaving the
opposing ends of the system.