# Persistent and anti-Persistent Motion in Bounded and Unbounded Space:
Resolution of the First-Passage Problem

Daniel Marris<EMAIL_ADDRESS>School of Engineering Mathematics
and Technology, University of Bristol, BS8 1TW  Luca Giuggioli
<EMAIL_ADDRESS>School of Engineering Mathematics and Technology,
University of Bristol, BS8 1TW

###### Abstract

The presence of temporal correlations in random movement trajectories is a
widespread phenomenon across biological, chemical and physical systems. The
ubiquity of persistent and anti-persistent motion in many natural and
synthetic systems has led to a large literature on the modelling of temporally
correlated movement paths. Despite the substantial body of work, little
progress has been made to determine the dynamical properties of various
transport related quantities, including the first-passage or first-hitting
probability to one or multiple absorbing targets when space is bounded. To
bridge this knowledge gap we generalise the renewal theory of first-passage
and splitting probabilities to correlated discrete variables. We do so in
arbitrary dimensions on a lattice for the so-called correlated or persistent
random walk, the one step non-Markovian extension of the simple lattice random
walk in bounded and unbounded space. We focus on bounded domains and consider
both persistent and anti-persistent motion in hypercubic lattices as well as
the hexagonal lattice. The discrete formalism allows us to extend the notion
of the first-passage to that of the directional first-passage, whereby the
walker must reach the target from a prescribed direction for a hitting event
to occur. As an application to spatio-temporal observations of correlated
moving cells that may be either repelled or attracted to hard surfaces, we
compare the first-passage statistics to a target within a reflecting domain
depending on whether an interaction with the reflective interface invokes a
reversal of the movement direction or not. With strong persistence we observe
multi-modality in the first-passage distribution in the former case, which
instead is greatly suppressed in the latter.

##  I Introduction

Advances in tracking sensors for animals [1], cells [2] and particles [3] have
made conspicuous the importance of quantifying the spatio-temporal dynamics of
movement. To this end, two fundamental frameworks, the lattice random walk
(LRW) [4, 5, 6], and its continuous counterpart, the Brownian walk [7], have
been widely employed to model and analyse Markovian trajectories, that is when
the dynamics arise from history independent random processes.

Much of the versatility of the two frameworks, comes from the analytic
description of diffusive processes, the quintessential example of a Markovian
system. The exact mathematical representation of diffusion has in fact
contributed to its widespread application across disciplines, from cell
biology [8] and ecology [9] to social sciences [10] and economics [11].
However, the increased resolution of modern tracking technologies has made
apparent the approximate nature of the diffusive paradigm and is highlighting
the need to model accurately the non-Markovian statistics within movement
paths.

One such example of non-Markovian features in spatio-temporal trajectories is
persistence or anti-persistence in the movement steps, that is by having,
respectively, a higher tendency to continue in the same direction or to
reverse it. Although both movement statistics have been observed, persistent
motion is more common and has been observed at small scales, e.g., bacteria
undergoing run and tumble dynamics [12], as well as being seen at large
scales, e.g., animals foraging [13] and pedestrians in dense crowds [14],
while anti-persistent motion has been reported in hemocytes [15], cell
extracts [16], telomeres of bone osteosarcoma cells [17] through to the
evolution of the score difference in a game of basketball [18]. Persistent
motion also plays a central role in the exciton coherence question [19, 20].

The ubiquity of persistent motion in many natural and artificial phenomena has
generated a wealth of literature on the subject. Within the space-time
continuous description much is known about the one-dimensional occupation
probability [21, 5, 22, 23, 24]. However, unlike the case of Brownian walks,
with correlation in the motion, the standard renewal theory approach [4, 25,
26] to find the first-passage is not sufficient to determine the temporal
dependence of transport related quantities such as the first-passage or first-
hitting statistics to one or multiple absorbing targets. While advances have
been made in bounded space, these studies are limited to one dimension [27,
28] or can only access the mean first-passage time [29].

The main issue with a continuous spatial description in higher dimensions
relates to the need to track the previous movement direction as an internal
degree of freedom. In two or higher dimensions the turning angle (relative to
the previous movement direction) is a continuous variable and it thus requires
an infinite number of internal states [30, 31]. Moreover, as the absorption
property of a target needs to be specified for each direction in which it can
be reached, formalising first-passage processes in continuous space in higher
dimensions remains a challenge.

By employing a discrete spatial description [4, 6] it is possible to overcome
this limitation. To do that we use the so-called correlated random walk (CRW),
proposed for the first time in the 1920s [32, 33], as the natural non-
Markovian extension to the LRW. For the CRW, instead of prescribing a
probability to move in a given direction as in the standard LRW, the movement
steps are determined by the probability that the walker continues in the same
direction as the previous step [5, 21]. In other words, with one parameter one
creates persistence, or positive correlation, if the walker is more likely to
continue and anti-persistence, or negative correlation, if backtracking is
more likely.

Despite the long modelling history of persistent motion in the spatially
discrete literature, the only major work on the CRW first-passage dynamics has
appeared very recently [34]. That study considers an unbounded lattice both in
one and two dimensions under the assumption that the walker moves its first
step in one predetermined direction. With this simplification the one-
dimensional problem may be reduced to a known first-passage quotient [35] and
the exact first-passage dynamics to a single target is derived, while in two
dimensions the dynamics are obtained asymptotically. Despite this recent
progress, determining the general first-passage time dependence to one or
multiple targets has remained an elusive task, with insights for the one
target case in the bounded space relying exclusively upon studies of the
global mean first-passage time in periodically bounded space [36].

In this paper, we generalise the discrete renewal equation and develop a
general theory to determine the first-passage probability to a finite number
of targets for a CRW in bounded and unbounded domains of arbitrary dimensions.
Focusing on bounded space, we first find the analytic representation of the
CRW occupation probability when space is bounded by periodic and reflective
boundaries in hypercubic lattices and by periodic boundaries in a six-sided
hexagonal lattice. We derive a general form of the first-passage probability
and elucidate that with persistence, even in one-dimensional domains, the
first-passage probability can be multi-modal. Furthermore, we extend the
notion of first-passage time to the directional first-passage time, which we
define as the probability that the walker reaches a target site from a
specific direction for the first time.

As an application to flagellated bacterium that may be attracted or repelled
by solid walls [37], we present an iterative process to extract occupation and
first-passage probabilities with arbitrary boundary conditions. With this
method we study the effect of two different reflecting boundary conditions,
namely whether the boundary invokes a change of direction on various transport
statistics, and show that is the former, which aids a searching walker at
intermediate times.

The contents of the rest of the paper is laid out as follows. Section II
treats the formulation of the unbounded problem. The bounded propagators or
Green’s function of the occupation probability in one dimensional domains are
derived in Sec. III, while their higher dimensional counterparts, together
with the case of the hexagonal lattice, which we use to model chiral peristent
motion, are treated in Sec. IV. In Sec. V we turn our attention to the
placement of absorbing traps and discuss the relationship between their
location and the imposition of absorbing boundary conditions on the lattice.
It is here that we generalise the discrete renewal equation to correlated
motion. In Sec. VI, we derive the first moments of these distributions, then
in Sec. VII, we display the iterative method for calculating transport
statistics for cases where the analytic solutions are not known. Concluding
remarks and possible future avenues are then discussed in Sec. VIII.

Figure 1:  (Colour Online). Schematics of the two internal state
representation of a persistent or anti-persistent movement process in one
dimension. A spatial site consists of two states, $m=1,2$ shown, respectively,
in green (upper circle) and blue (lower circle). Allowed steps are drawn with
arrows and the parameters of their corresponding jump probabilities are shown
next to the corresponding arrows. In (a), the periodic and unbounded lattices
are displayed together, because the transition probabilities are identical,
with the only difference arising as the wrapping between sites $1$ and $N$ in
the periodic case. In panel (b) the absorbing state is drawn in red, while
panels (c) and (d) represent the two types of reflection, which we label
‘bounced’ and ‘bunched’ due to the different interaction at the boundary. In
the former, as a walker tries to move out of the domain, it reverses its
movement direction, while in the latter it does not.

##  II Unbounded Propagator

Although the CRW is a one-step non-Markov process, its dynamics can be
represented by a set of coupled second-order Markov equations [38], as amply
shown across the literature both in continuous [39, 24, 40] and discrete time
[34]. For the latter, in one dimension it is given by

$\displaystyle
Q(n,1,t\\!+\\!1)\\!=\\!fQ(n\\!-\\!1,1,t)\\!+\\!bQ(n\\!-\\!1,2,t)\\!+\\!cQ(n,1,t),$
(1) $\displaystyle
Q(n,2,t\\!+\\!1)\\!=\\!bQ(n\\!+\\!1,1,t)\\!+\\!fQ(n\\!+\\!1,2,t)\\!+\\!cQ(n,2,t),$

where $Q(n,m,t)$ is the occupation probability at time $t$ to be at site $n$
and state $m$. The state (or internal degree of freedom) represents the
direction the walker may enter the site, namely a walker in state 1 must have
arrived from the left, while a walker in state 2 from the right (see schematic
in Fig. 1). The parameter $f$ is the probability of continuing forward in the
direction last travelled while $b$ is the probability of backtracking away
from the previous direction travelled. $c=1-(f+b)$ is the so-called lazy or
sojourn parameter whereby $0<f+b\leq 1$, with $f+b=1$ corresponding to a
walker that moves at every time step. When $f>b$ the walker is in the
persistent regime, while for $f<b$ the walker is subject to anti-persistence.
When $f=b$, one recovers the diffusive dynamics modelled by the LRW. The
effect of the lazy parameter is akin to a rescaling of the velocity in the
continuous limit of the CRW, the telegrapher’s equation (see supplementary
material [41] for details).

To solve Eq. (1), it is convenient to write it in matrix form as

$\bm{Q}(n,t+1)=\mathbb{A}\cdot\bm{Q}(n-1,t)+\mathbb{B}\cdot\bm{Q}(n+1,t)+\mathbb{C}\cdot\bm{Q}(n,t),$
(2)

where $\bm{Q}(n,t)=[Q(n,1,t),Q(n,2,t)]^{\intercal}$,
$\mathbb{A}=\bigl{[}\begin{smallmatrix}f&b\\\ 0&0\end{smallmatrix}\bigr{]}$,
$\mathbb{B}=\bigl{[}\begin{smallmatrix}0&0\\\ b&f\end{smallmatrix}\bigr{]}$,
and $\mathbb{C}=\bigl{[}\begin{smallmatrix}c&0\\\
0&c\end{smallmatrix}\bigr{]}$ with the dot operator representing matrix
multiplication.

We assume that the walker enters its initial site from the left with some
arbitrary probability $\alpha_{1}$ and from the right with
$\alpha_{2}=1-\alpha_{1}$, such that the initial condition is written as
$\bm{Q}(n,0)=\delta_{n,n_{0}}\bm{U_{\bm{m}_{0}}}$ where
$\bm{U_{\bm{m}_{0}}}=\left[\alpha_{1},\alpha_{2}\right]^{\intercal}$, i.e.,
the probability of starting in state $m_{0_{1}}=\alpha_{1}$ while the
probability of starting in state $m_{0_{2}}=\alpha_{2}$, and the propagation
is symmetric if $\alpha_{1}=\alpha_{2}$. Using known procedures for random
walks with internal degrees of freedom [42, 35] (outlined in Sec. II of the
supplementary Material [41] for convenience), the generating function,
$\widetilde{f}(z)=\sum_{t=0}^{\infty}f(t)z^{t}$, solution of Eq. (2), the so-
called lattice Green’s function or unbounded propagator, can be obtained
analytically as a column vector of the occupation probability in each internal
state:

$\widetilde{\bm{Q}}_{n_{0}}(n,z)=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{-i\xi(n-n_{0})}\big{[}\mathbb{I}-z\bm{\lambda}(\xi)\big{]}^{-1}\\!\cdot\bm{U_{\bm{m}_{0}}}\text{d}\xi,$
(3)

where $\mathbb{I}$ is the identity matrix and

$\bm{\lambda}(\xi)=\left(\begin{array}[]{cc}fe^{i\xi}+c&be^{i\xi}\\\
be^{-i\xi}&fe^{-i\xi}+c\end{array}\right),$ (4)

is the so-called structure function with $\xi$ the Fourier variable. The
$(i,j)$-th element of $\bm{\lambda}(\xi)$ governs the movement from state $j$
to state $i$, for example, the element $\bm{\lambda}_{1,2}(\xi)$ encodes
moving from state 2 to state 1, which is achieved only via a backtracking,
which increases the coordinate $n$. $\bm{\lambda}(0)$ allows to verify the
normalisation
$\left(\sum_{i=1}^{2}\bm{\lambda}_{i,j}(0)=1,j\in\\{1,2\\}\right)$, as
expected from a stochastic matrix. The associated expression in arbitrary
dimensions may be found in ref. [41], where we show a generalisation of the
known unbounded propagator [43].

The integral in Eq. (3) can be evaluated exactly (see [41]), and after summing
over the internal states
$\widetilde{Q}(n,z)=\sum_{i=1}^{2}\widetilde{\bm{Q}}(n,z)_{i}$, one obtains,
as a generating function, the unbounded occupation probability for the entire
site as

$\displaystyle\widetilde{Q}_{n_{0}}(n,z)=\frac{1}{(1-2cz)\sqrt{\left(1+\delta
z^{2}\right)^{2}-4\varepsilon^{2}z^{2}}}\left[\frac{1-cz}{r(z)^{|n-n_{0}|}}+\frac{\alpha_{1}(b-f)z}{r(z)^{|n-n_{0}+1|}}+\frac{\alpha_{2}(b-f)z}{r(z)^{|n-n_{0}-1|}}\right],$
(5)

where $r(z)=\frac{1+\delta z^{2}}{2\varepsilon
z}\left(1+\sqrt{1-\frac{4\varepsilon^{2}z^{2}}{[1+\delta z^{2}]^{2}}}\right)$,
$\varepsilon=f(1-cz)(1-2cz)^{-1}$ and
$\delta=(f^{2}-b^{2}+c^{2})(1-2cz)^{-1}$, and where the appearance of the
higher order $z$ terms in Eq. (5) makes conspicuous that the system contains
internal degrees of freedom. Using the relation $f+b+c=1$, the scaling term in
Eq. (5) can be rewritten as
$\left[(1-z)[1-z(1-2b)][1-z(1-2f)]\\{1-z[1-2(f+b)]\\}\right]^{-1/2}$. In the
limit $z\rightarrow 1$, as the remaining terms of Eq. (5) give the constant
$2b$, we obtain $\widetilde{Q}_{n_{0}}(n,z)\sim b/\sqrt{2fb(f+b)(1-z)}$,
alternatively the time dependence through Tauberian theorems is
$Q_{n_{0}}(n,t)\sim b/\sqrt{2\pi fb(f+b)t}$ for $t\rightarrow\infty$. As
expected, this scaling is the same as the one of the diffusive LRW, and to
which it reduces when $f=b=q/2$. It is in fact a straightforward exercise to
show that when $f=b$, irrespective of the values of $\alpha_{1}$ and
$\alpha_{2}$, Eq. (5) reduces to the known solution of the symmetric unbounded
walk [6].

##  III Bounded Propagators

We now consider a CRW confined to a bounded domain of size $N$ ($n\in[1,N]$)
and seek the dynamics of the occupation probability in the case of periodic
and reflecting boundaries. Below we use the superscript notation $p$ and $r$,
respectively, to distinguish them and will use the letter $\bm{P}$ to denote a
bounded propagator compared to $\bm{Q}$ for the infinite case.

###  III.1 Periodic Boundary

For the periodic domain, the boundary condition is imposed in the same way as
for the standard Pólya walk [4, 6] making it straightforward to find the
periodic propagator via the method of images (see [41] Sec. III) as

$\bm{P}^{(p)}_{n_{0}}(n,t)\\!=\\!\frac{1}{N}\\!\sum_{\kappa=0}^{N-1}e^{\frac{-2\pi
i\kappa(n-n_{0})}{N}}\bm{\lambda}\left(\frac{2\pi\kappa}{N}\right)^{t}\\!\cdot\bm{U_{\bm{m}_{0}}},$
(6)

where $\bm{\lambda}(\cdot)^{t}$ denotes the matrix power. Diagonalising
$\bm{\lambda}\left(\frac{2\pi\kappa}{N}\right)$ to take the matrix power
explicitly and performing a summation over the states, one finds the site
occupation probability
$P^{(p)}_{n_{0}}(n,t)=\sum_{i=1}^{2}\bm{P}^{(p)}_{n_{0}}(n,t)_{i}$, and when
$\alpha_{1}=\alpha_{2}=\frac{1}{2}$ we obtain

$P^{(p)}_{n_{0}}(n,t)\\!=\\!\frac{1}{2N}\\!\sum_{\kappa=0}^{N-1}\\!\cos\left(\frac{2\pi\kappa(n-n_{0})}{N}\right)\\!\Big{[}f(\kappa,t)+h(\kappa,t)\Big{]},$
(7)

where

$\displaystyle f(\kappa,t)\\!=\\!\lambda_{+}^{t}\\!+\\!\lambda_{-}^{t},\qquad
h(\kappa,t)\\!=\\!\frac{b\cos\\!\left(\frac{2\pi\kappa}{N}\right)\\!\left[\lambda_{+}^{t}\\!-\\!\lambda_{-}^{t}\right]}{\sqrt{f^{2}\\!\cos^{2}\\!\left(\frac{2\pi\kappa}{N}\right)\\!+\\!b^{2}\\!-\\!f^{2}}},$
(8)

and

$\lambda_{\pm}=c+f\cos\left(\frac{2\pi\kappa}{N}\right)\pm\sqrt{f^{2}\cos^{2}\left(\frac{2\pi\kappa}{N}\right)+b^{2}-f^{2}},$
(9)

are the eigenvalues of the matrix structure function
$\bm{\lambda}(\frac{2\pi\kappa}{N})$. The non-symmetric case is less compact
and we give it explicitly in Eq. (S23) of the supplementary material [41].

In the diffusive limit, where $f=b=\frac{q}{2}$, Eq. (6), reduces to the
expression for the propagator for the LRW on a one-dimensional periodic
lattice [6], while in the $f=1$, limit one has
$P^{(p)}_{n_{0}}(n,t)=N^{-1}\sum_{\kappa=0}^{N-1}\left\\{\alpha_{1}\cos\left[2\pi\kappa(n-n_{0}+t)/N\right]\right.+\left.\alpha_{2}\cos\left[2\pi\kappa(n-n_{0}-t)/N\right]\right\\}$,
which shows a superposition of two weighted ballistic waves travelling in
opposite directions. In the other limit $(b=1)$, corresponding to the
perfectly anti-persistent case, the walker constantly hops between the initial
condition, which it visits on even time steps, and the lattice sites on either
side at odd time steps, as one may see from the resulting expression:
$P^{(p)}_{n_{0}}(n,t)=(2N)^{-1}\sum_{\kappa=0}^{N-1}\exp[-2\pi
i\kappa\left(n-n_{0}\right)/N]\big{\\{}[1+(-1)^{t}]-[-1+(-1)^{t}][\alpha_{1}\exp(2\pi
i\kappa/N)+\alpha_{2}\exp(-2\pi i\kappa/N)]\big{\\}}$.

###  III.2 Reflective Barrier

When the correlated particle is confined within a reflecting box, there are
two classes of boundary interactions that one might consider as displayed in
panel (c) and (d) of Fig. 1. More specifically, one must choose whether the
direction of persistence is flipped upon hitting the boundary or not. When the
persistence direction is not flipped then an accumulation of particles appears
at the boundary, which is reflected in the increase in the magnitude of the
steady state probability at the boundary [24, 44]. We call this accumulation a
‘bunching’. In contrast, when the boundary induces a change of direction,
which we call a ‘bounce’, there is no accumulation. It is the ‘bouncing’ case
we consider here and we find the propagator via two different methods, viz.
the method of images [6] and the so-called squeezing method [45], which allow,
respectively, for two variations on this boundary condition. We treat the
‘bunched’ case later in Sec. VII and provide a comparison of the transport
statistics with the ‘bounced’ case.

#### The Method of Images

The method of images is a dimensionally scalable method that allows to include
the effects of interfaces or boundary conditions on the walker dynamics from
the knowledge of the dynamics in the absence of boundaries [46]. To account
for the change in movement direction caused by a reflection, we ensure that
the bounded walker changes state upon interacting with the boundary. To
achieve this, one must specify boundary conditions for each internal state.
More specifically, we require that
$P^{(r)}(1,\\{2,1\\},t)=P^{(r)}(0,\\{1,2\\},t)$, and
$P^{(r)}(N,\\{1,2\\},t)=P^{(r)}(N+1,\\{2,1\\},t)$, where the compact notation
$\\{i,j\\}$ implies pairing the left index $i$ (right index $j$), on the LHS
of the equation, to the left index (right index) on the RHS. In this way, by
allowing the image of state two to interact with the dynamics on state one,
and vice versa, the persistent walker trying to leave the boundary will be
reflected back and continue its journey in the opposite direction (see A for
further discussion). With the boundary conditions set, before applying the
method of images, we need to ensure that the transport process is spatially
symmetric, which we accomplish by imposing
$\alpha_{1}=\alpha_{2}=\frac{1}{2}$.

By taking the state level boundary conditions given above, and performing a
summation over the internal states, one finds the boundary condition for the
entire lattice site. Upon doing so, one find the boundary condition to be
$P^{(r)}_{n_{0}}(1,t)=P^{(r)}_{n_{0}}(0,t)$,
$P^{(r)}_{n_{0}}(N,t)=P^{(r)}_{n_{0}}(N+1,t)$, which, as expected [47], is the
same boundary condition as the diffusive random walker. Since the image set
for this boundary condition has been derived recently [6], one can find the
reflective propagator as
$P^{(r)}_{n_{0}}(n,t)=\sum_{i=1}^{2}\bm{P}^{(r)}_{n_{0}}(n,t)_{i}$ where

$\bm{P}^{(r)}_{n_{0}}(n,t)=\frac{1}{N}\sum_{\kappa=0}^{N-1}\beta_{\kappa}\cos\left(\frac{\pi\kappa(n-\frac{1}{2})}{N}\right)\cos\left(\frac{\pi\kappa(n_{0}-\frac{1}{2})}{N}\right)\bm{\lambda}\left(\frac{\pi\kappa}{N}\right)^{t}\cdot\bm{U_{\bm{m}_{0}}},$
(10)

with $\beta_{\kappa}=1$ when $\kappa=0$, else $\beta_{\kappa}=2$ (see Sec. IV
of the supplementary material [41] for a full derivation).

From Eq. (10), and similarly in Eq. (6), one may notice that, in contrast to
the unbounded dynamics in Eq. (5), within periodic and reflecting domains we
are able to keep separate the spatial and temporal components of the
propagator. In doing so, one realises that, as expected, the respective
spatial dependence in the symmetric $\alpha_{1}=\alpha_{2}$ case coincides
with the one obtained for the so-called lazy Pólya’s walk in finite domains
[6] and the structure function pertains in both solutions as the boundary
condition does not impact the individual jump probabilities in the bulk of the
domain. Since the structure of the two propagators is analogous, one can also
expand Eq. (10) via the same procedure as the periodic case to find explicitly
$P^{(r)}_{n_{0}}(n,t)$. The semi-bounded bounced propagator can also be found
via the method of images and, for completeness, we have presented it in the
supplementary material [41].

#### The Squeezing Method

We exploit an alternate procedure to bound the propagator in a reflective box,
which uses only the properties of the periodic propagator [45, 48, 49]. While
being less scalable to higher dimensions, when compared to the method of
images, it is simple to construct analytically in one-dimensional systems.
This technique is particularly convenient to analyse the case when
$\alpha_{1}\neq\alpha_{2}$, since the use of an image set in a periodic
geometry wraps upon itself, which avoids the interplay between the two
independent image sets used to derive Eq. (10).

We follow ref. [45] and squeeze a torus into a reflecting lattice. More
specifically, we transform a periodic domain with $M\in 2\mathbb{Z}$ sites
into a reflective one with $N=(M+2)/2$ sites ($N\geq 2$) by summing the
probabilities of the sites a distance $l$ on either side of the
$1^{\text{st}}$ and $N^{\text{th}}$ site (see Fig. 1b in [45] or Fig. 1 of the
present supplementary material [41] for a schematic representation). In doing
so, we impose the boundary condition as implemented orignally by Chandrasekhar
for the LRW [50], which states that the walker attempting to escape from $N$
is bounced back to $N-1$ with certainty, and similarly for site $n=1$. The
reflective propagator with this condition is [45]

$\bm{P}_{n_{0}}^{(r_{s})}(n,t)=\bm{P}^{(p)}_{n_{0}}(n,t)+\mu(n)\bm{P}^{(p)}_{n_{0}}(M+2-n,t),$
(11)

where we indicate with the superscript $r_{s}$ that the reflecting domain has
been accounted by squeezing a periodic one, and where $\mu(n)=0$ when $n=1$ or
$N$, and $\mu(n)=1$, otherwise. We note that by construction, in this
instance, one should not think of the states as being representative of the
direction of motion previously travelled, since the baseline propagator on the
torus, from which Eq. (11) is built, considers trajectories that may enter
each lattice site from either direction.

Figure 2: (Colour Online). The occupation probability in a one-dimensional
reflecting domain of length $N=51$ found via the method of images, Eq. (10),
and the squeezing technique, Eq. (12). Panel (a) shows the evolution of the
occupation probability of both propagators, with the method of images $(r)$
solution plotted in purple and the squeezed solution $(r_{s})$ displayed in
blue, for the symmetric $\alpha_{1}=\alpha_{2}=0.5$ case. For clarity, as
there is no differences between the dynamics of the two propagators before
boundary interaction, we omit in the legend the entries for the method of
images solution at short times. The dynamics seen are with $f=0.7$, $b=0.1$
and $n_{0}=35$. Panel (b) depicts the occupation probability of the squeezed
propagator at $t=23$ with $n_{0}=25$, $f=0.8$, $b=0.1$. Here we examine the
effect of changing $\alpha_{1}$ and $\alpha_{2}=1-\alpha_{1}$. The dots on
some curves across the two panels represent the ensemble average of $10^{6}$
stochastic simulations.

Using Eqs. (6) and (11), and after some algebra, the reflective propagator is
given explicitly by

$\displaystyle\bm{P}^{(r_{s})}_{n_{0}}(n,t)=\frac{1}{2N-2}$
$\displaystyle\sum_{\kappa=0}^{2N-3}\bigg{[}\exp\left(\frac{-\pi
i\kappa(n-n_{0})}{N-1}\right)+\mu(n)\exp\left(\frac{-\pi i\kappa(2N-n-
n_{0})}{N-1}\right)\bigg{]}\bm{\lambda}\left(\frac{\pi\kappa}{N-1}\right)^{t}\cdot\bm{U_{\bm{m}_{0}}}.$
(12)

In Fig. 2(a), we plot the reflecting propagator found via both methods (Eqs.
(10) and (12)). For times small enough that no boundary interaction has
occurred, one cannot distinguish between the two representations, as expected.
However, after boundary interaction, the occupation probability differs
slightly due to the no-waiting scheme at the boundary that Eq. (12) obeys.
Since these sites can be viewed as ‘slippy’, i.e., the chance of leaving the
site is higher than the rest of the domain, the steady state distribution
shows a bench-like structure, with these boundary sites having a lower
probability of being occupied. In fact, it can be shown (see B) that when $f$
and $b$ are not identically one,
$P^{(r_{s})}_{n_{0}}(n,t\to\infty)=\frac{1}{(N-1)}$ ($n\neq\\{1,N\\}$) while
$P^{(r_{s})}_{n_{0}}(1,t\to\infty)=P^{(r_{s})}_{n_{0}}(N,t\to\infty)=\frac{1}{2(N-1)}$.

In Fig. 2(b), we display the effect of the weightings over the internal states
in the initial condition. When $\alpha_{1}\neq\alpha_{2}$, the propagation is
clearly asymmetric, as the higher weighting of trajectories that started right
are likely to continue heading right.

##  IV Propagators in Higher Dimensions

As mentioned in the introduction, in continuous space modelling persistence in
higher dimensions requires the inclusion of an infinite number of internal
variables [30]. A discrete representation, on the other hand, requires only
$Z$ internal variables, where $Z$ is the coordination number of the lattice.
In dimensions larger than one, we require a third parameter, $\ell$, to govern
the probability of turning, which we define as a change of movement direction
in any direction except turning back. For simplicity, we first assume the
probability of moving in any lateral direction, to be uniform, although this
can be relaxed with no additional mathematical burden, as we will show below
with the hexagonal lattice. Moreover, in the $d$ dimensional hypercubic
lattice we demand $f+b+\mathfrak{d}\ell+c^{(\mathfrak{d})}=1$, where
$\mathfrak{d}$ is the number of permissible lateral movements, defined as
$\mathfrak{d}=Z-2$ and $c^{(\mathfrak{d})}$ is the sojourn probability. The
initial condition must also account for the higher number of movement
directions, meaning that
$\bm{Q}_{\bm{n}_{0}}(\bm{n},0)=\prod_{i=1}^{d}\left(\delta_{n_{i},n_{0_{i}}}\right)\bm{U}_{\bm{m}_{0}}$,
where $\bm{U}_{\bm{m}_{0}}=[\alpha_{1},...,\alpha_{2d}]^{\intercal}$.

###  IV.1 Hypercubic Lattices

We again start by writing a Master equation to govern the unbounded dynamics.
In $d$ dimensions, this takes the form of $2d$ coupled equations, which we
have written explicitly for the $d=2$ case in the supplementary material [41].
Akin to the $d=1$ case, it is convenient to re-write the coupled equations in
matrix form, and for arbitrary dimension $d$, we obtain

$Q(\bm{n},t\\!+\\!1)=\sum_{i=1}^{d}\left\\{\mathbb{A}_{2d}^{(i)}\cdot
Q(\bm{n}\\!-\\!\bm{e}_{i},t)\\!+\\!\mathbb{B}_{2d}^{(i)}\cdot
Q(\bm{n}\\!+\\!\bm{e}_{i},t)\right\\}+\mathbb{C}_{2d}\cdot Q(\bm{n},t),$ (13)

where $\bm{e}_{i}$ is a unit vector along dimension $i$.
$\mathbb{A}_{2d}^{(i)}$, $\mathbb{B}_{2d}^{(i)}$ are $2d\times 2d$ matrices
that govern the movement probabilities along either direction of the
$i^{\text{th}}$ axis, which is encoded via elements along $(2i-1)^{\text{th}}$
and $2i^{\text{th}}$ rows, respectively, while all other rows in
$\mathbb{A}_{2d}^{(i)}$ and $\mathbb{B}_{2d}^{(i)}$ are identically zero.
$\mathbb{C}_{2d}$ governs the sojourn probability and is of the form
$\mathbb{C}_{2d}=c^{(\mathfrak{d})}\mathbb{I}$.

The method of images allows us to find the periodic and reflective propagator
in arbitrary dimensions. Following the supplementary material [41], we find
the propagator in arbitrary dimensions again as a summation of column vector
of state probabilities
$P^{(\gamma)}_{n_{01},...,n_{0d}}(n_{1},...,n_{d},t)=\sum_{m=1}^{2d}\bm{P}^{(\gamma)}_{n_{01},...,n_{0d}}(n_{1},...,n_{d},t)_{m}$,
where

$\displaystyle\bm{P}^{(\gamma)}_{n_{01},...,n_{0d}}(n_{1},...,n_{d},t)\\!=\\!\frac{1}{N^{d}}\\!\sum_{\kappa_{1}=0}^{N_{1}-1}\\!...\\!\sum_{\kappa_{d}=0}^{N_{d}-1}\left[\prod_{i=1}^{d}g_{\kappa_{i}}^{(\gamma)}(n_{i},n_{0_{i}}\\!)\right]\bm{\lambda}\left(\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}},...,\pi\mathcal{N}^{(\gamma)}_{\kappa_{d}}\right)^{t}\cdot\bm{U}_{\bm{m}_{0}},$
(14)

and $\gamma=p$ or $r$. The function
$g_{\kappa_{i}}^{(\gamma)}(n_{i},n_{0_{i}})$ is the spatial dependence of the
periodic or reflective bounded domain, given by
$g_{\kappa_{j}}^{(p)}(n_{j},n_{0_{j}})=\exp[-2\pi
i\kappa_{j}(n_{j}-n_{0_{j}})/N_{j}]$ and
$g_{\kappa_{j}}^{(r)}(n_{j},n_{0_{j}})=\beta_{\kappa_{j}}\cos[\pi\kappa_{j}(n_{j}-1)/(N_{j})]\cos[\pi\kappa_{j}(n_{0_{j}}-1)/(N_{j})]$,
and where $\mathcal{N}^{(p)}_{i}=2\kappa_{i}/N_{i}$, while
$\mathcal{N}^{(r)}_{i}=\kappa_{i}/N_{i}$. We note that, similarly to the $d=1$
case, the image method demands that when $\gamma=r$ the initial weighting over
the states to be uniform.

The $2d\times 2d$ matrix
$\bm{\lambda}\left(\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}},...,\pi\mathcal{N}^{(\gamma)}_{\kappa_{d}}\right)$
maintains a simple structure regardless of the dimension and also allows
straight forward calculation of the uniform steady state for both $\gamma=p$
or $r$ (see B). Explicit analytical expressions of the scalar propagator for
the site occupation probability depends on the analytical determination of the
eigensystem of
$\bm{\lambda}\left(\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}},...,\pi\mathcal{N}^{(\gamma)}_{\kappa_{d}}\right)$.
We discuss some of these analytic cases starting from the fully correlated or
anti-correlated limits. In those cases many elements of
$\bm{\lambda}\left(\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}},...,\pi\mathcal{N}^{(\gamma)}_{\kappa_{d}}\right)$
reduce to zero. When $b,\ell=0$, the matrix is diagonal, so it is trivial to
take the matrix power explicitly and the propagator is simply the sum of the
matrix diagonal elements weighted according to $\bm{U}_{\bm{m}_{0}}$.
Furthermore, when $f,c,\ell=0$, the structure function reduces to a
generalised permutation matrix, which one may also take to arbitrary power in
closed form.

As
$\sum_{t=0}^{\infty}z^{t}\bm{\lambda}(\cdot)^{t}=[\mathbb{I}-z\bm{\lambda}(\cdot)]^{-1}$
in employing generating functions, it may become more convenient to take a
matrix inverse than calculating the eigensystem to determine the propagator in
time. To illustrate this aspect, we consider the simplified cases discussed in
[51, 36], whereby the authors make the assumption that the probability of a
lateral movement is equivalent to that of a backward movement, which
corresponds to setting $b=\ell$ in our notation. We display the generating
function of the occupation probability for this case in Sec. V of ref. [41].

Furthermore, carrying out the same procedure, but now setting $f=\ell$, one
obtains a structure function matrix where the eigenvalues can be calculated
explicitly (see [41]). Following this, the propagator for the entire site
probability can in fact be determined in the time domain as

$\displaystyle\widetilde{P}^{(\gamma)}_{n_{0_{1}},n_{0_{2}}}(n_{1},n_{2},t)=\frac{1}{2N^{2}}\sum_{\kappa_{1}=0}^{N-1}\sum_{\kappa_{2}=0}^{N-1}\left[\prod_{i=1}^{2}g_{\kappa_{i}}^{(\gamma)}(n_{i},n_{0_{i}})\right]\left[f(\kappa_{1},\kappa_{2},t)+\frac{1-2\ell}{2\ell}g(\kappa_{1},\kappa_{2},t)\right],$
(15)

where

$\displaystyle
f(\kappa_{1},\kappa_{2},t)=\psi_{+}(\kappa_{1},\kappa_{2})^{t}+\psi_{-}(\kappa_{1},\kappa_{2})^{t},$
(16) $\displaystyle g(\kappa_{1},\kappa_{2},t)$
$\displaystyle=\frac{\sigma(\kappa_{1},\kappa_{2})\left[\psi_{+}(\kappa_{1},\kappa_{2})^{t}-\psi_{-}(\kappa_{1},\kappa_{2})^{t}\right]}{\sqrt{\sigma(\kappa_{1},\kappa_{2})^{2}+1-4\ell}},$
(17)
$\psi_{\pm}(\kappa_{1},\kappa_{2})=\sigma(\kappa_{1},\kappa_{2})\pm\sqrt{\sigma(\kappa_{1},\kappa_{2})^{2}+1-4\ell},$
(18)

and

$\sigma(\kappa_{1},\kappa_{2})=\ell\left(\cos\left[\pi\mathcal{N}_{\kappa_{1}}^{(\gamma)}\right]+\cos\left[\pi\mathcal{N}_{\kappa_{2}}^{(\gamma)}\right]\right).$
(19)

###  IV.2 Hexagonal Lattice

There has been some attempts in the past to study correlated motion in a
hexagonal setting. The steady state of a run and tumble model which can change
direction by a rotation of $60^{\circ}$ was studied analytically in [52],
while the path of a persistent photon with varying angles of incidence is
studied as a $12^{\text{th}}$ order Markov chain in a reflecting hexagonal
geometry in [53], where long time diffusion constants and moments are derived.
Moreover, the CRW on a hexagonal lattice has been used in [13] to determine
optimal search strategies for central place foragers.

When considering the hexagonal lattice, it is convenient to use Her’s three
axis coordinate system [54] where each coordinate is represented by three
linearly dependent coordinates $(n_{1},n_{2},n_{3})$ such that
$n_{1}+n_{2}+n_{3}=0$. The finite domain is six-sided, with the number of
lattice points $\Omega=3R^{2}+3R+1$ governed by the circumradius of the
hexagon $R$. As the coordination number in this case is $Z=6$, we require
$f+b+\sum_{j=1}^{4}\ell_{j}+c^{(4)}=1$ and $\bm{U}_{{\bm{m}_{0}}}$ is a
$6\times 1$ column vector. Each lateral movement direction is defined in
relation to the previous movement step where $\ell_{1}$ governs the
probability of a $60^{\circ}$ rotation anticlockwise, $\ell_{2}$ a
$120^{\circ}$ anticlockwise rotation, $\ell_{3}$ a $120^{\circ}$ clockwise
rotation and $\ell_{4}$ that of a $60^{\circ}$ clockwise rotation (see [41]
(Fig. 2) for pictorial representation). With these parameters, we find the
propagator for correlated motion in a periodic hexagonal geometry as (see
[41])

$\displaystyle\bm{P}^{(\mathcal{H})}_{\bm{n}_{0}}(\bm{n},t)=\Bigg{\\{}\frac{\bm{\lambda}^{(\mathcal{H})}(0,0)^{t}}{\Omega}+\frac{1}{\Omega}\sum_{r=0}^{R-1}\sum_{s=0}^{3r+2}e^{\frac{2\pi
i\bm{\kappa}\cdot(\bm{n}-\bm{n}_{0})}{\Omega}}\bm{\lambda}^{(\mathcal{H})}\big{(}\bm{\kappa}(r,s)\big{)}^{t}\Bigg{.}+e^{\frac{-2\pi
i\bm{\kappa}\cdot(\bm{n}-\bm{n}_{0})}{\Omega}}\bm{\lambda}^{(\mathcal{H})}\big{(}-\bm{\kappa}(r,s)\big{)}^{t}\Bigg{\\}}\cdot\bm{U}_{{\bm{m}_{0}}},$
(20)

with the vector
$\bm{\kappa}(r,s)=2\pi[\kappa_{1}(r,s),\kappa_{2}(r,s)]/\Omega$ having
components $\kappa_{1}(r,s)=R(s+1)+s-r$ and $\kappa_{2}(r,s)=R(2-s+3r)+r+1$
[55] and with $\bm{\lambda}^{(\mathcal{H})}(\cdot,\cdot)$ defined in Eq. (S45)
of the supplementary material [41].

In a periodic hexagon the non-orthogonality of the coordinates leads to a
shift as the walker crosses the domain boundaries [55]. In the context of this
present work, if a ballistic walker traverses across the upper boundary it
will continue its path via the bottom edge, but shifted along one coordinate.
One can choose which way to shift the walker, and we have arbitrarily chosen
the right shift here. For such a case, in Fig. 3(a), we display the
probability when the dynamics are ballistic, where one observes pulses of
probability with magnitude 1/6, which originally travelled from the origin to
the six corners, have traversed the boundary and are no longer on their
original axis. Since each side of the hexagon is shifted in the same way, the
rotational symmetry remains, which is also seen in a less persistent case in
panel (b).

Figure 3: (Colour Online). Correlated motion in hexagonal lattices. Panels
(a)-(c) show the occupation probability from Eq. (20) at time $t=11$ in a
hexagonal domain with circumradius $R=8$ ($\Omega=217$) and initial condition
at the origin (shown via the red circles) with
$\bm{U_{\bm{m}_{0}}}=\left(1,1,1,1,1,1\right)^{\intercal}/6$. The boundary
condition is taken as the right-shift periodic boundary condition, and we
provide further details of how this impacts the dynamics in ref. [41]. In
panel (a) the walker is ballistic ($f=1$) and there is no steady state. There
are instead probability pulses of magnitude $\frac{1}{6}$ that travel around
the domain for all time. In panel (b) $f=0.6$, $b=0.06$ and
$\ell_{\\{1,2,3,4\\}}=0.085$, while in panel (c) $f=0.8$, $\ell_{1}=0.2$
($b=\ell_{\\{2,3,4\\}}=0$). Panel (d) depicts a sample trajectory, with the
same movement probabilities as panel (c). The trajectory is over $2000$
timesteps and the colour of the trajectory indicates the time, with light red
indicating short times and dark red indicating long times. The domain is of
size $R=700$, which is chosen so the walker does not interact with the
boundary, and we omit the boundary for visual ease. The initial condition is
shown by the orange circle at (0,0,0) and the coordinate axes for the
hexagonal lattice are shown via the black arrows.

An interesting feature of this boundary condition is that the ballistic walker
in the periodic hexagonal lattice is guaranteed to visit every site, which we
show mathematically in Sec. VI of the supplementary material [41]. In
contrast, in the square lattice, either periodic or reflecting, the ballistic
walker is confined to a genuine sublattice, i.e. only the coordinates
$[(1,n_{0_{2}}),...,(N,n_{0_{2}}),(n_{0_{1}},1),...,(n_{0_{2}},N)]$, which can
be evinced as when $b$ and $s$ are identically zero,
$\left[\mathbb{I}-\bm{\lambda}\left(\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}},\pi\mathcal{N}^{(\gamma)}_{\kappa_{2}}\right)\right]^{-1}$
is singular away from $\kappa_{1}=\kappa_{2}=0$, which occurs only when the
walk is confined to a genuine sublattice [56]. This has interesting
implications on the search time of the walkers in these lattices. Clearly,
when the target is away from the coordinates of the sublattice, the search
time of the walker in the square lattice becomes increasingly longer as
persistence is increased [36] until the probability of a first-passage event
ever taking place becomes zero. However, on the hexagonal lattice, the
ballistic walker will, with probability one, hit the target in $t\leq\Omega$
timesteps making an increase of persistence a good search strategy with the
periodic hexagonal boundary condition.

#### Chiral Motion

By relaxing the assumption that all lateral movement directions are uniformly
distributed, we introduce the possibility of looping patterns into the
movement. We present this as a discrete alternative to the model of ‘chiral’
persistent motion [30, 57, 58, 59], where an angular bias is introduced into
the turning angle, which is of interest as a model of a charged particle
subjected to soft scattering by aligned magnetic domains [30] and particles
driven by microscopic torques such as active colloids [60] and spermatozoa
[61], and chlamydomonas [62]. While chiral motion appears at first to lend
itself to a continuous space representation, it is known that even the two-
dimensional case leads to an analytical intractable Fokker-Planck equation
[30], meaning occupation probabilities and further transport quantities are
difficult to obtain. On the other hand, in the discrete paradigm, to create
looping, one merely alters the movement probabilities in
$\bm{\lambda}^{(\mathcal{H})}(\cdot,\cdot)$ allowing to obtain with ease the
occupation probability.

We create this motion, with anti-clockwise loops, by first setting $c^{(4)}=0$
and then $f,\ell_{1}>b,\ell_{\\{2,3,4\\}}$, i.e. the probability of either
continuing forward or turning $60^{\circ}$ is much larger than any other
movement direction. For this chiral example, using the parameters $f=0.8$,
$\ell_{1}=0.2$, we plot the occupation probability in Fig. 3(c) and a sample
trajectory in Fig. 3(d) in a large domain so that no boundary interaction
takes place. The single trajectory displays periods of persistent motion
before making loopy turns, which are anti-clockwise in most cases. To control
these loops, one may either employ a larger $\ell_{1}$, while decreasing $f$,
so that their frequency increases, or one may increase $\ell_{2}$ at the
expense of $\ell_{1}$, to tighten loops. To create clockwise loops, one
instead places the turning probability into $\ell_{3}$ or $\ell_{4}$ instead.

##  V Absorbing Traps and First-Passage Processes

In the presence of absorbing traps, the one-dimensional boundaries to be
imposed on the equation governing the occupation probability for the
persistent random walk are known [39]. As shown pictorially in Fig. 1(b) for
the one dimensional case, to constrain entirely the domain between sites $n=1$
and $n=N$, the Master equation needs to satisfy the boundary condition
$P^{(a)}(N,1,t)=P^{(a)}(1,2,t)=0$, where now the superscript $a$ indicates the
presence of absorbing boundaries.

In two and higher dimensions, the constraints to represent a fully absorbing
domain can be understood by considering, for simplicity, the square lattice.
Here, with $Z=4$ the states 1,2,3 and 4 represent, respectively, a walker that
has moved West, East, North and South in the last movement step. In this case,
the walker can reach the West of the domain only by increasing its $n_{1}$
coordinate, thus reaching any site with coordinate $n_{1}=N$ in state 1, the
north of the domain by increasing its $n_{2}$ coordinate, thus reaching site
$n_{2}=N$ in state 3, etc. As such, there is only one movement direction the
walker can undertake which results in an absorption at each boundary, leading
to the boundary condition
$P^{(a)}(N,n_{2},1,t)=P^{(a)}(1,n_{2},2,t)=P^{(a)}(n_{1},N,3,t)=P^{(a)}(n_{1},1,4,t)=0$.
One may extend this logic to any lattice, that is, to create a fully absorbing
boundary traps need to be placed in only the one state in each of the boundary
sites that corresponds to the movement direction into that boundary.

On the other hand, in the cases when the traps are not complete along the
domain boundaries, are partially absorbing, or more generally, when traps are
away from boundary sites, one needs to consider all the states in a site to
make it fully absorbing. To do that one must choose an appropriate absorption
free propagator. While the method of images solution in a periodic domain is
amenable to placing traps, the reflective method of images solution is hardly
useful when extracting first-passage distributions. To explain why, recall, as
shown in Sec. III.2 (and A), that the construction of the reflective
propagator is reliant on the interplay between the symmetric propagation in
state one (two) and its oppposite image in state two (one). Hence, once a trap
is placed on an individual state, which stops movement along one direction,
the dynamics with the trap in site one (two) is no longer symmetric with its
image counterpart representing movement along state two (one). To overcome
this issue, in this Section we use the squeezed propagator viz. Eq. (10), for
the reflective occupation probability, which requires no symmetry between
individual image sets, since it is built using periodic propagators.

We are interested in studying first-passage processes both when a single
target consists of a subset of the $i=1,...,M$ internal degrees of freedom or
all of them (the entire site), in the cases where the initial probability is
shared between any or all of the internal states in $\bm{n}_{0}$. Since the
set of trajectories that originate in state $\bm{n}_{0},m_{i}$ are independent
of those starting in $\bm{n}_{0},m_{j}$ and those terminating in
$\bm{s},m_{\bm{s}_{i}}$ are independent of those terminating in
$\bm{s},m_{\bm{s}_{j}}$ $(i\neq j)$, it suffices to consider them separately
and sum the result. Furthermore, when there are multiple targets, the
summation is carried out over both the internal states and the site location.

In the presence of an unordered set of absorbing traps
$\mathdutchcal{S}=\left\\{(\bm{s},m_{\bm{s}})_{i},...,(\bm{s},m_{\bm{s}})_{S}\right\\}$
($S$ is the cardinality of the set), where each element represents a trap in
the internal state at $(\bm{s}_{i},m_{\bm{s}_{i}})$, we obtain (see C for the
details) the following propagator generating function

$\displaystyle\widetilde{P}^{(a)}_{\bm{n}_{0}}(\bm{n},z)=\widetilde{P}_{\bm{n}_{0}}(\bm{n},z)-\sum_{m=1}^{M}\sum_{j=1}^{M}\sum_{i=1}^{S}\alpha_{m_{\bm{s}_{i}}}\widetilde{P}_{\bm{s},m_{\bm{s}_{i}}}(\bm{n},m,z)\frac{\det[\mathbb{H}^{(i)}(\bm{n}_{0},m_{0_{j}},z)]}{\det[\mathbb{H}(z)]},$
(21)

where
$\mathbb{H}(z)_{l,k}=\alpha_{m_{\bm{s}_{k}}}\widetilde{P}_{\bm{s}_{k},m_{\bm{s}_{k}}}(\bm{s}_{l},m_{\bm{s}_{l}},z)$,
$\mathbb{H}(z)_{k,k}=\alpha_{m_{\bm{s}_{k}}}\widetilde{P}_{\bm{s}_{k},m_{\bm{s}_{k}}}(\bm{s}_{k},m_{\bm{s}_{k}},z)$
and $\mathbb{H}^{(i)}(\bm{n}_{0},m_{0},z)$ is the same, but with the
$i^{\text{th}}$ column replaced with
$\alpha_{m_{0_{j}}}\left[\widetilde{P}_{\bm{n}_{0},m_{0_{j}}}(\bm{s}_{1},m_{\bm{s}_{1}},z),...,\widetilde{P}_{\bm{n}_{0},m_{0_{j}}}(\bm{s}_{S},m_{\bm{s}_{S}},z)\right]^{\intercal}$.
The notation
$\alpha_{m_{0_{j}}}P_{\bm{n}_{0},m_{0_{j}}}(\bm{n},m,t)=\alpha_{m_{0_{j}}}\bm{e}_{m}^{\intercal}\cdot\bm{P}_{\bm{n}_{0}}(\bm{n},t)\cdot\bm{e}_{m_{0_{j}}}$,
is the occupation probability that the walker occupies site and state
$(\bm{n},m)$ given its initial position was $(\bm{n}_{0},m_{0_{j}})$, which is
weighted by $\alpha_{m_{0_{j}}}$.

To help elucidate the notation we note that while the site-state pair defining
the trap location is unique, it is possible to have multiple states within one
site containing traps. For example, consider the case of a fully absorbing
trap in site $\bm{s}$ in a two-dimensional square lattice. Since a trap in
state $m_{\bm{s}}=1$ of site $\bm{s}$ absorbs a walker reaching $\bm{s}$ from
the left (movement along the West direction) only, traps must also be placed
in sites $m_{\bm{s}}=2,3$ and 4 to absorb a walker entering the site from any
direction. Hence, for this case
$\mathdutchcal{S}=\left\\{(\bm{s},1),(\bm{s},2),(\bm{s},3),(\bm{s},4)\right\\}$.

### First-Passage and Splitting Dynamics

Knowledge of the occupation probability for the absorbing propagator allows us
to study other transport statistics such as the first exit or first-passage
time [25, 26]. By performing a spatial summation, we find the generating
function of the survival probability, namely
$\widetilde{S}_{\bm{n}_{0}}(z)=\sum_{\bm{n}}\widetilde{P}^{(a)}_{\bm{n}_{0}}(\bm{n},z)$,
and from the known relation
$\widetilde{F}_{\bm{n}_{0}}(\bm{n},z)=1-(1-z)\widetilde{S}_{\bm{n}_{0}}(z)$,
we obtain the generating function for the first-passage probability to an
arbitrary number of target states as

$\widetilde{F}_{\bm{n}_{0}}(\mathdutchcal{S},z)=\sum_{j=1}^{M}\sum_{i=1}^{S}\alpha_{m_{\bm{s}_{i}}}\frac{\det[\mathbb{H}^{(i)}(\bm{n}_{0},m_{0_{j}},z)]}{\det[\mathbb{H}(z)]}.$
(22)

With multiple targets ($S>1$), the theory also allows the determination of the
so-called splitting probability, that is the probability of reaching
$\bm{s}_{k}$ and not any of the other targets in $\mathcal{S}$. Following ref.
[63], one finds

$\widetilde{F}_{\bm{n}_{0}}(\bm{s}_{k}|\bm{s}_{1};...;\bm{s}_{k}-1;\bm{s}_{k}+1;...;\bm{s}_{S},z)=\alpha_{m_{\bm{s}_{k}}}\frac{\sum_{j=1}^{M}\det[\mathbb{H}^{(k)}(\bm{n}_{0},m_{0_{j}},z)]}{\det[\mathbb{H}(z)]}.$
(23)

Equations (22) and (23) are the sought-after quantities. They are the solution
for the multi-target first-passage and splitting problems, respectively, for
random walks with internal degrees of freedom, of which the focus of this
study, the CRW, is one such example. The expressions are general and accounts
for an arbitrary distribution over $\bm{U}_{{\bm{m}_{0}}}$ with any number of
target states, which may be located in one or multiple lattice sites.

Depending on the problem Eq. (22) can be simplified to more compact
expressions. The simplest case is the one target state for the one-dimensional
unbounded CRW. In such a scenario, if the target state is located at site
$s>n_{0}$, one may put the target at $(s,1)$ since the walker may only arrive
at $s$ from the left, as one can evince by looking at Fig. 1(b). With a single
target state, Eq. (22) can be written as

$\widetilde{F}_{\bm{n}_{0},m_{0}}(\bm{s},m_{\bm{s}},z)=\frac{\sum_{j=1}^{M}\alpha_{m_{0_{j}}}\widetilde{P}_{n_{0},m_{0_{j}}}(\bm{s},m_{\bm{s}},z)}{\widetilde{P}_{\bm{s},m_{\bm{s}}}(\bm{s},m_{\bm{s}},z)},$
(24)

which can be further reduced to the known first-passage quotient for random
walks with internal degrees of freedom
$\widetilde{F}_{\bm{n}_{0},m_{0_{k}}}(\bm{s},m_{\bm{s}},z)=\frac{\widetilde{P}_{\bm{n}_{0},m_{0_{k}}}(\bm{s},m_{\bm{s}},z)}{\widetilde{P}_{\bm{s},m_{\bm{s}}}(\bm{s},m_{\bm{s}},z)}$
[35, 34] when the initial condition is localised at state $m_{0_{k}}$.

For the bounded one-dimensional lattice, where the walker may enter from
either direction, the case with two target states at the same spatial site,
e.g., $\mathcal{S}=\\{(s,1),(s,2)\\}$, is particularly relevant. Here, Eq.
(22) can be expanded out (see Sec. VII of ref. [41]) to find

$\displaystyle\widetilde{F}_{n_{0}}(s,z)$
$\displaystyle=\frac{\alpha_{1}\\!\left(\\!\widetilde{P}_{n_{0},1}(s,1,z)\left[\widetilde{P}_{s,2}(s,2,z)\\!-\\!\widetilde{P}_{s,1}(s,2,z)\right]\\!+\\!\widetilde{P}_{n_{0},1}(s,2,z)\\!\left[\widetilde{P}_{s,1}(s,1,z)\\!-\\!\widetilde{P}_{s,2}(s,1,z)\right]\\!\right)}{\widetilde{P}_{s,1}(s,1,z)\widetilde{P}_{s,2}(s,2,z)-\widetilde{P}_{s,1}(s,2,z)\widetilde{P}_{s,2}(s,1,z)\\!}$
$\displaystyle+\frac{\alpha_{2}\\!\left(\widetilde{P}_{n_{0},2}(s,1,z)\left[\widetilde{P}_{s,2}(s,2,z)\\!-\\!\widetilde{P}_{s,1}(s,2,z)\right]\\!+\\!\widetilde{P}_{n_{0},2}(s,2,z)\\!\left[\widetilde{P}_{s,1}(s,1,z)\\!-\\!\widetilde{P}_{s,2}(s,1,z)\right]\right)}{\widetilde{P}_{s,1}(s,1,z)\widetilde{P}_{s,2}(s,2,z)-\widetilde{P}_{s,1}(s,2,z)\widetilde{P}_{s,2}(s,1,z)\\!},$
(25)

Figure 4: (Colour Online). First-passage probability in a finite one
dimensional domain of size $N=22$ with initial condition $n_{0}=8$ and a
target at $s=N$ obtained via a numerical $z$-inverse transform [64, 65] of Eq.
(25), with the periodic propagator, Eq. (6), used for panel (a) and the
squeezed reflective propagator, Eq. (12), used for panel (b). We note here the
use of the two state first-passage in the squeezed reflective propagator is
needed as it is constructed via the periodic dynamics. In the main figure of
panel (a) and (b) we use $\alpha_{1}=\alpha_{2}=1/2$, $f=0.65$ and we vary
$b$, while in the insets we set $b=0$ and vary $f$. In (c) we fix $f=0.82$,
$b=0.08$ and study how the first-passage probability changes as we change
$\alpha_{1}$. The lines are again obtained from Eq. (25) with the appropriate
propagator. Across all three panels the dots represent the ensemble average of
$10^{6}$ stochastic simulations.

where we have lightened the notation with
$\widetilde{F}_{n_{0}}(s,z)=\widetilde{F}_{n_{0}}(\\{(s,1),(s,2)\\},z)$,
$\alpha_{m_{0_{1}}}=\alpha_{1}$ and similarly for the second state weighting.
The alternative derivation via a generalised renewal equation of Eqs. (24) and
(25) can be found in Sec. VII of ref. [41].

A prominent feature of the first-passage probability for a bounded CRW is the
appearance of multiple peaks. While in periodic domains bi-modal [66, 67, 68]
and multi-modal [69] first-passage distributions for Markov LRWs have been
reported, in those cases the appearance of multiple peaks is a consequence of
the periodicity of the boundaries and the presence of a bias for which there
are repeated occasions when a fraction of the trajectories may hit a target.
This is not the case for the CRW, as we observe multiple peaks with periodic
as well as reflecting domains.

With correlated motion we may in fact identify trajectories that travel both
towards and away from the target. With the existence of two peaks documented
in Fig. 4, the first peak corresponds to those trajectories, which travel
persistently towards the target and the second peak represents those that
travel away and later come back via, either, a bouncing reflection or the
periodic boundary. Since the peaks correspond to the two sets of trajectories,
the initial weighting controls which peak contains larger probability, while
the overall level of persistence governs their collective height. Furthermore,
for a given domain size, the path that reaches the target via the periodic
boundary is always shorter than the path that bounces back via the reflective
boundary. Hence, the time at which the second peak occurs is always lower in
periodic domains when compared to its reflective counterpart, which may also
be seen by comparing panels (a) and (b).

In panel (c), the effect of the individual weighting on the height of the
peaks is seen explicity, where we have placed the target at $s=22$ and the
walker intially at $n_{0}=8$. With a smaller proportion of trajectories set to
move in the right direction at $t=0$ ($\alpha_{1}=0.2$), the periodic case has
a large peak at short time, which correspond to the larger fraction of
trajectories that move towards the target via the left boundary. The smaller
peak at later times is due instead to those walker paths that first headed
right. In contrast, with a reflecting boundary at $n=1$ we do not see a peak
until those trajectories that initially headed right reach the target, with a
second smaller peak appearing later once the trajectories initially travelling
left are reflected back from $n=1$. In contrast, we observe a very different
dynamics when we weigh heavily the walker trajectories so that the majority
starts moving right ($\alpha_{1}=0.99$). In such a case we see far less
deviation between the two boundary conditions as most persistent trajectories
do not encounter the left boundary. We note also that the oscillatory
behaviour at short times seen in panels (a) and (c) is a simple consequence of
a low sojourn probability as in $c=0$ limit, one creates parity issues [70].
While here, we are away from this limit, with a low $c$ in short times the
chance of being on an even site at even times (given an even initial
condition) is still higher. The oscillations dampen with time as the
probability of undertaking a sojourn (which breaks the parity) increases.

Figure 5: (Colour Online). First-passage probabilities of correlated random
walks in two-dimensional periodic domains. In panel (a), through a numerical
inversion of Eq. (24), we plot the directional first-passage probability for
each of the four movement directions in a square periodic lattice of length
$N=8$ to a single site at $\bm{s}=(7,3)$. For example, in order to reach
$(\bm{s},m_{\bm{}{s}}=1)$ the walker has to travel West, to reach
$(\bm{s},m_{\bm{}{s}}=3)$ the walker has to travel North, etc. The initial
condition is at $\bm{n}_{0}=(4,2)$ with equal probability to move in any
direction and with $f=0.42$, $b=0.04$ and $\ell=0.18$. Panel (b) depicts
results on a hexagonal lattice and shows, via a numerical inversion of Eq.
(22), the first-passage probability to a set $\mathdutchcal{S}$ of one, two,
and three targets in a periodic hexagon of circumradius $R=4$ ($\Omega=61$) at
locations $\bm{s}_{1}=(-3,3,0)$, $\bm{s}_{2}=(4,-1,-3)$ and
$\bm{s}_{3}=(-3,-1,4)$, which we add sequentially starting from $\bm{s}_{1}$.
Since the walker can enter from any direction, all the internal states in
$\bm{s}_{i}$ are targets. The initial condition in this case is just off
centre at $\bm{n}_{0}=(0,1,-1)$ and $f=0.80$, $b=0.02$, $\ell=0.01$. The
approximate location within the domain of the initial condition and target
sites are shown via the schematic in the inset, with the initial condition a
red open circle and the targets shown with crosses. Dots in (a) and (b) are
the result of $5\cdot 10^{6}$ and $10^{6}$ stochastic simulations,
respectively.

We now turn our attention to the periodic square lattice and study the
probability that the walker reaches a site for the first time moving in a
specific direction, and we call this quantity the directional first-passage
probability. This may be useful in a context where the target is only visible
from a certain direction and is obtained by placing a target in only one
specific state in a lattice site. We show, in Fig. 5(a), the cases in which
the walker reaches the target travelling in each of the four directions. We
note that the directional first-passage is a normalised probability
distribution as it does not exclude trajectories that first reach the target
via an alternate direction. One may ask for the probability of reaching a site
for the first time travelling in a specific direction conditioned on the
walker having not first reached for the first time from any other direction
via the splitting probability in Eq. (23).

With the initial condition at site $\bm{n}_{0}=(4,2)$, to reach $\bm{s}=(7,3)$
for the first time at $t=4$ or $t=5$, the walker has to enter the site
travelling either North or West. However, travelling West gives a higher
directional first-passage at these times because the walker travelling North
into the target at $t=4$ must have made a left-hand turn from $\bm{n}=(7,2)$
at the preceding time. In contrast, to access the site travelling West at an
identical time, the walker has three possible trajectories. After $t=4$, the
West directional first-passage probability sees a sharp decrease, which is a
consequence of this timescale being dominated by unlikely trajectories that
have had to make a sojourn or an even number of backtracks between the intial
condition and the target. Thus, it is here that the Eastern directional first-
passage becomes more likely as this timescale coincides with the time in which
the Eastern travelling persistent trajectories are able to reach the target.

To analyse an example with multiple absorbing targets (Fig. 5(b)), we consider
the hexagonal lattice and place the walker initially in the site
$\bm{n}_{0}=(0,1,-1)$, just to the South-West of the centre of the domain,
with uniform initial weighting $\alpha_{i}=1/6$ for $i=1,...,6$. We place
three targets, the first to the West of the initial condition, the second to
the South-East and the third to the North (see the caption of Fig. 5 for the
exact coordinates). Since the walker can enter any target via all six
directions, the three target case requires 18 target states.

The sustained oscillations, which have only been seen in a small number of
systems, including quantum walks [71], biased random walks in hexagonal
domains [69] and other non-Markov systems [72] are the natural, higher-
dimensional analogue of the bimodality seen in the one-dimensional systems,
where again, the peaks correspond to the timescales for which the walker moves
away from the target as well as those trajectories that miss the target and
must travel round the entire domain. Naturally, as the indirect trajectories
take over and more and more trajectories reach the targets, the oscillations
dampen and decay to zero [69].

While the multi-modal dynamics is present with one, two and three targets, the
placement of targets has a strong effect on the height and location of the
peaks. In Fig. 5(c), since the initial condition is slightly off-centered, the
high persistence of the walker means the first target (placed directly along
the horizontal circumradius) is likely bypassed leading to low search success
until a trajectory, which first set off in the opposite direction traverses
the boundary. When a second and third target are added just away from the
circumradius, the shift imposed by the boundary condition aids the persistent
walker in exploring new sections of the domain, which leads to the emergence
of the large peak in the two and three target cases. Finally, the third target
then kills the second peak since the placement of the third target blocks the
route many persistent trajectories took to $\bm{s}_{1}$, subsequently
increasing the height of the first peak.

##  VI Summary Statistics and the Mean First-Passage Time

With the generating function of the first-passage time distribution to
multiple targets, it becomes an algebraic exercise to extract the
$n^{\text{th}}$ moments, and in particular its first moment, the mean first-
passage time (MFPT),
$\mathdutchcal{F}_{\bm{n_{0}}\to\\{\mathdutchcal{S}\\}}=\left.\frac{\partial\widetilde{F}_{\bm{n_{0}}}(\mathdutchcal{S},z)}{\partial
z}\right|_{z=1}$.

When the initial condition is uniformly distributed across the site, the
structure of Eq. (22) is very similar to the Markov first-passage expression
for multiple targets [63], and one differentiates
$\widetilde{F}_{n_{0}}(\mathdutchcal{S},z)$ to find

$\mathdutchcal{F}_{\bm{n_{0}}\to\\{\mathdutchcal{S}\\}}=\frac{\det[\mathbb{T}_{0}]}{\det[\mathbb{T}_{1}]-\det[\mathbb{T}]}.$
(26)

The matrices in Eq. (26) are made up of the MFPT between localised sites
$\mathdutchcal{F}_{\bm{n}_{0},m_{0}\to\bm{s}_{i},m_{\bm{s}_{i}}}$. More
precisely, $\mathbb{T}_{i,i}=0$,
$\mathbb{T}_{i,k}=\mathdutchcal{F}_{\bm{s}_{k},m_{\bm{s}_{k}}\to\bm{s}_{i},m_{\bm{s}_{i}}}$
($i\neq k$),
$\mathbb{T}_{{0}_{i,k}}=\mathbb{T}_{i,k}-(2d)^{-1}\sum_{j=1}^{2d}\mathdutchcal{F}_{\bm{n}_{0},m_{0_{j}}\to\bm{s}_{i},m_{\bm{s}_{i}}}$
and $\mathbb{T}_{1_{i,k}}=\mathbb{T}_{i,k}-1$. To find these quantities one
evaluates expressions of the form
$\mathdutchcal{F}_{\bm{n}_{0},m_{0}\to\bm{s}_{i},m_{\bm{s}_{i}}}=\frac{\partial}{\partial
z}\left.\frac{\widetilde{P}_{\bm{n}_{0},m_{0}}(\bm{s},m_{\bm{s}},z)}{\widetilde{P}_{\bm{s},m_{\bm{s}}}(\bm{s},m_{\bm{s}},z)}\right|_{z=1},$
and we give details and explicit quantities in Sec. VIII of ref. [41].
Performing a summation of Eq. (26) over $\bm{n}_{0}$ yields the global mean
first-passage time (GMFPT) [35], which we denote as
$\mathdutchcal{G}_{\\{\mathdutchcal{S}\\}}$. In such a case, one replaces
$\mathbb{T}_{0}$ with $\mathbb{G}$ where
$\mathbb{G}_{i,j}=\mathbb{T}_{i,j}-(2dN)^{-1}\sum_{\bm{n}_{0}}\sum_{j=1}^{2d}\mathdutchcal{F}_{\bm{n}_{0},m_{0_{j}}\to\bm{s}_{i},m_{\bm{s}_{i}}}$,
while $\mathbb{T}$ and $\mathbb{T}_{1}$ are unaffected by the summation. We
discuss this quantity further in ref. [41].

Figure 6: (Colour Online). First-passage probabilities to a single target site
and their related mean first-passage times in one-dimensional lattices with
periodic (orange) and reflecting (green) boundary conditions. The meaning of
the colours, given explicitly in the legend of panel (a), extends over all
three plots. In panel (a) we study the MFPT in an $N=35$ lattice for a walker
starting at $n_{0}=30$ with equal weighting of the states in the anti-
persistent regime ($b>f$) with the target at $N=s$, by fixing $b=0.6$ and
varying $f$. Dots represent Eq. (26) evaluated with Eq. (S74) of ref. [41],
for the periodic lattice, and Eq. (S75) of ref. [41], for the reflecting
lattice. Crosses represent the result of $10^{6}$ stochastic simulations. For
the same domain and initial conditions, panel (b) shows the full first-passage
probability when $f=0.22$ and $b=0.6$. The solid lines are a numerical
inversion of Eq. (24), the dots are the result of $10^{6}$ stochastic
simulations and the corresponding MFPTs are shown with the arrows. The inset
of panel (b) is a blown-up look at the first-passage probabilities in the
tails for $t\in[2000,5000]$. Panel (c) displays the first-passage probability
in the persistent regime ($f=0.86$, $b=0.06$) where $N=s=19$ and $n_{0}=14$.
Solid lines are generated from the numerical inversion of Eq. (24), the dots
are obtained from $10^{6}$ stochastic simulations and the arrows indicate the
MFPT.

In Fig. 6(a), we plot the MFPT for anti-persistent walkers $(f<b)$ as a
function of $f$, and a given $b$. The increase of $f$ has two effects on the
walker’s ability to spread, namely it decreases both the chance of a sojourn
and the probability of turning backwards. This leads to an exponential
decrease in the mean search time as the dynamics approaches the diffusive
regime. Even though the target is placed only five lattice sites away from the
initial condition, the walker takes, on average, 13651 time steps to complete
its search in the strong anti-persistent regime $(f=0.02$). This is due to, in
part, those trajectories that have moved far from the target and subsequently
have got ‘stuck’ in cycles of constant backtracking.

It is in panel (b) that one further sees the effect of those trajectories that
get stuck far from the target as we display the first-passage probability in
another anti-persistent regime where ($f=0.22$, $b=0.6$). Since, the initial
condition is close to the target and the exploration is extremely slow, there
is little difference between periodic and reflecting domains at short and
intermediate times. Differences surface at long times since in the periodic
domain those walkers that have explored the domain in its entirety may find a
short route back to the target, while in the reflective domain walkers are
forced to travel back the whole length of the domain. It is these trajectories
that contribute to longer tails in the first-passage probability and
subsequently the larger MFPT in the reflective case. We note that, similarly
to the persistent case, the local minima (and subsequent kinks) seen at short
times is again due to the low sojourn probability making the probability of a
first-passage event at even times lower.

Figure 6(c) draws attention to the importance of studying the entire first-
hitting distribution in the presence of persistent motion ($f>b$). When first-
passage processes display multiple peaks in their time-dependent probability,
the MFPT may in fact provide misleading values of the time it takes to hit a
target. An extreme example is the ballistic limit, which has two distinct and
narrowly defined hitting times given at, for example,
$\mathdutchcal{F}_{n_{0}\to s}^{(r_{s})}=s-1$ $(s>n_{0})$,
$\mathdutchcal{F}_{n_{0}\to s}^{(r_{s})}=N-s$ $(s<n_{0})$ (for
$\alpha_{1}=\alpha_{2}$). In such a case the MFPT would give a temporal value
in between the two during which the probability of a first-passage event is
identically zero. As shown in panel (c), this aspect may also appear away from
the ballistic limit, in the regime of high persistence.

##  VII Iterative Procedure and Comparison Between Reflective Interactions

In the preceding sections we have exploited the analytical expressions of
propagator generating functions. For the cases in which propagators are not
known in closed form we present an iterative procedure to find the first-
passage probability. While the procedure is valid in arbitrary dimensions, for
notational ease we present it here explicitly for the one-dimensional case,
and use it as a tool to make comparisons between the two forms of reflection
mentioned in Sec. III.2. We differentiate between the two reflections based on
whether an interaction with the boundary invokes a reversal in the direction
of movement, which we call the ‘bouncing’ case, or if the movement direction
persists, which we call the ‘bunching’ case (see panels (c) and (d) in Fig. 1
for pictorial comparison). For an in depth appreciation of the differences
between the two boundary conditions we also compare here the occupation
probability and the mean square displacement (MSD) in the two cases.

We begin by rewriting the Master equation, Eq. (2), for a finite system using
transition matrices that describe the possible movement between lattice sites
and link the different internal states. In other words we represent spatial
jumps for which the internal state does not change as well as those for which
changes between internal states occur. Noting that the probability across the
entire site is given by
$\bm{P}_{n_{0}}(t)=\bm{P}_{n_{0}}(1,t)+\bm{P}_{n_{0}}(2,t)$, where
$\bm{P}_{n_{0}}(1,t)$ is an $N\times 1$ column vector of the occupation
probabilities of state 1 in each lattice site at time $t$, we have the Master
equation

$\displaystyle\bm{P}_{n_{0}}(1,t+1)$
$\displaystyle=\mathbb{D}\cdot\bm{P}_{n_{0}}(1,t)+\mathbb{E}\cdot\bm{P}_{n_{0}}(2,t),$
(27) $\displaystyle\bm{P}_{n_{0}}(2,t+1)$
$\displaystyle=\mathbb{F}\cdot\bm{P}_{n_{0}}(1,t)+\mathbb{G}\cdot\bm{P}_{n_{0}}(2,t),$

where $\mathbb{D}$ represents the transition probabilities from state 1 to
state 1, $\mathbb{E}$ from state 2 to state 1, $\mathbb{F}$ from state 1 to
state 2 and $\mathbb{G}$ encoding the transitions from state 2 to state 2. Due
to the finite space, each matrix is of size $N\times N$ and as each matrix
only governs a subset of the movement possibilities, they are spars and have
the following forms: $\mathbb{D}_{i,i}=c$, for $i\neq 1,N$,
$\mathbb{D}_{i,i+1}=a$ and $\mathbb{D}_{1,1}=\delta+c$,
$\mathbb{D}_{N,N}=\sigma+c$; $\mathbb{E}_{i,i+1}=b$ and
$\mathbb{E}_{1,1}=\omega$; $\mathbb{F}_{i,i-1}=b$ and
$\mathbb{F}_{N,N}=\omega$ and $\mathbb{G}_{i,i-1}=a$ and
$\mathbb{G}_{1,1}=\sigma+c$, $\mathbb{G}_{N,N}=\delta+c$, where $\delta$,
$\sigma$, $\omega$ are dependent on the chosen boundary conditions. With
initial conditions $P_{n_{0}}(n,1,0)=\alpha_{1}$,
$P_{n_{0}}(n,2,0)=\alpha_{2}$, that is
$\bm{P}_{n_{0}}(1,0)=\alpha_{1}\bm{e}_{n_{0}}$ and
$\bm{P}_{n_{0}}(2,0)=\alpha_{2}\bm{e}_{n_{0}}$, it is an iterative task to
find the occupation probability at time $t$.

For the ‘bunched’ case, $\delta=b$, $\sigma=a$, $\omega=0$, while for the
‘bouncing’ case, $\delta=b$, $\omega=a$ and $\sigma=0$. To include absorbing
traps one sets the outgoing probabilities of a trap state $(s,m_{s})$ as
$(1-\rho_{s,m_{s}})A(s,m_{s},n,m)$ where $\rho_{s,m_{s}}$ is the partial
absorption probability at $(s,m_{s})$ and with $A(s,m_{s},n,m)$ representing
the elements of the transition matrix encoding movement out of $(s,m_{s})$. As
an illustration, consider a trap at site $s=N$ and state $m_{s}=1$. In such a
case we multiply the $N^{\text{th}}$ column of $\mathbb{D}$ and $\mathbb{F}$
by $(1-\rho_{N,1})$, and the first-passage is obtained via
$F_{n_{0}}(n,t)=S_{n_{0}}(t-1)-S_{n_{0}}(t)$ with $F_{n_{0}}(n,0)=0$, where
$S_{n_{0}}(t)=\sum_{n=1}^{N}P_{{n}_{0}}(n,t)$.

In Fig. 7(a), we plot the occupation probability for the bouncing and bunching
reflective boundary conditions at varying time steps with uniform initial
weighting $\alpha_{1}=\alpha_{2}=1/2$. Even at short times one sees higher
probability accumulating at the boundaries in the ‘bunched’ case, with the
asymmetry in heights due to the initial condition being placed slightly to the
right. At longer times, with the chance of more boundary interaction, the
differences become more stark, which pertains until we recover the expected
steady states, i.e., uniform probability for the ‘bouncing’ case and increased
probability at the boundaries for the ‘bunched’ [24].

Figure 7: (Colour Online). A comparison of the CRW dynamics between two forms
of reflective boundaries, the ‘bunching’ and ‘bouncing’ interaction. Across
all four panels, the blue-green colours represent the ‘bunching’ boundary
condition, while the red-purple colours represent the ‘bouncing’ condition.
All lines are found via iteratively solving Eq. (27), while crosses are
obtained using the analytic propagator in Eq. (10) and all dots are the result
of $10^{6}$ stochastic simulations. Panel (a) displays the occupation
probability at varying time steps with $f=0.75$, $b=0.12$ with $n_{0}=16$ and
$N=30$, and panel (b) shows the MSD for a variety of persistence levels with
($b=1-f$) in an $N=19$ lattice with $n_{0}=10$, that is with the initial
condition in the middle of the domain. Panel (c) shows the first passage
distribution to site $s=N=30$, while its inset is a close up of the region of
discrepancy between the two boundary conditions when $f=0.75$, $b=0.12$. Panel
(d) shows the first-passage probability to a partially absorbing target placed
at $s=N-5$ ($N=30$), with initial condition $n_{0}=12$ and $f=0.75$, $b=0.12$.
The inset displays the same type of information, but the movement
probabilities are now $f=0.8$ and $b=0.04$.

In panel (b), we compare the MSD, $\Delta(t)=\langle(n(t)-n_{0})^{2}\rangle$,
with $\langle\cdot\rangle$ denoting the ensemble average, of a walker confined
between either boundary. In the diffusive limit, both boundary conditions are
equivalent, which is confirmed by the identical MSD in this case. However,
with increasing persistence one sees striking differences between the MSD
after boundary interactions have occurred. As predicted from the continuous
space-time equivalent analysis for the bouncing case [73], with high positive
correlation the MSD shows oscillatory dynamics. In the ballistic limit, these
oscillations persist for all times as the walker is guaranteed to return to
its initial condition periodically every $N$ timesteps. However, with non-zero
backtracking probability, the MSD displays oscillations at short times, which
subsequently decay as the result of the different times at which trajectories
return to the starting location. Then, due to the uniform steady state when
$b\neq 0$ the MSD eventually saturates to $N^{2}/12$. In contrast, the
bunching case sees a clear maxima before quickly lowering to a saturation
value, which is dependent on the level of persistence in the system. Since the
shape of the steady state is dependent on the level of persistence in the
system (where higher persistence levels leads to larger values of the steady
state at the boundary [24]), the MSD saturates at higher value for larger
$f-b$ with the maximal displacement remaining for all times in the ballistic
limit.

In panels (c) and (d), we study the first-passage and as expected, as we
increase the possibility of the walker interacting with the boundary, more
differences emerge in the first-passage dynamics. To illustrate, in the fully
absorbing case, seen in panel (c), differences in the dynamics emerge only at
intermediate times when the ‘bouncing’ walkers display a much higher second
peak compared to the ‘bunching’ ones as they are allowed to reverse their
trajectory back towards the target. Instead, in panel (d), where the target is
partially absorbing and away from the boundary, the walker may now interact
with both the left and right boundary and two timescales emerge where
differences occur. The first timescale is due to the trajectories overcoming
the target and subsequently being reflected from the right boundary, while the
second results from the trajectories that are reflected from the left. Thus,
it is the second timescale that is present for all values of $\rho$ (including
$\rho=1$), while the first only emerges for a smaller absorption probability.
As we lower $\rho$ further and enter the reaction limited regime (when the
timescale for a reaction to occur is much higher than the transport time to
reach the target [74]) the two timescales lengthen and begin to merge. For
example, when $\rho=0.05$, one sees a sustained region where the ‘bouncing’
case has a higher first-passage probability as the walker subject to this
boundary condition continuously passes over a target with low adsorption many
times.

##  VIII Conclusion

First-passage processes in Markovian LRWs have been studied extensively over
the last few years in a multitude of settings [75, 76, 77, 66, 55, 6, 78, 79,
80, 81]. On the other hand, for non-Markovian LRWs a general formalism to
derive the first-passage dynamics has been lacking. The simplest non-Markovian
LRW extension, the CRW, has received attention, but much of the analysis on
first-passage statistics has relied either on time-consuming computational
methods [13] or on the global MFPT [36], with very few results known
analytically. Here, we have constructed a general formalism to determine the
CRW first-passage probability with partially or fully absorbing sites.
Furthermore, by deriving closed form expressions for the occupation
probability in some instances, we have been able to obtain the first-passage
probability generating function analytically. For those cases for which the
occupation probability is not known in closed form, we have constructed an
iterative procedure to obtain the first-passage probability avoiding the use
of computationally burdening stochastic simulations.

While there have been generalisations of renewal processes with continuous
variables in the past, e.g., to non-Markovian dynamics [29] and ageing
processes [82], they are limited by their analytic tractability, for example,
obtaining only moments or asymptotic expansions. Our study instead generalises
the renewal formulation in discrete space and time presented by Erdös, Feller
and Pollard in 1949 and allows us to link analytically the occupation
probability to the first-passage probability for persistent and anti-
persistent walks with one or multiple targets. Our formalisms have highlighted
the natural appearance of multiple modes in the first-passage dynamics when
CRWs with strong persistence occur within finite periodic and reflecting
domains. We have also quantified the directional first-passage statistics,
that is the first-hitting dynamics based on the direction from which a walker
reach a target.

We have obtained analytically the MFPT to one or multiple targets in periodic
and reflecting domains, and found a general expression for the GMFPT, which
was previously only known for periodic domains [36]. By comparing the first-
passage dynamics with its respective MFPT, we find that the multi-modality of
the first-passage may lead to cases where the mean actually falls between two
modes in a temporal range with a very low chance for a first-passage event to
occur.

The theory presented here is rather flexible. It can be extended to biased-
correlated random walks [83, 84, 85], may cope with the inclusion of resetting
[75, 68] and inert spatial heterogeneities [77] such as permeable barriers
[86, 87, 88] or different media [76], and may also be used to account for the
dynamics in random environments [89]. We wish also to draw the readers’
attention to its applicability to other processes that are conveniently
modelled using random walks with internal states such as walks over non-
Bravais lattices [55, 35], double-diffusivity [90, 91], dimer migration on a
crystalline surface [92, 40] and chromatographic processes [5].

Our findings may also open up the possibility to analyse, the dynamics of
active particles in high spatial dimensions [51], encounter statistics [63] of
persistent walkers, exclusion processes [93, 94, 95], cover times [96], record
statistics [97] and walks with longer range memory such as the so-called
alzheimers [98] and the elephant random walks [99, 100]. Finally, owing to the
connection between the two-step CRW and quantum random walks [101], this work
in the classical paradigm could help provide tools for its quantum
counterpart.

## Acknowledgements

The authors acknowledge useful discussions with Toby Kay, Hernán Larralde,
Ralf Metzler, and Seeralan Sarvaharman and would like to thank the Isaac
Newton Institute for Mathematical Sciences for support and hospitality during
the programme Mathematics of movement: an interdisciplinary approach to mutual
challenges in animal ecology and cell biology when part of the work on this
paper was undertaken. This work was supported by: EPSRC grant number
EP/R014604/1 and was carried out using the computational facilities of the
Advanced Computing Research Centre, University of Bristol -
http://www.bris.ac.uk/acrc/. DM acknowledges funding from an Engineering and
Physical Sciences Research Council (EPSRC) DTP studentship, while LG
acknowledges funding from Biotechnology and Biological Sciences Research
Council (BBSRC) Grant No. BB/T012196/1 and the Natural Environment Research
Council (NERC) Grant No. NE/W00545X/1.

##  Appendix A Reflective Boundary Conditions

To employ the method of images, we note the following. Since a walker moving
right (left) must enter into state $m=1$, ($m=2$), at the right (left)
boundary, it is only the contribution in state $m=1$ ($m=2$) that needs to be
matched by an image. Then to model the ‘bouncing’, we invoke a change of state
to allow propagation back in the opposing direction by allowing the image of
the other state to take over. In Fig. 8, we show this interaction where in the
bottom panel, in the succeeding time step, the trajectories that continue left
over boundary will be replaced by those in state $m=1$, indicating the change
in direction.

Figure 8: (Colour Online). The interplay between the bounded domain (shown
between the two red dotted lines) and the images in three neighbouring domains
at timestep $t=1$ (left) and $t=4$ (right) found by plotting the scalar
expansion of Eq. (10). The occupation probability of state one is shown in
green while that of state two in blue. In both panels $f=0.8$ and $b=0.2$.

This boundary conditions are generalisable to higher dimensions, e.g. the case
$d=2$ requires

$\displaystyle P^{(r)}(1,n_{2},\\{2,1\\},t)$
$\displaystyle=P^{(r)}(0,n_{2},\\{1,2\\},t),$ (28) $\displaystyle
P^{(r)}(N_{1},n_{2},\\{1,2\\},t)$
$\displaystyle=P^{(r)}(N_{1}+1,n_{2},\\{2,1\\},t),$ $\displaystyle
P^{(r)}(n_{1},1,\\{4,3\\},t)$ $\displaystyle=P^{(r)}(n_{1},0,\\{3,4\\},t),$
$\displaystyle P^{(r)}(n_{1},N_{2},\\{3,4\\},t)$
$\displaystyle=P^{(r)}(n_{1},N_{2}+1,\\{4,3\\},t),$

whose sum over the entire site gives,
$P^{(r)}(n_{1},\\{1,N_{2}\\},t)=P^{(r)}(n_{1},\\{0,N_{2}+1\\},t)$ and
$P^{(r)}(\\{1,N_{1}\\},n_{2},t)=P^{(r)}(\\{0,N_{1}+1\\},n_{2},t)$. This allows
one to see that the boundary conditions along either axis are independent of
the other, which allows the easy use of the methods of images (see details in
the supplementary Material [41]).

##  Appendix B The Structure Function

The structure (or characteristic) function is the discrete Fourier transform
of the individual step probabilities of the walk [102, 4], which allows for
easy inspection of the transport properties of the walk [56]. In systems with
no internal degrees of freedom the structure function is scalar, while the
inclusion of internal states makes necessary the generalisation to matricial
functions of size $\mathfrak{n}\times\mathfrak{n}$, where $\mathfrak{n}$ is
the number of internal states.

Presently, the matrix
$\bm{\lambda}\left(\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}},...,\pi\mathcal{N}^{(\gamma)}_{\kappa_{d}}\right)$
is of size $2d\times 2d$ and has a common structure for all the lattices we
have considered here. Namely, the probability of remaining at the same state,
i.e. persisting in the direction of the previous movement step or sojourning,
is encoded in the diagonal elements, the probability of backtracking is given
on the upper and lower diagonal in alternating rows, while all the remaining
elements relate to the option of turning laterally. These features are evident
in Eq. (4) for the $d=1$ case, and are seen explicitly also for the two-
dimensional hypercubic lattice when we write out the matrix

$\bm{\lambda}\left(\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}},\pi\mathcal{N}^{(\gamma)}_{\kappa_{2}}\right)\\!=\\!\begin{bmatrix}fe^{i\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}}}+c^{(2)}&be^{i\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}}}&\ell
e^{i\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}}}&\ell
e^{i\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}}}\\\
be^{-i\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}}}&fe^{-i\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}}}+c^{(2)}&\ell
e^{-i\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}}}&\ell
e^{-i\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}}}\\\ \ell
e^{i\pi\mathcal{N}^{(\gamma)}_{\kappa_{2}}}&\ell
e^{i\pi\mathcal{N}^{(\gamma)}_{\kappa_{2}}}&fe^{i\pi\mathcal{N}^{(\gamma)}_{\kappa_{2}}}+c^{(2)}&be^{i\pi\mathcal{N}^{(\gamma)}_{\kappa_{2}}}\\\
\ell e^{-i\pi\mathcal{N}^{(\gamma)}_{\kappa_{2}}}&\ell
e^{-i\pi\mathcal{N}^{(\gamma)}_{\kappa_{2}}}&be^{-i\pi\mathcal{N}^{(\gamma)}_{\kappa_{2}}}&fe^{-i\pi\mathcal{N}^{(\gamma)}_{\kappa_{2}}}+c^{(2)}\end{bmatrix}\\!.$
(29)

As mentioned in the main text (Sec II), the $(i,j)$-th element of
$\bm{\lambda}\left(\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}},...,\pi\mathcal{N}^{(\gamma)}_{\kappa_{d}}\right)$
represents the movement that may occur from state $m_{j}$ to $m_{i}$. The
elements of one row, which have identical exponents in the exponentials,
represent the one (unique) permissible movement direction of the walker
entering the state $m_{i}$. For example,
$i\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}}$ in the first row indicates that to
enter state $m_{1}$, the walker must increase its $n_{1}$ coordinate, that is
by travelling West, $-i\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}}$ in the second
row shows that to enter $m_{2}$ a decrease in the $n_{1}$ coordinate occurs by
moving East, etc.

From the structure function one may also gain information about the steady
state of the system. When $f,b\neq 1$, for hypercubic lattices and hexagonal
lattices, in the limit $t\to\infty$, one has, respectively,
$\bm{\lambda}\left(\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}},...,\pi\mathcal{N}^{(\gamma)}_{\kappa_{d}}\right)^{t}=\bm{0}$
and $\bm{\lambda}^{(\mathcal{H})}\\!\big{(}\bm{\kappa}\big{)}^{t}=\bm{0}$ when
$\bm{\kappa}\not=0$, while
$\lim_{t\to\infty}\bm{\lambda}\left(0,...,0\right)^{t}=(2d)^{-1}\mathbb{J}$,
and $\bm{\lambda}^{(\mathcal{H})}\\!\big{(}0,0\big{)}^{t}=(6)^{-1}\mathbb{J}$,
where $\mathbb{J}$ is the all ones matrix. Therefore,
$\bm{P}^{(\gamma)}_{n_{01},...,n_{0d}}(n_{1},...,n_{d},t\to\infty)=(2d)^{-1}N^{-d}[1,...,1]^{\intercal}$
and
$\bm{P}^{(\mathcal{H})}_{\bm{n}_{0}}(\bm{n},t\to\infty)=(6\Omega)^{-1}[1,...,1]^{\intercal}$,
which means that
$P^{(\gamma)}_{n_{01},...,n_{0d}}(n_{1},...,n_{d},t\to\infty)=N^{-d}$ for
$\gamma\in\\{p,r\\}$ and
$P^{(\mathcal{H})}_{n_{01},n_{02}}(n_{1},n_{2},t\to\infty)=\Omega^{-1}$ as
expected from an irreducible (ergodic) aperiodic finite Markov chain with
doubly stochastic transition matrices. When $\gamma=r_{s}$, that is the
squeezed reflecting propagator, in the same limit
$\bm{\lambda}\left(\frac{\pi\kappa}{N-1}\right)^{t}=\bm{0}$ for $\kappa\neq 0$
and $\bm{\lambda}\left(0\right)^{t}=2^{-1}\mathbb{J}$. Hence, the steady state
in this case can be found as
$\bm{P}^{(r_{s})}_{n_{0}}(n,t\to\infty)=(N-1)^{-1}[\alpha_{1},\alpha_{2}]^{\intercal}$
($n\neq\\{1,N\\}$) while
$\bm{P}^{(r_{s})}_{n_{0}}(1,t\to\infty)=\bm{P}^{(r_{s})}_{n_{0}}(N,t\to\infty)=[2(N-1)]^{-1}[\alpha_{1},\alpha_{2}]^{\intercal}$.
In contrast, when either $f$ or $b$ equals 1,
$\bm{\lambda}\left(\pi\mathcal{N}^{(\gamma)}_{\kappa_{1}},...,\pi\mathcal{N}^{(\gamma)}_{\kappa_{d}}\right)$
and $\bm{\lambda}^{(\mathcal{H})}\\!\big{(}\bm{\kappa}\big{)}$ reduce to
generalised permutation matrices, which makes clear that no steady state will
be reached.

When $f=1$, in the periodically bounded domain, one can show that the walker
will be at $n_{0}$ with certainty every $t=\mathfrak{c}N$ steps
($\mathfrak{c}\in\mathbb{N}\cup\\{0\\}$), that is the system is a Markov chain
with periodicity $N$. To do so, let us consider the one-dimensional case (the
analysis is readily extended to arbitrary dimensions), where
$\bm{\lambda}\left(\frac{2\pi\kappa}{N}\right)=\left(\begin{smallmatrix}e^{\frac{2\pi
i\kappa}{N}}&0\\\ 0&e^{\frac{-2\pi i\kappa}{N}}\end{smallmatrix}\right)$.
Consequently,
$\bm{P}_{n_{0}}^{(p)}(n,t)=N^{-1}\sum_{\kappa=0}^{N-1}\exp\left(-\frac{2\pi
i\kappa(n-n_{0})}{N}\right)\left(\begin{smallmatrix}e^{\frac{2\pi i\kappa
t}{N}}&0\\\ 0&e^{\frac{-2\pi i\kappa
t}{N}}\end{smallmatrix}\right)\cdot\left[\begin{matrix}\alpha_{1}\\\
\alpha_{2}\end{matrix}\right]$. Upon inspection, one clearly sees that when
$t=\mathfrak{c}N$ and $n=n_{0}$, the propagator reduces to
$\bm{P}_{n_{0}}^{(p)}(n,t)=\left[\begin{matrix}\alpha_{1}&\alpha_{2}\end{matrix}\right]^{\intercal}$,
such that $P^{(p)}_{n_{0}}(n_{0},\mathfrak{c}N)=1$.

##  Appendix C Defect Technique in Correlated Random Walks

We consider the (unordered) set
$\mathdutchcal{S}=\left\\{(\bm{s}_{i},m_{\bm{s}_{i}}),...,(\bm{s}_{S},m_{\bm{s}_{S}})\right\\}$
of $S$ localised targets. To begin, we write, in arbitrary dimensions and with
arbitrary numbers of internal degrees of freedom, the Master equation for the
state level dynamics

$\displaystyle P(\bm{n},m,t+1)=\sum_{\bm{n}^{\prime}}\sum_{m^{\prime}=1}^{M}$
$\displaystyle\bigg{[}A(\bm{n},m,\bm{n}^{\prime},m^{\prime})P(\bm{n}^{\prime},m^{\prime},t)+\sum_{i=1}^{S}\delta_{\bm{n},\bm{s}_{i}}\delta_{m,m_{\bm{s}_{i}}}(1-\rho_{\bm{s}_{i},m_{\bm{s}_{i}}})A(\bm{s}_{i},m_{\bm{s}_{i}},\bm{n}^{\prime},m^{\prime})P(\bm{n}^{\prime},m^{\prime},t)\bigg{]},$
(30)

where $A(\bm{n},m,\bm{n^{\prime}},m^{\prime})$ is the transition probability
from state $\bm{n}^{\prime},m^{\prime}$ to state $\bm{n},m$, and
$\rho_{\bm{b}_{i},m_{\bm{b}_{i}}}$, $(0<\rho_{\bm{s}_{i},m_{\bm{s}_{i}}}\leq
1)$ governs the probability of getting absorbed at defect
$\bm{b}_{i},m_{\bm{b}_{i}}$ where $\rho_{\bm{s}_{i},m_{\bm{s}_{i}}}=1$
represents perfect trapping efficiency at that site. To proceed, one first
considers the $\rho_{\bm{s}_{i},m_{\bm{s}_{i}}}\neq 1$ case.

Assuming
$\sum_{i=1}^{S}\delta_{\bm{n},\bm{s}_{i}}\delta_{m,m_{\bm{s}_{i}}}(1-\rho_{\bm{s}_{i},m_{\bm{s}_{i}}})A(\bm{s}_{i},m_{\bm{s}_{i}},\bm{n}^{\prime},m^{\prime})P(\bm{n}^{\prime},m^{\prime},t)$,
for any $m$, as a given known function the formal solution of Eq. (30) is the
convolution in space and time of the absorbing propagator with the known terms
[63] and proceeding in the $z$-domain (i.e., taking the generating function)
we find

$\displaystyle\widetilde{P}^{(a)}(\bm{n},m,$ $\displaystyle
z)=\sum_{\bm{n}_{0}^{\prime}}\sum_{m_{0}^{\prime}=1}^{M}\widetilde{P}_{\bm{n}_{0}^{\prime},m_{0}^{\prime}}(\bm{n},m,z)P(\bm{n},m,0)$
$\displaystyle-z\sum_{\bm{n}_{0}^{\prime}}\sum_{m_{0}^{\prime}=1}^{M}\widetilde{P}_{\bm{n}_{0}^{\prime},m_{0}^{\prime}}(\bm{n},m,z)\sum_{i=1}^{S}\rho_{\bm{s}_{i},m_{\bm{s}_{i}}}\sum_{\bm{n}^{\prime}}\sum_{m^{\prime}=1}^{M}A(\bm{s}_{i},m_{\bm{s}_{i}},\bm{n}^{\prime},m^{\prime})\widetilde{P}(\bm{n}^{\prime},m^{\prime},z).$
(31)

By rearranging, and substituting the $z$-transform of the second term on the
RHS of Eq. (30) into the second term of Eq. (31) the generating function of
the formal solution of Eq. (30) is found as

$\displaystyle\widetilde{P}^{(a)}(\bm{n},m,z)$
$\displaystyle=\\!\sum_{\bm{n}_{0}^{\prime}}\sum_{m_{0}^{\prime}=1}^{M}\\!\bigg{[}\widetilde{P}_{\bm{n}_{0}^{\prime},m_{0}^{\prime}}(\bm{n},m,z)P(\bm{n},m,0)\\!+\\!\sum_{i=1}^{S}\frac{\rho_{\bm{s}_{i},m_{\bm{s}_{i}}}}{\rho_{\bm{s}_{i},m_{\bm{s}_{i}}}\\!-\\!1}\widetilde{P}_{\bm{s}_{i},m_{\bm{s}_{i}}}(\bm{n},m,\\!z)\\!\left[\widetilde{P}^{(a)}(\bm{s},m_{\bm{s}_{i}},\\!z)\\!-\\!P^{(a)}(\bm{s},m_{\bm{s}_{i}},0)\right]\\!\bigg{]}\\!.$
(32)

where we have used $\widetilde{P}_{\bm{n}_{0},m_{0}}(\bm{n},m,z)$ to denote
any valid (see the main text) defect free propagator.

An initial condition spatially localised over site $\bm{n}_{0}$ is given by
$P^{(a)}(\bm{n},0)=\delta_{\bm{n},\bm{n}_{0}}\sum_{m=1}^{M}P^{(a)}(\bm{n},m,0)$,
where the contribution by each state is
$P^{(a)}(\bm{n},m,0)=\delta_{m,m_{0}}\alpha_{m_{0}}\left[(1-\rho_{\bm{s}_{i},m_{\bm{s}_{i}}})\delta_{(\bm{n}_{0},m_{0})\in\mathdutchcal{S}}+\delta_{(\bm{n}_{0},m_{0})\notin\mathdutchcal{S}}\right]$.
Substitution of this initial condition into Eq. (32) leads directly to

$\displaystyle\widetilde{P}^{(a)}_{\bm{n}_{0},m_{0}}(\bm{n},m,z)=\alpha_{m_{0}}\widetilde{P}_{\bm{n}_{0},m_{0}}(\bm{n},m,z)+\sum_{i=1}^{S}\frac{\rho_{\bm{s}_{i},m_{\bm{s}_{i}}}}{\rho_{\bm{s}_{i},m_{\bm{s}_{i}}}-1}\alpha_{m_{\bm{s}_{i}}}\widetilde{P}_{\bm{s}_{i},m_{\bm{s}_{i}}}(\bm{n},m,z)\widetilde{P}^{(a)}_{\bm{n}_{0},m_{0}}(\bm{s},m_{\bm{s}_{i}},z).$
(33)

and following the standard defect technique procedure (see the Supplementary
Material [41] for further discussion) i.e., solving via Cramer’s rule [63] and
taking $\rho_{\bm{s}_{i},m_{\bm{s}_{i}}}\to 1$ for all
$(\bm{s}_{i},m_{\bm{s}_{i}})$, we find

$\widetilde{P}^{(a)}_{\bm{n}_{0},m_{0}}(\bm{n},m,z)=\alpha_{m_{0}}\widetilde{P}_{\bm{n}_{0},m_{0}}(\bm{n},m,z)-\sum_{i=1}^{S}\alpha_{m_{\bm{s}_{i}}}\widetilde{P}_{\bm{s},m_{\bm{s}_{i}}}(\bm{n},m,z)\frac{\det[\mathbb{H}^{(i)}(\bm{n}_{0},m_{0},z)]}{\det[\mathbb{H}(z)]},$
(34)

where
$\mathbb{H}(z)_{l,k}=\alpha_{m_{\bm{s}_{k}}}\widetilde{P}_{\bm{s}_{k},m_{\bm{s}_{k}}}(\bm{s}_{l},m_{\bm{s}_{l}},z)$,
$\mathbb{H}(z)_{k,k}=\alpha_{m_{\bm{s}_{k}}}\widetilde{P}_{\bm{s}_{k},m_{\bm{s}_{k}}}(\bm{s}_{k},m_{\bm{s}_{k}},z)$
and $\mathbb{H}^{(j)}(\bm{n}_{0},m_{0},z)$ is the same but with the
$j^{\text{th}}$ column replaced with
$\alpha_{m_{0}}\left[\widetilde{P}_{\bm{n}_{0},m_{0}}(\bm{s}_{1},m_{\bm{s}_{1}},z),...,\widetilde{P}_{\bm{n}_{0},m_{0}}(\bm{s}_{S},m_{\bm{s}_{S}},z)\right]^{\intercal}$,
i.e, they are all known propagators via
$\alpha_{m_{0}}\widetilde{P}_{\bm{n}_{0},m_{0}}(\bm{n},m,z)=\alpha_{m_{0}}\bm{e}_{m}^{\intercal}\cdot\bm{\widetilde{P}}_{\bm{n}_{0}}(\bm{n},z)\cdot\bm{e}_{m_{0}}$.

Since
$\widetilde{P}^{(a)}_{\bm{n}_{0}}(\bm{n},z)=\sum_{m=1}^{M}\sum_{j=1}^{M}\widetilde{P}_{\bm{n}_{0},m_{0_{j}}}^{(a)}(\bm{n},m,z)$,
and
$\widetilde{P}^{(\gamma)}_{\bm{n}_{0}}(\bm{n},z)=\sum_{m=1}^{M}\sum_{j=1}^{M}\alpha_{m_{0_{j}}}\widetilde{P}_{\bm{n}_{0},m_{0_{j}}}^{(\gamma)}(\bm{n},m,z)$,
where there is no $\alpha_{m_{0}}$ multiplier on the absorbing case as it is
already implemented via the initial condition in Eq. (33), one may perform the
double summation over both sides of Eq. (34) to arrive at Eq. (21) of the main
text.

## References

  * [1] R. Nathan, C. T. Monk, R. Arlinghaus, T. Adam, J. Alós, M. Assaf, H. Baktoft, C. E. Beardsworth, M. G. Bertram, A. I. Bijleveld, et al.  Big-data approaches lead to an increased understanding of the ecology of animal movement.  Science, 375(6582):eabg1780, 2022.
  * [2] E. Meijering, O. Dzyubachyk, and I. Smal.  Methods for cell and particle tracking.  Method Enzymol., 504:183–200, 2012.
  * [3] N. Chenouard, I. Smal, F. De Chaumont, M. Maška, I. F. Sbalzarini, Y. Gong, J. Cardinale, C. Carthel, S. Coraluppi, M. Winter, et al.  Objective comparison of particle tracking methods.  Nature methods, 11(3):281–289, 2014.
  * [4] E. W. Montroll and G. H. Weiss.  Random walks on lattices. II.  J. Math. Phys., 6(2):167–181, 1965.
  * [5] G. H. Weiss.  Aspects and applications of the random walk.  Elsevier Science & Technology, 1994.
  * [6] L. Giuggioli.  Exact spatiotemporal dynamics of confined lattice random walks in arbitrary dimensions: A century after Smoluchowski and Pólya.  Phys. Rev. X, 10:021045, 2020.
  * [7] P. Mörters and Y. Peres.  Brownian motion, volume 30.  Cambridge University Press, 2010.
  * [8] P. C. Bressloff.  Stochastic processes in cell biology, volume 41.  Springer, 2014.
  * [9] V.M. Kenkre and L. Giuggioli.  Theory of the spread of epidemics and movement ecology of animals: an interdisciplinary approach using methodologies of physics and mathematics.  Cambridge University Press, 2021.
  * [10] A. Bejan and G. W. Merkx.  Constructal theory of social dynamics.  Springer, 2007.
  * [11] P. Embrechts, C. Klüppelberg, and T. Mikosch.  Modelling extremal events: for insurance and finance, volume 33.  Springer, 2013.
  * [12] D. Selmeczi, S. Mosler, P. H. Hagedorn, N. B. Larsen, and H. Flyvbjerg.  Cell motility as persistent random motion: theories from experiments.  Biophys. J., 89(2):912–931, 2005.
  * [13] B. R. G. Prasad and R. M. Borges.  Searching on patch networks using correlated random walks: space usage and optimal foraging predictions using Markov chain models.  J. Theor. Biol., 240(2):241–249, 2006.
  * [14] I. Echeverría-Huarte, A. Garcimartín, R.C. Hidalgo, C. Martín-Gómez, and I. Zuriguel.  Estimating density limits for walking pedestrians keeping a safe interpersonal distancing.  Sci. Rep., 11(1):1534, 2021.
  * [15] N. Korabel, G. D. Clemente, D. Han, F. Feldman, T. H. Millard, and T. A. Waigh.  Hemocytes in drosophila melanogaster embryos move via heterogeneous anomalous diffusion.  Commun. Phys., 5(1):269, 2022.
  * [16] K. Speckner and M. Weiss.  Single-particle tracking reveals anti-persistent subdiffusion in cell extracts.  Entropy, 23(7):892, 2021.
  * [17] R. Benelli and M. Weiss.  Probing local chromatin dynamics by tracking telomeres.  Biophys. J., 121(14):2684–2692, 2022.
  * [18] A. Gabel and S. Redner.  Random walk picture of basketball scoring.  J. Quant. Anal. Sports, 8(1), 2012.
  * [19] V.M. Kenkre, V. Ern, and A. Fort.  Coherence effects in triplet-exciton transport via time-dependent delayed fluorescence.  Phys. Rev. B, 28(2):598, 1983.
  * [20] V.M. Kenkre and Y.M. Wong.  Effect of transport coherence on trapping: Quantum-yield calculations for excitons in molecular crystals.  Phys. Rev. B, 23(8):3748, 1981.
  * [21] J. Rudnick and G. Gaspari.  Elements of the Random Walk, volume 1.  Cambridge University Press, 2004.
  * [22] S. Jose, D. Mandal, M. Barma, and K. Ramola.  Active random walks in one and two dimensions.  Phys. Rev. E, 105(6):064103, 2022.
  * [23] K. Martens, L. Angelani, R. Di Leonardo, and L. Bocquet.  Probability distributions for the run-and-tumble bacterial dynamics: An analogy to the lorentz model.  Eur. Phys. J. E, 35:1–6, 2012.
  * [24] K. Malakar, V. Jemseena, A. Kundu, K. V. Kumar, S. Sabhapandit, S. N. Majumdar, S. Redner, and A. Dhar.  Steady state, relaxation and first-passage properties of a run-and-tumble particle in one-dimension.  J. Stat. Mech-Theory E., 2018(4):043215, 2018.
  * [25] S. Redner.  A guide to first-passage processes.  Cambridge university press, 2001.
  * [26] R. Metzler, S. Redner, and G. Oshanin.  First-passage phenomena and their applications, volume 35.  World Scientific, 2014.
  * [27] P. C. Bressloff.  Encounter-based model of a run-and-tumble particle.  J. of Stat. Mech. Theory Exp., 2022(11):113206, 2022.
  * [28] L. Angelani.  Run-and-tumble particles, telegrapher’s equation and absorption problems with partially reflecting boundaries.  J. Phys. A-Math. Theor., 48(49):495003, 2015.
  * [29] T. Guérin, N. Levernier, O. Bénichou, and R. Voituriez.  Mean first-passage times of non-Markovian random walkers in confinement.  Nature, 534(7607):356–359, 2016.
  * [30] H. Larralde.  Transport properties of a two-dimensional “chiral” persistent random walk.  Phys. Rev. E, 56(5):5004, 1997.
  * [31] J. Masoliver.  Three-dimensional telegrapher’s equation and its fractional generalization.  Phys. Rev. E, 96(2):022101, 2017.
  * [32] R. Furth.  The brownian motion with consideration of the longevity of the direction of movement.  Z. Phys., 2:244–256, 1920.
  * [33] G. I. Taylor.  Diffusion by continuous movements.  Proc. London Math. Soc., 2(1):196–212, 1922.
  * [34] H. Larralde.  First-passage probabilities and mean number of sites visited by a persistent random walker in one-and two-dimensional lattices.  Phys. Rev. E, 102(6):062129, 2020.
  * [35] E. W. Montroll.  Random walks on lattices. III. Calculation of first-passage times with application to exciton trapping on photosynthetic units.  J. Math. Phys., 10(4):753–765, 1969.
  * [36] V. Tejedor, R. Voituriez, and O. Bénichou.  Optimizing persistent random searches.  Phys. Rev. Lett., 108(8):088103, 2012.
  * [37] H. Shum and E. A. Gaffney.  Hydrodynamic analysis of flagellated bacteria swimming in corners of rectangular channels.  Phys. Rev. E, 92(6):063016, 2015.
  * [38] N. G. van Kampen.  Remarks on non-Markov processes.  Braz. J. Phys., 28:90–96, 1998.
  * [39] J. Masoliver, J. M. Porra, and G. H. Weiss.  Solutions of the telegrapher’s equation in the presence of traps.  Phys. Rev. A, 45(4):2222, 1992.
  * [40] U. Landman and M. F. Shlesinger.  Stochastic theory of multistate diffusion in perfect and defective systems. ii. case studies.  Phys. Rev. B, 19(12):6220, 1979.
  * [41] See supplementary material for further details on the derivations.
  * [42] B. D. Hughes.  Random walks and random environments Volume 1: Random Walks.  Clarendon Press Oxford ; New York, 1995.
  * [43] M. H. Ernst.  Random walks with short memory.  J. Stat. Phys., 53:191–201, 1988.
  * [44] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe.  Active particles in complex and crowded environments.  Rev. Mod. Phys., 88(4):045006, 2016.
  * [45] Z. Kalay.  Effects of confinement on the statistics of encounter times: exact analytical results for random walks in a partitioned lattice.  J. Phys. A-Math. Theor., 45(21):215001, 2012.
  * [46] J. B. Keller.  The scope of the image method.  Commun. Pur. Appl. Math., 6(4):505–512, 1953.
  * [47] J. Masoliver, J. M. Porra, and G. H. Weiss.  Solution to the telegrapher’s equation in the presence of reflecting and partly reflecting boundaries.  Phys. Rev. E, 48(2):939, 1993.
  * [48] E. W. Montroll.  Random walks in multidimensional spaces, especially on periodic lattices.  J. Soc. Ind. Appl. Math., 4(4):241–260, 1956.
  * [49] G. Beylkin, C. Kurcz, and L. Monzón.  Fast algorithms for Helmholtz Green’s functions.  P. Roy. Soc. A-Math. Phy., 464(2100):3301–3326, 2008.
  * [50] S. Chandrasekhar.  Stochastic problems in physics and astronomy.  Rev. Mod. Phys., 15(1):1, 1943.
  * [51] K. Proesmans, R. Toral, and C. Van den Broeck.  Phase transitions in persistent and run-and-tumble walks.  Physica A Stat., 552:121934, 2020.
  * [52] N. R. Smith, P. Le Doussal, S. N. Majumdar, and G. Schehr.  Exact position distribution of a harmonically confined run-and-tumble particle in two dimensions.  Phys. Rev. E, 106(5):054133, 2022.
  * [53] M. Miri and H. Stark.  Persistent random walk in a honeycomb structure: Light transport in foams.  Phys. Rev. E, 68(3):031102, 2003.
  * [54] I. Her.  Geometric transformations on the hexagonal grid.  IEEE Trans. Imag. Process., 4(9):1213–1222, 1995.
  * [55] D. Marris, S. Sarvaharman, and L. Giuggioli.  Exact spatiotemporal dynamics of lattice random walks in hexagonal and honeycomb domains.  Phys. Rev. E, 107(5):054139, 2023.
  * [56] F. Spitzer.  Principles of random walk, volume 34.  SpringerVerlag, second edition, 1976.
  * [57] H. Larralde and F. Leyvraz.  Three-dimensional diffusion with helical persistence.  J. Phys. A-Math. Theor., 48(26):265001, 2015.
  * [58] J. Masoliver, J. M. Porra, and G. H. Weiss.  Some two and three-dimensional persistent random walks.  Physica A: Statistical Mechanics and its Applications, 193(3-4):469–482, 1993.
  * [59] C. Hargus, J. M. Epstein, and K. K. Mandadapu.  Odd diffusivity of chiral random motion.  Phys. Rev. Lett., 127(17):178001, 2021.
  * [60] F. Kümmel, B. Ten Hagen, R. Wittkowski, I. Buttinoni, R. Eichhorn, G. Volpe, H. Löwen, and C. Bechinger.  Circular motion of asymmetric self-propelling particles.  Phys. Rev. Lett., 110(19):198302, 2013.
  * [61] I. H. Riedel, K. Kruse, and J. Howard.  A self-organized vortex array of hydrodynamically entrained sperm cells.  Science, 309(5732):300–303, 2005.
  * [62] P. Fuchter and H. Bloomfield-Gadêlha.  The three-dimensional coarse-graining formulation of interacting elastohydrodynamic filaments and multi-body microhydrodynamics.  J. R. Soc. Interface, 20(202):20230021, 2023.
  * [63] L. Giuggioli and S. Sarvaharman.  Spatio-temporal dynamics of random transmission events: from information sharing to epidemic spread.  J. Phys. A-Math. Theor., 55(37):375005, 2022.
  * [64] J. Abate and W. Whitt.  Numerical inversion of probability generating functions.  Oper. Res. Lett., 12(4):245–251, 1992.
  * [65] J. Abate, G. L. Choudhury, and W. Whitt.  An introduction to numerical transform inversion and its application to probability models.  In Computational probability, pages 257–323. Springer, 2000.
  * [66] S. Sarvaharman and L. Giuggioli.  Closed-form solutions to the dynamics of confined biased lattice random walks in arbitrary dimensions.  Phys. Rev. E, 102(6):062124, 2020.
  * [67] J. Holhouse and S. Redner.  First-passage on disordered intervals.  arXiv preprint arXiv:2307.08879, 2023.
  * [68] O. L. Bonomo and A. Pal.  First passage under restart for discrete space and time: application to one-dimensional confined lattice random walks.  Phys. Rev. E, 103(5):052129, 2021.
  * [69] L. Giuggioli, S. Sarvaharman, D. Das, D. Marris, and T. Kay.  Multi-target search in bounded and heterogeneous environments: a lattice random walk perspective.  arXiv preprint arXiv:2311.00464, 2023.
  * [70] M. Kac.  Random walk and the theory of brownian motion.  Am. Math. Mon., 54(7):369–391, 1947.
  * [71] M. Kulkarni and S. N. Majumdar.  First detection probability in quantum resetting via random projective measurements.  arXiv preprint arXiv:2305.15123, 2023.
  * [72] T. Verechtchaguina, I.M. Sokolov, and L. Schimansky-Geier.  First passage time densities in non-Markovian models with subthreshold oscillations.  Europhys. Lett., 73(5):691, 2006.
  * [73] L. Giuggioli, J.R. Potts, and S. Harris.  Predicting oscillatory dynamics in the movement of territorial animals.  J. R. Soc. Interface., 9(72):1529–1543, 2012.
  * [74] V. M. Kenkre and P. Reineker.  Exciton dynamics in molecular crystals and aggregates.  Springer, 2006.
  * [75] D. Das and L. Giuggioli.  Discrete space-time resetting model: Application to first-passage and transmission statistics.  J. Phys. A-Math. Theor., 55(42):424004, 2022.
  * [76] D. Das and L. Giuggioli.  Dynamics of lattice random walk within regions composed of different media and interfaces.  J. Stat. Mech-Theory E., 2023(1):013201, 2023.
  * [77] S. Sarvaharman and L. Giuggioli.  Particle-environment interactions in arbitrary dimensions: a unifying analytic framework to model diffusion with inert spatial heterogeneities.  Phys. Rev. Res., 5(4):043281, 2023.
  * [78] O. Bénichou and R. Voituriez.  From first-passage times of random walks in confinement to geometry-controlled kinetics.  Phys. Rep., 539(4):225–284, 2014.
  * [79] S. Condamin, O. Bénichou, and M. Moreau.  First-passage times for random walks in bounded domains.  Phys. Rev. Lett., 95(26):260601, 2005.
  * [80] A. Godec and R. Metzler.  First passage time distribution in heterogeneity controlled kinetics: going beyond the mean first passage time.  Sci. Rep-UK, 6(1):20349, 2016.
  * [81] A. Godec and R. Metzler.  Universal proximity effect in target search kinetics in the few-encounter limit.  Phys. Rev. X, 6(4):041037, 2016.
  * [82] J. H. P. Schulz, E. Barkai, and R. Metzler.  Aging renewal theory and application to random walks.  Phys. Rev. X, 4:011028, 2014.
  * [83] S. Godoy and S. Fujita.  Reflection principles for biased correlated walks. simple applications.  J. Math. Phys., 33(9):2998–3003, 1992.
  * [84] V. Rossetto.  The one-dimensional asymmetric persistent random walk.  J. Stat. Mech-Theory E., 2018(4):043204, 2018.
  * [85] R. García-Pelayo.  Solution of the persistent, biased random walk.  Physica A, 384(2):143–149, 2007.
  * [86] T. Kay and L. Giuggioli.  Diffusion through permeable interfaces: fundamental equations and their application to first-passage and local time statistics.  Phys. Rev. Res., 4(3):L032039, 2022.
  * [87] P.C. Bressloff.  A probabilistic model of diffusion through a semi-permeable barrier.  Proc. R. Soc. A, 478(2268):20220615, 2022.
  * [88] T. Kosztołowicz.  Random walk in a discrete and continuous system with a thin membrane.  Phys. A: Stat. Mech. Appl., 298(3-4):285–296, 2001.
  * [89] D. Szász and B. Tóth.  Persistent random walks in a one-dimensional random environment.  J. Stat. Phys., 37(1-2):27–38, 1984.
  * [90] J. M. Hill.  A discrete random walk model for diffusion in media with double diffusivity.  ANZIAM J., 22(1):58–74, 1980.
  * [91] J. M. Hill and B.D. Hughes.  On the general random walk formulation for diffusion in media with diffusivities.  ANZIAM J., 27(1):73–87, 1985.
  * [92] U. Landman and M. F. Shlesinger.  Stochastic theory of multistate diffusion in perfect and defective systems. i. mathematical formalism.  Phys. Rev. B, 19(12):6207, 1979.
  * [93] S. Zhang, A. Chong, and B. D. Hughes.  Persistent exclusion processes: Inertia, drift, mixing, and correlation.  Phys. Rev. E, 100(4):042415, 2019.
  * [94] E. Teomy and R. Metzler.  Transport in exclusion processes with one-step memory: density dependence and optimal acceleration.  J. Phys. A-Math. Theor., 52(38):385001, 2019.
  * [95] E. Gavagnin and C. A. Yates.  Modeling persistence of motion in a crowded environment: The diffusive limit of excluding velocity-jump processes.  Phys. Rev. E, 97(3):032416, 2018.
  * [96] L. Régnier, M. Dolgushev, S. Redner, and O. Bénichou.  Universal exploration dynamics of random walks.  Nat. Comms., 14(1):618, 2023.
  * [97] L. Régnier, M. Dolgushev, and O. Bénichou.  Record ages of non-Markovian scale-invariant random walks.  Nat. Commun., 14(1):6288, 2023.
  * [98] J. C. Cressoni, M. A. A. da Silva, and G.M. Viswanathan.  Amnestically induced persistence in random walks.  Phys, Rev. Lett., 98(7):070603, 2007.
  * [99] G. M. Schütz and S. Trimper.  Elephants can always remember: Exact long-range memory effects in a non-Markovian random walk.  Phys. Rev. E, 70(4):045101, 2004.
  * [100] V.M. Kenkre.  Analytic formulation, exact solutions, and generalizations of the elephant and the alzheimer random walks.  arXiv preprint arXiv:0708.0034, 2007.
  * [101] S. Boettcher, S. Falkner, and R. Portugal.  Renormalization group for quantum walks.  In J. Phys. Conf. Ser., volume 473, page 012018. IOP Publishing, 2013.
  * [102] J. Klafter and I. M. Sokolov.  First steps in random walks: from tools to applications.  OUP Oxford, 2011.
  * [103] G. H. Weiss.  Some applications of persistent random walks and the telegrapher’s equation.  Physica A, 311(3-4):381–410, 2002.
  * [104] D. S. Bernstein.  Matrix mathematics: theory, facts, and formulas.  Princeton University Press, 2009.

Supplemental Materials: Persistent and anti-Persistent Motion in Bounded
Space: Resolution of the First-Passage Problem

## Continuum Limit of the Master Equation

It is straightforward to reduce the two state Master equation i.e., Eq. (1) in
the main text to the telegrapher’s equation when $c=0$ [103], that is when the
walker has no option to stay. Here we show the limiting procedure when
$c\not=0$.

We define the time between successive steps as $\Delta t$ and the spacing
between the lattice sites as $\Delta x$ and following standard procedures, see
e.g., [103], we expand Eq. (1) to first order in $\Delta t$ and $\Delta x$

$\displaystyle Q_{1}+\Delta t\frac{\partial Q_{1}}{\partial t}$
$\displaystyle=f\left(Q_{1}-\Delta x\frac{\partial Q_{1}}{\partial
x}\right)+b\left(Q_{2}-\Delta x\frac{\partial Q_{2}}{\partial
x}\right)+cQ_{1},$ (S1) $\displaystyle Q_{2}+\Delta t\frac{\partial
Q_{2}}{\partial t}$ $\displaystyle=f\left(Q_{2}+\Delta x\frac{\partial
Q_{2}}{\partial x}\right)+b\left(Q_{1}+\Delta x\frac{\partial Q_{1}}{\partial
x}\right)+cQ_{2},$

where $Q_{1}$ and $Q_{2}$ represent the two states. We define a scaling limit
for the probabilities, namely

$f=\alpha-\frac{\Delta t}{2T},\quad b=\frac{\Delta t}{2T}\quad\text{and}\quad
c=1-\alpha$ (S2)

$0\leq\alpha\leq 1$ where $T$ has the dimensions of time and $\lim_{\Delta
x,\Delta t\to 0}\frac{\Delta x}{\Delta t}=v$. After some algebra, dividing
through by $\Delta t$ and letting $\Delta t$ and $\Delta x$ tend to zero, we
find

$\displaystyle\frac{\partial Q_{1}}{\partial t}$
$\displaystyle=\frac{1}{2T}(Q_{2}-Q_{1})-\alpha v\frac{\partial
Q_{1}}{\partial x}$ (S3) $\displaystyle\frac{\partial Q_{2}}{\partial t}$
$\displaystyle=\frac{1}{2T}(Q_{1}-Q_{2})+\alpha v\frac{\partial
Q_{2}}{\partial x}.$

Defining $U=Q_{1}+Q_{2}$ and $W=Q_{1}-Q_{2}$ and summing and subtracting the
two lines of Eq. (S3), one obtains

$\frac{\partial^{2}U}{\partial t^{2}}=-\frac{1}{T}\frac{\partial U}{\partial
t}+\alpha^{2}v^{2}\frac{\partial^{2}U}{\partial x^{2}},$ (S4)

which shows that $\alpha=1-c$ is simply a rescaling of the velocity $v$.

## Derivation of the Unbounded Propagator

Analogously to the $d=1$ case, in higher dimensions one writes the Master
equation for the dynamics into each state. As an example we show explicitly
the $d=2$ ($Z=4$) case, which gives the following set of four coupled
equations

$\displaystyle Q(\bm{n},1,t\\!+\\!1)$
$\displaystyle=fQ(n_{1}-1,n_{2},1,t)\\!+\\!bQ(n_{1}-1,n_{2},2,t)\\!+\\!\ell
Q(n_{1}-1,n_{2},3,t)\\!+\\!\ell
Q(n_{1}-1,n_{2},4,t)\\!+\\!c^{(2)}Q(n_{1},n_{2},1,t),$ (S5) $\displaystyle
Q(\bm{n},2,t\\!+\\!1)$
$\displaystyle=bQ(n_{1}+1,n_{2},1,t)\\!+\\!fQ(n_{1}+1,n_{2},2,t)\\!+\\!\ell
Q(n_{1}+1,n_{2},3,t)\\!+\\!\ell
Q(n_{1}+1,n_{2},4,t)\\!+\\!c^{(2)}Q(n_{1},n_{2},2,t),$ $\displaystyle
Q(\bm{n},3,t\\!+\\!1)$ $\displaystyle=\ell Q(n_{1},n_{2}-1,1,t)\\!+\\!\ell
Q(n_{1},n_{2}-1,2,t)\\!+\\!fQ(n_{1},n_{2}-1,3,t)\\!+\\!bQ(n_{1},n_{2}-1,4,t)\\!+\\!c^{(2)}Q(n_{1},n_{2},3,t),$
$\displaystyle Q(\bm{n},4,t\\!+\\!1)$ $\displaystyle=\ell
Q(n_{1},n_{2}+1,1,t)\\!+\\!\ell
Q(n_{1},n_{2}+1,2,t)\\!+\\!bQ(n_{1},n_{2}+1,3,t)\\!+\\!fQ(n_{1},n_{2}+1,4,t)\\!+\\!c^{(2)}Q(n_{1},n_{2},4,t),$

where the movement parameters are defined in the first paragraph of Sec. IV of
the main text.

From Eq. (S5) it is straightforward to write a Master equation in matrix form
as

$\displaystyle\bm{Q}(\bm{n},t+1)=$ (S6)
$\displaystyle\mathbb{A}^{(1)}_{4}\cdot\bm{Q}(n_{1}-1,n_{2},t)+\mathbb{B}^{(1)}_{4}\cdot\bm{Q}(n_{1}+1,n_{2},t)+\mathbb{A}^{(2)}_{4}\cdot\bm{Q}(n_{1},n_{2}-1,t)+\mathbb{B}^{(2)}_{4}\cdot\bm{Q}(n_{1},n_{2}+1,t)+\mathbb{C}_{4}\cdot\bm{Q}(n_{1},n_{2},t),$

where

$\mathbb{A}^{(1)}_{4}=\begin{bmatrix}f&b&\ell&\ell\\\ 0&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ \end{bmatrix},\;\mathbb{B}^{(1)}_{4}=\begin{bmatrix}0&0&0&0\\\
b&f&\ell&\ell\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{bmatrix},\;\mathbb{A}^{(2)}_{4}=\begin{bmatrix}0&0&0&0\\\ 0&0&0&0\\\
\ell&\ell&f&b\\\ 0&0&0&0\\\
\end{bmatrix},\;\mathbb{B}^{(2)}_{4}=\begin{bmatrix}0&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ \ell&\ell&b&f\\\
\end{bmatrix}\;\text{and}\;\mathbb{C}_{4}=c^{(2)}\mathbb{I}.$ (S7)

In fact, this procedure can be generalised to arbitrary dimension $d$, and one
can write a matricial Master equation, also reported in Eq. (13) of the main
text, of the form

$Q(\bm{n},t+1)=\sum_{i=1}^{d}\left\\{\mathbb{A}_{2d}^{(i)}\cdot
Q(\bm{n}-\bm{e}_{i},t)+\mathbb{B}_{2d}^{(i)}\cdot
Q(\bm{n}+\bm{e}_{i},t)\right\\}+\mathbb{C}_{2d}\cdot Q(\bm{n},t),$ (S8)

where $\bm{e}_{i}$ is a unit vector along dimension $i$.

While the dynamics across the states is dependent on the internal degree of
freedom, the site dynamics are translationally invariant. This allows one to
apply the discrete Fourier transform
$\widehat{f}(\xi)=\sum_{n=-\infty}^{\infty}e^{in\xi}f(n)$, and we also
introduce the $z$-transform $\widetilde{f}(z)=\sum_{t=0}^{\infty}z^{t}f(t)$,
to find the unbounded solution in Fourier and $z$-domain as [42]

$\widehat{\widetilde{\bm{Q}}}(\xi_{1},...,\xi_{d},z)=\big{[}\mathbb{I}-z\bm{\lambda}(\xi_{1},...,\xi_{d})\big{]}^{-1}\cdot\widehat{\bm{Q}}(\xi_{1},...,\xi_{d},0).$
(S9)

Introducing the Fourier transform of the initial condition, namely
$\sum_{n_{1}=-\infty}^{\infty}...\sum_{n_{d}=-\infty}^{\infty}e^{i\bm{\xi}^{\intercal}\cdot\bm{n}}\left(\prod_{i=1}^{d}\delta_{n,n_{0_{i}}}\right)\bm{U_{m_{0}}}=e^{i\bm{\xi}^{\intercal}\cdot\bm{n}_{0}}\bm{U_{m_{0}}}$,
with $\bm{\xi}=(\xi_{1},...,\xi_{d})$, into Eq. (S9), we find

$\widehat{\widetilde{\bm{Q}}}_{\bm{n}_{0}}(\xi_{1},...,\xi_{d},z)=e^{i\bm{\xi}^{\intercal}\cdot\bm{n}_{0}}\left[\mathbb{I}-z\bm{\lambda}(\xi_{1},...,\xi_{d})\right]^{-1}\cdot\bm{U_{m_{0}}},$
(S10)

which gives through an inverse Fourier transformation the formal lattice
Green’s function

$\widetilde{\bm{Q}}_{\bm{n}_{0}}(\bm{n}_{1},...,\bm{n}_{d},z)=\frac{1}{(2\pi)^{d}}\int_{-\pi}^{\pi}...\int_{-\pi}^{\pi}e^{-i\bm{\xi}^{\intercal}\cdot(\bm{n}-\bm{n}_{0})}\left[\mathbb{I}-z\bm{\lambda}(\xi_{1},...,\xi_{d})\right]^{-1}\cdot\bm{U_{m_{0}}}\text{d}^{d}\xi,$
(S11)

where the structure of the matrix $\bm{\lambda}(\xi_{1},...,\xi_{d})$ is
discussed in Appendix B of the main text. As noted in ref. [42], to obtain the
state level probability one needs only to integrate over the specific element
in the matrix in Eq. (S11), e.g.,
$\widetilde{Q}_{\bm{n}_{0},m_{0}}(\bm{n}_{1},...,\bm{n}_{d},m,z)=(2\pi)^{-d}\int_{-\pi}^{\pi}...\int_{-\pi}^{\pi}\bm{e}_{m}^{\intercal}\cdot\widetilde{\bm{Q}}_{\bm{n}_{0}}(\bm{n}_{1},...,\bm{n}_{d},z)\cdot\bm{e}_{m_{0}}\text{d}^{d}\xi$.

We note here that other representations of the unbounded occupation
probability viz. Eqs. (S10) and (S11) are known. For example, in ref. [43], a
form of the propagator is given where the walk is decomposed into the movement
steps along each spatial dimension and the full occupation probability is
obtained via a summation over each spatial direction. Such a summation yields
the same outcome as performing
$\bm{Q}_{\bm{n}_{0}}(\bm{n}_{1},...,\bm{n}_{d},z)=\sum_{i=1}^{2d}\widetilde{\bm{Q}}_{\bm{n}_{0}}(\bm{n}_{1},...,\bm{n}_{d},z)_{i}$
in the present context. However, the present approach is more flexible as it
is able to discern exactly the state level occupation probability, easily
include a sojourn probability $c$ and an arbitrary distribution over
$\bm{U_{m_{0}}}$, with very little algebraic burden.

### Closed-form solution for $d=1$ case:

In one dimension it is possible to calculate explicitly the integral in Eq.
(S11). With $\bm{U_{m_{0}}}=[\alpha_{1},\alpha_{2}]^{\intercal}$ and
$\bm{\lambda}(\xi)$ defined in Eq. (4) of the main text one has

$\widetilde{\bm{Q}}_{n_{0}}(n,z)=\int_{-\pi}^{\pi}\frac{e^{-i\xi(n-n_{0})}}{1-2z\big{[}c-f\cos(\xi)\big{]}+z^{2}\big{[}f^{2}-b^{2}+c^{2}+2fc\cos(\xi)\big{]}}\begin{bmatrix}1-z\left(c+fe^{-i\xi}\right)&bze^{i\xi}\\\
bze^{-i\xi}&1-z\left(c+fe^{i\xi}\right)\end{bmatrix}\cdot\begin{bmatrix}\alpha_{1}\\\
\alpha_{2}\end{bmatrix}\text{d}\xi.$ (S12)

which is the explicit version Eq. (3) in the main text. By rewriting Eq. (S12)
as

$\widetilde{\bm{Q}}_{n_{0}}(n,z)=\int_{-\pi}^{\pi}\frac{e^{-i\xi(n-n_{0})}}{(1-2cz)\big{[}1-z\varepsilon\left(e^{i\xi}+e^{-i\xi}\right)+\delta
z^{2}\big{]}}\begin{bmatrix}1-z\left(c+fe^{-i\xi}\right)&bze^{i\xi}\\\
bze^{-i\xi}&1-z\left(c+fe^{i\xi}\right)\end{bmatrix}\cdot\begin{bmatrix}\alpha_{1}\\\
\alpha_{2}\end{bmatrix}\text{d}\xi,$ (S13)

and defining

$\displaystyle r_{\pm}(z)=\frac{1+\delta z^{2}}{2\varepsilon
z}\left(1\pm\sqrt{1-\frac{4\varepsilon^{2}z^{2}}{[1+\delta
z^{2}]^{2}}}\right),$ (S14)

where $r_{+}=r_{-}^{-1}$, $|r_{-}|<1<|r_{+}|$,
$\varepsilon=f(1-cz)(1-2cz)^{-1}$ and $\delta=(f^{2}-b^{2}+c^{2})(1-2cz)^{-1}$
one can follow closely to the solution of a similar integral [83] and find the
explicitly generating functions for $\widetilde{Q}_{n_{0}}(n,1,z)$ and
$\widetilde{Q}_{n_{0}}(n,2,z)$. Namely,

$\begin{bmatrix}\widetilde{Q}_{n_{0},1}(n,1,z)\\\
\widetilde{Q}_{n_{0},1}(n,2,z)\end{bmatrix}=\frac{\alpha_{1}}{(1-2cz)\sqrt{\left(1+\delta
z^{2}\right)^{2}-4\varepsilon^{2}z^{2}}}\begin{bmatrix}\frac{1-cz}{r_{+}(z)^{|n-n_{0}|}}-\frac{az}{r_{+}(z)^{|n-n_{0}+1|}}\\\
\frac{bz}{r_{+}(z)^{|n-n_{0}+1|}}\end{bmatrix}$ (S15)

when the walker is assumed to have reached $n_{0}$ from the left at $t=0$, and

$\begin{bmatrix}\widetilde{Q}_{n_{0},2}(n,1,z)\\\
\widetilde{Q}_{n_{0},2}(n,2,z)\end{bmatrix}=\frac{\alpha_{2}}{(1-2cz)\sqrt{\left(1+\delta
z^{2}\right)^{2}-4\varepsilon^{2}z^{2}}}\begin{bmatrix}\frac{bz}{r_{+}(z)^{|n-n_{0}-1|}}\\\
\frac{1-cz}{r_{+}(z)^{|n-n_{0}|}}-\frac{az}{r_{+}(z)^{|n-n_{0}-1|}},\end{bmatrix}.$
(S16)

when it reached from the right.

The solution at each state is irrespective of how the walker arrived at
$n_{0}$ at time $t=0$ and we sum the rows in Eqs. (S15) and (S16) to find

$\displaystyle\begin{bmatrix}\widetilde{Q}_{n_{0}}(n,1,z)\\\
\widetilde{Q}_{n_{0}}(n,2,z)\end{bmatrix}=\frac{1}{(1-2cz)\sqrt{\left(1+\delta
z^{2}\right)^{2}-4\varepsilon^{2}z^{2}}}\begin{bmatrix}\frac{\alpha_{1}(1-cz)}{r(z)^{|n-n_{0}|}}-\frac{\alpha_{1}fz}{r(z)^{|n-n_{0}+1|}}+\frac{\alpha_{2}bz}{r(z)^{|n-n_{0}-1|}}\\\
\frac{\alpha_{2}(1-cz)}{r(z)^{|n-n_{0}|}}-\frac{\alpha_{2}fz}{r(z)^{|n-n_{0}-1|}}+\frac{\alpha_{1}bz}{r(z)^{|n-n_{0}+1|}}\end{bmatrix},$
(S17)

where for simplicity we have dropped the subscript + from the definition of
$r_{+}$. A summation of the two states then gives Eq. (5) of the main text.

## Periodic Boundary Conditions

The periodic images in arbitrary dimensions are written as [42]

$\displaystyle\widetilde{\bm{P}}^{(p)}_{n_{0_{1}},...,n_{0_{d}}}(n_{1},...,n_{d},z)=\sum_{\kappa_{1}=-\infty}^{\infty}...\sum_{\kappa_{d}=-\infty}^{\infty}\widetilde{\bm{Q}}_{n_{0_{1}}+\kappa_{1}N_{1},...,n_{0_{d}}+\kappa_{d}N_{d}}(n_{1},...,n_{d},z),$
(S18)

inserting this into Eq. (S11) we find that

$\widetilde{\bm{P}}^{(p)}_{n_{0_{1}},...,n_{0_{d}}}(n_{1},...,n_{d},z)=\frac{1}{(2\pi)^{d}}\int_{-\pi}^{\pi}...\int_{-\pi}^{\pi}\sum_{\kappa_{1}=-\infty}^{\infty}...\sum_{\kappa_{d}=-\infty}^{\infty}e^{-i(\bm{n}^{\intercal}-\bm{n}_{0}^{\intercal}+N\bm{\kappa}^{\intercal})\cdot\bm{\xi}}\left[\mathbb{I}-z\bm{\lambda}(\xi_{1},...,\xi_{d})\right]^{-1}\cdot\bm{U_{m_{0}}}\text{d}^{d}\xi.$
(S19)

Using the $d$-dimensional Poisson summation formula

$\sum_{\kappa_{1}=-\infty}^{\infty}...\sum_{\kappa_{d}=-\infty}^{\infty}\mu_{d}\delta(\bm{\xi}-\bm{\kappa}^{\intercal}\cdot\bm{\mu})=\sum_{\kappa_{1}=-\infty}^{\infty}...\sum_{\kappa_{d}=-\infty}^{\infty}e^{-i\bm{\kappa}^{\intercal}\cdot\bm{N}\bm{\xi}}$
(S20)

where $\mu_{d}=\frac{2\pi}{N_{d}}$, one finds

$\widetilde{\bm{P}}^{(p)}_{n_{0_{1}},...,n_{0_{d}}}(n_{1},...,n_{d},z)=\frac{1}{N^{d}}\int_{-\pi}^{\pi}...\int_{-\pi}^{\pi}\sum_{\kappa_{1}=-\infty}^{\infty}...\sum_{\kappa_{d}=-\infty}^{\infty}e^{-i(\bm{n}^{\intercal}-\bm{n_{0}}^{\intercal})\cdot\bm{\xi}}\delta(\bm{\xi}-\bm{\kappa}^{\intercal}\cdot\bm{\mu})\left[\mathbb{I}-z\bm{\lambda}(\xi_{1},...,\xi_{d})\right]^{-1}\cdot\bm{U_{m_{0}}}\text{d}^{d}\xi.$
(S21)

Since the region of integration is over the Brillouin zone of the lattice
sites, the singularities of
$\delta(\bm{\xi}-\bm{\kappa}^{\intercal}\cdot\bm{\mu})$ lie in the integral
region if and only if $\bm{\kappa}$ is a site in the lattice (for a general
discussion of solving integrals of this type see [42]). Solving the integral
in Eq. (S21) using this property leads straightforwardly to

$\widetilde{\bm{P}}^{(p)}_{n_{0_{1}},...,n_{0_{d}}}(n_{1},...,n_{d},z)=\frac{1}{N^{d}}\sum_{\kappa_{1}=0}^{N-1}...\sum_{\kappa_{d}=0}^{N-1}e^{-i(\bm{n}^{\intercal}-\bm{n_{0}}^{\intercal})\cdot\frac{2\pi}{N}\bm{\kappa}}\left[\mathbb{I}-z\bm{\lambda}\left(\frac{2\pi\kappa_{1}}{N},...,\frac{2\pi\kappa_{d}}{N}\right)\right]^{-1}\cdot\bm{U_{m_{0}}},$
(S22)

which after taking the inverse $z$-transform via the identity
$\sum_{t=0}^{\infty}z^{t}\bm{\lambda}(\bm{\xi})^{t}=\left[\mathbb{I}-z\bm{\lambda}(\bm{\xi})\right]$,
one finds Eq. (14) of the main text.

In the $1d$ case, one finds, explicitly, the probability over the entire site
as

$\displaystyle P^{(p)}_{n_{0}}(n,t)=$
$\displaystyle\frac{1}{N}\sum_{\kappa=0}^{N-1}\exp\left(\frac{-2\pi
i\kappa(n-n_{0})}{N}\right)\frac{be^{-\frac{2\pi
i\kappa}{N}}}{2\sqrt{b^{2}-f^{2}+f^{2}\cos^{2}\left(\frac{2\pi\kappa}{N}\right)}}\times$
(S23)
$\displaystyle\left[(\lambda_{-})^{t}\left(1+v_{-}\left(\frac{2\pi\kappa}{N}\right)\right)\left(\alpha_{2}v_{+}\left(\frac{2\pi\kappa}{N}\right)-\alpha_{1}\right)+(\lambda_{+})^{t}\left(1+v_{+}\left(\frac{2\pi\kappa}{N}\right)\right)\left(\alpha_{2}v_{-}\left(\frac{2\pi\kappa}{N}\right)-\alpha_{1}\right)\right],$

where
$v_{\pm}\left(\frac{2\pi\kappa}{N}\right)=\frac{1}{2b}\left[-f-fe^{\frac{4\pi\kappa}{N}}\pm
2e^{\frac{2\pi\kappa}{N}}\sqrt{b^{2}-f^{2}+f^{2}\cos^{2}\left(\frac{2\pi\kappa}{N}\right)}\right]$,
and $\lambda_{\pm}$ are defined in Eq. (9) of the main text. Equation (S23)
reduces to Eq. (7) of the main text when $\alpha_{1}=\alpha_{2}$.

## Reflective Boundary Conditions

We detail here the corresponding procedure for the reflective boundary
condition.

### Derivation of the $1d$ Semi-bounded Propagator

The image set used, for the semi-infinite propagator, is dependent on whether
the reflection is to the left or the right of the initial condition. We
assume, without loss of generality that the barrier is placed between sites
$n=0$ and $n=1$ meaning $P(1,\\{2,1\\},t)=P(0,\\{1,2\\},t)$.

Therefore, in the symmetric $\alpha_{1}=\alpha_{2}=\frac{1}{2}$ case, the
semi-bounded propagator,
$\widetilde{\bm{P}}^{(\mathdutchcal{r})}_{n_{0}}(n,z)$, can be written as

$\widetilde{\bm{P}}^{(\mathdutchcal{r})}_{n_{0}}(n,z)=\widetilde{\bm{Q}}_{n_{0}}(n,z)+\widetilde{\bm{Q}}_{-n_{0}+1}(n,z),$
(S24)

which is written, after summing the components, as

$\displaystyle\widetilde{P}^{(\mathdutchcal{r})}_{n_{0}}(n,z)=\frac{1}{(1-2cz)\sqrt{\left(1+\delta
z^{2}\right)^{2}-4\varepsilon^{2}z^{2}}}\bigg{[}$
$\displaystyle\frac{1-cz}{r(z)^{|n-n_{0}|}}+\frac{(b-f)z}{2r(z)^{|n-n_{0}+1|}}+\frac{(b-f)z}{2r(z)^{|n-n_{0}-1|}}+$
(S25)
$\displaystyle\frac{1-cz}{r(z)^{|n+n_{0}-1|}}+\frac{(b-f)z}{2r(z)^{|n+n_{0}|}}+\frac{(b-f)z}{2r(z)^{|n+n_{0}-2|}}\bigg{]}.$

In the ($b=0$) ballistic limit, the boundary condition can be simplified to
$P(0,1,t+1)=P(1,2,t)$ as there will never be an instance of the walker,
entering state $m=1$ and then crossing the boundary at a later time. In fact,
the trajectories in state $m=1$, will remain in $m=1$ for all time, meaning
the image of this part of the dynamics will deterministically travel towards
minus infinity and never interact with the domain of interest. As such, one
can write the propagator in this case as

$\displaystyle\widetilde{P}^{(\mathdutchcal{r})}_{n_{0}}(n,z)$
$\displaystyle=\widetilde{Q}_{n_{0}}(n,1,z)+\widetilde{Q}_{n_{0}}(n,2,z)+\widetilde{Q}_{-n_{0}+1}(n,2,z)=\frac{1}{(1-2cz)\sqrt{\left(1+\delta
z^{2}\right)^{2}-4\varepsilon^{2}z^{2}}}\bigg{[}\frac{1-cz}{r(z)^{|n-n_{0}|}}$
$\displaystyle+\frac{1-cz}{r(z)^{|n+n_{0}-1|}}-\frac{zf}{2}\bigg{(}\frac{1}{r(z)^{|n-n_{0}+1|}}+\frac{1}{r(z)^{|n-n_{0}-1|}}+\frac{1}{r(z)^{|n+n_{0}|}}\bigg{)}\bigg{]}.$
(S26)

### Fully Bounded Propagator

We derive here the propagator in arbitrary dimensions. Employing the method of
images for the reflective case [6],

$\displaystyle\widetilde{\bm{P}}^{(r)}_{n_{0_{1}},...,n_{0_{d}}}(n_{1},...,n_{d},z)=\sum_{\kappa_{1}=-\infty}^{\infty}...\sum_{\kappa_{d}=-\infty}^{\infty}$
$\displaystyle\bigg{\\{}\widetilde{\bm{Q}}_{n_{0_{1}}+2\kappa_{1}N_{1},...,n_{0_{d}}+2\kappa_{d}N_{d}}(n_{1},...,n_{d},z)$
(S27)