# The spectral boundary of the Asymmetric Simple Exclusion Process (ASEP) -
free fermions, Bethe ansatz and random matrix theory

Goran Nakerst  Institut für Theoretische Physik, Technische Universität
Dresden, D-01062 Dresden, Germany  Department of Theoretical Physics, Maynooth
University, Co. Kildare, Ireland  Tomaž Prosen  University of Ljubljana,
Faculty for Mathematics and Physics, Jadranska 19, Ljubljana, Slovenia
Masudul Haque  Institut für Theoretische Physik, Technische Universität
Dresden, D-01062 Dresden, Germany  Department of Theoretical Physics, Maynooth
University, Co. Kildare, Ireland  Max-Planck-Institut für Physik komplexer
Systeme, D-01187 Dresden, Germany

###### Abstract

In non-equilibrium statistical mechanics, the Asymmetric Simple Exclusion
Process (ASEP) serves as a paradigmatic example. We investigate the spectral
characteristics of the ASEP, focusing on the spectral boundary of its
generator matrix. We examine finite ASEP chains of length $L$, under periodic
(pbc) and open boundary conditions (obc). Notably, the spectral boundary
exhibits $L$ spikes for pbc and $L+1$ spikes for obc. Treating the ASEP
generator as an interacting non-Hermitian fermionic model, we extend the model
to have tunable interaction. In the non-interacting case, the analytically
computed many-body spectrum shows a spectral boundary with prominent spikes.
For pbc, we use the coordinate Bethe ansatz to interpolate between the
noninteracting case to the ASEP limit, and show that these spikes stem from
clustering of Bethe roots. The robustness of the spikes in the spectral
boundary is demonstrated by linking the ASEP generator to random matrices with
trace correlations or, equivalently, random graphs with distinct cycle
structures, both displaying similar spiked spectral boundaries.

##  I Introduction

The Asymmetric Simple Exclusion Process (ASEP) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
11] is a well-studied paradigmatic stochastic many-body model that has been
used to understand a wide range of non-equilibrium phenomena. This paper
explores the spectral boundary of the markov matrix (the generator of ASEP),
with a focus on a characteristic spiky formation, by establishing connections
between the ASEP, non-interacting fermions, and random matrices featuring
trace correlations.

The ASEP model has proven instrumental in shedding light on phenomena like
non-equilibrium phase transitions [12, 5, 13, 14, 15, 16], and shock formation
[17, 18, 19, 16, 20, 21], among others. Its versatility extends across various
domains, such as protein synthesis [22, 23, 24], intracellular transport [24,
25, 26], traffic flows [27], and quantum dots [28]. Another major incentive
for its study is the association of the ASEP with interface dynamics and its
connection to the Kardar-Parisi-Zhang equation in 1D (or equivalent noisy
Burgers’ equation) [29, 30, 31, 32].

The ASEP is a model where particles move stochastically on a one-dimensional
lattice, adhering to exclusion interactions that restrict each site to a
single particle, mirroring volume exclusion in real systems. Particles move to
adjacent sites only if these sites are unoccupied. The process is termed
’asymmetric’ due to the unequal probabilities for particle movement to the
left or right, leading to directional bias. In cases where movement is limited
to one direction, the model is referred to as Totally Asymmetric Simple
Exclusion Process (TASEP).

A probability vector $P$ of particle configurations evolves according to the
equation

$\frac{d}{dt}P(t)=HP(t),$ (1)

where $H$ is the generator matrix that governs the dynamics of the system.
This markov (stochastic) matrix is a cornerstone of our study as it
encapsulates all the dynamical information of the ASEP. The spectrum of $H$ is
particularly insightful: it informs us about the various rates at which
different states of the system evolve, which is crucial for understanding how
the system approaches its steady state.

The asymmetry of the ASEP implies that the matrix $H$ is non-hermitian and its
eigenvalues are generally complex. The real part of these eigenvalues relates
to the relaxation times of eigenmodes, indicating how quickly the system
returns to the steady state after a disturbance. The imaginary part, on the
other hand, determines the oscillatory behavior of the system, setting the
time scales of periodic or quasi-periodic patterns in the system evolution.

In this paper, we focus on finite chains of length $L$ and either periodic
(pbc) or open boundary conditions (obc). The finite-dimensional nature of $H$
in these cases leads to a discrete and bounded spectrum. Analyzing this
spectrum, especially establishing tight bounds on it, provides valuable
insights into the aforementioned time scales and the overall dynamical
properties of the system.

Our primary objective is to investigate and explain an intriguing feature of
the shape of the spectral boundary, namely, the prominent spikes clearly seen
in Figs. 1(a,b) and also in previous studies [33, 34]. The formation of these
spikes – $L$ spikes for pbc and $L+1$ for obc – present a fascinating aspect
of the spectral characteristics of the ASEP. Unraveling the mechanisms behind
the formation of these spikes in the spectral boundary is a major focus of
this work. We elucidate the emergence of spectral spikes through three
approaches.

Firstly, the generator matrix $H$ is modeled as an interacting, non-Hermitian,
spinless fermion system with interaction strength $U=1$. For $U=0$, $H$
reduces to a non-interacting fermion model. Although this is not a Markov
matrix, it is instructive to study the $U=0$ case as it is solvable as a non-
Hermitian free-fermion Hamiltonian. (We refer to this as the “non-interacting
ASEP”.) The many-body spectrum of $H$ in this case, expressible as sums of
single-particle eigenvalues on ellipses (circles for TASEP) in the complex
plane, exhibits $L$ spikes ($L+1$ for obc) at its spectral boundary.

In the second approach, we extend the coordinate Bethe ansatz method,
traditionally used for calculating the spectrum of $U=1$ with pbc [32], to
encompass arbitrary interaction strengths $U$. For TASEP, the many-body
spectrum is constituted by sums of Bethe roots, which exhibit an elliptical
clustering in the complex plane within the range $0\leq U\leq 1$. By focusing
on the cluster sizes and disregarding finer Bethe root details, we demonstrate
that the spectral boundary, akin to the $U=0$ case, is defined by sums of
Bethe roots from neighboring clusters, resulting in a prominent display of $L$
spikes.

Lastly, we underscore the resilience of these spiky spectral boundaries by
relating the TASEP to a random graph ensemble. In TASEP, the number of updates
required to revert to a specific configuration is a multiple of $L$ ($L+1$ for
obc) [33]. We examine random graphs wherein all cycle lengths are divisible by
$L$ ($L+1$ for obc). Our findings reveal that the spectral boundaries of both
the adjacency matrix (analogous to $U=0$ in TASEP) and the Laplacian matrices
(corresponding to $U=1$ in TASEP) of this random graph ensemble are
characterized by the presence of $L$ ($L+1$) spikes.

The resilience of the spiky spectral boundary is noteworthy. This feature,
inherent in the non-interacting fermion model, remarkably withstands the
reintroduction of interactions. Furthermore, it prevails even when all aspects
of $H$ are disregarded, except for the cycle lengths in the many-body graph.

The article is organized as follows: In Sec. II we introduce the generator
matrix of ASEP with pbc and obc. In Sec. III and Sec. IV we present results of
the non-interacting ASEP ($U=0$) with pbc and obc, respectively. In Sec. V we
investigate the interacting TASEP ($0\leq U\leq 1$) with pbc by Bethe ansatz.
In Sec. VI we compare TASEP to random graphs with the aforementioned cycle
structure. We conclude in Sec. VII. Appendices A and B provide additional
information on solving the non-interacting TASEP with obc. Appendix C details
the derivation of Bethe equations for any $U$ with pbc, and Appendix D
presents numerical specifics for solving these equations to determine the full
spectrum of the generator matrix $H$.

##  II Generator matrix of ASEP

Figure 1: Spectrum of the generator matrix $H$ of TASEP (a,b) and the non-
interacting TASEP (c,d) on $L=11$ sites. The spectrum shows $L$ spikes in
(a,c) for pbc with $N=5$ particles and $L+1$ spikes in (b,d) for obc. Red
solid lines in (c,d) denote the spectral boundary according to Eq. (24).

In this section, we will introduce the generator matrix $H$ of the ASEP for
pbc and obc as non-Hermitian fermion models, along with essential notation.

We consider ASEP chains of length $L$. The number of particles in the chain is
denoted by $N$ and the particle density by $\rho=N/L$. The probability for a
particle to hop right or left in time $dt$ is $p\,dt$ or $q\,dt$,
respectively, with the convention $p+q=1$ unless specified otherwise.

Let us introduce

$H=H_{U}=H_{0}+U\mathcal{I},$ (2)

where $H_{0}$ is a matrix with non-negative off-diagonal elements and zero
diagonal entries and $\mathcal{I}$ is a diagonal matrix. The term $U$ denotes
the interaction strength. In the ASEP context, $H_{0}$ represents a non-
interacting spinless fermion model, and $\mathcal{I}$ a 4-point (2-body)
fermion interaction.

The generator of the ASEP is $H=H_{1}$ with interaction strength $U=1$. Here,
$H_{1}$ is the generator of a stochastic Markov process and a stochastic
matrix, where the sums of all columns of $H_{1}$ equal zero. This property is
ensured by the diagonal elements of $\mathcal{I}$ equalling the sums of the
corresponding columns of $H_{0}$,

$\mathcal{I}_{jj}=\sum_{k}\left(H_{0}\right)_{kj}.$ (3)

Whenever $U\neq 1$, $H_{U}$ ceases to be a Markov matrix and does not generate
the ASEP or any other stochastic process.

Studying $H_{U}$ with $U\neq 1$ could elucidate the $U=1$ case for two
reasons. First, the analyticity of $H_{U}$ in $U$ suggests that its properties
at $U\neq 1$ could be extrapolated to $U=1$. Second, the diagonal matrix
$\mathcal{I}$ exists only to ensure the Markov property of $H_{1}$ and,
according to Eq. (3), is entirely determined by $H_{0}$. Therefore, ignoring
$\mathcal{I}$ in the $U=0$ case likely retains some features of the Markov
matrix $H_{1}$.

###  II.1 Periodic Boundary Conditions

For pbc the matrices $H_{0}$ and $\mathcal{I}$ are given by

$\displaystyle H_{0}$
$\displaystyle=\sum_{j=1}^{L}\left(p\sigma_{j+1}^{+}\sigma_{j}^{-}+q\sigma_{j}^{+}\sigma_{j+1}^{-}\right)$
(4) $\displaystyle\mathcal{I}$
$\displaystyle=\frac{1}{4}\sum_{j=1}^{L}\left(\sigma_{j}^{z}\sigma_{j+1}^{z}-1\right).$
(5)

The symbols $\sigma^{\pm}$ denote spin raising and lowering operators, while
$\sigma^{z}$ denotes the $z$-component of the spin. The spin-up state is
interpreted as a particle present, while the spin-down state is interpreted as
a particle absent.

Without loss of generality we can assume $q\leq p$. For $p,q\neq 0$ the matrix
$H$ can be mapped to an XXZ spin $1/2$ chain with non-Hermitian, twisted
boundary conditions [9]. For $p=q$ the matrix $H$ is Hermitian and for $U=1$
reduces to the Heisenberg spin chain.

The matrix $H$ can be written in terms of fermions by a Jordan-Wigner
transformation

$c_{j}^{(\dagger)}=e^{i\pi\sum_{k<j}\sigma_{k}^{+}\sigma_{k}^{-}}\sigma_{j}^{-(+)},$
(6)

where $c^{(\dagger)}_{j}$ are fermionic annihiliation (creation) operators.
The corresponding fermionic operator $H$ is then given by

$\displaystyle H_{0}$
$\displaystyle=\sum_{j=1}^{L-1}\left(pc_{j+1}^{\dagger}c_{j}+qc_{j}^{\dagger}c_{j+1}\right)+(-1)^{N+1}(pc_{1}^{\dagger}c_{L}+qc_{L}^{\dagger}c_{1})$
(7) $\displaystyle\mathcal{I}$
$\displaystyle=\sum_{j=1}^{L}c_{j}^{\dagger}c_{j}c_{j+1}^{\dagger}c_{j+1}-N.$
(8)

$H_{0}$ is the Hamiltoninian of non-hermitian free spinless fermions, while
$\mathcal{I}$ denotes a fermionic quartic interaction. The off-diagonal
elements of $H$ given by $H_{0}$ are non-negative, while the diagonal of the
diagonal matrix $\mathcal{I}$ consists of non-positive values.

The number of particles $N$ (spin-up states) is conserved by $H$ for all
interaction strengts $U$.

###  II.2 Open Boundary Conditions

For obc, the matrix $H_{0}$ is given by

$H_{0}=\sum_{j=1}^{L-1}\left(p\sigma_{j+1}^{+}\sigma_{j}^{-}+q\sigma_{j}^{+}\sigma_{j+1}^{-}\right)+\alpha\sigma_{1}^{+}+\gamma\sigma_{1}^{-}+\beta\sigma_{L}^{-}+\delta\sigma_{L}^{+},$
(9)

while the diagonal $\mathcal{I}$ is given by

$\displaystyle\mathcal{I}$
$\displaystyle=\frac{1}{4}\sum_{j=1}^{L-1}\left(\sigma_{j}^{z}\sigma_{j+1}^{z}-1\right)$
$\displaystyle+\frac{1}{2}\left[(p-q-\alpha+\gamma)\sigma^{z}_{1}+(q-p-\delta+\beta)\sigma^{z}_{L}\right]$
$\displaystyle-\frac{1}{2}\left[\alpha+\beta+\gamma+\delta\right].$ (10)

The bulk term of $H_{0}$ for obc is the same as for pbc. The terms at the
edges of the chain on site 1 and $L$ with parameters
$\alpha,\beta,\gamma,\delta$ denote particles hopping in and out of the chain
from an infinite reservoir of particles. Similar to pbc and $p,q\neq 0$, $H$
can be mapped to an XXZ chain with non-Hermitian, twisted boundary conditions
[35].

As in the pbc case, the operator $H$ can be written in terms of fermions. The
single spin operators at the end of the chain on site 1 and $L$ hinder a
straightforward application of a Jordan-Wigner transformation. Instead, we
treat the infinite reservoir as an additional site. We enlarge the chain of
length $L$ to a ring of length $L+1$ and change the terms connecting to site
$L+1$ accordingly. This is formally done by application of the well-known
Kramers-Wannier duality transformation [36]
$\sigma_{j}^{x}\to\prod_{l=1}^{j}\sigma_{l}^{z}$ and
$\sigma_{j}^{z}\to\sigma_{j}^{x}\sigma_{j+1}^{x}$. The details are in Appendix
A. Adding a site to the chain comes with the caveat that the multiplicity of
every eigenvalue of the so-transformed $H_{0}$ is doubled.

To keep the algebra simpler we restrict to the TASEP case $p=1$ and
$q=\gamma=\delta=0$, leaving $\alpha$ and $\beta$ as free parameters. The
following results can be straightforwardly generalized to arbitrary
$p,q,\gamma,\delta$. As outlined in Appendix A the Hamiltonian $H_{0}$ is
expressible in terms of spinless fermions $c,c^{\dagger}$ as

$\displaystyle H_{0}=$
$\displaystyle\alpha(c_{L+1}-c_{L+1}^{\dagger})c_{1}^{\dagger}+\sum_{j=1}^{L-1}\left[c_{j}c_{j+1}^{\dagger}\right]$
$\displaystyle+(-1)^{L}\mathcal{P}_{c}\beta c_{L}(c_{L+1}+c_{L+1}^{\dagger}),$
(11)

where $\mathcal{P}_{c}$ denotes the parity of the fermion number

$\displaystyle\mathcal{P}_{c}=(-1)^{\sum_{j=1}^{L+1}c^{\dagger}_{j}c_{j}}=(-1)^{N},$
(12)

which is conserved by $H_{0}$. Restricted to a fixed parity sector, $H_{0}$ is
a quadratic Hamiltonian. The corresponding spectrum is the same for each
parity sector leading to the aforementioned doubling of the spectral
multiplicity. This will be shown in detail in Sec. IV.2.

In summary, the non-interacting TASEP $H_{0}$ on $L$ sites with obc can be
written as a free fermion model on $L+1$ sites, with twisted pbc and
‘superconducting’ terms $c_{L,1}^{(\dagger)}c_{L+1}^{(\dagger)}$ connecting to
the additional site $L+1$.

###  II.3 Spectrum

All eigenvalues of $H$ are either real or come in complex conjugate pairs.
This characteristic stems from the fact that $H$ can be represented as a real
matrix. Specifically, for the case where $U=1$, the stochastic nature of $H$
dictates that its spectrum is situated in the left half of the complex plane.

Fig. 1 presents the spectrum of TASEP on a lattice with $L=11$ sites. The
spectral boundary shows $L$ spikes for pbc ($N=5$ particles) for $U=1$ in (a)
and $U=0$ in (c) and $L+1$ spikes for obc and $U=1$ in (b) and $U=0$ in (d).
For obc the parameters corresponding to the reservoirs are chosen as
$\alpha=\beta=1$ and $\gamma=\delta=0$. The subsequent sections primarily aim
to derive the mechanism responsible for the spikes in the spectral boundary.

Panels (c) and (d) of Fig. 1 reveal a highly structured spectrum for the non-
interacting TASEP $H_{0}$, exhibiting rotational invariance at angles $2\pi/L$
for pbc and $2\pi/(L+1)$ for obc. This characteristic stems from a “quasi-
symmetry” of $H_{0}$, which is investigated in detail in Sections III and IV.

For TASEP with obc, the spectral boundary spikes are always prominent, as
illustrated for the non-interacting TASEP in Sec. IV. However, this is not the
case for pbc. In Fig. 3(a) the spectrum of the pbc TASEP ($U=1$) and in (b)
its non-interacting variant ($U=0$) are presented for $L=40$ sites and $N=2$
particles, without any noticeable spikes in the spectral boundary. Sec. III
will demonstrate that, technically, the spectral boundary of the non-
interacting TASEP has $L=40$ spikes, but their distinctiveness fades in the
dilute limit where $\rho\to 0$.

##  III “Non-interacting” ASEP with pbc

In this section, we investigate the spectrum of the non-interacting ASEP
$H_{0}$ for pbc given by Eq. (4) and Eq. (7), respectively. Sec. III.1 is
devoted to the calculation of the single-body eigenvalues of $H_{0}$. In Sec.
III.2 we show the rotational invariance of the many-body spectrum of TASEP and
in Sec. III.3 we combine the results from the preceding subsections and show
how the spiky spectral boundary emerges. We quantify the prominence of the
spikes in Sec. III.4 and comment on whether they survive in the limit of large
$L$.

###  III.1 Single-body spectrum

Figure 2: Spectrum of the non-interacting TASEP $H_{0}$ on $L=11$ sites with
pbc. Single-body eigenvalues with $p=1$ and $q=0$ in (a) and $p=0.7$ and
$q=0.3$ in (b). In (c) we show part of the many-body spectrum with $N=5$
particles highlighting the tips of the spikes (red) and other boundary
eigenvalues (blue). All boundary eigenvalues are located on circles of radius
$1$, with crosses marking the midpoints.

Let us focus on the totally asymmetric case $p=1$ and $q=0$ first. Considering
the single-body sector of $H_{0}$ as given in Eq. (7), we see that the single-
body spectrum $\lambda$ is given by roots of the polynomial

$\lambda^{L}+(-1)^{N+1}.$ (13)

The roots are given by $\lambda=\omega^{j}$, where $\omega=e^{i\pi/L}$ and
$0\leq j<2L$ runs over all even (odd) integers when $N$ is odd (even). Thus
the single-body spectrum lies on the unit circle. In Fig. 2(a) the single-body
spectrum for $p=1$ and $q=0$ and $L=11$ and odd $N$ is shown together with the
unit circle.

For arbirtrary values of $p$ and $q$, the single-body spectrum is represented
as

$\lambda=p\omega^{j}+q\omega^{-j},$ (14)

with $j$ defined as previously. This spectrum lies on an ellipse with foci at
$\pm 2\sqrt{pq}$ and semi-major axis $p+q$ and semi-minor axis $p-q$,

$\\{(p+q)\cos(t)+i(p-q)\sin(t):0\leq t\leq 2\pi\\}.$ (15)

Figure 2 (b) illustrates the single-body spectrum for $p=0.7$ and $q=0.3$,
alongside the ellipse defined by Eq. (15).

The structure of the single-body spectrum for any $p,q$ suggests a
straightforward relation with the totally asymmetric scenario $q=0$. By
modifying the imaginary component while maintaining the real part constant,

$z\to\operatorname{Re}z+i\frac{p+q}{p-q}\operatorname{Im}z,$ (16)

we can convert the single-body eigenvalues for general $p,q$ values to those
corresponding to the $q=0$ case. This transforms the ellipse into a circle of
radius $p+q$. Without loss of generality, we restrict ourselves to $p=1$ and
$q=0$ for the remainder of this section.

###  III.2 Rotational invariance

With $p=1$ and $q=0$, the single-body spectrum remains unchanged under complex
plane rotations of $2\pi/L$. This rotational invariance also applies to the
many-body spectrum, which comprises sums of single-body eigenvalues.

Furthermore, this symmetry is evident in $H_{0}$ when transforming spin and
fermionic operators. Transforming $c_{j}\to e^{-i2\pi j/L}c_{j}=\tilde{c}_{j}$
and $c_{j}^{\dagger}\to e^{i2\pi j/L}c_{j}^{\dagger}=\tilde{c}_{j}^{\dagger}$,
or in terms of spin operators $\sigma_{j}^{\pm}\to e^{\pm i2\pi
j/L}\sigma_{j}^{\pm}=\tilde{\sigma}_{j}^{\pm}$, results in
$e^{i2\pi/L}H_{0}=\tilde{H}_{0}$. Here, $\tilde{H}_{0}$ is constructed like
$H_{0}$, but using the modified operators $\tilde{c},\tilde{c}^{\dagger}$
($\tilde{\sigma}^{\pm}$). Since these altered operators maintain their
respective (anti-)commutation relations, the spectra of $H_{0}$ and
$\tilde{H}_{0}$ are identical. Therefore, the spectrum of $H_{0}$ is invariant
under $2\pi/L$ rotations.

###  III.3 Spectral boundary

The structure of the many-body spectrum as observed in Fig. 1 is now a
consequence of the relation of single-body to many-body eigenvalues and the
rotational symmetry.

For ease of notation, we define $\lambda_{j}=\omega^{2j}$ when $N$ is odd, and
$\lambda_{j}=\omega^{2j+1}$ for even $N$. The many-body eigenvalues are
obtained by adding $N$ of these $L$ single-body eigenvalues. More precisely,
the many-body eigenvalues $E$ correspond uniquely to configurations
$s=(s_{1},\dots,s_{L})\in\\{0,1\\}^{L}$, where $\sum_{j}s_{j}=N$, and are
given by

$E=\sum_{j=1}^{L}s_{j}\lambda_{j}.$ (17)

The many-body eigenvalues $E_{t}$ which appear at the spike tips, have the
highest absolute values and are derived from configurations $s$ with
contiguous non-zero $s_{j}$ entries. Specifically, each of the $L$ tips
$E_{t}(j_{0})$ is linked to an index $1\leq j_{0}\leq L$ and a configuration
$s=s_{t}(j_{0})$ with

$s_{j}=\begin{cases}1&j_{0}\leq j\leq j_{0}+N-1\\\
0&\text{otherwise.}\end{cases}$ (18)

Here, $j\equiv j-L$ is applied for $j>L$. The eigenvalues $E_{t}(j_{0})$ are
calculated as

$E_{t}(j_{0})=\sum_{j=j_{0}}^{j_{0}+N-1}\lambda_{j}.$ (19)

Configurations $s$ that lead to spike tips are termed ’domain wall
configurations’. The many-body eigenvalues $E_{t}$ are depicted as red circles
(light colored in print) in Fig. 2(c).

Boundary eigenvalues in the many-body spectrum arise from ’interpolating’
between configurations of adjacent spike tips. In these configurations, the
domain walls differ by a shift of one site. The interpolation process between
these two domain walls involves moving a single particle (or executing a
single spin flip). As a result, the configurations formed contain a maximum of
two domain walls, each separated by one site. Specifically, boundary
configurations $s=s_{b}(j_{0},l_{0})$ are associated with indices $1\leq
j_{0}\leq L$ and $j_{0}\leq l_{0}\leq j_{0}+N$, defined as

$s_{j}=\begin{cases}1&j_{0}\leq j\leq j_{0}+N\text{ and }j\neq l_{0}\\\
0&j=l_{0}\\\ 0&\text{otherwise.}\end{cases}$ (20)

Again, $j\equiv j-L$ is used for $j>L$. The corresponding boundary eigenvalues
$E_{b}(j_{0},l_{0})$ are computed by

$E_{b}(j_{0},l)=\sum_{j=j_{0};j\neq l}^{j_{0}+N}\lambda_{j}.$ (21)

When $l_{0}=j_{0}$ or $l_{0}=j_{0}+N$ (indicating a single domain wall), the
boundary eigenvalue matches a spike tip, $E_{b}(j_{0},j_{0})=E_{t}(j_{0}+1)$
or $E_{b}(j_{0},j_{0}+N)=E_{t}(j_{0})$, respectively. The boundary eigenvalues
$E_{b}(j_{0},l)$ for $j_{0}<l<j_{0}+N$ are those many-body eigenvalues located
’between’ the spike tips $E_{t}(j_{0})$ and $E_{t}(j_{0}+1)$, depicted as blue
circles in Fig. 2(c).

Eq.(21) can be reformulated as

$E_{b}(j_{0},l)=\sum_{j=j_{0}}^{j_{0}+N}\lambda_{j}-\lambda_{l}.$ (22)

Given $|\lambda_{l}|=1$ and the independence of the sum from $l$, all boundary
eigenvalues are on $L$ circles of radius 1. For $N\leq L/2$, the circle
midpoints are the many-body spectrum tips $E_{t}^{(N+1)}(j_{0})$ with $N+1$
particles. The tips $E_{t}^{(N)}$ intersect two adjacent circles. This is
illustrated in Fig. 2(c) with circles as black lines and midpoints as gray
crosses.

According to Eq. (19), all tips reside on a circle with radius $R$, defined as

$R=\left|\frac{1-e^{i2\pi N/L}}{1-e^{i2\pi/L}}\right|=\frac{\sin(\pi
N/L)}{\sin(\pi/L)}.$ (23)

This radius, combined with the circular pattern of the boundary eigenvalues,
enables us to establish a continuous boundary for the many-body spectrum. It
is formed by the intersection of all circles of radius 1 with the disc of
radius $R$ from Eq. (23). The boundary is parameterized by

$z_{B}(t)=e^{-if(t)}\left(\gamma_{1}+\gamma_{2}e^{ig(kt)}\right),$ (24)

with $\gamma_{1}=\frac{\sin(\pi\rho)}{\sin(\pi/L)}$ and $\gamma_{2}=1$, with
piece-wise constant $f$,

$f(t)=\frac{\pi}{L}\left(2\left\lfloor\frac{Lt}{2\pi}\right\rfloor-1\right),$
(25)

and $g$ is piece-wise the identity,

$g(t)=\pi(1-\rho)+\rho(t\operatorname{mod}\ 2\pi).$ (26)

The continuous boundary $z_{B}(t)$ is illustrated as a red (gray in print)
curve in Fig. 1(c) for $L=11$ and $N=5$ and in Fig. 3(b) for $L=40$ and $N=2$.
As expected, all boundary eigenvalues reside on the continuous boundary
parametrized by $z_{B}(t)$.

Eq. (24) is related to the spectral boundary of random matrices with higher-
order cyclic correlations between $L$-tuples of matrix elements, akin to
random graphs with a dominant cycle structure [37]. Their spectral boundary
forms a hypotrochoidic curve, which is recovered from Eq. (24) by letting
$f(t)=g(t)=t$. This relation hints at the connection between the spectral
boundary of the non-interacting TASEP and random matrices; we explore this
connection in Section VI.

###  III.4 Quantification of spikes

This subsection aims to measure the sharpness of the spectral boundary in the
non-interacting TASEP, particularly focusing on whether spikes persist in
large system sizes and, if so, how. For simplicity, we consider particle
densities $0\leq\rho\leq 1/2$. As the ASEP spectrum is invariant under
changing $\rho\to 1-\rho$ this comes with no loss of generality.

To assess the spikiness of the spectral boundary, we examine the ratio between
two distances: $d_{t}$, the distance between spike tips, and $d_{b}$, the
maximum extension of the spectral boundary beyond a circle of radius $R$.
Recall, this circle of radius $R$ represents the smallest enclosing disk for
the TASEP spectrum. $d_{b}$ measures how far the radius 1 circles, carrying
the boundary eigenvalues, reach into the enclosing circle. A larger $d_{b}$
relative to $d_{t}$ indicates that these radius 1 circles extend more into the
enclosing spectrum. Therefore, the ratio $2d_{b}/d_{t}$ quantifies the
spikiness of the boundary. A value close to 1 suggests a spiky boundary, while
a significantly smaller ratio implies a less spiky boundary. This factor of
two arises because $d_{t}$ pertains to the diameter of the boundary circles,
whereas $d_{b}$ is compared to their radius.

Following some simple trigonometry one finds that the distances $d_{t}$ and
$d_{b}$ are given by

$d_{t}=2\sin(\pi\rho)$ (27)

and

$d_{b}=1-\frac{\cos(\pi\rho+\pi/(2L))}{\cos(\pi/(2L))}.$ (28)

The fraction $2d_{b}/d_{t}$ then simplifies to

$\frac{2d_{b}}{d_{t}}=\tan(\pi\rho/2)+\tan(\pi/(2L)).$ (29)

Eq. (29) shows a monotonic increase with $\rho$, indicating that the spectral
boundary becomes more pronouncedly spiky at higher $\rho$ values. Due to the
invariance of the spectrum under the transformation $\rho\to 1-\rho$, the
boundary reaches its maximum spikiness at $\rho=1/2$.

The analytical findings are confirmed by panels (c) and (b) in Fig. 1 and Fig.
3, respectively. In Fig. 1(c), the many-body spectrum of $H_{0}$ is markedly
spiky for $\rho=5/11\approx 0.45$, whereas in Fig. 3(b), the spectral boundary
is nearly circular, aligning with the low $\rho$ value of 2/40=0.05.

In examining the large $L$ limit, we will explore two scenarios: the
“thermodynamic” limit, where both $N$ and $L$ increase to infinity while
maintaining a fixed $\rho$, and the few-particle (dilute) limit, where $N$
remains constant and only $L$ approaches infinity.

####  III.4.1 “Thermodynamic” limit

In the thermodynamic limit, the distance $d_{t}$ remains constant, whereas
$d_{b}$ approaches $1-\cos(\pi/\rho)$. Consequently, the ratio $2d_{b}/d_{t}$
tends towards $\tan(\pi\rho/2)$. This implies that for any non-zero $\rho$,
the spiky structure of the spectral boundary is preserved in the thermodynamic
limit, becoming more pronounced with increasing $\rho$.

Fig. 1(c) presents the many-body spectrum of the non-interacting TASEP for
$L=11$ and $N=5$, with Fig. 2(c) offering a closer view of the spectral
boundary. Here, $\rho\approx 0.45$ and $2d_{b}/d_{t}\approx 1.01$ indicate
pronounced spikes of the spectral boundary, as evident.

Regarding the length scales at which these spikes are observable, consider the
following: The radius $R$ of the spectrum scales as $O(L)$, necessitating a
rescaling of the spectrum by $1/L$ to ensure a well-defined spectral density
in the thermodynamic limit. At an infinite $L$, this rescaled spectrum densely
fills the unit circle. For finite $L$, the tips of the spikes are spaced at a
distance of $d_{t}=O(1/L)$, and the distance $d_{b}$ of the spectral boundary
from the unit circle is also $O(1/L)$. Therefore, at the length scale of
$1/L$, the spiky nature of the spectral boundary is distinctly visible.

####  III.4.2 Dilute limit (large $L$, constant $N$)

Figure 3: Spectrum of the generator matrix $H$ of TASEP (a) and the “non-
interacting” TASEP (b) on $L=40$ sites with $N=2$ particles (dilute limit).
The red solid line in (b) denotes the spectral boundary according to Eq. (24).
The spectral boundary appears smooth and non-spiky in both panels.

In the scenario where $N$ is fixed and $L$ increases, both distances $d_{t}$
and $d_{b}$ decrease, scaling as $O(1/L)$ and $O(1/L^{2})$, respectively.
Consequently, the ratio $2d_{b}/d_{t}$ tends towards 0, as indicated by Eq.
(29). Therefore, in this limit, the spiky structure of the spectral boundary
does not persist.

Fig. 3 shows the many-body spectrum of the TASEP for $L=40$ and $N=2$,
representative of the dilute limit. We show both a TASEP case ($U=1$) and a
non-interacting TASEP case ($U=0$). With a $2d_{b}/d_{t}$ ratio of $\approx
0.01$ it reveals a non-spiky spectral boundary, barely distinguishable from a
circle, as shown by the red curve (gray in print) in Fig. 3(b).

##  IV “Non-interacting” TASEP with obc

In this section, we will present the analytical derivation of the many-body
spectrum of the non-interacting TASEP $H_{0}$ with obc, specifically for $p=1$
and $q=\gamma=\delta=0$. Generalizations to arbitrary $p,q,\gamma,\delta$ are
straightforward.

In Sec. IV.1 we establish the rotational invariance of the spectrum of
$H_{0}$. In Sec. IV.2 we derive its single-particle spectrum and demonstrate
its relation to the many-body eigenvalues. Sec. IV.3 demonstrates that the
spectral boundary of $H_{0}$, similar to the pbc case, is defined by the
intersection of circles with a disk, featuring $L+1$ spikes. In the limit of
large $L$, this boundary is akin to the pbc case with density $\rho=1/2$,
highlighted in Sec. IV.4.

###  IV.1 Rotational symmetry

The spectrum of the non-interacting TASEP $H_{0}$ is invariant under rotations
of angle $\frac{2\pi}{L+1}$. Similar to the pbc case, consider the change of
operators $c_{j}^{\dagger}\to
e^{i\frac{2\pi}{L+1}j}c_{j}^{\dagger}=\tilde{c}_{j}^{\dagger}$ and $c_{j}\to
e^{-i\frac{2\pi}{L+1}j}c_{j}=\tilde{c}_{j}$ or, equivalently,
$\sigma_{j}^{\pm}\to e^{\pm
i\frac{2\pi}{L+1}j}\sigma_{j}^{\pm}=\tilde{\sigma}^{\pm}$. This change implies
that $e^{i\frac{2\pi}{L+1}}H_{0}=\tilde{H}_{0}$, where $\tilde{H}_{0}$ is
$H_{0}$ with $c,c^{\dagger}$ ($\sigma$) replaced by the tilde operators. As
the tilde operators fulfill the canonical (anti-)commutation relations of
fermion operators (Pauli matrices), the spectrum of the non-interacting TASEP
is invariant under rotations of angle $\frac{2\pi}{L+1}$.

###  IV.2 Single- and many-body spectrum

Before we diagonalize $H_{0}$, let us specify the parity sector as
$s=(-1)^{L}\mathcal{P}_{c}$. To simplify the following arguments, we will
abuse notation and not distinguish between $H_{0}$ and $H_{0}$ restricted to a
subspace of constant parity. At the end of this subsection, we will take the
difference into account properly.

Let us collect the Dirac fermion operators $c,c^{\dagger}$ into a
$(2L+2)$-dimensional vector
$\boldsymbol{c}=(c_{1},\dots,c_{L+1},c_{1}^{\dagger},\dots,c_{L+1}^{\dagger})^{t}$.
We express $H_{0}$ given by Eq. (11) as

$H_{0}=\frac{1}{2}\boldsymbol{c}^{\dagger}\begin{pmatrix}A&B\\\
C&-A^{t}\end{pmatrix}\boldsymbol{c}=\frac{1}{2}\boldsymbol{c}^{\dagger}M_{c}\boldsymbol{c}$
(30)

where the $(L+1)\times(L+1)$-matrices $A,B$ and $C$ are given by

$\displaystyle A_{ij}$ $\displaystyle=-\delta_{i,j+1\mathrm{mod}(L+1)}$
$\displaystyle+(1-\beta
s)\delta_{i,L+1}\delta_{j,L}+(1-\alpha)\delta_{i,1}\delta_{j,L+1},$ (31)
$\displaystyle B_{ij}$
$\displaystyle=\alpha(\delta_{i,1}\delta_{j,L+1}-\delta_{i,L+1}\delta_{j,1}),$
(32) $\displaystyle C_{ij}$ $\displaystyle=\beta
s(\delta_{i,L}\delta_{j,L+1}-\delta_{i,L+1}\delta_{j,L}),$ (33)

and $\delta$ denotes the Kronecker-delta symbol.

The matrix $A$ is, up to deformations in the $(1,L+1)$th and $(L+1,L)$th
entries, a circulant matrix with only one non-zero off-diagonal. The matrices
$B$ and $C$ only contain two non-zero entries. Thus, the solutions $\lambda$
and $u$ to the eigenvalue problem

$M_{c}u=\lambda u,$ (34)

are closely related to the eigen-decomposition of circulant matrices, which in
turn are given by Fourier transforms. As shown in detail in Appendix B, the
eigenvalues $\lambda$ are solutions of

$\lambda^{2L+2}=4(\alpha\beta)^{2}(-1)^{L},$ (35)

and are independent of the parity sector $s$. Since the polynomial in Eq. (35)
is of even degree, its roots appear in pairs of $\pm\lambda$.

The Hamiltonian $H_{0}$ in Eq. (11) is non-Hermitian, preventing the direct
use of the (Hermitian) Bogoliubov-de-Gennes formalism for linking the
eigenvalues of $M_{c}$ to the many-body spectrum of $H_{0}$. Hence, we will
pursue an alternative method. We proceed as in [38] and express
$c,c^{\dagger}$ in terms of Majorana fermions

$\phi_{j,1}=\frac{1}{\sqrt{2}}(c_{j}+c_{j}^{\dagger}),\quad\phi_{j,2}=\frac{1}{i\sqrt{2}}(c_{j}-c_{j}^{\dagger}).$
(36)

After collecting the Majorana fermions $\phi_{j,l}$ into a column vector
$\boldsymbol{\phi}=(\phi_{1,1},\phi_{1,2},\dots,\phi_{L+1,1},\phi_{L+1,2})^{t}$,
$H_{0}$ can be written as

$\displaystyle
H_{0}=\frac{1}{2}\boldsymbol{\phi}^{t}M_{\phi}\boldsymbol{\phi},$ (37)

where the matrix $M_{\phi}$ is a complex and anti-symmetric
$(2L+2)\times(2L+2)$-matrix. The transformation of Majorana fermions $\phi$ to
Dirac fermions $c$ via Eq. (36) is unitary, making $M_{\phi}$ and $M_{c}$
unitarily equivalent and hence sharing the same eigenvalues.

As $M_{\phi}$ is anti-symmetric, it can be factorized [38] as

$M_{\phi}=\frac{1}{2}V\Lambda JV^{t}$ (38)

where

$V^{t}V=J=\operatorname{Id}_{L+1}\otimes\begin{pmatrix}0&1\\\
1&0\end{pmatrix},$ (39)

$\operatorname{Id}_{L+1}$ denotes the $(L+1)\times(L+1)$ identity matrix and
$\Lambda$ is a diagonal matrix containing the eigenvalues of $M_{\phi}$
($M_{c}$). The anti-symmetry of $M_{\phi}$ implies that its eigenvalues come
in pairs $\pm\lambda$, which is consistent with the solutions of Eq. (35). The
diagonal of $\Lambda$ is ordered as
$\lambda_{1},-\lambda_{1},\dots\lambda_{L+1},-\lambda_{L+1}$. We fix the
choice between $\lambda_{j}$ and $-\lambda_{j}$ by requiring
$\operatorname{Re}\lambda_{j}\geq 0$.

Let us define another type of Dirac fermions $b,b^{\prime}$ as

$(b_{1},b_{1}^{\prime},\dots,b_{L+1},b_{L+1}^{\prime})^{t}=\left(V^{t}\boldsymbol{\phi}\right).$
(40)

These fulfill the usual anti-commutation relations of Dirac fermions [38], but
$b^{\prime}$ is in general not the Hermitian adjoint of $b$. Nevertheless, the
Hamiltonian $H_{0}$ becomes diagonal in terms of $b,b^{\prime}$,

$H_{0}=\sum_{j=1}^{L+1}\lambda_{j}b_{j}^{\prime}b_{j}-\frac{1}{2}\sum_{j=1}^{L+1}\lambda_{j}.$
(41)

The eigenstates of $H_{0}$ are given by creation operators $b_{j}^{\prime}$
acting on the vacuum $\left|0\right\rangle_{b}$, which are $2^{L+1}$ in total.
But not all eigenstates correspond to an eigenvalue of $H_{0}$ given by Eq.
(11). We have to take into account that the Dirac fermions $b,b^{\prime}$ are
only defined on fixed parity subspaces.

We numerically find that the parity operator $\mathcal{P}_{b}$ of the
$b,b^{\prime}$ fermions obeys

$\mathcal{P}_{b}=-s\mathcal{P}_{c},$ (42)

where $P_{c}$ denotes the parity operator of the $c$ fermions. Recall that we
let $s=(-1)^{L}\mathcal{P}_{c}$ at the beginning of this subsection. Thus, the
admissible $b^{\prime}$-fermion states must have $b$-parity
$\mathcal{P}_{b}=-(-1)^{L}=(-1)^{L+1}$. Especially, the parity of the
admissible $b$-states does not depend on $s$. Thus both parity sectors give
rise to the same many-body spectrum of $H_{0}$ in Eq. (41), as required.

In summary, the many-body spectrum of the non-interacting TASEP, subject to a
global shift in the complex plane, is represented by the sums of the $L+1$
roots from Eq. (35) with positive real parts. These are scaled roots of $\pm
1$ with magnitude proportional to $(\alpha\beta)^{1/(L+1)}$. Depending on
whether $L$ is odd or even, an even or odd number of summands, respectively,
are included in the sums.

###  IV.3 Spectral boundary

Figure 4: Many-body spectrum of the non-interacting TASEP with obc on (a)
$L=6$ and (b) $L=7$ sites. Similar to pbc in Fig. 2, all boundary eigenvalues
lie on circles, with midpoints denoted by crosses.

The emergence of the many-body spectrum of the non-interacting TASEP with obc
follows a similar principle than for pbc discussed in Sec. III: the many-body
spectrum consists of sums of (scaled) roots of $\pm 1$. In the following, we
describe how the spiky spectral boundary emerges for obc. Especially, we will
demonstrate that, akin to the pbc case, the spectral boundary resides on $L+1$
circles, each with a radius of $(2\alpha\beta)^{1/(L+1)}$, and provide a
comparable parametrization for this boundary.

In the following, we focus exclusively on the spectral boundary associated
with the most negative real parts. This is illustrated in Fig. 4, parts (a)
and (b), where the eigenvalues of the relevant sectors are marked with blue
and red circles. The rotational symmetry of the spectrum means that the
structure of the boundary is a repetitive pattern reflecting the shape of
sectors with the smallest real parts. Hence, restriction to sectors with the
most negative real part eigenvalues comes with no loss of generality.

Let us first consider even $L$. Recall that the many-body spectrum is given by
sums of an odd number of positive real part roots of the polynomial in Eq.
(35). Let us denote the $L+1$ roots with non-negative real part by
$\lambda_{1},\dots,\lambda_{L+1}$. Then the $L+1$ many-body eigenvalues with
the smallest real parts are eigenvalues lying on the spectral boundary and
given by

$\lambda_{j}-\frac{1}{2}\sum_{l=1}^{L+1}\lambda_{l}.$ (43)

If we label $\lambda_{l}$ by increasing angle with branch-cut on the negative
imaginary axis then the tips of the spectrum are given by the indices $j=1$
and $j=L+1$.

In Fig. 4(a) we show the spectrum of the non-interacting TASEP with obc on
$L=6$ sites. The spectrum shows $L+1=7$ spikes. The boundary and tips
according to Eq. (43) are shown as blue and red markers, respectively. The
markers lie on a circle with midpoint
$-\frac{1}{2}\sum_{l=1}^{L+1}\lambda_{l}$ and radius
$|\lambda_{j}|=(2\alpha\beta)^{1/(L+1)}$.

Let us now consider the slightly more complicated case of odd $L$. In Fig.
4(b) we show the many-body spectrum on $L=7$ sites. The tip of the spectral
edge with the smallest real part is given by an ‘empty’ sum of $\lambda_{l}$’s
and thus is $-\frac{1}{2}\sum_{l=1}^{L+1}\lambda_{l}$. The boundary
eigenvalues are given by the following (shifted) sum of two single-particle
eigenvalues:

$\lambda_{j}+\lambda_{1,L}-\frac{1}{2}\sum_{l=1}^{L+1}\lambda_{l},$ (44)

where $2\leq j\leq L$ and $\lambda_{1}$ corresponds to the lower spectral
boundary in Fig. 4(b) while $\lambda_{L}$ corresponds to the upper part. The
midpoints of the circles are given by
$\lambda_{1,L}-\frac{1}{2}\sum_{l=1}^{L+1}\lambda_{l}$ and the radius again by
$|\lambda_{j}|=(2\alpha\beta)^{1/(L+1)}$.

Similar to the pbc case, we can establish a continuous boundary for the many-
body spectrum, parametrized by Eq. (24). In the obc case the constants
$\gamma_{1,2}$ are given by

$\displaystyle\gamma_{1}$
$\displaystyle=(2\alpha\beta)^{1/(L+1)}\frac{1}{2\sin(\pi/(2L+2))}$ (45)
$\displaystyle\gamma_{2}$ $\displaystyle=(2\alpha\beta)^{1/(L+1)},$ (46)

while the piece-wise constant $f$ and the piece-wise identify function $g$ are
given by

$\displaystyle f(t)$
$\displaystyle=\frac{\pi}{L}\left(2\left\lfloor\frac{Lt}{2\pi}\right\rfloor-1\right)$
(47) $\displaystyle g(t)$
$\displaystyle=\pi\frac{L+2}{2L+2}+\frac{L}{2L+2}(t\operatorname{mod}\ 2\pi).$
(48)

The continuous boundary $z_{B}(t)$ with the above parameters is illustrated in
Fig. 1(d) as a red (gray in print) curve for $L=11$.

###  IV.4 Spikes in the large $L$ limit

The parametrization of the spectral boundary for obc shows a clear link to the
spectral boundary for pbc. Specifically, in the large $L$ limit with constant
$\alpha,\beta$, the obc spectral boundary aligns with the pbc case at
$\rho=1/2$. This relation is immediately evident for $\gamma_{2}$, $f$, and
$g$. Further, a series expansion of $\gamma_{1}$ for large $L$ reveals that
its leading term, $\gamma_{1}=L/\pi+O(1)$, is identical in both cases, with
differences emerging only at $O(1)$.

Consequently, in the large $L$ limit, the spiky spectral boundary in the obc
case remains pronounced. Rescaling the spectrum by $1/L$, the spectral density
approaches filling the unit disk as $L\to\infty$. For finite $L$, the tips are
spaced by $O(1/L)$, and the maximum deviation of the boundary from the unit
circle is also $O(1/L)$.

##  V Pbc TASEP by Bethe ansatz

In Sec. III, we showed that in the non-interacting TASEP ($U=0$) with pbc, the
spiky boundary of the many-body spectrum emerges essentially as sums of evenly
spaced single-body eigenvalues $\lambda_{1},\dots,\lambda_{L}$. This section
expands that concept to interaction strengths $0<U$. Employing the coordinate
Bethe ansatz, we generalize the single-body framework to Bethe roots, which
tend to cluster close to $\lambda_{1},\dots,\lambda_{L}$. This clustering,
combined with TASEP many-body eigenvalues being sums of Bethe roots, results
in a spiky spectrum boundary for any interaction strength $0\leq U\leq 1$.

This section focuses on $\rho\approx 1/2$, where the most prominent spectral
boundary spikes in the non-interacting ASEP were observed. In the low-density
limit ($\rho$ approaching zero), we anticipate a spectral boundary for the
usual ASEP similar to the non-interacting case, characterized by a smooth,
circular boundary without spikes. Fig. 3 partly supports this, showing similar
many-body spectra for TASEP with $U=1$ (a) and $U=0$ (b), both featuring
smooth, non-spiky spectral boundaries.

Sec. V.1 generalizes the coordinate Bethe ansatz to arbitrary $U$, with
derivation details and numerical solution methods detailed in Appendices C and
D. In Sec. V.2, we demonstrate the clustering of solutions to the Bethe
equations and in Sec. V.3, we establish how this clustering results in a spiky
spectral boundary.

###  V.1 Coordinate Bethe ansatz

We start by determining the many-body spectrum of $H$, as described in Eqs.
(4) and (5), for arbitrary $U$. We closely follow the application of the
coordinate Bethe ansatz to the ASEP in [32], which dealt with $U=1$. The
coordinate Bethe ansatz has since been used extensively for ASEP [39, 32, 40,
6, 41, 42, 43, 44, 45, 9, 46, 47, 48, 49, 10, 50, 51, 33, 52, 53, 34, 11, 54,
55, 56, 57, 58, 59, 60] and for various variants/extensions of ASEP [61, 62,
63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 9, 78, 79, 80, 81,
82, 83, 84, 85, 86, 87, 58, 88, 89, 90, 91, 92]. However, to the best of our
knowledge an extension to general $U$ has not been presented before; we
present this extension in this work. The details of the derivation of the
Bethe equations are provided in Appendix C; in this section we present the
results and their application to understanding the spectral boundary.

The Bethe ansatz eigenvalues $E$ for arbitrary $U$ are given by

$E=\sum_{j=1}^{N}\big{(}pz_{j}+qz_{j}^{-1}-U\big{)},$ (49)

where $z_{j}$ are complex numbers, the so-called Bethe roots, which in turn
are solutions of the following recurrent relations

$z_{j}^{L}=\prod_{k=1;k\neq
j}^{N}\left(-\frac{q+pz_{j}z_{k}-Uz_{j}}{q+pz_{j}z_{k}-Uz_{k}}\right).$ (50)

The solutions of Eq. (50) are $N$-tuples $(z_{1},\dots,z_{N})$ and each
$N$-tuple gives rise to an eigenvalue $E$ of the TASEP $H$ via Eq. (49).

Numerical data indicates that in small systems, each eigenvalue is a sum of
Bethe roots, although a formal proof of the completeness of the Bethe ansatz
is lacking [93, 94, 95]. In our finite ASEP system investigations, all
eigenvalues conformed to the Bethe ansatz.

For $p=1$ and $q=0$, the eigenvalue equation simplifies to

$E=\frac{1}{2}\sum_{j=1}^{N}\big{(}Z_{j}-U\big{)},$ (51)

and the Bethe equations transform into

$\left(U+Z_{j}\right)^{L-N}\left(U-Z_{j}\right)^{N}=-2^{L}\prod_{k=1}^{N}\frac{Z_{k}-U}{Z_{k}+U},$
(52)

with

$Z_{k}=2z_{k}-U$ (53)

representing scaled, shifted Bethe roots. We refer to either the $z_{k}$ or
the $Z_{k}$ as the Bethe roots, depending on the context. In Eqs. (52), the
main simplification from the general $p,q$ case is the independence of the
right-hand side from $j$, which makes the solutions $Z_{j}$ roots of the
polynomial $P(Z)=(U+Z)^{L-N}(U-Z)^{N}-Y$, with $Y$ given by the right-hand
side of Eq. (52). This not only simplifies numerical computation of the Bethe
roots $Z_{j}$, but also ensures their continuity in $U$ [96]. Consequently, we
will focus on the specific case of $p=1$ and $q=0$ for the rest of this
section.

Appendix D details the numerical solution process for the Bethe Eqs. (52) and
the systematic retrieval of all Bethe roots.

Eq. (49) and (51) establish that many-body eigenvalues are sums of Bethe
roots, up to a global shift. To demonstrate a spiky spectral boundary, we will
show numerically a sufficient clustering of Bethe roots, which is the focus of
the rest of this section.

###  V.2 Clustering of the Bethe roots

Figure 5: All $N\cdot\binom{L}{N}$ Bethe roots $Z_{j}$ of the TASEP. (a)-(d)
$L=8$, $N=4$, for different values of $U$. (e) $L=14$, $N=7$, $U=1$.

To examine the spectral boundary in terms of the Bethe roots, we will consider
in the complex plane the Bethe roots ($z_{j}$ or $Z_{j}$) corresponding to
each of the $\binom{L}{N}$ eigenstates. There are thus $N\times\binom{L}{N}$
Bethe roots in total, for any value of $U$. Such plots are shown in Figure 5.

For $U=0$, the Bethe roots $z_{j}$ satisfy the equation
$z_{j}^{L}=(Z_{j}/2)^{L}=(-1)^{N+1}$, and agree with the single-body
eigenvalues of $H_{0}$ as stated in Eq. (14). Therefore, the many-body
spectrum derived via the Bethe ansatz for $U=0$ aligns with that of the non-
interacting ASEP model discussed in Section III, as expected. An illustrative
example of the Bethe roots $Z_{j}=2z_{j}$ for $U=0$ is provided in Figure 5
(a) for $L=8$ and $N=4$. Here, each solution of the Bethe equations
contributes $N=4$ roots, which together describe one of the $\binom{8}{4}$
eigenstates. We plot all the $4\times\binom{8}{4}$ roots together in a single
plot. Since for $U=0$ every solution to the Bethe equations is a subset of the
8 single-body eigenvalues of $H_{0}$, the union of all solutions is highly
degenerate and only 8 unique markers show up in Fig. 5(a).

For $U>0$ the degeneracy of the $U=0$ case is lifted and the
$4\times\binom{8}{4}$ Bethe roots $Z_{j}$ become distinct, as observed in Fig.
5(b-d) for $U=0.33$, $0.66$, and $U=1$, respectively. The continuity of Bethe
roots $z_{j}$ in $U$ suggests that for small $U$, these roots should be
proximate to the $L$th roots of $(-1)^{N+1}$. Numerically, this is confirmed
as the Bethe roots $z_{j}$ tend to cluster around the $L$th roots of
$(-1)^{N+1}$ for small $U$. As depicted in Figure 5 (b) and (c) for $U=0.33$
and $U=0.66$ respectively, the $Z_{j}$’s distinctly form $L=8$ clusters around
the Bethe roots for $U=0$. This clustering is even discernible for $U=1$, as
shown in Figure 5 (d), where the $L=8$ clusters remain identifiable.

For larger $L$, the Bethe root clusters overlap at $U=1$, evident from Fig.
5(e) for $L=14$ and $N=7$. However, the statistical width of these clusters
diminishes with larger $L$. This is demonstrated in Fig. 7, where the average
cluster width decreases as $L^{-1/2}$ in the thermodynamic limit with
$\rho=N/L=1/2$ and $N,L\to\infty$.

We define the locations and widths of these clusters by fitting a Gaussian
mixture model of $L$ independent Gaussians $\mathcal{N}$ with complex means to
the Bethe roots. The Bethe roots distribution is approximated as
$\frac{1}{L}\sum_{j=1}^{L}f_{j}$, with $f_{j}$ representing Gaussian
densities. We label the Gaussians of the optimal fit as $\mathcal{N}_{j}$,
each characterized by its mean $\mu_{j}$ and standard deviation $\sigma_{j}$.

###  V.3 Structure of the many-body spectrum

Figure 6: (a) The many-body spectrum of TASEP with $L=14$ and $N=7$
(multiplied by 2 and shifted by $N$). (b) Probability density function of the
many-body spectrum of the random Bethe roots $Z$ for $L=14$ and $N=7$ capped
at $10^{-3}$. Red (gray in print) dots are the means of the complex Gaussians.

In the following, we will show that by considering only the centers and widths
of Bethe root clusters, and not their specific structure, we can approximate a
many-body spectrum that mirrors key characteristics of the TASEP many-body
spectrum, particularly its spiky boundary.

Recall that for $U=0$ each many-body eigenvalue $E$ is a sum of $N$ out of $L$
single-body eigenvalues. Specifically, $E$ is given by

$E=\sum_{j=1}^{L}s_{j}\lambda_{j}=\sum_{s_{j}\neq 0}\lambda_{j},$ (54)

where $s\in\\{0,1\\}^{L}$ is a configuration with $\sum_{j}s_{j}=N$ and
$\lambda_{j}$ are the single-particle eigenvalues determined in Sec. III.3. By
Eq. (51) every many-body eigenvalue of the TASEP ($U=1$) corresponds to a sum
of $N$ Bethe roots $(Z_{1},\dots,Z_{N})$ and by the continuation from $U=1$ to
$U=0$ each Bethe root $Z_{j}$ belongs to one of the $L$ clusters. Instead of
summing solutions of the Bethe Eqs. (52) we employ a statistical ansatz and
consider random many-body eigenvalues of the form

$E_{rand}=\frac{1}{L}\sum_{j=1}^{L}s_{j}\mathcal{N}_{j}=\mathcal{N}_{s},$ (55)

where $\mathcal{N}_{s}$ denotes a Gaussian with mean
$\sum_{j=1}^{L}s_{j}\mu_{j}$ and variance $\sum_{j=1}^{L}s_{j}\sigma_{j}^{2}$.
We refer to $\mathcal{N}_{s}$ as many-body Gaussians. The full random many-
body spectrum is then given by

$\frac{1}{\mathcal{Z}}\sum_{\begin{subarray}{c}s\in\\{0,1\\}^{L}\\\
s_{1}+\dots+s_{L}=N\end{subarray}}\mathcal{N}_{s}$ (56)

where $\mathcal{Z}=L\binom{L}{N}$ is a normalization constant. Keep in mind
that the Gaussians $\mathcal{N}_{s}$ for different configurations $s$ are
independent. The many-body spectrum of the TASEP is a specific sample of the
distribution in Eq. (56). For $U=0$ the random spectrum becomes deterministic
and agrees with the non-interacting many-body spectrum presented in Sec. III.

In Fig. 6(b), we present the probability density from Eq.(56) for $L=14$,
$N=7$, and $U=1$, with the density capped at $10^{-3}$ for clarity. The red
markers indicate the means $\sum_{j=1}^{L}s_{j}\mu_{j}$ of the many-body
Gaussians $\mathcal{N}_{s}$. Both the discrete means and the continuous
density exhibit pronounced spikes at the boundary. When these means are
compared to the TASEP many-body spectrum shown in Fig. 6(a), even finer
details of the spectrum are discernible in the structure of the means.

The boundary of the random many-body spectrum is mainly determined by
Gaussians $\mathcal{N}_{s}$, associated with domain wall configurations of one
or two domain walls, separated by at most one empty site, due to the
exponential decay of the Gaussian probability density function. These
configurations are identical to those defining the spectral boundary in the
non-interacting case.

The random Bethe spectrum and the TASEP spectrum share a remarkably similar
overall shape. However, differences do exist, e.g., the boundary of the random
Bethe spectrum is not skewed leftwards in the complex plane. This is
attributed to the additional structure in the Bethe root clusters seen in Fig.
5, not represented by rotationally invariant Gaussians.

###  V.4 “Thermodynamic limit”

Figure 7: The width $\sigma$ of the complex Gaussians fitted to the clusters
of the Bethe roots for $U=1$ at half-filling $N=L/2$. The solid line denotes
the average $\langle\sigma_{j}\rangle_{j}=\frac{1}{L}\sum_{j}\sigma_{j}$ of
the cluster widths and the dotted line guides the eye to $L^{-1/2}$. The inset
shows the absolute value of the centers of the complex Gaussians $|\mu|$.
Black solid line indicates the average.

Similar to the non-interacting case with $U=0$, we demonstrate that the spiky
boundary persists in the thermodynamic limit as $L$ and $N$ increase while
maintaining a fixed density $\rho=N/L$.

Let us first focus on the centers $\sum_{j=1}^{L}s_{j}\mu_{j}$ of the many-
body Gaussians $\mathcal{N_{s}}$, depicted as red dots in Fig. 6. According to
the inset of Fig. 7, the absolute values of $|\mu_{j}|$ appear to be
independent of $L$. This independence suggests that the non-interacting case
scenario also applies to the many-body Gaussian centers. For boundary
configurations $s$, these centers, being sums of $N=\rho L$ nearby $\mu_{j}$,
scale with $L$. Given that both the tip distance ($d_{t}$ from Sec.III.4) and
boundary depth ($d_{b}$ from Sec.III.4) are proportional to 1, the spiky
structure of the boundary Gaussian centers is maintained in the thermodynamic
limit.

However, this does not automatically mean that the spiky spectral boundary of
the random spectrum, as defined in Eq. (56), persists in the thermodynamic
limit. For this to hold true, the widths of the Gaussians $\mathcal{N}_{j}$ in
the mixture model must decrease sufficiently fast.

Fig. 7 displays the widths $\sigma_{j}$ of $\mathcal{N}_{j}$ for the TASEP
case ($U=1$) at half-filling ($N=L/2$), with $L$ ranging from 8 to 22. The
cluster widths $\sigma_{j}$ vary, being larger for clusters with smaller
$|\operatorname{Re}Z|$ and smaller for those with larger
$|\operatorname{Re}Z|$, as also observed in Fig. 5(e). Despite this variation,
the widths $\sigma_{j}$ are centered around their average
$\langle\sigma_{j}\rangle_{j}=\frac{1}{L}\sum_{j=1}^{L}\sigma_{j}$, which
decreases approximately as $\propto L^{-1/2}$, as shown by the dashed line in
Fig. 7. Consequently, the variance
$\sigma_{s}=\sum_{j=1}^{L}s_{j}\sigma_{j}^{2}$ of the Gaussians
$\mathcal{N}_{s}$ scales as $\propto 1$. This indicates that the standard
deviation of the boundary Gaussians $\mathcal{N}_{s}$ remains on the order of
$\propto 1$ even as $L$ increases, aligning with the scale of both the tip
distance and spike depths. Therefore, the spiky structure of the statistical
many-body spectrum for $U=1$ is preserved in the thermodynamic limit, as in
the $U=0$ case presented in Sec. III.

##  VI The random matrix picture

Figure 8: TASEP spectrum (pbc) with $L=12$ and $N=6$ for (a) $U=0$ and (c)
$U=1$. In (b,d) spectral density of random graphs with cycle length divisible
by $L$; in (b) of the adjacency matrix and in (d) of the negative
(combinatorial) Laplacian. In (e) traces of powers of the non-interacting
TASEP generator $H_{0}$ (squares) and random graph adjacency matrix (circles).

In the previous sections, we showed that the spikes of the spectral boundary
of the TASEP are a consequence of the many-body spectrum being generated by
summing single-particle-like clusters.

This section demonstrates that the spiky spectral boundary is a prevalent
characteristic in a broad range of systems, extending beyond free fermions or
those solvable by the Bethe ansatz. Specifically, this feature is typical in
systems where the many-body graph exhibits a particular cycle structure, with
cycle lengths being integer multiples of the spike count.

###  VI.1 From TASEP to graphs

The matrix elements of the generator of the non-interacting TASEP $H_{0}$ are
either zero or one. Thus the generator matrix is naturally interpreted as the
adjacency matrix of a directed graph. This graph, which we will call the many-
body graph of TASEP, has vertices representing particle configurations in the
chain and edges indicating permissible transitions. For TASEP with $U=1$, its
generator matrix $H$ is the negative combinatorial Laplacian of this graph.

###  VI.2 Cycles of TASEP

The permissible transitions between particle configurations impose constraints
on the structure of the many-body graph. Our focus is on the nature of cycles
in the many-body graph, which are closed walks with only the start and end
vertices being the same.

The cycle lengths in the TASEP many-body graph are divisible by $L$ for pbc
and by $L+1$ for obc [33]. This is evident in cycles among configurations,
which only contain a single particle. These cycles consist of $L$ particle
movements ($L+1$ for obc) such that the particle arrives at its original
position.

The number of closed walks with length $k$ is related to entries of the $k$th
power of the adjacency matrix $A$ ($A=H_{0}$ in the case of ASEP). The element
$(A^{k})_{ij}$ denotes the number of distinct walks of length $k$ from vertex
$i$ to $j$. Thus $(A^{k})_{ii}$ counts the number of distinct closed walks
with length $k$ starting and ending at vertex $i$ and
$\operatorname{tr}(A^{k})$ aggregates the total number of closed walks with
length $k$. Especially, if $\operatorname{tr}(A^{k})=0$ then the graph does
not contain any closed walks, thus any cycle, of length $k$.

In Fig. 8(e), we depict $\operatorname{tr}(A^{k})+1$ as blue squares, where
$A=H_{0}$, plotted against $k=1,\dots,2L$ for a system of $L=12$ sites and pbc
with $N=6$ particles. The addition of $+1$ facilitates a logarithmic scale on
the y-axis. Here, $\operatorname{tr}(A^{k})$ equals zero for all values of $k$
not divisible by $L$, indicating the absence of cycles in the graph with
length $k\ \mathrm{mod}\ L\neq 0$. Similarly, for obc,
$\operatorname{tr}(A)^{k}=0$ if and only if $k\ \operatorname{mod}\ L+1=0$
(not shown).

###  VI.3 Random graph model

To demonstrate the robustness of the spiky spectral boundary, we compare the
TASEP spectrum with the spectral density of a random graph ensemble
characterized only by cycles whose lengths are divisible by $L$. This
comparison is focused on the TASEP with pbc, noting that the obc scenario can
be similarly analyzed by simply adjusting $L$ to $L+1$.

We sample the random graph by initially forming a directed cycle with $D$
vertices. Next, we randomly choose a vertex and traverse the graph randomly
for $L-1$ steps. The vertex reached after $L-1$ steps is connected back to the
starting vertex, creating a cycle of length $L$. This process is repeated
until the graph contains a predetermined total number $n$ of edges.

Typically, the longest closed walk in the graph is the initial directed cycle
linking all $D$ vertices. When the number of vertices $D$ is divisible by $L$,
the construction of the graph ensures that all cycle lengths in the random
graph are also divisible by $L$.

Fig. 8 contrasts the random graph ensemble to the TASEP with $L=12$ sites and
pbc with $N=6$ particles. Quantities of the random graph ensemble are averaged
over $2,000$ samples, with the cycle length set to $L$ and the number of
vertices $D=924$, matching the Hilbert space dimension of the TASEP.

In Fig. 8(e) we present $\operatorname{tr}(A^{k})+1$ for the random graph
ensemble, shown as red circles. In this ensemble, $\operatorname{tr}(A^{k})$
is zero for all $k$ that are not integer multiples of $L$. Whenever $k$ is an
integer multiple of $L$, $\operatorname{tr}(A^{k})$ for the adjacency matrix
$A$ of the random graph ensemble is comparable in magnitude to
$\operatorname{tr}(A^{k})$ for $A=H_{0}$, the generator matrix of TASEP. This
similarity suggests that the number of closed walks in the random graph
ensemble is on par with that in the TASEP many-body graph.

Fig. 8(a-d) displays a comparison between the random graph ensemble and TASEP,
matching the parameters used in (e). In (a) and (b), we show the non-
interacting TASEP spectrum alongside the estimated spectral density of the
graph ensemble - both featuring $L$ distinct spikes.

In Fig. 8(c) and (d), the focus is on the spectrum of TASEP ($U=1$) and the
spectral density of the negative graph Laplacian for the random graph
ensemble. Notably, the random graph Laplacian also presents $L$ pronounced
spikes. The spike patterns, particularly their ”bending” towards the left,
show a resemblance to the TASEP spikes. The overall shape of the spectral
density (ignoring the spikes) takes on a spindle-like form, characteristic of
(sparse) random Markov matrices [97, 98, 99, 100].

##  VII Conclusion and Discussion

In this work, we explored the connections among the spectral problems for
ASEP, free fermion models, and random matrix theory, focusing particularly on
the distinctive spiky shape of the ASEP spectral boundary. We reformulated the
ASEP generator matrices as non-Hermitian fermionic models with a variable
interaction parameter $U$, where $U=1$ corresponds to the standard ASEP. We
analytically demonstrated that in the non-interacting ASEP ($U=0$), this spiky
spectral boundary arises from aggregating single-particle eigenvalues
positioned on ellipses (circles for TASEP). For pbc, we extended this concept
to interacting TASEP, showing that the spiky boundary remains and originates
from the summation of clustered Bethe roots. Lastly, we confirmed the
robustness of this spiky boundary by considering only the cycle structure in
the many-body graph, revealing that corresponding random graphs exhibit a
similar spiky spectral boundary.

This research opens up several questions for further exploration. We
demonstrated the spiky spectral boundary in TASEP, largely attributed to Bethe
roots clustering. It is intriguing to consider whether such clustering also
occurs in ASEP. The straightforward connection between TASEP and ASEP in their
non-interacting forms suggests that the spiky spectral boundary might extend
to standard ASEP (with $U=1$) as well. However, it remains to be seen how
introducing interactions influences Bethe roots clustering and the potential
emergence of a spiky spectral boundary.

In this study, we concentrated on the Bethe ansatz for pbc. The ASEP with obc
is also solvable via the Bethe ansatz, though the equations are more complex,
as detailed in various studies [101, 102, 82, 103, 104]. One might ask whether
the spiky spectral boundary in the obc case is also associated with a
clustering of Bethe roots similar to the pbc case.

The spectral boundary of random graphs with dominant cycle lengths typically
follows a hypotrochoidic curve, as noted in [37]. These graphs usually lack
cycles shorter than $L$ but can have cycles longer than $L$. However, the
random graph ensemble we introduced deviates from this standard hypotrochoidic
pattern, likely due to its more restricted cycle structure, where all cycles
are of lengths divisible by $L$. Extending the hypotrochoidic law to encompass
this specific graph ensemble would be a valuable advancement.

This study concentrated on the spiky spectral boundary of the ASEP. Formation
of spikes has as well been observed in the off-diagonals of reduced density
matrices in the symmetric simple exclusion process (XXX model) [105] and the
observable representation of Ising chain Glauber dynamics [106]. These
observations together with the robustness of the spiky spectral boundary to
perturbations make the the investigation of other models, both classical and
quantum, that possess a similar cycle structure in their many-body graphs or
comparable trace correlations in their generator matrices, an intriguing
direction for future research.

###### Acknowledgements.

GN and MH thank R. Taggart and P.C. Burke and GN thanks T. Giamarchi and I.
Lobaskin for helpful discussions. This research is supported by the Deutsche
Forschungsgemeinschaft through SFB No. 1143 (Project ID No. 247310070) (GN and
MH) and the Irish Research Council Government of Ireland Postgraduate
Scholarship Scheme (GOIPG/2019/58) (GN). TP is supported by the Grants
N1-0219, N1-0334 as well as Program P1-0402 of Slovenian Research and
Innovation Agency (ARIS).

##  Appendix A Quadratic Fermion Model for obc

In this section, we will show that the non-interacting ASEP $H_{0}$ with obc
is a quadratic fermion model. Especially, we will prove Eq. (11).

Recall

$H_{0}=\sum_{j=1}^{L-1}\left(p\sigma_{j+1}^{+}\sigma_{j}^{-}+q\sigma_{j}^{+}\sigma_{j+1}^{-}\right)+\alpha\sigma_{1}^{+}+\gamma\sigma_{1}^{-}+\beta\sigma_{L}^{-}+\delta\sigma_{L}^{+},$
(57)

For that, we first apply, as mentioned in the main text, the Kramers-Wannier
duality transformation [36]

$\sigma_{j}^{x}\to\prod_{l=1}^{j}\sigma_{l}^{z},\quad\sigma_{j}^{z}\to\sigma_{j}^{x}\sigma_{j+1}^{x},$
(58)

where we implicitly have enlarged the chain of length $L$ by one additional
site to a chain of length $L+1$. Thus the multiplicity of every eigenvalue of
the so-transformed $H_{0}$ is doubled. Applying a Jordan-Wigner transformation

$w_{j}=\left(\prod_{l=1}^{j-1}\sigma_{l}^{z}\right)\sigma_{j}^{-},\quad
w_{j}^{\dagger}=\left(\prod_{l=1}^{j-1}\sigma_{l}^{z}\right)\sigma_{j}^{+},$
(59)

and rewriting in terms of Majorana “real” and “imaginary” parts of the Dirac
fermions $w,w^{\dagger}$,

$\gamma_{j,1}=w_{j}^{\dagger}+w_{j},\quad\gamma_{j,2}=i(w_{j}^{\dagger}-w_{j}),$
(60)

the Hamiltonian $H_{0}$ is given by

$\displaystyle H_{0}=$
$\displaystyle\sum_{j=1}^{L-1}\left[\frac{p+q}{4}\left(i\gamma_{j+1,1}\gamma_{j+1,2}-i\gamma_{j,2}\gamma_{j+2,1}\right)+\frac{p-q}{4}\left(\gamma_{j+1,1}\gamma_{j+2,1}+\gamma_{j,2}\gamma_{j+1,2}\right)\right]$
$\displaystyle+\frac{1}{2}\left[(\alpha+\gamma)i\gamma_{1,1}\gamma_{1,2}+(\alpha-\gamma)\gamma_{1,1}\gamma_{2,1}\right]+\frac{1}{2}\left(\prod_{j=1}^{L+1}i\gamma_{j,1}\gamma_{j,2}\right)\left[(\delta+\beta)i\gamma_{L+1,1}\gamma_{L+1,2}-(\delta-\beta)\gamma_{L,2}\gamma_{L+1,2}\right].$
(61)

The string of Majoranas
$\prod_{j=1}^{L+1}\left(i\gamma_{j,1}\gamma_{j,2}\right)=(-1)^{L+1}\mathcal{P}_{w}$
equals, up to a sign, the parity operator $\mathcal{P}_{w}$ of Dirac fermions
$w,w^{\dagger}$, which commutes with $H_{0}$. Thus, restricted to the sub-
spaces of constant parity, the Hamiltonian $H_{0}$ becomes quadratic.

Note that $H_{0}$ in terms of the Majorana fermions $\gamma_{j,l}$ is acting
non-trivially on the additional site $L+1$.

To keep the algebra simpler let us consider from now on the case $p=1$ and
$q=\gamma=\delta=0$. The following calculations can be straightforwardly
generalized to arbitrary $p,q,\gamma,\delta$. Thus $H_{0}$ in terms of the
Majorana fermions $\gamma$ simplifies to

$\displaystyle H_{0}$
$\displaystyle=\sum_{j=1}^{L-1}\left(\sigma_{j+1}^{+}\sigma_{j}^{-}\right)+\alpha\sigma_{1}^{+}+\beta\sigma_{L}^{-}$
$\displaystyle=\frac{1}{2}\alpha[i\gamma_{1,1}\gamma_{1,2}+\gamma_{1,1}\gamma_{2,1}]$
$\displaystyle+\frac{1}{2}(-1)^{L+1}\mathcal{P}_{w}\beta[i\gamma_{L+1,1}\gamma_{L+1,2}+\gamma_{L,2}\gamma_{L+1,2}]$
$\displaystyle+\frac{1}{4}\sum_{j=1}^{L-1}\bigg{[}(\gamma_{j,2},\
\gamma_{j+1,1})\begin{pmatrix}1&-i\\\
i&1\end{pmatrix}\begin{pmatrix}\gamma_{j+1,2}\\\
\gamma_{j+2,1}\end{pmatrix}\bigg{]}.$ (62)

The eigenvalues of the $2\times 2$-matrix are 0 and 2, while the eigenvectors
are $(1,-i)^{t}$ and $(1,i)^{t}$, respectively. Thus the following pairing of
Majorana fermions

$c_{j}^{\dagger}=\frac{1}{2}(\gamma_{j,2}-i\gamma_{j+1,1}),\quad
c_{j}=\frac{1}{2}(\gamma_{j,2}+i\gamma_{j+1,1}),$ (63)

into Dirac fermions $c,c^{\dagger}$ drastically simplifies the bulk term. By
identifying $\gamma_{L+2,1}=\gamma_{1,1}$ the pairing given by Eq. (63) turns
the chain on sites 1 to $L+1$ into a ring, connecting site 1 and $L+1$. The
Hamiltonian $H_{0}$ is given in terms of $c,c^{\dagger}$ as

$\displaystyle H_{0}=$
$\displaystyle\alpha(c_{L+1}-c_{L+1}^{\dagger})c_{1}^{\dagger}+\sum_{j=1}^{L-1}\left[c_{j}c_{j+1}^{\dagger}\right]$
$\displaystyle+(-1)^{L}\mathcal{P}_{c}\beta c_{L}(c_{L+1}+c_{L+1}^{\dagger}),$

where $\mathcal{P}_{c}$ denotes the parity of the Dirac fermions
$c,c^{\dagger}$. This is Eq.(11).

##  Appendix B Diagonalizing $M_{c}$

In this section we calculate the eigenvalues and eigenvectors of $M_{c}$ given
by Eq. (30) and Eqs.(31-33) thereafter. We denote the eigenvalue equation by
$M_{c}\boldsymbol{u}=\lambda\boldsymbol{u}$ with the $2L+2$ dimensional vector
$\boldsymbol{u}$. In terms of
$\boldsymbol{u}=(u_{1},\dots,u_{L+1},u_{1}^{\prime},\dots,u_{L+1}^{\prime})$
the eigenvalue equation reads

$\displaystyle\lambda u_{1}$ $\displaystyle=-\alpha(u_{L+1}-u_{L+1}^{\prime})$
(64) $\displaystyle\lambda u_{2}$ $\displaystyle=-u_{1}$ (65)
$\displaystyle\dots$ $\displaystyle\lambda u_{L-1}$ $\displaystyle=-u_{L-2}$
(66) $\displaystyle\lambda u_{L}$ $\displaystyle=-u_{L-1}$ (67)
$\displaystyle\lambda u_{L+1}$ $\displaystyle=-s\beta u_{L}-\alpha
u_{1}^{\prime}$ (68)

and

$\displaystyle\lambda u_{1}^{\prime}$ $\displaystyle=u_{2}^{\prime}$ (69)
$\displaystyle\lambda u_{2}^{\prime}$ $\displaystyle=u_{3}^{\prime}$ (70)
$\displaystyle\dots$ $\displaystyle\lambda u_{L-1}^{\prime}$
$\displaystyle=u_{L}^{\prime}$ (71) $\displaystyle\lambda u_{L}^{\prime}$
$\displaystyle=\beta s(u_{L+1}+u_{L+1}^{\prime})$ (72) $\displaystyle\lambda
u_{L+1}^{\prime}$ $\displaystyle=-\beta su_{L}+\alpha u_{1}^{\prime}.$ (73)

Combining the Eqs. (65)-(67) with $u_{1},\dots,u_{L}$ and Eqs. (69)-(71) with
$u_{1}^{\prime},\dots,u_{L}^{\prime}$ recursively we get for $2\leq j\leq L$

$u_{j}=-\lambda^{-1}u_{j-1}=\dots=(-\lambda)^{-j+1}u_{1}$ (74)

and

$u_{j}^{\prime}=\lambda u_{j-1}^{\prime}=\dots=\lambda^{j-1}u_{1}^{\prime}.$
(75)

Especially, the following holds

$\displaystyle u_{L}$ $\displaystyle=(-\lambda)^{-L+1}u_{1}$ (76)
$\displaystyle u_{L}^{\prime}$ $\displaystyle=\lambda^{L-1}u_{1}^{\prime}.$
(77)

By substituting Eq. (76) and Eq. (77) into Eq. (68) and Eq. (73),
respectively, we get

$\displaystyle u_{1}$
$\displaystyle=\alpha\lambda^{-1}(-u_{L+1}+u^{\prime}_{L+1})$ (78)
$\displaystyle u_{L+1}$ $\displaystyle=\beta
s(-\lambda)^{-L}u_{1}-\alpha\lambda^{-1}u_{1}^{\prime}$ (79) $\displaystyle
u_{1}^{\prime}$ $\displaystyle=\lambda^{-L}\beta s(u_{L+1}+u_{L+1}^{\prime})$
(80) $\displaystyle u_{L+1}^{\prime}$ $\displaystyle=\beta
s(-\lambda)^{-L}u_{1}+\alpha\lambda^{-1}u_{1}^{\prime}.$ (81)

Adding and subtracting Eq. (79) and Eq. (81), respectively, leads to

$\displaystyle u_{L+1}+u_{L+1}^{\prime}$ $\displaystyle=2\beta
s(-\lambda)^{-L}u_{1}$ (82) $\displaystyle-u_{L+1}+u_{L+1}^{\prime}$
$\displaystyle=2\alpha\lambda^{-1}u_{1}^{\prime},$ (83)

which in turn implies that

$\displaystyle u_{1}^{\prime}$
$\displaystyle=2(-1)^{L}\lambda^{-2L}\beta^{2}u_{1}$ (84) $\displaystyle
u_{1}$ $\displaystyle=2\alpha^{2}\lambda^{-2}u_{1}^{\prime},$ (85)

by using Eqs. (78) and (80). Combining the last two equations leads to

$u_{1}=4(\alpha\beta)^{2}(-1)^{L}\lambda^{-2L-2}u_{1},$ (86)

which implies, for $u_{1}\neq 0$, the eigenvalue Eq. (35)

$\lambda^{2(L+1)}=(-1)^{L}4(\alpha\beta)^{2}.$ (87)

The roots of this polynomial are given by

$\lambda=(2\alpha\beta)^{\frac{1}{L+1}}\begin{cases}\exp\left(\frac{i\pi}{2L+2}2k\right)&L\text{
even,}\\\ \exp\left(\frac{i\pi}{2L+2}(2k-1)\right)&L\text{ odd,}\end{cases}$

where $k=1,\dots,2L+2$.

##  Appendix C Bethe equations for ASEP with pbc

In this section, we derive the Bethe equations presented in Sec. V.1. These
results extend the usual $U=1$ ASEP Bethe ansatz [32] to the case of arbitrary
$U$.

By $\left|x_{1},\dots,x_{N}\right\rangle$ we denote the state of $N$ particles
at position $x_{1},\dots,x_{N}$. In the following, we let $x_{1}<\dots<x_{N}$
up to an overall shift in the indices. The wavefunction
$\left|\psi\right\rangle$ in the basis of
$\left|x_{1},\dots,x_{N}\right\rangle$ is given by

$\left|\psi\right\rangle=\sum_{x_{1}<\dots<x_{N}}\psi(x_{1},\dots,x_{N})\left|x_{1},\dots,x_{N}\right\rangle,$
(88)

where $\psi(x_{1},\dots,x_{N})$ denotes the coefficient of
$\left|\psi\right\rangle$ with respect to
$\left|x_{1},\dots,x_{N}\right\rangle$. Now, let $\left|\psi\right\rangle$ be
an eigenstate of the generalized Markov matrix $H$ with eigenvalue $E$, i.e.
$H\left|\psi\right\rangle=E\left|\psi\right\rangle$. Recall that we can write
the generator matrix $H$ as

$H=\sum_{i=1}^{L}\left(p\sigma_{i}^{-}\sigma_{i+1}^{+}+q\sigma_{i}^{+}\sigma_{i+1}^{-}\right)+\frac{U}{4}\sum_{i=1}^{L}\left(\sigma_{i}^{z}\sigma_{i+1}^{z}-1\right).$
(89)

Let us first focus on the action of the off-diagonal term in eq. (89) on
$\left|x_{1},\dots,x_{N}\right\rangle$. It is easy to see that

$\sum_{i=1}^{L-1}\sigma_{i}^{-}\sigma_{i+1}^{+}\left|x_{1},\dots,x_{N}\right\rangle=\sum_{j=1}^{N-1}(1-\delta(x_{j+1}-x_{j},1))\left|x_{1},\dots,x_{j}+1,\dots,x_{N}\right\rangle,$
(90)

and

$\sum_{i=1}^{L-1}\sigma_{i}^{+}\sigma_{i+1}^{-}\left|x_{1},\dots,x_{N}\right\rangle=\sum_{j=2}^{N}(1-\delta(x_{j}-x_{j-1},1))\left|x_{1},\dots,x_{j}-1,\dots,x_{N}\right\rangle,$
(91)

where $\delta(x,y)$ equals one whenever $x=y$ and is zero otherwise. The
remaining boundary terms are determined as follows. If $x_{N}\neq L$ then
$\sigma_{L}^{-}\sigma_{1}^{+}\left|x_{1},\dots,x_{N}\right\rangle=0$, so let
$x_{N}=L$. Then

$\displaystyle\sigma_{L}^{-}\sigma_{1}^{+}\left|x_{1},\dots,x_{N}\right\rangle$
$\displaystyle=(1-\delta(x_{1},1))\left|1,x_{1},\dots,x_{N-1}\right\rangle$
(92)
$\displaystyle=(1-\delta(x_{1}-x_{N}\pmod{L},1))\left|x_{1},\dots,x_{N-1},X_{N}+1\right\rangle$
(93)

by identifying
$\left|x_{1},\dots,x_{N-1},L+1\right\rangle=\left|1,x_{1},\dots,x_{N-1}\right\rangle$.
On the other hand, whenever $x_{1}\neq 1$ we have
$\sigma_{L}^{+}\sigma_{1}^{-}\left|x_{1},\dots,x_{N}\right\rangle=0$, while
for $x_{1}=1$ we get

$\displaystyle\sigma_{L}^{+}\sigma_{1}^{-}\left|x_{1},\dots,x_{N}\right\rangle$
$\displaystyle=(1-\delta(x_{N},L))\left|x_{2},\dots,x_{N},L\right\rangle$ (94)
$\displaystyle=(1-\delta(x_{1}-x_{N}\pmod{L},1))\left|x_{1}-1,x_{2},\dots,x_{N}\right\rangle,$
(95)

where we identified
$\left|x_{2},\dots,x_{N},L\right\rangle=\left|0,x_{2},\dots,x_{N}\right\rangle$.
Taking everything together we have

$\displaystyle\sum_{i=1}^{L}\sigma_{i}^{-}\sigma_{i+1}^{+}\left|x_{1},\dots,x_{N}\right\rangle$
$\displaystyle=\sum_{j=1}^{N}(1-\delta(x_{j+1}-x_{j}\pmod{L},1))\left|x_{1},\dots,x_{j}+1,\dots,x_{N}\right\rangle,$
(96)
$\displaystyle\sum_{i=1}^{L}\sigma_{i}^{+}\sigma_{i+1}^{-}\left|x_{1},\dots,x_{N}\right\rangle$
$\displaystyle=\sum_{j=1}^{N}(1-\delta(x_{j}-x_{j-1}\pmod{L},1))\left|x_{1},\dots,x_{j}-1,\dots,x_{N}\right\rangle.$
(97)

The diagonal term in Eq. (89) acts on $\left|x_{1},\dots,x_{N}\right\rangle$
as

$\displaystyle\sum_{i=1}^{L}\left(\sigma_{i}^{z}\sigma_{i+1}^{z}-1\right)\left|x_{1},\dots,x_{N}\right\rangle$
$\displaystyle=\sum_{j=1}^{N-1}\delta(x_{j+1}-x_{j},1)+\delta(x_{1}-x_{N},1-L)-N$
$\displaystyle=\sum_{j=1}^{N}\delta(x_{j+1}-x_{j}\pmod{L},1)-N,$ (98)

where we note that $\sigma^{z}_{i}=2n_{i}-1$ and thus

$\displaystyle\frac{1}{4}\sum_{i=1}^{L}\left(\sigma_{i}^{z}\sigma_{i+1}^{z}-1\right)$
$\displaystyle=\frac{4}{4}\sum_{i=1}^{L}n_{i}n_{i+1}-\frac{2}{4}\sum_{i=1}^{L}n_{i}-\frac{2}{4}\sum_{i=1}^{L}n_{i+1}=\left[\sum_{i=1}^{L}n_{i}n_{i+1}\right]-N.$
(99)

Summarizing, the action of $H$ on $\left|x_{1},\dots,x_{N}\right\rangle$ is

$\displaystyle H\left|x_{1},\dots,x_{N}\right\rangle$
$\displaystyle=p\sum_{j=1}^{N}(1-\delta(x_{j+1}-x_{j}\pmod{L},1))\left|x_{1},\dots,x_{j}+1,\dots,x_{N}\right\rangle$
$\displaystyle+q\sum_{j=1}^{N}(1-\delta(x_{j}-x_{j-1}\pmod{L},1))\left|x_{1},\dots,x_{j}-1,\dots,x_{N}\right\rangle$
$\displaystyle-U\sum_{j=1}^{N}(1-\delta(x_{j+1}-x_{j}\pmod{L},1))\left|x_{1},\dots,x_{N}\right\rangle.$
(100)

Now, consider the eigenvalue equation
$H\left|\psi\right\rangle=E\left|\psi\right\rangle$,

$\displaystyle H\left|\psi\right\rangle$
$\displaystyle=\sum_{x_{1}<\dots<x_{N}}\psi(x_{1},\dots,x_{N})H\left|x_{1},\dots,x_{N}\right\rangle=\sum_{x_{1}<\dots<x_{N}}\psi(x_{1},\dots,x_{N})E\left|x_{1},\dots,x_{N}\right\rangle.$
(101)

Let us concentrate on the term in Eq. (C) proportional to $p$

$\displaystyle
p\sum_{j=1}^{N}\sum_{x_{1}<\dots<x_{N}}\psi(x_{1},\dots,x_{N})(1-\delta(x_{j+1}-x_{j}\pmod{L},1))\left|x_{1},\dots,x_{j}+1,\dots,x_{N}\right\rangle.$
(102)

After a change of variables $\tilde{x}_{i}=x_{i}$ for $i\neq j$ and
$\tilde{x}_{j}=x_{j}+1$ the above equation reads

$\displaystyle
p\sum_{j=1}^{N}\sum_{\tilde{x}_{1}<\dots<\tilde{x}_{N}}\psi(\tilde{x}_{1},\dots,\tilde{x}_{j}-1,\dots,\tilde{x}_{N})(1-\delta(\tilde{x}_{j}-\tilde{x}_{j-1}\pmod{L},1))\left|\tilde{x}_{1},\dots,\tilde{x}_{N}\right\rangle.$
(103)

Let us now focus on the term in Eq. (C) proportional to $q$. One finds with
the change $\tilde{x}_{j}=x_{j}-1<x_{j+1}-1=\tilde{x}_{j+1}-1$, thus
$\tilde{x}_{j+1}-\tilde{x}_{j}>1$ and
$\tilde{x}_{j-1}=x_{j-1}<x_{j}-1=\tilde{x}_{j}$, that this term equals

$\displaystyle
q\sum_{\tilde{x}_{1}<\dots<\tilde{x}_{N}}\psi(\tilde{x}_{1},\dots,\tilde{x}_{j}+1,\dots,\tilde{x}_{N})(1-\delta(\tilde{x}_{j+1}-\tilde{x}_{j}\pmod{L},1))\left|\tilde{x}_{1},\dots,\tilde{x}_{N}\right\rangle,$
(104)

where the first constraint is realized via the delta term and the second
constraint by the summation.

By orthogonality of $\left|x_{1},\dots,x_{N}\right\rangle$ the eigenvalue
equation $H\left|\psi\right\rangle=E\left|\psi\right\rangle$ turns into
$\binom{L}{N}$ equations for the wavefunction coefficients

$\displaystyle p\sum_{j=1}^{N}\big{(}1-\delta(x_{j}-x_{j-1}\mod
L,1))(\psi(x_{1},\dots,x_{j}-1,\dots,x_{N})-U\psi(x_{1},\dots,x_{N})\big{)}$
$\displaystyle+$ $\displaystyle
q\sum_{j=1}^{N}\big{(}1-\delta(x_{j+1}-x_{j}\mod
L,1))(\psi(x_{1},\dots,x_{j}+1,\dots,x_{N})-U\psi(x_{1},\dots,x_{N})\big{)}$
$\displaystyle=E\psi(x_{1},\dots,x_{N}).$ (105)

Here we additionally used that $p+q=1$. Now, we make the ansatz for the
wavefunction coefficient

$\psi(x_{1},\dots,x_{N})=\sum_{\tau\in
S_{N}}A(\tau)\prod_{j=1}^{N}z_{\tau(j)}^{x_{j}},$ (106)

where the summation runs over all elements of the symmetric group $S_{N}$ and
the $z_{j}$’s and $A(\tau)$’s are complex numbers. Let us consider a
configuration $x_{1}<\dots<x_{N}$ where all particles have at least distance
1, i.e. no consecutive particles. Plugging the ansatz into the term
proportional to $p$ results in

$\displaystyle\sum_{j=1}^{N}\psi(x_{1},\dots,x_{j}-1,\dots,x_{N})-U\psi(x_{1},\dots,x_{N})$
$\displaystyle=\sum_{j=1}^{N}\sum_{\tau\in
S_{N}}A(\tau)\left(z_{\tau(j)}^{x_{j}-1}\prod_{l=1;l\neq
j}^{N}z_{\tau(l)}^{x_{l}}-U\prod_{l=1}^{N}z_{\tau(l)}^{x_{l}}\right)$ (107)
$\displaystyle=\sum_{\tau\in
S_{N}}A(\tau)\prod_{l=1}^{N}z_{\tau(l)}^{x_{l}}\sum_{j=1}^{N}\left(z_{\tau(j)}^{-1}-U\right).$
(108)

Similarly, one gets the analogous expression for the term proportional to $q$
with the change $z_{\tau(j)}^{-1}\to z_{\tau(j)}$. Thus Eq. (105) in terms of
the Bethe ansatz reads

$E=\sum_{j=1}\left(pz_{j}+qz_{j}^{-1}-U\right),$ (109)

which is Eq. (49) in Sec. V.1. Now, consider a configuration
$\left|x_{1},\dots,x_{N}\right\rangle$ with two particles adjacent to each
other. Then

$A(\dots,l,\dots,k,\dots)=-\frac{p+qz_{l}z_{k}-Uz_{l}}{p+qz_{l}z_{k}-Uz_{k}}A(\dots,k,\dots,l,\dots).$
(110)

The periodic boundary condition enforces
$\psi(x_{1},\dots,x_{N-1},L+1)=\psi(1,x_{1},\dots,x_{N-1})$, which implies

$A(\tau(1),\dots,\tau(N))z_{\tau(N)}^{L}=A(\tau(N),\tau(1),\dots,\tau(N-1)).$
(111)

Combining both boundary constraints leads to the Bethe equations

$z_{j}^{L}=\prod_{k=1;k\neq
j}^{N}\left(-\frac{p+qz_{j}z_{k}-Uz_{j}}{p+qz_{j}z_{k}-Uz_{k}}\right).$ (112)

In the case of TASEP with $q=1$ and $p=0$ the Bethe equations reduce to

$\displaystyle z_{j}^{L}$ $\displaystyle=\prod_{k=1;k\neq
j}^{N}\left(-\frac{z_{j}z_{k}-Uz_{j}}{z_{j}z_{k}-Uz_{k}}\right)=\frac{z_{j}^{N}}{(z_{j}-U)^{N}}(-1)^{N-1}\prod_{k=1}^{N}\frac{z_{k}-U}{z_{k}},$
(113)

so

$z_{j}^{L-N}(z_{j}-U)^{N}=(-1)^{N-1}\prod_{k=1}^{N}\frac{z_{k}-U}{z_{k}}.$
(114)

Denoting $Z_{k}=2z_{k}-U$ we get

$(U+Z_{j})^{L-N}(U-Z_{j})^{N}=-2^{L}\prod_{k=1}^{N}\frac{Z_{k}-U}{Z_{k}+U}.$
(115)

##  Appendix D Solving the Bethe equations numerically

Figure 9: Visualization of solving the Bethe Eqs. (52) of TASEP (pbc) for
$L=6$ and $N=3$. All markers are roots of the polynomial $P$ (Eq. (116)) for
different $Y$. The outer (square) markers are the roots for initial
$Y^{(1)}=10\times 2^{L}$, the inner (triangles) markers for $Y$ converged, and
the circles denote roots of $P$ for intermediate $Y$. Red markers (upper
complex plane) are chosen to calculate the next $Y$. Gray circle has radius
$|Y^{(1)}|^{1/L}$.

In this section, we will describe how to self-consistently solve the Bethe
equations numerically. We will mostly follow the approach in [9] with some
additional tweaks.

Restricting to $p=1$ and $q=0$ reduces the difficulty of solving the Bethe
equations considerably because the right-hand side of Eq. (52) does not depend
on $j$, as does the right-hand side of Eq. (50) for general $p,q$.

Consider the polynomial $P(z)$,

$P(z)=\left(U+z\right)^{L-N}\left(U-z\right)^{N}-Y,$ (116)

where $Y$ denotes an arbitrary complex number and let us denote the right-hand
side of Eq. (52) by

$\tilde{Y}(Z_{1},\dots,Z_{N})=2^{L}\prod_{k=1}^{N}\frac{Z_{k}-U}{Z_{k}+U}.$
(117)

Then every solution $Z_{1},\dots,Z_{N}$ of Eq. (52) are roots of the
polynomial $P$ with $Y=\tilde{Y}(Z_{1},\dots,Z_{N})$. To find a solution to
the Bethe equations one first calculates the roots
$Z_{1}^{(1)},\dots,Z_{L}^{(1)}$ of $P$ for an initial $Y^{(1)}$. Of these $L$
roots of $P$ one chooses $N$ roots, $Z_{1}^{(1)},\dots,Z_{N}^{(1)}$, and
evaluates the next $Y^{(2)}=\tilde{Y}(Z_{1}^{(1)},\dots,Z_{N}^{(1)})$. Again,
the roots $Z_{1}^{(2)},\dots,Z_{L}^{(2)}$ of $P$ with $Y=Y^{(2)}$ are
calculated and $N$ roots $Z_{1}^{(1)},\dots,Z_{N}^{(1)}$ are chosen to
evaluate the next $Y^{(3)}=\tilde{Y}(Z_{1}^{(2)},\dots,Z_{N}^{(2)})$. This
procedure is then iterated until convergence all of the $N$ chosen roots is
reached, $Z_{j}^{(l)}\approx Z_{j}^{(l+1)}$ for all $1\leq j\leq N$.

The convergence of this procedure presupposes consistency of the choice of the
$N$ roots out of $L$ roots of the polynomial $P$ [51, 34]. The first choice of
$Z_{1}^{(1)},\dots,Z_{N}^{(1)}$ out of $Z_{1}^{(1)},\dots,Z_{L}^{(1)}$ is
arbitrary. Subsequent roots $Z_{1}^{(l)},\dots,Z_{N}^{(l)}$ are chosen to be
closest to the previous roots

$Z_{j}^{(1)}=\operatorname*{argmin}_{Z_{k}^{(l)}:1\leq k\leq
L}|Z_{k}^{(l)}-Z_{j}^{(l-1)}|,$ (118)

where the minimum runs over all roots $Z_{1}^{(l)},\dots,Z_{L}^{(l)}$ of $P$
with $Y=Y^{(l)}$. If multiple $Z_{k}^{(l)}$ are close to $Z_{j}^{(l-1)}$ we do
not update $Y^{(l+1)}$ with $Z_{j}^{(l)}$ but with a linear combination of
$Z_{j}^{(l)}$ and $Z_{j}^{(l-1)}$, i.e. $Y^{(l+1)}=\tilde{Y}(\dots,dY\
Z_{j}^{(l)}+(1-dY)Z_{j}^{(l-1)},\dots)$ where $0<dY\leq 1$ denotes the
fraction of interpolation between $Z_{j}^{(l)}$ and $Z_{j}^{(l-1)}$.

The above-described procedure typically leads to convergence of
$Z_{1}^{(l)},\dots,Z_{N}^{(l)}$ and thus to a solution of the Bethe Eqs. (52).
In Fig. 9 we show the roots $Z_{1}^{(l)},\dots,Z_{6}^{(l)}$ obtained during
the above algorithm for $L=6$ and $N=3$. The square markers denote the initial
$Z_{1}^{(1)},\dots,Z_{6}^{(1)}$ with $Y^{(1)}=10\times 2^{L}$, while the
triangles denote the final and converged $Z_{1}^{(end)},\dots,Z_{6}^{(end)}$
(relative or absolute error of Eq. (52) $<10^{-3}$). The circles indicate
intermediate roots. Initially, the 3 red squares (upper half-plane) are chosen
as $Z_{1}^{(1)},\dots,Z_{3}^{(1)}$, and subsequent roots (upper half-plane in
red) according to their previous closest roots. For visualization purposes,
$dY$ was chosen to be $dY=0.5$.

To find all solutions to the Bethe Eqs. (52) systematically we use different
combinations of initial $Y^{(1)}$ and initial root choices. Namely, we
typically choose $Y^{(1)}$ with $|Y^{(1)}|^{1/L}\gg U$. This ensures that the
roots of $P$ with $Y=Y^{(1)}$ are close to the circle with radius
$|Y^{(1)}|^{1/L}$. In Fig. 9 the roots of $P$ for $Y=Y^{(1)}=10\times 2^{6}$
denoted by the square markers are close to the circle with radius $2\times
10^{1/6}\approx 2.9$. Then we solve the Bethe equations for every combination
of $N$ roots out of $L$. This typically gives us almost all solutions of the
Bethe Eqs. (52). By iterating this procedure for a handful of initial
$Y^{(1)}$ we found all Bethe roots for the systems we investigated (up to
$L=22$).

## References

  * Spitzer [1970] F. Spitzer, Interaction of Markov processes, Advances in Mathematics 5, 246 (1970).
  * Liggett [1985] T. M. Liggett,  _Interacting Particle Systems_ , Vol. 2 (Springer, 1985).
  * Spohn [1991] H. Spohn,  _Large Scale Dynamics of Interacting Particles_ (Springer, 1991).
  * Schütz and Domany [1993] G. Schütz and E. Domany, Phase transitions in an exactly soluble one-dimensional exclusion process, Journal of Statistical Physics 72, 277 (1993).
  * Derrida  _et al._ [1993a] B. Derrida, M. R. Evans, V. Hakim, and V. Pasquier, Exact solution of a 1D asymmetric exclusion model using a matrix formulation, Journal of Physics A: Mathematical and General 26, 1493 (1993a).
  * Derrida [1998] B. Derrida, An exactly soluble non-equilibrium system: The asymmetric simple exclusion process, Physics Reports 301, 65 (1998).
  * Liggett [1999] T. M. Liggett,  _Stochastic Interacting Systems: Contact, Voter and Exclusion Processes_, Vol. 324 (Springer, 1999).
  * Schütz [2001] G. M. Schütz, Exactly solvable models for many-body systems far from equilibrium, Phase Transitions and Critical Phenomena 19, 1 (2001).
  * Golinelli and Mallick [2006] O. Golinelli and K. Mallick, The asymmetric simple exclusion process: an integrable model for non-equilibrium statistical mechanics, Journal of Physics A: Mathematical and General 39, 12679 (2006).
  * Chou  _et al._ [2011] T. Chou, K. Mallick, and R. K. P. Zia, Non-equilibrium statistical mechanics: from a paradigmatic model to biological transport, Reports on Progress in Physics 74, 116601 (2011).
  * Mallick [2015] K. Mallick, The exclusion process: A paradigm for non-equilibrium behaviour, Physica A: Statistical Mechanics and its Applications 418, 17 (2015).
  * Krug [1991] J. Krug, Boundary-induced phase transitions in driven diffusive systems, Physical Review Letters 67, 1882 (1991).
  * Sasamoto [1999] T. Sasamoto, One-dimensional partially asymmetric simple exclusion process with open boundaries: orthogonal polynomials approach, Journal of Physics A: Mathematical and General 32, 7109 (1999).
  * Blythe  _et al._ [2000] R. A. Blythe, M. R. Evans, F. Colaiori, and F. H. L. Essler, Exact solution of a partially asymmetric exclusion model using a deformed oscillator algebra, Journal of Physics A: Mathematical and General 33, 2313 (2000).
  * Parmeggiani  _et al._ [2003] A. Parmeggiani, T. Franosch, and E. Frey, Phase coexistence in driven one-dimensional transport, Physical Review Letters 90, 086601 (2003).
  * Blythe and Evans [2007] R. A. Blythe and M. R. Evans, Nonequilibrium steady states of matrix-product form: a solver’s guide, Journal of Physics A: Mathematical and Theoretical 40, R333 (2007).
  * Derrida  _et al._ [1993b] B. Derrida, S. A. Janowsky, J. L. Lebowitz, and E. R. Speer, Exact solution of the totally asymmetric simple exclusion process: Shock profiles, Journal of Statistical Physics 73, 813 (1993b).
  * Mallick [1996] K. Mallick, Shocks in the asymmetry exclusion model with an impurity, Journal of Physics A: Mathematical and General 29, 5375 (1996).
  * Jafarpour [2005] F. H. Jafarpour, Multiple shocks in a driven diffusive system with two species of particles, Physica A: Statistical Mechanics and its Applications 358, 413 (2005).
  * Kim  _et al._ [2011] M. Kim, L. Santen, and J. D. Noh, Asymmetric simple exclusion process in one-dimensional chains with long-range links, Journal of Statistical Mechanics: Theory and Experiment 2011, P04003 (2011).
  * Arita [2015] C. Arita, Synchronized shocks in an inhomogeneous exclusion process, Europhysics Letters 112, 40001 (2015).
  * MacDonald and Gibbs [1969] C. T. MacDonald and J. H. Gibbs, Concerning the kinetics of polypeptide synthesis on polyribosomes, Biopolymers 7, 707 (1969).
  * MacDonald  _et al._ [1968] C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin, Kinetics of biopolymerization on nucleic acid templates, Biopolymers 6, 1 (1968).
  * Bressloff and Newby [2013] P. C. Bressloff and J. M. Newby, Stochastic models of intracellular transport, Reviews of Modern Physics 85, 135 (2013).
  * Neri  _et al._ [2013] I. Neri, N. Kern, and A. Parmeggiani, Modeling cytoskeletal traffic: An interplay between passive diffusion and active transport, Physical Review Letters 110, 098102 (2013).
  * Fang  _et al._ [2019] X. Fang, K. Kruse, T. Lu, and J. Wang, Nonequilibrium physics in biology, Reviews of Modern Physics 91, 045004 (2019).
  * Helbing [2001] D. Helbing, Traffic and related self-driven many-particle systems, Reviews of Modern Physics 73, 1067 (2001).
  * Karzig and Von Oppen [2010] T. Karzig and F. Von Oppen, Signatures of critical full counting statistics in a quantum-dot chain, Physical Review B 81, 045317 (2010).
  * Meakin  _et al._ [1986] P. Meakin, P. Ramanlal, L. M. Sander, and R. C. Ball, Ballistic deposition on surfaces, Physical Review A 34, 5091 (1986).
  * Kardar  _et al._ [1986] M. Kardar, G. Parisi, and Y.-C. Zhang, Dynamic scaling of growing interfaces, Phys. Rev. Lett. 56, 889 (1986).
  * Krug and Spohn [1991] J. Krug and H. Spohn, Kinetic roughening of growing surfaces, in  _Solids far from equilibrium_ (Cambridge University Press, 1991) p. 479.
  * Gwa and Spohn [1992] L. H. Gwa and H. Spohn, Bethe solution for the dynamical-scaling exponent of the noisy Burgers equation, Physical Review A 46, 844 (1992).
  * Prolhac [2013] S. Prolhac, Spectrum of the totally asymmetric simple exclusion process on a periodic lattice—bulk eigenvalues, Journal of Physics A: Mathematical and Theoretical 46, 415001 (2013).
  * Prolhac [2015] S. Prolhac, Current fluctuations and large deviations for periodic TASEP on the relaxation scale, Journal of Statistical Mechanics: Theory and Experiment 2015, P11028 (2015).
  * Essler and Rittenberg [1996] F. H. L. Essler and V. Rittenberg, Representations of the quadratic algebra and partially asymmetric diffusion with open boundaries, Journal of Physics A: Mathematical and General 29, 3375 (1996).
  * Kogut [1979] J. B. Kogut, An introduction to lattice gauge theory and spin systems, Reviews of Modern Physics 51, 659 (1979).
  * Aceituno  _et al._ [2019] P. V. Aceituno, T. Rogers, and H. Schomerus, Universal hypotrochoidic law for random matrices with cyclic correlations, Physical Review E 100, 010302 (2019).
  * Prosen [2008] T. Prosen, Third quantization: a general method to solve master equations for quadratic open Fermi systems, New Journal of Physics 10, 043026 (2008).
  * Dhar [1987] D. Dhar, An exactly solved model for interfacial growth, Phase Transitions 9, 51 (1987).
  * Kim [1995] D. Kim, Bethe ansatz solution for crossover scaling functions of the asymmetric XXZ chain and the Kardar-Parisi-Zhang-type growth model, Physical Review E 52, 3512 (1995).
  * Derrida and Lebowitz [1998] B. Derrida and J. L. Lebowitz, Exact large deviation function in the asymmetric exclusion process, Physical Review Letters 80, 209 (1998).
  * Golinelli and Mallick [2004a] O. Golinelli and K. Mallick, Bethe ansatz calculation of the spectral gap of the asymmetric exclusion process, Journal of Physics A: Mathematical and General 37, 3321 (2004a).
  * Golinelli and Mallick [2004b] O. Golinelli and K. Mallick, Hidden symmetries in the asymmetric exclusion process, Journal of Statistical Mechanics: Theory and Experiment 2004, P12001 (2004b).
  * Golinelli and Mallick [2005a] O. Golinelli and K. Mallick, Spectral gap of the totally asymmetric exclusion process at arbitrary filling, Journal of Physics A: Mathematical and General 38, 1419 (2005a).
  * Golinelli and Mallick [2005b] O. Golinelli and K. Mallick, Spectral degeneracies in the totally asymmetric exclusion process, Journal of Statistical Physics 120, 779 (2005b).
  * Appert-Rolland  _et al._ [2008] C. Appert-Rolland, B. Derrida, V. Lecomte, and F. van Wijland, Universal cumulants of the current in diffusive systems on a ring, Physical Review E 78, 021122 (2008).
  * Prolhac and Mallick [2008] S. Prolhac and K. Mallick, Current fluctuations in the exclusion process and Bethe ansatz, Journal of Physics A: Mathematical and Theoretical 41, 175002 (2008).
  * Prolhac [2010] S. Prolhac, Tree structures for the current fluctuations in the exclusion process, Journal of Physics A: Mathematical and Theoretical 43, 105002 (2010).
  * Mallick [2011] K. Mallick, Some exact results for the exclusion process, Journal of Statistical Mechanics: Theory and Experiment 2011, P01024 (2011).
  * Simon [2011] D. Simon, Bethe ansatz for the weakly asymmetric simple exclusion process and phase transition in the current distribution, Journal of Statistical Physics 142, 931 (2011).
  * Motegi  _et al._ [2012] K. Motegi, K. Sakai, and J. Sato, Exact relaxation dynamics in the totally asymmetric simple exclusion process, Physical Review E 85, 042105 (2012).
  * Corwin [2014] I. Corwin, Two ways to solve ASEP, in  _Topics in Percolative and Disordered Systems_, Springer Proceedings in Mathematics & Statistics (Springer, New York, 2014) pp. 1–13.
  * Prolhac [2014] S. Prolhac, Spectrum of the totally asymmetric simple exclusion process on a periodic lattice-first excited states, Journal of Physics A: Mathematical and Theoretical 47, 375001 (2014).
  * Prolhac [2016a] S. Prolhac, Finite-time fluctuations for the totally asymmetric exclusion process, Phys. Rev. Lett. 116, 090601 (2016a).
  * Prolhac [2016b] S. Prolhac, Extrapolation methods and Bethe ansatz for the asymmetric exclusion process, Journal of Physics A: Mathematical and Theoretical 49, 454002 (2016b).
  * Prolhac [2017] S. Prolhac, Perturbative solution for the spectral gap of the weakly asymmetric exclusion process, Journal of Physics A: Mathematical and Theoretical 50, 315001 (2017).
  * Prolhac [2020] S. Prolhac, Riemann surface for TASEP with periodic boundaries, Journal of Physics A: Mathematical and Theoretical 53, 445003 (2020).
  * Liu  _et al._ [2020] Z. Liu, A. Saenz, and D. Wang, Integral formulas of ASEP and q-TAZRP on a ring, Communications in Mathematical Physics 379, 261 (2020).
  * Prolhac [2021] S. Prolhac, From the riemann surface of TASEP to ASEP, Journal of Physics A: Mathematical and Theoretical 54, 395002 (2021).
  * Baik and Liu [2021] J. Baik and Z. Liu, Periodic TASEP with general initial conditions, Probability Theory and Related Fields 179, 1047 (2021).
  * Noh and Kim [1994] J. D. Noh and D. Kim, Interacting domain walls and the five-vertex model, Physical Review E 49, 1943 (1994).
  * Schütz [1997] G. M. Schütz, Exact solution of the master equation for the asymmetric exclusion process, Journal of Statistical Physics 88, 427 (1997).
  * Alimohammadi  _et al._ [1998] M. Alimohammadi, V. Karimipour, and M. Khorrami, Exact solution of a one-parameter family of asymmetric exclusion processes, Physical Review E 57, 6370 (1998).
  * Sasamoto and Wadati [1998] T. Sasamoto and M. Wadati, One-dimensional asymmetric diffusion model without exclusion, Physical Review E 58, 4181 (1998).
  * Derrida and Evans [1999] B. Derrida and M. R. Evans, Bethe ansatz solution for a defect particle in the asymmetric exclusion process, Journal of Physics A: Mathematical and General 32, 4833 (1999).
  * Alcaraz and Bariev [1999] F. C. Alcaraz and R. Z. Bariev, Exact solution of the asymmetric exclusion model with particles of arbitrary size, Physical Review E 60, 79 (1999).
  * Karimipour [1999] V. Karimipour, A remark on integrability of stochastic systems solvable by matrix product ansatz, Europhysics Letters 47, 501 (1999).
  * Roshani and Khorrami [1999] F. Roshani and M. Khorrami, Asymmetric one-dimensional exclusion processes: A two-parameter exactly solvable example, Physical Review E 60, 3393 (1999).
  * Roshani and Khorrami [2001] F. Roshani and M. Khorrami, Solvable multispecies reaction-diffusion processes, Physical Review E 64, 011101 (2001).
  * Alcaraz and Bariev [2002] F. C. Alcaraz and R. Z. Bariev, Exact integrable spin chains and transfer matrices related to models with stochastic dynamics, Physica A: Statistical Mechanics and its Applications 306, 51 (2002).
  * Ferreira and Alcaraz [2002] A. A. Ferreira and F. C. Alcaraz, Anomalous tag diffusion in the asymmetric exclusion model with particles of arbitrary sizes, Physical Review E 65, 052102 (2002).
  * Roshani and Khorrami [2002] F. Roshani and M. Khorrami, Annihilation-diffusion processes: An exactly solvable model, Journal of Mathematical Physics 43, 2627 (2002).
  * Priezzhev [2003a] V. B. Priezzhev, Exact nonstationary probabilities in the asymmetric exclusion process on a ring, Physical Review Letters 91, 050601 (2003a).
  * Povolotsky  _et al._ [2003] A. M. Povolotsky, V. B. Priezzhev, and C.-K. Hu, Finite size behavior of the asymmetric avalanche process, Physica A: Statistical Mechanics and its Applications 321, 280 (2003).
  * Priezzhev [2003b] V. B. Priezzhev, Exact nonstationary probabilities in the asymmetric exclusion process on a ring, Physical Review Letters 91, 050601 (2003b).
  * Povolotsky [2004] A. M. Povolotsky, Bethe ansatz solution of zero-range process with nonuniform stationary state, Physical Review E 69, 061109 (2004).
  * Roshani and Khorrami [2005] F. Roshani and M. Khorrami, A family of discrete-time exactly-solvable exclusion processes on a one-dimensional lattice, The European Physical Journal B - Condensed Matter and Complex Systems 45, 399 (2005).
  * Povolotsky and Mendes [2006] A. M. Povolotsky and J. F. F. Mendes, Bethe ansatz solution of discrete time stochastic processes with fully parallel update, Journal of Statistical Physics 123, 125 (2006).
  * Povolotsky and Priezzhev [2006] A. M. Povolotsky and V. B. Priezzhev, Determinant solution for the totally asymmetric exclusion process with parallel update, Journal of Statistical Mechanics: Theory and Experiment 2006, P07002 (2006).
  * Povolotsky and Priezzhev [2007] A. M. Povolotsky and V. B. Priezzhev, Determinant solution for the totally asymmetric exclusion process with parallel update: II. Ring geometry, Journal of Statistical Mechanics: Theory and Experiment 2007, P08018 (2007).
  * Majumdar  _et al._ [2008] S. N. Majumdar, K. Mallick, and S. Nechaev, Bethe ansatz in the Bernoulli matching model of random sequence alignment, Physical Review E 77, 011110 (2008).
  * Simon [2009] D. Simon, Construction of a coordinate Bethe ansatz for the asymmetric simple exclusion process with open boundaries, Journal of Statistical Mechanics: Theory and Experiment 2009, P07017 (2009).
  * Lazo and Ferreira [2010] M. J. Lazo and A. A. Ferreira, Asymmetric exclusion model with impurities, Physical Review E 81, 050104 (2010).
  * Lazo and Ferreira [2012] M. J. Lazo and A. A. Ferreira, Asymmetric exclusion model with several kinds of impurities, Journal of Statistical Mechanics: Theory and Experiment 2012, P05017 (2012).
  * Lee [2012] E. Lee, The current distribution of the multiparticle hopping asymmetric diffusion model, Journal of Statistical Physics 149, 50 (2012).
  * Borodin  _et al._ [2015] A. Borodin, I. Corwin, L. Petrov, and T. Sasamoto, Spectral theory for interacting particle systems solvable by coordinate Bethe ansatz, Communications in Mathematical Physics 339, 1167 (2015).
  * Derbyshev  _et al._ [2015] A. E. Derbyshev, A. M. Povolotsky, and V. B. Priezzhev, Emergence of jams in the generalized totally asymmetric simple exclusion process, Physical Review E 91, 022125 (2015).
  * Chen  _et al._ [2022] Z. Chen, J. de Gier, I. Hiki, T. Sasamoto, and M. Usui, Limiting current distribution for a two species asymmetric exclusion process, Communications in Mathematical Physics 395, 59 (2022).
  * Saenz  _et al._ [2022] A. Saenz, C. A. Tracy, and H. Widom, Domain walls in the Heisenberg-Ising spin-1/2 chain, in  _Toeplitz Operators and Random Matrices: In Memory of Harold Widom_, Operator Theory: Advances and Applications (Springer International Publishing, Cham, 2022) pp. 9–47.
  * Gier  _et al._ [2023] J. d. Gier, W. Mead, and M. Wheeler, Transition probability and total crossing events in the multi-species asymmetric exclusion process, Journal of Physics A: Mathematical and Theoretical 56, 255204 (2023).
  * Ishiguro  _et al._ [2023] Y. Ishiguro, J. Sato, and K. Nishinari, Asymmetry-induced delocalization transition in the integrable non-Hermitian spin chain, Physical Review Research 5, 033102 (2023).
  * Lobaskin  _et al._ [2023] I. Lobaskin, M. R. Evans, and K. Mallick, Integrability of two-species partially asymmetric exclusion processes, Journal of Physics A: Mathematical and Theoretical 56, 165003 (2023).
  * Dorlas [1993] T. C. Dorlas, Orthogonality and completeness of the Bethe ansatz eigenstates of the nonlinear Schroedinger model, Communications in Mathematical Physics 154, 347 (1993).
  * Langlands and Saint-Aubin [1995] R. P. Langlands and Y. Saint-Aubin, Algebro-geometric aspects of the Bethe equations, in  _Strings and Symmetries_, Vol. 447 (Springer, Berlin, Heidelberg, 1995) pp. 40–53.
  * Brattain  _et al._ [2017] E. Brattain, N. Do, and A. Saenz, The completeness of the Bethe ansatz for the periodic ASEP (2017).
  * Brillinger [1966] D. R. Brillinger, The Analyticity of the Roots of a Polynomial as Functions of the Coefficients, Mathematics Magazine 39, 145 (1966).
  * Timm [2009] C. Timm, Random transition-rate matrices for the master equation, Physical Review E 80, 021140 (2009).
  * Denisov  _et al._ [2019] S. Denisov, T. Laptyeva, W. Tarnowski, D. Chruściński, and K. Życzkowski, Universal spectra of random Lindblad operators, Physical Review Letters 123, 140403 (2019).
  * Tarnowski  _et al._ [2021] W. Tarnowski, I. Yusipov, T. Laptyeva, S. Denisov, D. Chruściński, and K. Życzkowski, Random generators of Markovian evolution: A quantum-classical transition by superdecoherence., Physical review. E 104, 034118 (2021).
  * Nakerst  _et al._ [2023] G. Nakerst, S. Denisov, and M. Haque, Random sparse generators of Markovian evolution and their spectral properties, Phys. Rev. E 108, 014102 (2023).
  * de Gier and Essler [2005] J. de Gier and F. H. L. Essler, Bethe ansatz solution of the asymmetric exclusion process with open boundaries, Physical Review Letters 95, 240601 (2005).
  * de Gier and Essler [2006] J. de Gier and F. H. L. Essler, Exact spectral gaps of the asymmetric exclusion process with open boundaries, Journal of Statistical Mechanics: Theory and Experiment 2006, P12011 (2006).
  * de Gier  _et al._ [2011] J. de Gier, C. Finn, and M. Sorrell, The relaxation rate of the reverse-biased asymmetric exclusion process, Journal of Physics A: Mathematical and Theoretical 44, 405002 (2011).
  * Crampe  _et al._ [2011] N. Crampe, E. Ragoucy, and D. Simon, Matrix coordinate Bethe ansatz: applications to XXZ and ASEP models, Journal of Physics A: Mathematical and Theoretical 44, 405003 (2011).
  * Alba [2015] V. Alba, Eigenstate thermalization hypothesis and integrability in quantum spin chains, Phys. Rev. B 91, 155123 (2015).
  * Gaveau  _et al._ [2006] B. Gaveau, L. S. Schulman, and L. J. Schulman, Imaging geometry through dynamics: the observable representation, Journal of Physics A: Mathematical and General 39, 10307 (2006).