# Lifted TASEP: A Solvable Paradigm for Speeding up Many-Particle Markov Chains

Fabian H L Essler, Werner Krauth

► To cite this version:

Fabian H L Essler, Werner Krauth. Lifted TASEP: A Solvable Paradigm for Speeding up Many-Particle Markov Chains. Physical Review X, 2024, 14 (4), pp.041035. 10.1103/physrevx.14.041035. hal-04991145

HAL Id: hal-04991145

https://hal.science/hal-04991145v1

Submitted on 14 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Lifted TASEP: A Solvable Paradigm for Speeding up Many-Particle Markov Chains

Fabian H. L. Essler<sup>1,2</sup> and Werner Krauth<sup>2,1</sup>

<sup>1</sup>Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory, Oxford OX1 3PU, United Kingdom

<sup>2</sup>Laboratoire de Physique de l'Ecole normale supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université de Paris Cité, 24 rue Lhomond, 75005 Paris, France

<figure><img src="image_2.png" /></figure>

(Received 17 October 2023; revised 2 May 2024; accepted 5 September 2024; published 6 November 2024)

Virtually all Markov-chain Monte Carlo algorithms used for sampling a given distribution are reversible, and they satisfy the detailed-balance condition. For local chains, this leads to a slow, diffusive exploration of sample space. Significant speedups can be achieved through nonreversible algorithms with the given distribution as a targeted steady state. However, nonreversible algorithms for sampling are difficult to set up and to analyze, and exact speedup results for interacting many-particle systems are very rare. Here, we introduce the “lifted” totally asymmetric simple exclusion process (TASEP) as an exactly solvable paradigm for nonreversible many-particle Markov chains. It samples the same hard-sphere distribution as the Metropolis algorithm for symmetrically diffusing hard-core particles on a one-dimensional lattice. We solve the lifted TASEP by an unusual kind of coordinate Bethe ansatz and show that it exhibits polynomial (in particle number) speedups in the relaxation time for the asymptotic approach of the steady state, as well as the nonasymptotic mixing time, compared to both Metropolis and Kardar-Parisi-Zhang-based dynamics. The lifted TASEP is the reduction onto the one-dimensional lattice of the successful hard-sphere event-chain Monte Carlo algorithm, and we discuss that it can likewise be generalized to soft interaction potentials.

DOI: 10.1103/PhysRevX.14.041035

Subject Areas: Computational Physics, Statistical Physics

## I. INTRODUCTION

Markov-chain Monte Carlo (MCMC), the field of stochastic algorithms built on the concept of sampling, has countless applications in physics, chemistry, biology, statistics, and other fields. Although Markov chains and Monte Carlo methods are often used to describe explicitly time-dependent and nonequilibrium phenomena, as in nuclear reactions and elementary-particle collisions, the major objective of MCMC is the sampling of a given probability distribution. In classical physics, this given probability distribution is normally the Boltzmann distribution, and MCMC sampling has provided crucial insights into a range of problems [1,2] including phase transitions in particle systems and spin models, the behavior of polymers, and the crumpling of interfaces. In quantum physics, MCMC has elucidated the behavior of helium at the λ transition [3], has played a key role in our understanding of strongly correlated electron and spin systems [4], and lies at the origin of the vast field of lattice gauge theory [5].

All the above applications can be formulated as the sampling of a given distribution $\pi$ of configurations $x$ in some high-dimensional sample space $\Omega$. The sampling is, in principle, achieved by initializing the system under consideration in an initial configuration $x_0$, thus, an initial distribution $\pi^{\{t=0\}}(x) = \delta(x - x_0)$, and then employing a Markovian stochastic update rule to generate a sequence $x_0, x_1, \dots, x_t, x_{t+1}, \dots$, with $x_t$ distributed with probability $\pi^{\{t\}}(x_t)$. For sufficiently late times $t$, the set $\{x_t\}$ obtained by repeating the procedure $t$ times practically samples the target distribution, which is nothing but $\pi \approx \pi^{\{t\}}$. Ever since the inception of the field in 1953 [6], the chosen stochastic update rule has mostly been reversible (that is, invariant under time reversal). This implies that the transition matrix, the probability $P(x_t, x_{t+1})$ to move from $x_t$ to $x_{t+1}$, satisfies the celebrated detailed-balance condition

$$ \pi_{x_t} P(x_t, x_{t+1}) = \pi_{x_{t+1}} P(x_{t+1}, x_t). \quad (1) $$

Together with mild conditions (irreducibility and aperiodicity), detailed balance is a sufficient condition for convergence of the sequence $\pi^{\{t=0\}}, \pi^{\{t=1\}}, \dots$ toward $\pi$, whereas a necessary condition is global balance, as discussed in the following. The workhorses of MCMC such as the Metropolis [6] and heat-bath [7–9] algorithms are all reversible; that is, they satisfy the detailed-balance condition of Eq. (1). As a

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
consequence of detailed balance, all net flows vanish at late times. For local Markov chains, which move particles only in their immediate neighborhoods, the resulting dynamics is typically slow and diffusive, which is highly undesirable.

In dimensions greater than one, the relaxation time (the inverse gap of the transition matrix that governs the asymptotic approach to the steady state as well as the steady-state dynamics) and the mixing time (that governs the initial, nonasymptotic approach of the steady state from a most unfavorable initial configuration at $t=0$) can become astronomical (see Appendix A4 for rigorous definitions). For the two-dimensional hard-disk model, for example, no practically converging Monte Carlo algorithm was known for decades, until the advent of nonreversible Markov chains [10,11]. This is all the more problematic as nonlocal cluster methods, which have revolutionized Monte Carlo simulations of spin models [12,13], have not been as successful in particle models at high densities [2,14,15], and, consequently, “local” methods have remained at the forefront of research. In recent years, massively parallel MCMC algorithms have been implemented, often on graphic processing units [16–19], but these algorithms do not reduce the overall operation count of their single-core implementations. The development of fast single-core algorithms (that reduce this count and that can themselves often be parallelized), thus, retains all its interest.

### A. Speedups through “lifting”

It is well known that Markov chains designed to describe nonequilibrium processes do not satisfy detailed balance and are, thus, nonreversible. Their nonvanishing flows, even in the steady state, allow them to relax faster than diffusively. However, in standard applications, their steady states are not a priori known, and characterizing them is, in fact, a key objective [20,21]. The fast relaxation of nonreversible Markov chains can, in principle, be harnessed to speed up the sampling from a target distribution $\pi$ by breaking the detailed balance of Eq. (1) and satisfying only the global-balance condition

$$ \pi_x = \sum_{y \in \Omega} \pi_y P(y, x) \quad \forall x \in \Omega. \qquad (2) $$

Because of $\pi^{\{t+1\}} = \pi^{\{t\}} P$, for any Markov chain, Eq. (2) imposes that the target distribution $\pi$ is the steady state of the Markov-chain dynamics. Under the conditions of irreversibility and aperiodicity, any initial distribution $\pi^{\{t=0\}}$ again converges toward $\pi$ (see Appendix A for more details on the connection between reversible and nonreversible Markov chains). We then want to construct a nonreversible transition matrix $P$, in other words, a Monte Carlo algorithm, corresponding to an out-of-equilibrium process which samples the same $\pi$, but at much earlier times. In real-life sampling problems such as interacting many-particle systems, this program could not be put into practice for many decades, and

local reversible algorithms continued to be the norm. The situation changed only with the advent of “lifted” Markov chains [22,23], that allow one to systematically set up nonreversible algorithms with the required characteristics for sampling a target distribution $\pi$ [24,25].

To introduce to the concept of lifting, we follow Ref. [23] and consider a single particle conducting a random walk on a one-dimensional lattice of $L$ sites with open (that is, hard-wall) boundary conditions. The sample space $\Omega^{RW} = \{1, ..., L\}$ is the set of possible configurations of the particle, and $\mathcal{L}^{RW} = \{-1, +1\}$ denotes the set of moves (see Appendix A3 for precise definitions). At each time $t=0, 1, ..., $, a random move $\sigma$ (forward or backward) is applied to a configuration $x_t$, resulting in a configuration $x_{t+1}$ at the subsequent time step:

$$ \underbrace{\begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix}}_{x_t \in \Omega^{RW}} \rightarrow \underbrace{\begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix}}_{x_t \in \Omega^{RW}, \sigma \in \mathcal{L}^{RW}} \rightarrow \underbrace{\left\{ \begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix} \quad \begin{matrix} p = \frac{1}{2} \\ p = \frac{1}{2} \end{matrix} \right\}}_{x_{t+1} \in \Omega^{RW}} . \qquad (3) $$

Moves across the boundaries at the end points are rejected:

$$ \underbrace{\begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix}}_{x_t=1 \text{ or } x_t=L} \rightarrow \underbrace{\begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix}}_{\text{rejected back move}} \rightarrow \underbrace{\left\{ \begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix} \quad \begin{matrix} p = \frac{1}{2} \\ p = \frac{1}{2} \end{matrix} \right\}}_{x_{t+1} \in \Omega^{RW}}. \qquad (4) $$

The random walk is reversible with a flat steady-state distribution $\pi^{RW}(x) = 1/L$ for all $x$. It takes a time $\sim L^2$ for the particle to move through the lattice and explore the entire sample space, signaling diffusive behavior.

The lifted random walk duplicates (pictorially, “lifts”) each $x_t \in \Omega^{RW}$ into an enlarged sample space $\Omega^{L-RW} = \Omega^{RW} \times \mathcal{L}^{RW}$:

$$ \underbrace{\begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix}}_{x_t \in \Omega^{RW}} \rightarrow \underbrace{\left\{ \begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix} \quad \begin{matrix} (i, \sigma) \\ \sigma = \frac{\alpha}{1-\alpha} \end{matrix} \right\}}_{x_t \in \Omega^{L-RW}}, \qquad (5) $$

which now contains $2L$ elements. In the bulk of the system, the lifted random walk first moves from $(i, \sigma)$ to $(i+\sigma, \sigma)$ and then, with probability $\alpha = 1 - \bar{\alpha}$, flips the direction from $\sigma$ to $-\sigma$:

$$ \underbrace{\begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix}}_{x_t \in \Omega^{L-RW}} \xrightarrow{p=1} \underbrace{\begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix}}_{i+\sigma} \rightarrow \underbrace{\left\{ \begin{bmatrix} \bullet & \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet & \bullet \end{bmatrix} \quad \begin{matrix} p = \alpha \\ p = \bar{\alpha} \end{matrix} \right\}}_{x_{t+1} \in \Omega^{L-RW}, \alpha \ll 1, \bar{\alpha} = 1-\alpha}, \qquad (6) $$
while at a boundary it mostly flips $\sigma$

$$ \underbrace{\bullet\bullet\bullet}_{x_t \text{ at end point}} \xrightarrow{p=1} \underbrace{\bullet\bullet\bullet}_{\sigma\to-\sigma} \to \underbrace{\left\{ \begin{array}{ll} \bullet\bullet\bullet & p = \alpha \\ \bullet\bullet\bullet & p = \bar{\alpha} \end{array} \right.}_{x_{t+1} \in \Omega^{L-RW}} \quad (7) $$

The lifted random walk is nonreversible [the move from $(i, \sigma)$ to $(i + \sigma, \sigma)$ cannot be reversed], but its steady state still corresponds to a flat distribution (i.e., the equal-probability mixture of all $2L$ lifted configurations). If we choose a small probability $\alpha = 1/L$, the particle moves in $\sim L$ steps through the lifted lattice of $2L$ sites, much faster than the original random walk. This speedup of the sampling in the augmented space $\Omega^{L-RW}$ is fully reflected in the exploration of the original $\Omega^{RW}$. Finally, while the random walk has mixing and relaxation times $\sim N^2$, the lifted random walk, in spite of its larger sample space, has much reduced times $\sim N = \sqrt{N^2}$ [23]. The square-root reduction from $\sim N^2$ to $\sim N$ is the optimal speedup that can be achieved for a general lifting [22]. For a precise definition of lifted Markov chains, see Appendix A 3.

## B. Interacting systems: Event-chain Monte Carlo and other algorithms

The lifting approach always involves, as in Eq. (5), that a configuration of the original sample space $\Omega$ corresponds to a number of lifted configurations. The resulting enlarged space is often constructed as a product of $\Omega$ with a set of lifting variables obtained from the move set $\mathcal{L}$. This setting has provided a blueprint for constructing fast, nonreversible Markov chains for real-life sampling problems. In the field of statistics, for example, the “zigzag” algorithm generalizes the single-particle diffusion and has allowed the faster estimation of Bayesian probabilities [26]. In the context of interacting particle models, the lifted random walk of Eqs. (6) and (7) generalizes to the event-chain Monte Carlo algorithm, both for hard spheres and related systems [24,27,28] and for particle and spin models with general potentials [25,29]. This has resulted in vastly improved algorithms and has lead to breakthrough applications in statistical mechanics [10,30] and in chemical physics [31,32]. The event-chain Monte Carlo algorithm has proven particularly powerful for hard-sphere systems, with a reported three-orders-of-magnitude speedup compared to its parent Metropolis algorithm for large system sizes [11]. These simulations remain very time-consuming, and it has not yet been clarified whether the large speedup corresponds to an improved scaling of relaxation times with system size. The model we discuss in the following, the lifted totally asymmetric simple exclusion process (TASEP), is the reduction of hard-sphere event-chain Monte Carlo algorithm to a one-dimensional lattice, and

we demonstrate such an improved scaling for it, with respect to its reversible parent MCMC algorithm.

As the name suggests, the lifted TASEP is a lifting of the celebrated symmetric simple exclusion process (SSEP) [33–35] and its asymmetric generalization. The latter are paradigmatic generalizations of the random walk to interacting many-particle systems and have been applied to the description of a wide range of phenomena [21], including the dynamics of ribosomes along RNA [36,37] and traffic flow [38].

## C. Overview of results and plan

In the present work, we start from the hard-sphere Metropolis algorithm on a one-dimensional lattice, which is equivalent to the SSEP [33–35]. We interpret the much-studied TASEP as a lifting of the SSEP. More precisely, what we call the “bidirectional” TASEP is a lifting of the SSEP as defined in Ref. [22], with the TASEP being “one-half” of the lifting. The focus of this work, the lifted TASEP, is, thus, a second-generation lifting of the SSEP. By construction, the SSEP, the TASEP, and the lifted TASEP all have the same steady state. Numerical simulations show that the lifted TASEP features radically reduced time correlations. We then demonstrate that it can be solved exactly by a coordinate Bethe ansatz. We derive the Bethe equations characterizing the eigenvalues of the transition matrix and show that they fully describe its spectrum. The lifted TASEP has a free parameter, the “pullback” $\alpha$, which does not influence its steady state. For a critical value of the pullback that we obtain analytically, numerical simulations and Bethe ansatz both indicate that the relaxation time of the algorithm scales with an exceptionally small exponent with the system size. This exponent is smaller than for other known models of one-dimensional particle systems and, in particular, smaller than the exponent of the Kardar-Parisi-Zhang (KPZ) class. We, thus, suggest that the lifted TASEP belongs to a new class of Bethe-ansatz solvable models that, as we have argued above, may also explain the exceptionally fast convergence of the event-chain dynamics in higher-dimensional models.

The plan of this work is as follows: In Sec. II, we define the SSEP and its liftings, the TASEP and the lifted TASEP, and establish basic properties such as the steady state that they share, namely, the flat Boltzmann distribution typical of hard-sphere models. We discuss the pullback, a dynamic parameter of the lifted TASEP, which does not affect the steady state but affects how quickly the process relaxes. We also provide a representation of the lifted TASEP as a non-Hermitian quantum spin chain. In Sec. III, we analytically compute the critical pullback, for which the process exhibits a polynomial speedup (in particle number) in its relaxation to the steady state. We also show, using Monte Carlo simulations, the influence of the pullback on the autocorrelation time of the structure factor and estimate
the mixing time of the lifted TASEP through exact enumeration and through large-scale simulations where we follow the evolution of compact initial states toward equilibrium. This gives a crucial upper bound for the relaxation time of our system. In Sec. IV, we derive the Bethe equations and then, in Sec. V, derive bounds on the spectral gap of the lifted TASEP, both for the critical pullback $\alpha_{crit}$ and for generic values of $\alpha$. The bounds combine the results of our numerical solutions with the Bethe-ansatz solution. Section VI contains our conclusions. In a series of appendixes, we present mathematical background on Markov chains, lifting, and the convergence in the presence of nonreversibility (Appendix A) and showcase transition matrices for small systems for the SSEP, the TASEP, and the lifted TASEP (Appendix B). We furthermore present a step-by-step derivation of the Bethe equations for three particles (Appendix C).

## II. EXCLUSION MODELS: NONREVERSIBILITY AND LIFTING

The SSEP implements the local reversible Metropolis algorithm for $N$ indistinguishable hard spheres on an $L$-site lattice with periodic boundary conditions (Sec. II A). We interpret the much-studied TASEP [21,39–44] as a nonreversible lifting of the SSEP (Sec. II B) and then introduce a “second-generation” lifting of the SSEP, the lifted TASEP (Sec. II C). It generalizes the lifted random walk of a single particle (discussed in Sec. I) to the case of $N$ interacting particles and achieves considerable speedups, as we show in Secs. III and V.

### A. Symmetric simple exclusion process

In the SSEP, configurations are described by $N$ occupied sites $\{\dots, j, k, l\dots\}$ on a one-dimensional $L$-site lattice. At each time step, a randomly chosen active particle

$$ \underbrace{\left( \begin{smallmatrix} j & k & l \\ \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet \end{smallmatrix} \right)}_{x_t \in \Omega^{\text{SSEP}}} \rightarrow \left\{ \begin{array}{ll} \underbrace{\left( \begin{smallmatrix} j & k & l \\ \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet \end{smallmatrix} \right)}_{\text{choice of active particle}} & p = \frac{1}{N} \\ \underbrace{\left( \begin{smallmatrix} j & k & l \\ \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet \end{smallmatrix} \right)}_{\text{choice of active particle}} & p = \frac{1}{N} \\ \underbrace{\left( \begin{smallmatrix} j & k & l \\ \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet \end{smallmatrix} \right)}_{\text{choice of active particle}} & p = \frac{1}{N} \end{array} \right. \rightarrow \dots, \quad (8) $$

attempts a random forward or backward move:

$$ \dots \rightarrow \underbrace{\left( \begin{smallmatrix} j & k & l \\ \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet \end{smallmatrix} \right)}_{\text{choice of move}} \rightarrow \left\{ \begin{array}{ll} \underbrace{\left( \begin{smallmatrix} j & k & l \\ \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet \end{smallmatrix} \right)}_{x_{t+1} \in \Omega^{\text{SSEP}}} & p = \frac{1}{2} \\ \underbrace{\left( \begin{smallmatrix} j & k & l \\ \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet \end{smallmatrix} \right)}_{x_{t+1} \in \Omega^{\text{SSEP}}} & p = \frac{1}{2} \end{array} \right. . \quad (9) $$

Moves that would violate the exclusion condition (at most one particle per site) are rejected, resulting in $x_{t+1} = x_t$. The SSEP is reversible, and its steady state is a flat distribution of allowed configurations, $\pi_x = \text{const} \ \forall \ x \in \Omega^{\text{SSEP}}$ (see Appendix B). Its relaxation time scales as $t_{\text{rel}}^{\text{SSEP}} \sim N^3$.

In contrast, the initial approach toward equilibrium from a worst-case initial configuration takes place on a non-asymptotic [45] mixing timescale $t_{\text{mix}} \sim N^3 \log N$. The initial approach to equilibrium, thus, takes longer than the decorrelation in equilibrium, giving rise to the “cutoff phenomenon” [34,35,45], a sudden convergence of $\pi^{(t)}$ toward $\pi$ at $t \sim t_{\text{mix}}$.

Given the simplicity of the steady state of the SSEP (and its generalizations), some comments are in order. The flat distribution $\pi$, with its equal probabilities for all configurations of $N$ nonoverlapping particles on $L$ sites, can clearly be sampled directly as a random combination in $O(N)$ steps, much more rapidly than through any Markov chain. The related problem of $N$ hard spheres on a continuous line can likewise be sampled in $O(N)$ single steps (see Sec. VI. 1 in Ref. [2]). Hence, for this particular state, there is no practical need to employ a nonreversible lifted Markov chain to address the general computational sampling problem we are concerned with. The motivation for considering of the SSEP is as follows: (i) As noted above, the SSEP serves as a model for the sampling of the much-studied hard-sphere distributions in more than one dimensions, where the sampling algorithms are highly nontrivial because of the regions with “illegal” zero-weight configurations in between legal configurations with equal (that is, “flat”) weights. Crucially, efficient direct-sampling or nonlocal cluster algorithms are unavailable at the relatively high densities relevant for the study of phase transitions. (ii) Even though the SSEP has a very simple steady state, its dynamics is highly complex, as it describes an interacting many-particle system. Hence, studying it can be expected to provide key insights on how liftings speed up the dynamics even though the steady state is simple.

### B. Totally asymmetric simple exclusion process

A lifting of the SSEP, with an augmented sample space $\Omega^{\text{TASEP}} = \Omega^{\text{SSEP}} \times \{-1, +1\}$ and a reduced move set which now consists in only the choice of the active particle but not its direction of motion, is provided by the “bidirectional TASEP.” A bidirectional-TASEP configuration, thus, consists in the general direction of motion (either forward or backward) in addition to the particle configuration of the SSEP. At each time step, the direction of the moves is as indicated in the lifted configurations, but the particle to move is chosen randomly. The direction of the moves changes under certain conditions (see Appendix B 3 for an explicit construction of the transition matrix). For periodic boundary conditions, the bidirectional TASEP separates into forward and backward moving sectors, and we then can restrict our attention without loss of generality to the former. In the resulting TASEP, at each time step $t = 0, 1, \dots$, a randomly chosen active particle
$$E\psi^r_{\{\dots,j,\vec{k},l,\dots\}} = \alpha\psi^r_{\{\dots j,k+1,l,\dots\}} + \bar{\alpha}\psi^r_{\{\dots j,k+1,l,\dots\}} \quad (15)$$

and for $k+1=l$:

$$E\psi^r_{\{\dots j, \vec{k}, l, \dots\}} = \alpha\psi^r_{\{\dots j, \vec{k}, k+1, \dots\}} + \bar{\alpha}\psi^r_{\{\dots j, k, k+1, \dots\}} \quad (16)$$

Clearly, these are solved for $E=1$ by $\psi^r = \text{const}$, so that $P^{\text{l-T}}$ is a valid stochastic matrix with unit row sums. The transition matrix has nonzero diagonal elements [see Eq. (14)], so that it is aperiodic for any pullback $0 < \alpha < 1$. The nonzero column elements of the transition matrix are as follows.

For $j < k - 1$ and $l > k + 1$,

$$\begin{aligned} P^{\text{l-T}}(\{\dots j, \vec{k-1}, l, \dots\}, \{\dots j, \vec{k}, l, \dots\}) &= \bar{\alpha}, \\ P^{\text{l-T}}(\{\dots j, k, \vec{l-1}, \dots\}, \{\dots j, \vec{k}, l, \dots\}) &= \alpha. \end{aligned} \quad (17)$$

For $j = k - 1$ and $l > k + 1$,

$$\begin{aligned} P^{\text{l-T}}(\{\dots \vec{k-1}, k, l, \dots\}, \{\dots j, \vec{k}, l, \dots\}) &= \bar{\alpha}, \\ P^{\text{l-T}}(\{\dots k-1, k, \vec{l-1}, \dots\}, \{\dots j, \vec{k}, l, \dots\}) &= \alpha. \end{aligned} \quad (18)$$

For $j < k - 1$ and $l = k + 1$,

$$\begin{aligned} P^{\text{l-T}}(\{\dots j, \vec{k-1}, k+1, \dots\}, \{\dots j, \vec{k}, l, \dots\}) &= \bar{\alpha}, \\ P^{\text{l-T}}(\{\dots j, \vec{k}, k+1, \dots\}, \{\dots j, \vec{k}, l, \dots\}) &= \alpha. \end{aligned} \quad (19)$$

For $j = k - 1$ and $l = k + 1$,

$$\begin{aligned} P^{\text{l-T}}(\{\dots \vec{j}, k, l, \dots\}, \{\dots j, \vec{k}, l, \dots\}) &= \bar{\alpha}, \\ P^{\text{l-T}}(\{\dots j, \vec{k}, l, \dots\}, \{\dots j, \vec{k}, l, \dots\}) &= \alpha. \end{aligned} \quad (20)$$

The two nonzero elements of each column of the transition matrix thus add up to one, so that $P^{\text{l-T}}$ is doubly stochastic. This implies that the steady state is constant for $0 < \alpha < 1$ and shows that the model is indeed a lifting of the TASEP. Furthermore, from any of these cases, one may directly write down eigenvalue equations $E\psi = \psi P^{\text{l-T}}$ for the left eigenvectors $\psi$ of $P^{\text{l-T}}$, which we further analyze using the Bethe ansatz. Again, $\psi = \pi = \text{const}$ is the left eigenvector with eigenvalue 1, as already follows from the double stochasticity.

The lifted TASEP differs from the integrable defect particle in the asymmetric exclusion process studied in Ref. [20] in two important ways. First, in our stochastic process, only the defect particle (that is, the active particle) moves. Second, the lifted TASEP dynamics is nonlocal in the space of the lifting variables, because the pullback acts between subsequent particles (without moving them), which may be far separated in space.

# 1. Lifted TASEP as a non-Hermitian quantum spin chain

We may consider the transition matrix to induce a continuous time evolution, which then can be cast in the form of an imaginary-time Schrödinger equation [51]. To do so, we define three basis states $|a\rangle_j$ with $a = 0, 1, 2$ on each site of the lattice and a basis of operators acting on them:

$$E_j^{ab} = |a\rangle_{jj}\langle b|, \quad a, b \in \{0, 1, 2\}. \quad (21)$$

The probability distributions of the stochastic process give rise to states by the map

$$|\psi\rangle = \sum_{j_1<j_2<\dots<j_N} \sum_{a=1}^{N} \psi_a(j) E_{j_1}^{\sigma_1 0} E_{j_2}^{\sigma_2 0} \dots E_{j_N}^{\sigma_N 0} |0\rangle. \quad (22)$$

The master equation takes the form of an imaginary-time Schrödinger equation

$$\frac{d}{dt} |P(t)\rangle = H |P(t)\rangle, \quad (23)$$

where the “Hamiltonian” is given by

$$H = \sum_{j=1}^{N} [\bar{\alpha} + \alpha P_{12}] [E_j^{02} E_{j+1}^{20} + E_j^{12} E_{j+1}^{21}]. \quad (24)$$

Here, $P_{12}$ is a nonlocal operator that permutes the active particle with the regular particle preceding it:

$$\begin{aligned} P_{12} &= \sum_{k<j} E_k^{21} \Pi_{k+1}^{j-1} E_j^{12} + \sum_{j<k} \Pi_1^{j-1} E_j^{12} E_k^{21} \Pi_{k+1}^L, \\ \Pi_m^n &= \prod_{\ell=m}^{n} E_\ell^{00}. \end{aligned} \quad (25)$$

# 2. Generalizations of the lifted TASEP

The lifted TASEP can be taken to the continuum-space limit of $N$ particles on $L$ sites, with $L \to \infty$ [52]. The model is then conveniently represented as a Markov process in continuous time with the active particle moving with fixed constant positive velocity $v = 1$, rendering it equivalent to the hard-sphere event-chain Monte Carlo algorithm in one dimension. One may show that the critical pullback $\alpha_{\text{crit}}$ exactly compensates the system pressure, connecting to Ref. [49] and to the generalizations presented there to arbitrary potentials between particles.

On the lattice, the lifted TASEP furthermore generalizes to repulsive long-range interactions between subsequent particles:

$$\pi_{\dots j,k,l,\dots} = \dots \pi_{k-j} \pi_{l-k,\dots} \quad (26)$$
where the single-pair term $\pi_0 = 0$ is the vanishing weight of having two particles on the same site, $\pi_1$ the weight for having two particles at distance 1, and $\pi_2 \ge \pi_1$ for a distance of two, etc. In the generalized lifted TASEP, the deterministic first part of the lifted-TASEP move becomes itself stochastic:

$$ \underbrace{\left( \begin{array}{c|c|c} \bullet & \bullet & \bullet \\ \hline \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet \end{array} \right)}_{x_t \in \Omega^{L-T}} \rightarrow \left\{ \begin{array}{ll} \begin{array}{c} \bullet \\ \hline \bullet \\ \bullet \end{array} & \begin{array}{c} \bullet \\ \hline \bullet \\ \bullet \end{array} \\ \begin{array}{c} \bullet \\ \hline \bullet \\ \bullet \end{array} & \begin{array}{c} \bullet \\ \hline \bullet \\ \bullet \end{array} \end{array} \right. \begin{array}{l} p = \frac{\pi_{l-k+1}}{\pi_{l-k}} \dots, \\ 1-p \end{array} \quad (27) $$

so that the attempted move of particle k can be rejected by particle l already at a distance [as in the lower case in Eq. (27)] and transformed into a lifting move. The move in Eq. (27), which implements the factorized Metropolis algorithm [25], is followed by a pullback move, and the steady state is again the (nonflat) Boltzmann distribution in Eq. (26) [52].

### III. LIFTED TASEP: CRITICAL PULLBACK AND NUMERICAL COMPUTATIONS

We now turn to the role of the pullback $\alpha$. Vanishing pullback ($\alpha = 0$) makes the activity advance with positive velocity, while unit pullback ($\alpha = 1$) makes it recede on average, that is, acquire a negative velocity. We observe, using numerical simulations of large systems, that the lifted TASEP equilibrates much faster for a particular value $\alpha_{crit} = N/L$. This reproduces behavior previously observed for the event-chain Monte Carlo, in the continuous-space version of the lifted TASEP [49]. We then analyze the spectrum of the model for small system sizes (Sec. III B) and show that it is aperiodic and irreducible for $0 < \alpha < 1$ but not for $\alpha = 0, 1$. For all $0 < \alpha < 1$, the lifted TASEP thus converges toward the steady state without the need of restarts, which were required in earlier versions of event-chain Monte Carlo [49,53,54].

#### A. Critical pullback and MCMC simulations

While the value of the pullback $0 < \alpha < 1$ does not affect the steady state that is always flat, it plays a remarkable role in the dynamics of the lifted TASEP. Any lifted configuration can evolve either through the $\alpha$ move or the $\bar{\alpha}$ move [see Eqs. (13) and (14)]. An $\bar{\alpha}$ move advances the activity by one site. A series of $L-1$ sequential $\bar{\alpha}$ moves translates the entire configuration backward by one site. An $\alpha$ move, on average over all lifted configurations, pulls the activity back by $(L-N)/N$ sites. Averaging over $\alpha$ moves and $\bar{\alpha}$ moves shows that the activity drifts with velocity $-\alpha L/N + 1$, so that the drift vanishes for $\alpha_{crit} = N/L$. In event-chain Monte Carlo, the drift velocity is proportional to the system pressure [25], and in the one-dimensional harmonic model, the zero-pressure point

exhibits fast integrated autocorrelation times [55]. This behavior was numerically confirmed in a continuous-space variant of the lifted TASEP [49] (see also Ref. [27]). The pullback similarly generates a zero-pressure point at $\alpha_{crit}$, with the connection between the pressure and the pullback most visible in the continuum-space limit mentioned in Sec. II C 2.

Steady-state MCMC simulations for the lifted TASEP confirm the special role of the pullback. These simulations start from an equilibrium configuration $x_0$ that we obtain by direct sampling of a random combination and that gives rise to an equilibrium chain of configurations $\{x_0, x_1, x_2, \dots\}$ whose correlations depend only on the difference of times. For individual trajectories $x_t = \{j_1(t), \dots, j_N(t)\}$, where $j_n$ denote the particle positions, we analyze the integrated autocorrelation time $\tau_{S(q)}$ [56] of the equilibrium structure factor at time $t$ (see Appendix A 4 for the definition of $\tau_{S(q)}$):

$$ S(q, t) = \frac{1}{N} \left| \sum_{r=1}^{L} e^{iqr} \rho(r, t) \right|^2, \quad \rho(r, t) = \sum_{n=1}^{N} \delta_{r,j_n(t)}. \quad (28) $$

Here, $q = 2\pi/L$ is the smallest nonzero wave number. The structure factor is sensitive to long-range density fluctuations, which are expected to relax slowly in equilibrium. As a function of $\alpha$, $\tau_{S(q)}$ exhibits a pronounced dip at $\alpha = \alpha_{crit}$ [see Fig. 1, and see Appendix A 4 or Eq. (3.4) in Ref. [56] for a definition of the integrated autocorrelation function]. This means that the density autocorrelations in equilibrium decay much more quickly at $\alpha = \alpha_{crit}$.

We now turn to the $N$ dependence of $\tau_{S(q)}$ at fixed $\alpha$. We observe a power-law behavior with different exponents at $\alpha_{crit}$ and away from it (see also Ref. [49]). In Fig. 2, we show results for $L = 2N$ and several difference values of $\alpha$.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Integrated autocorrelation time $\tau_{S(q)}$ for the structure factor $S(q = 2\pi/L)$ [see Eq. (28)], as a function of the pullback $\alpha$, for different $N = L/2$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Integrated autocorrelation time $\tau_{S(q)}$ of the lifted TASEP for the structure factor $S(q = 2\pi/L)$, as a function of $N = L/2$, for different values of the pullback. Inset: $\tau_{S(q)}$ for the SSEP (scaling as $\sim N^3$), the TASEP ($\sim N^{5/2}$), and the lifted TASEP ($\sim N^{3/2}$), evidencing significant speedup through lifting.</figcaption></figure>

The behavior for large values of $N$ can be summarized as follows:

$$ \tau_{S(q)} \sim \begin{cases} N^{3/2} & \text{for } \alpha = \alpha_{\text{crit}}, \\ N^{5/2} & \text{for } \alpha \neq \alpha_{\text{crit}}. \end{cases} \qquad (29) $$

The speedup of the lifted TASEP at $\alpha = \alpha_{\text{crit}}$ with respect to the SSEP and the TASEP is impressive (see the inset in Fig. 2). At $\alpha_{\text{crit}}$, the motion of the activity is driftless but superdiffusive (see also Ref. [49]): It takes only $\sim L^{3/2}$ steps for the activity to move through the system of $L$ sites. All the results that we quote here for the integrated autocorrelation time of the structure factor for $\alpha = \alpha_{\text{crit}}$ are also found for its autocorrelation function, as well as for the autocorrelation of the active-particle position.

To estimate the mixing time $t_{\text{mix}}$ for $\alpha = \alpha_{\text{crit}} = \frac{1}{2}$, we take the compact state $x_{t=0} = \{\vec{1}, 2, \dots, N\}$ as our trial worst-case initial configuration and compute in a numerically exact way the total variation distance from the steady state $\pi$ (see Appendix A 4). At a difference with the simulations for the structure factor, these simulations, thus, have an explicit time dependence, and the compact initial configuration, at time $t = 0$, is the conjectured most-out-of-equilibrium possible configuration (see Chap. 4 in Ref. [57] for a discussion of mixing times). Our data (shown in Fig. 3) are compatible with

$$ t_{\text{mix}} \sim N^2. \qquad (30) $$

<figure><img src="image_5.png" /><figcaption>FIG. 3. Worst-case total variation distance $d_{\text{TV}}$ for the lifted TASEP as a function of rescaled time $t$ for small system sizes with $L = 2N$, at pullback $\alpha = \alpha_{\text{crit}} = \frac{1}{2}$. Inset: diameter of the lifted TASEP, compatible with a scaling $d = N^2 + O(N)$.</figcaption></figure>

This indicates that the lifted TASEP mixes much faster than the TASEP, for which both the relaxation time [40,46] and the mixing time [47] scale as $N^{5/2}$ in our time units. Given that the mixing time cannot be faster than the relaxation time, Eq. (30) indicates that the relaxation times of the lifted TASEP cannot be slower than $\sim N^2$. As the TASEP is in the KPZ universality class, this already shows that the lifted TASEP is in a different dynamical universality class, with a smaller dynamical exponent. Also, it cannot scale faster than $\sim N^2$ because of a recent diameter bound [58]. We also find that the diameter $d$ (the minimal time to connect any two configuration $x, y \in \Omega^{l-T}$) numerically scales as $\sim N^2$ with a coefficient equal to one, that is, $d = N^2 + O(N)$ (see the inset in Fig. 3). The diameter gives an $\alpha$-independent lower bound on the mixing time, $t_{\text{mix}} \ge d/2 \sim N^2$, which is compatible with Eq. (30) (see Ref. [57]).

## B. Lifted TASEP: Transition-matrix spectrum

To illustrate the spectrum of the transition matrix $P^{l-T}$ and its dependence on the pullback $\alpha$, we show in Fig. 4 the complex eigenvalues for $N = 4, L = 8$ as a function of $\alpha$.

<figure><img src="image_6.png" /><figcaption>FIG. 4. Eigenvalue spectrum of the lifted TASEP. (a) Spectrum for $N = 4, L = 8$ as a function of the pullback $\alpha$. For $\alpha = 0$ and 1, all eigenvalues $E$ are equally spaced on the complex unit circle. (b) Spectrum for $N = 7, L = 14$ for $\alpha = \alpha_{\text{crit}} = \frac{1}{2}$, with a single eigenvalue on the complex unit circle.</figcaption></figure>
For $L = 2N$ and $\alpha = 0$, the $N(2N)$ eigenvalues lie on $2N(2N-1)$ equally spaced points on the unit circle. Each of these eigenvalues has degeneracy

$$C_{N-1} = \frac{1}{N} \binom{2N-2}{N-1},$$

where $C_{N-1}$ is the Catalan number. For $\alpha = 1$, the eigenvalues lie on $2N^2$ equally spaced points on the complex unit circle, and all the eigenvalues except $E = (1,0)$ are again $C(N-1)$ times degenerate, while the eigenvalue $E = (1,0)$ has a degeneracy of $(2N)$ - $(2N^2-1)C(N-1)$ [see Fig. 4(a)]. Clearly, in these two cases, the transition matrix $P^{l-T}$ is not irreducible (except for $N=2, \alpha=0$), and it is periodic for all $N$. In contrast, for $0 < \alpha < 1$, the transition matrix $P^{l-T}$ has only a single eigenvalue with $|E|=1$ and is irreducible and aperiodic [see Fig. 4(b)].

## IV. LIFTED TASEP: EXACT SOLUTION BY A BETHE ANSATZ

In order to solve the eigenvalue equation for the transition matrix in the $N$-particle sector, we introduce the following compact notations:

$$\{j; a\} = \{j_1, \dots, j_{a-1}, \vec{j}_a, j_{a+1}, \dots, j_N\}. \quad (31)$$

Here, $j_1 < \dots < j_N$ are the positions of the $N$ particles, and $a$ specifies the active particle. The transition matrix of the lifted TASEP in Eq. (20) then takes the form

$$P^{l-T}(\{k; b\}, \{j; a\}) = \bar{\alpha}\delta_{a,b+1}\delta_{k_b+1,k_{b+1}}\delta_{j,k} \\ + \bar{\alpha}\delta_{b,a}(1 - \delta_{k_b+1,k_{b+1}})\delta_{j,k+e_b} \\ + \alpha\delta_{b,a+1}(1 - \delta_{k_b+1,k_{b+1}})\delta_{j,k+e_b} \\ + \alpha\delta_{b,a}\delta_{k_b+1,k_{b+1}}\delta_{k,j}. \quad (32)$$

Here, we define $e_a = (0, \dots, 0, 1, 0, \dots)$, where the nonzero entry is at position $a$ and $\delta_{j,k} = \prod_{n=1}^N \delta_{j_n, k_n}$. The left eigenvalue equation for $P^{l-T}$ then takes the form

$$E\psi_a(\mathbf{j}) = \bar{\alpha}\psi_{a-1}(\mathbf{j})\delta_{j_a, j_{a-1}+1} \\ + \bar{\alpha}\psi_a(\dots, j_a-1, \dots)(1-\delta_{j_a, j_{a-1}+1}) \\ + \alpha\psi_{a+1}(\dots, j_{a+1}-1, \dots)(1-\delta_{j_{a+1}, j_a+1}) \\ + \alpha\psi_a(\mathbf{j})\delta_{j_{a+1}, j_a+1}. \quad (33)$$

As already noted above, when viewed as a non-Hermitian quantum spin chain [51], this corresponds to a three-state model with infinite-range interactions (see Sec. II C 1). We solve the eigenvalue equation (33) by applying a nested Bethe ansatz

$$\psi_a(\mathbf{j}) = \sum_{Q \in S_N} A_a(Q) \prod_{j=1}^{N} (z_j)^{j_{Q_j}}, \quad (34)$$

where $S_N$ is the set of all permutations of $(1, \dots, N)$ and $\{z_1, \dots, z_N\}$ are complex spectral parameters that specify a given solution. We note in passing that we have not attempted to use the Bethe ansatz to determine the right eigenvectors of the transition matrix, which, in general, differ from the left eigenvectors but have the same eigenvalues. Substituting the ansatz of Eq. (34) into Eq. (33) for the probabilities and assuming that $|j_n - j_{n-1}| > 1 \ \forall \ n$, we obtain

$$A_a(Q) \left[ E - \frac{\bar{\alpha}}{z_{Q_a^{-1}}} \right] = \frac{\alpha}{z_{Q_{a+1}^{-1}}} A_{a+1}(Q). \quad (35)$$

These relations express all $A_a(Q)$ in terms of only $A_1(Q)$:

$$A_a(Q) = \left[ \prod_{b=1}^{a-1} \frac{z_{Q_{b+1}^{-1}}}{\alpha} \left( E - \frac{\bar{\alpha}}{z_{Q_b^{-1}}} \right) \right] A_1(Q). \quad (36)$$

Additional relations are obtained by considering “sector boundaries.”

### A. Sector boundary I: $j_{a-1}+1 < j_a = j_{a+1}-1$

Here, the eigenvalue equation

$$E\psi_a(\mathbf{j}) = \alpha\psi_a(\mathbf{j}) + \bar{\alpha}\psi_a(\dots, j_a-1, \dots) \quad (37)$$

[see Eq. (19)] implies the following relations between amplitudes:

$$\frac{A_a[C_a(Q)]}{A_a(Q)} = \frac{z_{Q_{a+1}^{-1}}}{z_{Q_a^{-1}}} T(z_{Q_{a+1}^{-1}}, z_{Q_a^{-1}}), \\ T(z_1, z_2) = -\frac{E - \alpha - \bar{\alpha}/z_2}{E - \alpha - \bar{\alpha}/z_1}, \quad (38)$$

where $C_a(Q) = (\dots Q_{a+1}^{-1}, Q_a^{-1}, \dots)^{-1} = (a, a+1)Q$ is the permutation obtained from $Q$ by swapping $a$ with $a+1$. These can be rewritten as conditions on $A_1(Q)$ using Eq. (36):

$$\frac{A_1[C_a(Q)]}{A_1(Q)} = T(z_{Q_{a+1}^{-1}}, z_{Q_a^{-1}}) \begin{cases} 1 & \text{if } a > 1, \\ \frac{z_{Q_{a+1}^{-1}}}{z_{Q_a^{-1}}} & \text{if } a = 1. \end{cases} \quad (39)$$

### B. Sector boundary II: $j_{a-1}+1 = j_a < j_{a+1}-1$

Here, the eigenvalue equation

$$E\psi_a(\mathbf{j}) = \bar{\alpha}\psi_{a-1}(\mathbf{j}) + \alpha\psi_{a+1}(\dots, j_{a+1}-1, \dots) \quad (40)$$
[see Eq. (18)] implies the following relations between amplitudes:

$$ \frac{A_a[C_{a-1}(Q)]}{A_a(Q)} = T(z_{Q_a^{-1}}, z_{Q_{a-1}^{-1}}) \frac{E - \bar{\alpha}/z_{Q_a^{-1}}}{E - \bar{\alpha}/z_{Q_{a-1}^{-1}}}. \quad (41) $$

These can be rewritten as conditions on $A_1(Q)$ using Eq. (36):

$$ \frac{A_1[C_{a-1}(Q)]}{A_1(Q)} = T(z_{Q_a^{-1}}, z_{Q_{a-1}^{-1}}), \quad a \ge 2. \quad (42) $$

As an arbitrary permutation can be generated from the transpositions $(a, a+1)$, the conditions of Eqs. (39) and (42) allow us to relate all amplitudes $A_1(Q)$ to $A_1(1, \dots, N)$. The various relations are equivalent to the following: Let $Q$ and $Q'$ be two permutations related by a single transposition:

$$ Q' = \Pi_{t_1, t_2} Q, \quad (43) $$

where $\Pi$ exchanges the entries $Q_{t_1}$ and $Q_{t_2}$ and we choose conventions such that we always have $Q_{t_1} < Q_{t_2}$. Then we have

$$ \frac{A_1(Q)}{A_1(Q')} = T(z_{t_1}, z_{t_2}) \times \begin{cases} \frac{z_{t_1}}{z_{t_2}} & \text{if } 1 \in \{Q_{t_1}, Q_{t_2}\}, \\ 1 & \text{else.} \end{cases} \quad (44) $$

The first case can be viewed as a scattering process between the active particle and an ordinary particle, while the latter corresponds to a scattering process involving two ordinary particles.

All conditions arising from other sector boundaries (e.g., three neighboring particles) reduce to Eqs. (36) and (44).

### C. Bethe equations

Imposing periodic boundary conditions gives the requirements

$$ \begin{aligned} \psi_N(j_1, \dots, j_{N-1}, L+1) &= \psi_1(1, j_1, \dots, j_{N-1}), \\ \psi_1(0, j_2, \dots, j_N) &= \psi_N(j_2, \dots, j_N, L). \end{aligned} \quad (45) $$

These conditions give rise to the Bethe equations that determine the spectral parameters $z_j$ and the eigenvalue $E$:

$$ \begin{aligned} z_a^{L-1} &= \frac{\bar{\alpha}/z_a - E}{\alpha} \prod_{b=1}^{N} T(z_a, z_b), \quad a = 1, \dots, N, \\ \prod_{j=1}^{N} (E - \bar{\alpha}/z_j) &= \alpha^N \prod_{k=1}^{N} \frac{1}{z_k}. \end{aligned} \quad (46) $$

A step-by-step derivation of Eq. (46) for $N = 3$ is given in Appendix C.

Some comments are in order. As always in integrable models, it is not a priori clear that the solutions of the Bethe equations provide a complete set of eigenstates of the transition matrix. To gain some insight into the question of completeness, we consider specific values of $\alpha$ for $L = 6 = 2N$ and $L = 8 = 2N$ and verify that the solutions of Eq. (46) indeed provide a complete set of eigenvalues $P^{l-T}$. In addition, we check for $L = 2N \le 14$ and several values of $\alpha$ that the eigenvalues $P^{l-T}$ with the largest real parts are all obtained from Eq. (46) to machine precision.

### D. Equivalent form of the Bethe equations for $L=2N$

In the next section, we consider in some detail the case $L = 2N$ (so that $\alpha_{crit} = 1/2$), in which it is convenient to reparametrize the Bethe-ansatz equations using

$$ \begin{aligned} \beta &= \frac{\bar{\alpha}}{E - \alpha}, & u_a &= \frac{2z_a}{\beta} - 1, & \delta &= \frac{\beta - \bar{\alpha}/\alpha}{\beta + \bar{\alpha}/\alpha}, \\ \mu &= \left(\frac{2}{\beta}\right)^L \frac{1 - \alpha(1 - \beta)}{2\alpha} \prod_{b=1}^{N} \frac{u_b - 1}{u_b + 1}. \end{aligned} \quad (47) $$

This maps Eq. (46) onto

$$ \begin{aligned} (1 - u_a^2)^{L/2} &= -\mu(u_a + \delta), & a &= 1, \dots, \frac{L}{2}, \\ \left( \frac{2\alpha}{1 - \alpha(1 - \beta)} \right)^{L/2} &= \prod_{b=1}^{L/2} (u_b + \delta). \end{aligned} \quad (48) $$

These equations are similar to the ones of the TASEP with periodic boundary conditions [39,40,59] in that the equations for $u_a$ can be written as simple polynomial equations involving a self-consistently determined constant. However, there are two key differences.

(1) The factor $u_a + \delta$ in the equation for $u_a$ makes the roots lie on a nontrivial contour, whereas for the TASEP they lie on a circle.

(2) The eigenvalue $E$ enters like a root of a nested Bethe equation rather than being a simple additive function of the roots.

### V. SPECTRAL PROPERTIES OF THE LIFTED TASEP FROM THE BETHE ANSATZ FOR $L=2N$

We now turn to an analysis of the spectrum of the transition matrix $P^{l-T}$ for large system sizes using the Bethe-ansatz equations (48). We first show that the steady state is indeed the fully mixed state at fixed particle number and identify the corresponding solution of Eq. (48). We then present an analysis of the spectral gap

$$ \Delta = \operatorname{Re} \ln[E^*], \quad (49) $$

where $E^*$ is the eigenvalue with magnitude closest to one. Our results are summarized as follows:
$$ \Delta(L) \ge \begin{cases} A(\frac{1}{2})L^{-2} + o(L^{-2}) & \text{if } \alpha = \alpha_{\text{crit}} = \frac{1}{2}, \\ A(\alpha)L^{-5/2} + o(L^{-5/2}) & \text{if } \alpha \neq \alpha_{\text{crit}}. \end{cases} \quad (50) $$

These bounds on the gap translate into bounds for the relaxation times, which scale at least as $L^2$ for $\alpha = \alpha_{\text{crit}}$ and at least as $L^{5/2}$ otherwise. Together with our bounds for the mixing time obtained from our numerical simulations, this, in turn, leads us to conjecture that the bounds in Eq. (50) are, in fact, saturated; i.e., we have an equality rather than an inequality in Eq. (50). At $\alpha = \alpha_{\text{crit}}$, this then implies a polynomial speedup for the relaxation time with respect to the SSEP and the TASEP.

### A. Steady state

By virtue of the double stochasticity of $P^{1-T}$, its steady state $\pi$ is the fully mixed state at fixed particle number and has $E = 1$. The corresponding Bethe equations read

$$ \beta = 1, \quad \delta = 2\alpha - 1. \quad (51) $$

Our equations for $L = 2N$ then become

$$ \frac{(1 - u_a^2)^{L/2}}{u_a - 1 + 2\alpha} = -\frac{2^L}{2\alpha} \prod_{b=1}^{L/2} \frac{u_b - 1}{u_b + 1}, \quad a = 1, \dots, \frac{L}{2}, \\ (2\alpha)^{L/2} = \prod_{b=1}^{L/2} (u_b - 1 + 2\alpha). \quad (52) $$

Clearly, taking $u_a \to 1$ provides a solution to these equations. The limit $u_a \to 1$ corresponds to $z_a \to 1$, and substituting this into our expression in Eq. (34) for $\psi_a(j)$ correctly gives $\pi$.

### B. Spectral gap for $\alpha = \alpha_{\text{crit}}$

For concreteness, we now fix $L = 4m + 2$ for integer $m$. In order to study the spectral gap, we first identify the structure of solutions of the Bethe equations that give rise to low-lying excitations by considering small system sizes $L \le 14$. We then follow such solutions for larger values of $L$. In particular, this leads us to consider the following family of solutions of the Bethe equations (48). Out of the $L + 2$ roots of

$$ (1 - u^2)^{(L+2)/2} + \mu_L(u + \delta_L) = 0, \quad (53) $$

we select a set $\{\bar{u}_1, \dots, \bar{u}_2\}$ that fulfils

$$ \frac{L}{2} \ln[\bar{u}_n - 1] + \frac{L}{2} \ln[-\bar{u}_n - 1] - \ln(-\mu_L) - \ln(\bar{u}_n + \delta_L) = 2\pi i I_n, \quad (54) $$

where

<figure><img src="image_7.png" /><figcaption>FIG. 5. Sets of roots $\{u_j|1 \le j \le L/2\}$ corresponding to the excited state of Eq. (55) for $L = 2N$ and $\alpha = \alpha_{\text{crit}} = \frac{1}{2}$, with $L = 214$ (blue) and $L = 494$ (yellow). For large $L$, the roots approach a nontrivial contour in the complex plane.</figcaption></figure>

$$ \{I_n\} = \left\{ -\frac{L-2}{4}, \dots, -1, 0, 1, \dots, \frac{L-2}{4} \right\}. \quad (55) $$

We find that the root distribution of the corresponding excited state is of the form shown in Fig. 5. We observe that as the system size is increased the roots lie on a closed contour in the complex plane. The logarithm of the eigenvalue of these states is shown as a function of system size in Fig. 6. The numerical results for $18 \le L \le 494$ are well described by fits of the form

$$ \operatorname{Im} \ln[E_1(L)] = \frac{a_0}{L} + \frac{a_1}{L^{3/2}} + \frac{a_2}{L^2}, \\ \operatorname{Re} \ln[E_1(L)] = \frac{a_3}{L^2} + \frac{a_4}{L^{5/2}} + \frac{a_5}{L^3}, \quad (56) $$

where

$$ \begin{align} a_0 &= 6.2847, & a_1 &= -0.04746, & a_2 &= 6.673, \nonumber \\ a_3 &= -20.092, & a_3 &= 6.0446, & a_4 &= -29.342. \tag{57} \end{align} $$

<figure><img src="image_8.png" /><figcaption>FIG. 6. Real and imaginary parts of the logarithm of the eigenvalue of the excited state of Eq. (55) as functions of system size for $L = 2N$ and $\alpha = \alpha_{\text{crit}} = \frac{1}{2}$. The solid lines are the fits Eqs. (56) and (57). The dashed red line is a one-parameter fit $\operatorname{Re} \ln[E_1(L)] = -19.7448L^{-2}$ and is indistinguishable from Eqs. (56) and (57) on the scale of the plot.</figcaption></figure>
Interestingly, the real part of $\ln[E_1(L)]$ scales as $L^{-2}$, while the imaginary part scales as $L^{-1}$. This is in contrast to the gap for the TASEP at half filling, which is real and scales— in our units of time—as $L^{-5/2}$ [40]. On the other hand, the TASEP gap below half filling has an imaginary part that scales (in our units of time) as $L^{-2}$ [41].

We now provide more details on how we select the state Eq. (55). We first consider small system sizes and identify solutions to the Bethe equations that correspond to the lowest few excited states. We then construct solutions of the Bethe equations for the corresponding distribution of integers $I_n$ for a range of system sizes $L$ and compute their eigenvalues. For sufficiently large values of $L$, the state Eq. (55) is the lowest excited state among those we have considered. For short and intermediate system sizes, other states have larger values of $\text{Re}\ln[E]$, but, as we increase $L$, “level crossings” of the real part occur. The same happens for other values of $\alpha$, and, importantly, level crossings can occur at large values of $L$. An example for $\alpha = 0.4$ is shown in Fig. 10. This is unusual for Bethe-ansatz solvable models and constitutes a significant complication. Therefore, $\text{Re}(\ln[E_1(L)])$ is (presently) only an upper bound for the spectral gap. The Bethe-ansatz result for the scaling of the spectral gap at $\alpha = 1/2$ is, thus,

$$ \text{Re}[\ln(E_1(L))] \geq \frac{-20.0924}{L^2} + o(L^{-2}). \quad (58) $$

This, in turn, implies that the relaxation time is larger than or equal to $\sim N^2$. On the other hand, it is a theorem that $t_{\text{mix}} \geq t_{\text{rel}}$ [see Eq. (12.37) in Ref. [57] and Appendix A 4]. Using the results of our numerical simulations in Sec. III A for the mixing time, we conclude that the relaxation time of the lifted TASEP at $\alpha = \alpha_{\text{crit}}$ indeed scales as $N^2$.

### C. Spectral gap for $\alpha \neq \alpha_{\text{crit}}$

For $L = 2N$, we have determined the scaling with system size of the eigenvalues of the solution specified by the same set of integers $\{I_n\}$ in Eq. (55) for several values of $\alpha \neq 1/2$ where the activity drifts with finite velocity. This again provides an upper bound of the spectral gap for large values of $L$. In Fig. 7, we show the root distributions for $\alpha = 0.4$ and $\alpha = 0.9$, respectively. The logarithm of the eigenvalue of these states is shown as a function of system size in Figs. 8 and 9. Here, we fit the spectral gap to the functional form

$$ \begin{aligned} \text{Re}[\ln(E_1(L))] &= \frac{a_0}{L^{5/2}} + \frac{a_1}{L^3} + \frac{a_2}{L^{7/2}} + \mathcal{O}(L^4), \\ \text{Im}[\ln(E_1(L))] &= \frac{b_0}{L} + \frac{b_1}{L^{3/2}} + \frac{b_2}{L^2} + \mathcal{O}(L^{5/2}), \end{aligned} \quad (59) $$

where, for $\alpha = 0.4$, we find

<figure><img src="image_9.png" /><figcaption>FIG. 7. Sets of roots $\{u_j|1 \le j \le L/2\}$ corresponding to the excited state of Eq. (55) for $\alpha = 0.9$, $N = L/2$, and $L = 272$ (yellow) and $\alpha = 0.4$, $N = L/2$, and $L = 290$ (blue).</figcaption></figure>

$$ \begin{aligned} a_0 &= -573.913, & a_1 &= 5651.58, & a_2 &= -18183.4, \\ b_0 &= -7.54557, & b_1 &= 0.167695, & b_2 &= -7.61005. \end{aligned} \quad (60) $$

For $\alpha = 0.9$, the eigenvalue is real for all $N$, and the fit parameters are

$$ a_0 = -19.1347, \quad a_1 = 24.5024, \quad a_2 = -1.94545. \quad (61) $$

As we have already mentioned, the determination of the spectral gap as a function of $L$ is complicated by the presence of level crossings between Bethe ansatz eigenvalues at large values of $L$. An example is shown for $\alpha = 0.4$ in Fig. 10. We see that in that case a Bethe state with a gap that scales as $L^{-3/2}$, in fact, gives a larger eigenvalue than the state in Eq. (55) for $L \le 82$ and only then gets overtaken. As we cannot rule out such level crossings between the state described by Eq. (55) and other Bethe states, our result for the eigenvalue provides only an upper bound of the gap.

Our Bethe-ansatz prediction for the scaling of the spectral gap at $\alpha \neq \alpha_{\text{crit}} = 1/2$ is summarized as

<figure><img src="image_10.png" /><figcaption>FIG. 8. Real part of the logarithm of the eigenvalue of the excited state described in the text for $\alpha = 0.4$ as functions of system size. The solid lines are fits to Eq. (59) with the parameters of Eq. (60).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 9. Logarithm of the eigenvalue of the excited state described in the text for $\alpha = 0.9$ as functions of system size. The solid lines are fits to Eq. (59) with the parameters of Eq. (61).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 10. Level crossing as a function of $L$ between a family of solutions to the Bethe equations with $L^{-3/2}$ scaling (blue) and a family with $L^{-5/2}$ scaling (orange) for $\alpha = 0.4$ and $L = 2N$. The crossing occurs around $L \approx 82$.</figcaption></figure>

$$ \mathrm{Re}[\ln(E_1(L))] \geq \frac{A(\alpha)}{L^{5/2}} + o(L^{-5/2}). \quad (62) $$

## VI. DISCUSSION AND OUTLOOK

Lifting is a powerful method with huge potential for speeding up MCMC algorithms, but exact results were so far restricted to noninteracting systems. In this work, we have introduced an exactly solvable lifting of an interacting many-particle system, the hard-sphere Metropolis algorithm on the one-dimensional lattice, which is equivalent to the SSEP. The latter is a key paradigm for the stochastic dynamics of interacting particle systems. Mixing and relaxation times of the SSEP are rigorously known to scale differently with system size [34,35]. Its $t_{rel} \sim N^3$ scaling can be viewed as a particle diffusion in the Edwards-Wilkinson universality class [60], and it also follows from the proof [61] of the Aldous conjecture which identifies the dynamics of exclusion models with single-particle diffusion [57]. Because liftings are by design local in space, it is known [22] that the maximal possible speedup (i.e., the reduction of mixing and relaxation timescales) scales at

best with an exponent that is half the one of the base model. This corresponds to a diffusive-to-ballistic speedup [50]. The best-case inverse-gap scaling for liftings of the SSEP is $\sim N^{3/2}$ (up to multiplicative logarithms). Our MCMC simulations of the lifted TASEP are indeed compatible with this optimal scaling behavior for pullback $\alpha = \alpha_{crit}$: Both the integrated autocorrelation time $\tau_{S(q)}$ of the structure factor and the motion of the activity appear to scale as $\sim N^{3/2}$, whereas for the SSEP, the integrated autocorrelation time scales as $\sim N^3$. Clearly, in our model, lifting achieves a spectacular speedup. For $\alpha \neq \alpha_{crit}$, the same analysis suggests $\tau_{S(q)} \sim N^{5/2}$. In contrast, our simulations of the mixing time for a class of compact initial configurations suggests a scaling of $t_{mix} \sim N^2$ for $\alpha = \alpha_{crit}$.

As we have shown, the lifted TASEP is exactly solvable by a nested Bethe ansatz in spite of its nonlocal structure in the space of the particle labels. Our numerical solution of the Bethe equations derived in Sec. IV yields exact eigenvalues of its transition matrix for system sizes which are totally out of reach of any other method. By analyzing specific families of solutions with $L = 2N$, we find from the Bethe ansatz that the relaxation time (that is, the inverse of the spectral gap) cannot scale faster than $\sim N^2$ at $\alpha = \alpha_{crit}$ and $\sim N^{5/2}$ for $\alpha \neq \alpha_{crit}$. On the other hand, the relaxation time cannot scale slower than the mixing time. At $\alpha = \alpha_{crit}$, we have excellent numerical data indicating that $t_{mix} \sim N^2$, which then fixes the relaxation time. The calculation of the mixing time should, of course, be made rigorous, as has been done for the SSEP and the TASEP.

An $N^2$ scaling of the relaxation time corresponds to a dynamical critical exponent $z = 1$, which is smaller than the exponents for both the KPZ ($z = 3/2$) and the Edwards-Wilkinson ($z = 2$) universality classes, and indicates a significant polynomial (in particle number) speedup in the approach to the steady state. The $\sim N^{5/2}$ behavior for $\alpha \neq \alpha_{crit}$ agrees with our simulations of $\tau_{S(q)}$, whereas at $\alpha = \alpha_{crit}$ our simulations suggest $\sim N^{3/2}$. While different scalings of integrated autocorrelation times and relaxation times are known in the literature [62], the bound obtained from the Bethe ansatz, in fact, rules out a $\sim N^{3/2}$ scaling of the relaxation time. A possible explanation of the simulation result is that for large $N$ the structure factor $S(q)$ has negligible overlap with the mode that causes the $\sim N^2$ scaling of the inverse gap seen in the Bethe-ansatz analysis at $\alpha = \alpha_{crit}$. We, thus, once again witness the exceptional difficulty of estimating the timescales $t_{mix}$ and $t_{rel}$ in MCMC, even in simplified models, and the utility of techniques such as the Bethe ansatz, which directly access the spectrum of the transition matrix. We find it of particular interest that the scaling of the relaxation time of the lifted TASEP can be improved by tuning the external parameter $\alpha$ of the nonreversible lifting of the SSEP. This illustrates a theorem presented in Ref. [22], which says that a sizable speedup of a reversible Markov chain can be
achieved only through a nonreversible lifting. This, in turn, implies that the reversible local Markov chains (at least within the lifting framework) are all more or less equivalent and that the nonreversibility at the heart of this paper is an essential ingredient.

Our work raises a number of questions that should be addressed in future research. First and foremost, a more complete analysis of the spectrum of the transition matrix in the limit of large system sizes needs to be carried out, and the scaling behavior of the low-lying eigenvalues should be derived analytically by recasting the Bethe equations in terms of as nonlinear integral equations; cf. Refs. [41,63–65]. Arguably the most important open question from the point of view of integrability is to determine the integrable structure of the lifted TASEP (its R matrix and L operator [66]). Given the spatially nonlocal nature of its transition matrix, it is not obvious how the lifted TASEP fits into the existing classification.

Knowing the integrable structure would be very helpful in addressing other pressing questions, like whether the partially asymmetric simple exclusion process has an integrable lifting as well and which boundary conditions are compatible with integrability. This question is natural from the point of view of Monte Carlo algorithms and also arises in the context of boundary-induced phase transitions in the TASEP. In particular, one could ask whether there is a lifting of the open TASEP that reproduces its highly nontrivial steady state [67] but has a significantly shorter relaxation time. We also expect that the lifted TASEP with several active particles remains integrable and expect the continuum limit of the lifted TASEP at vanishing density [49], that is, the event-chain algorithm with factor fields, to be integrable. We need to understand the role of the pullback $\alpha$ and the underlying reason why a vanishing drift of the activity, at $\alpha = \alpha_{crit}$, provides the optimal speedup of the relaxation of the lifted TASEP and its generalizations toward equilibrium. For the lifted TASEP itself, it appears necessary and possible to render our findings rigorous, especially to obtain mathematically rigorous proofs of the relaxation and especially the mixing times as they have been obtained for the SSEP [34,35] and the TASEP [47]. We also require better understanding of which observables actually show the $\sim N^2$ correlations at $\alpha = \alpha_{crit}$.

Recent decades have witnessed the convergence of rigorous mathematical research on Markov chains and of the application-driven development of Monte Carlo algorithms. The present work takes a significant step forward in this direction. On the one hand, the lifted TASEP invites a mathematically rigorous analysis following in the footsteps of recent work on the TASEP and the SSEP. On the other hand, this new model is precisely the one-dimensional lattice reduction of the event-chain Monte Carlo algorithm, which aims at bringing the spectacular speedups discussed in this work to real-life applications in statistical and

chemical physics. In this context, the Bethe ansatz has now entered the game as a powerful tool. We expect that the lifted TASEP is only the first member of a larger class of integrable lifted Markov chains and will be an important benchmark for nonreversible Monte Carlo algorithms relevant for applications in the sciences.

## ACKNOWLEDGMENTS

This work was supported in part by the Engineering and Physical Sciences Research Council under Grant No. EP/ S020527/1 (F. H. L. E.) and the Alexander-von-Humboldt Foundation (W. K.). W. K. gratefully acknowledges the hospitality of the Simons Center for Computational Physical Chemistry at New York University, where this research was carried out in part, and was supported by a grant from the Simons Foundation (Grant No. 839534, MET). We are grateful to A. C. Maggs, C. Monthus, G. Robichon, M. Staudacher, and X. Zhang for helpful discussions and D. Mugita for a careful reading of the manuscript.

## APPENDIX A: MARKOV CHAINS, LIFTINGS, AND CONVERGENCE

In this appendix, we collect some basic definitions and results about Markov chains, lifted Markov chains, and Markov processes and discuss the convergence toward the steady state.

### 1. Basic definition

We consider finite sample spaces $\Omega$ and Markov chains defined by (i) a transition matrix $P$, where $P(x, y)$ gives the conditional probability to move in one time step to $y$ if at $x$; (ii) an initial distribution $\pi^{\{t=0\}}$, generally concentrated on a fixed initial configuration $x_0$. The probability distribution at time $t = 1, 2, \dots$, starting from $x_0$, is denoted by $\pi^{\{t\}} = P^t(x_0, \cdot)$. It satisfies $\pi_x^{\{t\}} = \sum_y \pi_y^{\{t-1\}} P(y, x)$. An irreducible Markov chain has a unique steady state that satisfies

$$ \pi_x = \sum_{y \in \Omega} \pi_y P(y, x) \quad \forall x \in \Omega. \quad (A1) $$

This equation becomes a condition on $\pi$, the “global balance” condition, in our case where the target distribution $\pi$ is the imposed steady state.

The steady state $\pi$ is a row vector, and we refer to it as a left eigenvector of $P$ with eigenvalue 1 (that is, $E\pi = \pi P$, with $E = 1$). We call $\pi$ an equal-probability mixture if all its elements are identical. For an irreducible Markov chain that is aperiodic (that is, has no periodicity in its return to a given configuration $x$), $\pi$ is the limit of $\pi^{\{t\}}$ for $t \to \infty$. In the main text, we mention right eigenvectors, that is, column vectors $\psi^r$ that satisfy $\psi^r = P\psi^r$, and, in particular,
the column vector (1, ..., 1), which indicates that the sum of the possibilities to move away from a configuration x (including back to x) must sum up to one. This means that the matrix is stochastic. The eigenvalue equation for left eigenvectors reads

$$E\psi_x = \sum_{y \in \Omega} \psi_y P(y, x) \quad \forall x \in \Omega. \qquad (A2)$$

For an irreducible and aperiodic transition matrix, the eigenvalue $E = 1$ is the only one with unit norm. The detailed-balance condition

$$\pi_x P(x, y) = \pi_y P(y, x) \quad \forall x \in \Omega, y \in \Omega \qquad (A3)$$

implies the global-balance condition, as follows from summing over y and using the conservation of probabilities $\sum_y P(x, y) = 1$. As discussed in the main text, detailed balance is equivalent to reversibility, and nonreversible and reversible Markov chains both satisfy the global-balance condition Eq. (A1).

## 2. Markov chains vs Markov processes

In the main text, we consider discrete-time Markov chains, which move single particles at times $t = 0, 1, 2, ...$. This models the workings of a serial computer. The SSEP and the TASEP are usually defined as continuous-time Markov processes, where each of the elements of the respective move sets is sampled on average once per unit time. Thus, in the continuous-time SSEP, an average of $2N$ moves are attempted per unit time interval ($N$ attempted moves for the continuous-time TASEP). The discrete-time and the continuous-time versions are bijections of each other, with a rescaling of timescales by $2N$ for the SSEP and $N$ for the TASEP and, in addition, the replacement of equidistant discrete time steps by Poisson-distributed ones. As a consequence, the $\sim N^{5/2}$ relaxation and mixing timescales of the discrete-time TASEP become $\sim N^{3/2}$ timescales for continuous times.

## 3. Lifted Markov chains

Lifting [22,23] connects a Markov chain with sample space $\Omega$, transition matrix $P$, and steady state $\pi$ with a lifted Markov chain $\hat{\Pi}$ with $\hat{\Omega}$, $\hat{P}$, $\hat{\pi}$ through a mapping $f$ from $\hat{\Omega}$ to $\Omega$. Under this transformation, $\pi$ must be preserved:

$$\pi_v = \sum_{i \in f^{-1}(v)} \hat{\pi}_i, \qquad (A4)$$

where $f^{-1}(v)$ denotes the set of all $i \in \hat{\Omega}$ such that $f(i) = v$. In addition, the flows must be unchanged:

$$\underbrace{\pi_v P(v, u)}_{\text{collapsed flow}} = \sum_{i \in f^{-1}(v), j \in f^{-1}(u)} \overbrace{\hat{\pi}_i \hat{P}(i, j)}^{\text{lifted flow}}.$$

In this work, the lifted sample space $\hat{\Omega}$ is simply the Cartesian product of the “collapsed” sample space $\Omega$ with a space $\tilde{\mathcal{L}}$ of lifting variables, which correspond to the move set $\mathcal{L}$ or parts of it ($\tilde{\mathcal{L}} \subset \mathcal{L}$). While the lifted sample space is thus enhanced, the set of moves is smaller, and the lifted transition matrix is sparser. Although the transition matrix and its lifted version are closely related, the mixing and relaxation times may be decreased (at most) to their square roots [22,50].

Furthermore, in our case all the steady-state probabilities are constant, which implies that

$$P(v, u) = \frac{|\hat{\Omega}|}{|\Omega|} \sum_{i \in f^{-1}(v), j \in f^{-1}(u)} \hat{P}(i, j), \qquad (A5)$$

where $|\hat{\Omega}|$ denotes the size of the sample space $\hat{\Omega}$.

## 4. Convergence

For an irreducible and aperiodic Markov chain, $\pi^{(t)}$ converges toward $\pi$ for $t \to \infty$ from any initial configuration $x_0$. This convergence is commonly quantified by the worst-case distance $d_{\text{TV}}(t)$:

$$d_{\text{TV}}(t) = \max_{x_0} \|P^t(x_0, \cdot) - \pi\|_{\text{TV}}, \qquad (A6)$$

where $\|\mu - \nu\|_{\text{TV}}$ denotes the total variation distance (TVD) between two distributions $\mu$ and $\nu$:

$$\|\mu - \nu\|_{\text{TV}} = \frac{1}{2} \sum_{i \in \Omega} |\mu_i - \nu_i|. \qquad (A7)$$

The mixing time $t_{\text{mix}}(\epsilon)$ (which is known for the SSEP and the TASEP and which we determine numerically for the lifted TASEP) is defined as the earliest time at which $d_{\text{TV}}(t) = \epsilon$. It satisfies $d_{\text{TV}}[\ell t_{\text{mix}}(\epsilon)] \le (2\epsilon)^\ell$ for $\ell = 1, 2, ...$. This yields a bound for all times larger than the mixing time:

$$d_{\text{TV}}(t) \le \exp(-t/t_{\text{mix}}) \quad (\text{for } t > t_{\text{mix}}), \qquad (A8)$$

where we define $t_{\text{mix}} = t_{\text{mix}}(1/(2\epsilon))$. The mixing time $t_{\text{mix}}$ is a nonasymptotic timescale with a strict inequality for finite times given by Eq. (A8). Importantly, the mixing time can be larger than the inverse gap (that is, the relaxation time). In the SSEP, for example, $t_{\text{mix}} \sim N^3 \log N$, whereas the relaxation time (the inverse gap) scales as $t_{\text{rel}} \sim N^3$. This insight is at the basis of the modern theory of Markov chains [45,57].

However, for asymptotically large times, the worst-case distance approaches zero as

$$\lim_{t \to \infty} d_{\text{TV}}(t)^{1/t} = E^*, \qquad (A9)$$

where $E^*$ is the largest absolute value of transition matrix eigenvalues different from $E_1 = 1$ [see Eq. (12.37) in
Ref. [57] ]. Defining the spectral gap by $\gamma^* = -\ln(E^*)$ and the relaxation time as the inverse gap $t_{rel} = 1/\gamma^*$, one then has an asymptotic expression

$$d_{\text{TV}} \sim \exp(-\gamma^* t) \quad (\text{for } t \to \infty). \quad (\text{A10})$$

In the SSEP, for example, the mixing time scales more slowly than the relaxation time. This leads to a sudden collapse of $d_{\text{TV}}$ (on a timescale $t_{\text{mix}}$), as it switches over from a slow decay on the mixing timescale [see Eq. (A8)] to a faster decay on the relaxation timescale $t_{\text{rel}} = 1/\gamma^*$. We note that Eqs. (A8) and (A10) establish together that the relaxation time cannot be larger than the mixing time.

The relaxation time remains of relevance in equilibrium (rather than only for the asymptotic approach to equilibrium). For finite reversible Markov chains, one has

$$\text{Var}_{\pi} P^t g \le \text{Var}_{\pi} g \exp(-E_1 t), \quad (\text{A11})$$

where $g$ is an arbitrary function on the sample space.

In the main text, we discuss integrated autocorrelation times and autocorrelation functions. Following Ref. [56], we define the unnormalized autocorrelation function of an observable $f$ as the equilibrium average:

$$C(t) = \langle f_s f_{s+t} \rangle - \mu^2, \quad (\text{A12})$$

where $\mu$ is the sample mean value $\langle f \rangle$ of $f$. The normalized autocorrelation function is then $\rho(t) = C(t)/C(0)$, and the integrated autocorrelation time is given by

$$\tau_f = \frac{1}{2} \sum_{t=-\infty}^{\infty} \rho(t). \quad (\text{A13})$$

In our case, the function $f$ is the structure factor, with $f_t = S(q, t)$ defined in Eq. (28), and the given, smallest possible, nonzero wave number $q$.

## APPENDIX B: EXAMPLE TRANSITION MATRICES

For concreteness, we provide examples for transition matrices for the standard and the lifted random walks (Appendix B 1), as well as the SSEP with periodic and hard-wall boundary conditions (Appendix B 2). We further present example transition matrices for liftings of the SSEP, namely, the TASEP with both boundary conditions (Appendix B 3) and the lifted TASEP (Appendix B 4).

### 1. Random walk and lifted random walk

Here, we give the transition matrices of the random walk for $L = 4$ sites, in the basis:

$$1 = \Box \Box \Box \quad 2 = \Box \Box \Box \quad 3 = \Box \Box \Box \quad 4 = \Box \Box \Box \Box.$$

With hard-wall boundary conditions, the case discussed in the main text, we have

$$P_{\text{walls}}^{\text{RW}} = \frac{1}{2} \begin{bmatrix} 1 & 1 & \cdot & \cdot \\ 1 & \cdot & 1 & \cdot \\ \cdot & 1 & \cdot & 1 \\ \cdot & \cdot & 1 & 1 \end{bmatrix}. \quad (\text{B1})$$

Here and in the following, the “.” stands for zero. This transition matrix is irreducible and aperiodic by virtue of the presence of diagonal terms. It is doubly stochastic (each of its rows and of its columns sum to one), so that the steady state is constant. As $P_{\text{walls}}^{\text{RW}}$ is symmetric and doubly stochastic, the associated process is reversible.

In the basis

$$\begin{align*} 1 &= \Box \Box \Box \quad 3 &= \Box \Box \Box \quad 5 &= \Box \Box \Box \quad 7 &= \Box \Box \Box \Box \\ 2 &= \Box \Box \Box \quad 4 &= \Box \Box \Box \quad 6 &= \Box \Box \Box \quad 8 &= \Box \Box \Box \Box \end{align*} \quad (\text{B2})$$

the transition matrix of the lifted random walk with hard-wall boundary conditions reads

$$P_{\text{walls}}^{\text{RW}} = \begin{bmatrix} \cdot & \cdot & \bar{\alpha} & \alpha & \cdot & \cdot & \cdot \\ \bar{\alpha} & \alpha & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \bar{\alpha} & \alpha & \cdot & \cdot \\ \alpha & \bar{\alpha} & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \bar{\alpha} & \alpha & \cdot \\ \cdot & \cdot & \alpha & \bar{\alpha} & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \alpha & \bar{\alpha} \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \alpha \end{bmatrix},$$

where $\bar{\alpha} = 1 - \alpha$. It is again doubly stochastic, so that the steady state is constant. As required by Eq. (A5), its $2 \times 2$ blocks sum up to the corresponding single entries in $P_{\text{walls}}^{\text{RW}}$.

### 2. SSEP example transition matrices

For the SSEP with two particles on four sites, we use the basis

$$\begin{align*} 1 &= \Box \Box \Box \quad 2 &= \Box \Box \Box \quad 3 &= \Box \Box \Box \Box \tag{B3} \\ 4 &= \Box \Box \Box \quad 5 &= \Box \Box \Box \Box \quad 6 &= \Box \Box \Box \Box \end{align*}$$

and with periodic boundary conditions, its transition matrix is

$$P_{\text{pbc}}^{\text{SSEP}} = \frac{1}{4} \begin{bmatrix} 2 & 1 & \cdot & \cdot & 1 & \cdot \\ 1 & \cdot & 1 & 1 & \cdot & 1 \\ \cdot & 1 & 2 & \cdot & 1 & \cdot \\ \cdot & 1 & \cdot & 2 & 1 & \cdot \\ 1 & \cdot & 1 & 1 & \cdot & 1 \\ \cdot & 1 & \cdot & \cdot & 1 & 2 \end{bmatrix}. \quad (\text{B4})$$
For hard-wall boundary conditions, one has instead

$$P_{\text{walls}}^{\text{SSEP}} = \frac{1}{4} \begin{bmatrix} 3 & 1 & \cdots & \cdots & \cdots \\ 1 & 1 & 1 & 1 & \cdots \\ \vdots & \vdots & 1 & 2 & 1 \\ \vdots & \vdots & \vdots & 1 & 1 \\ \vdots & \vdots & \vdots & \vdots & 1 & 3 \end{bmatrix}. \quad (B5)$$

Both these transition matrices are doubly stochastic, so that their steady state is the equal-probability mixture. As the transition matrices are symmetric, these transition matrices are reversible.

### 3. Example transition matrices for the TASEP

We first discuss the bidirectional TASEP for two particles on four sites in the following basis of $2(\frac{4}{2}) = 12$ configurations, twice as many as for the SSEP:

$$ \begin{array}{l} 1 = \overline{\bullet \bullet \bullet \bullet} \\ 2 = \overline{\bullet \bullet \bullet \bullet} \\ 3 = \overline{\bullet \bullet \bullet \bullet} \\ \vdots \\ 11 = \overline{\bullet \bullet \bullet \bullet} \end{array} \qquad 4 = \overline{\bullet \bullet \bullet \bullet} \qquad (B6) $$

The bidirectional TASEP, which is a bona fide lifting of the SSEP, as defined in Ref. [22], rejects moves that violate the exclusion condition [see the first row in Eq. (10)]. It changes the overall direction upon collision with a wall (concretely, attempting to move the right particle in configuration “2” leads to 2, while moving the left particle leads to “1”). This yields the following transition matrix:

$$P_{\text{walls}}^{\text{TASEP}} = \frac{1}{2} \begin{bmatrix} 1 & 1 & \cdots & \cdots & \cdots \\ 1 & 1 & \cdots & \cdots & \cdots \\ \vdots & \vdots & 1 & 1 & \cdots \\ \vdots & \vdots & \vdots & 1 & 1 \\ \vdots & \vdots & \vdots & \vdots & 1 & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & 1 & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & 1 & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & 1 & 1 \end{bmatrix}. \quad (B7)$$

The transition matrix $P_{\text{walls}}^{\text{TASEP}}$ is doubly stochastic and concomitantly the steady state is flat (all configurations having the same weight), but it is nonsymmetric and, therefore, clearly nonreversible. The $2 \times 2$ blocks in Eq. (B7) again sum up to the corresponding entries of $P_{\text{walls}}^{\text{SSEP}}$ [see Eq. (B5)], and the bidirectional TASEP is, thus, a lifting of the SSEP.

With periodic boundary conditions, the transition matrix of the bidirectional TASEP can again be written in the basis of Eq. (B6), with appropriate flipping probabilities between the forward and backward sectors. The periodic boundary conditions allow one to separate the two sectors and to restrict one’s attention to the forward sector (see Ref. [50] for more details). For two particles on four sites, in the “forward” basis equivalent to that of Eq. (B3), this yields

$$P_{\text{pbc}}^{\text{TASEP}} = \frac{1}{2} \begin{bmatrix} 1 & 1 & \cdots & \cdots & \cdots \\ \vdots & \vdots & 1 & 1 & \cdots \\ \vdots & \vdots & \vdots & 1 & 1 \\ \vdots & \vdots & \vdots & \vdots & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & 1 & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & 1 & 1 \end{bmatrix}. \quad (B8)$$

This transition matrix is irreducible, and it is aperiodic because of its nonzero diagonal elements. As mentioned, it is an incomplete lifting of the SSEP, because of the absence of backward moves, but this could be easily amended. Nevertheless, as a doubly stochastic matrix, its steady state is flat (all configurations having the same weight).

### 4. Lifted TASEP transition matrix (example)

For simplicity, we consider the transition matrix of the lifted TASEP with periodic boundary conditions only in the forward sector. This leaves us with $N(\frac{L}{N})$ configurations. For two particles on four sites, we have

$$ \begin{align*} 1 &= \overline{\bullet \bullet \bullet \bullet} & 2 &= \overline{\bullet \bullet \bullet \bullet} & 3 &= \overline{\bullet \bullet \bullet \bullet} & 4 &= \overline{\bullet \bullet \bullet \bullet} \\ 5 &= \overline{\bullet \bullet \bullet \bullet} & 6 &= \overline{\bullet \bullet \bullet \bullet} & 7 &= \overline{\bullet \bullet \bullet \bullet} & 8 &= \overline{\bullet \bullet \bullet \bullet} \\ 9 &= \overline{\bullet \bullet \bullet \bullet} & 10 &= \overline{\bullet \bullet \bullet \bullet} & 11 &= \overline{\bullet \bullet \bullet \bullet} & 12 &= \overline{\bullet \bullet \bullet \bullet} \end{align*} $$

In this basis, the transition matrix of the lifted TASEP follows from Eqs. (17) and (18) to be
$$ P^{\text{l-T}} = \begin{bmatrix} \alpha & \bar{\alpha} & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdot & \alpha & \bar{\alpha} & \cdots & \cdots & \cdots & \cdots \\ \cdot & \cdot & \alpha & \bar{\alpha} & \cdots & \cdots & \cdots \\ \cdot & \cdot & \cdot & \alpha & \bar{\alpha} & \cdots & \cdots \\ \cdot & \cdot & \cdot & \cdot & \alpha & \bar{\alpha} & \cdots \\ \cdot & \cdot & \cdot & \cdot & \cdot & \alpha & \bar{\alpha} \\ \bar{\alpha} & \alpha & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bar{\alpha} \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bar{\alpha} \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \bar{\alpha} \end{bmatrix}, \quad (B9) $$

where $\bar{\alpha} = \alpha - 1$. $P^{\text{l-T}}$ is doubly stochastic [so that its steady state is flat (all configurations having the same weight)], irreducible, and aperiodic because of its nonzero diagonal elements. The transition matrix $P^{\text{l-T}}$ can be partitioned into 2 × 2 blocks that sum up to the analogous elements of $P_{\text{pbc}}^{\text{TASEP}}$, as required for a lifting for a constant $\pi$.

## APPENDIX C: THREE-PARTICLE BETHE ANSATZ

For concreteness, we present a step-by-step derivation of the Bethe-ansatz equations of the lifted TASEP for $N=3$ particles on an L-site lattice with periodic boundary conditions. The left eigenvalue equations for the transition matrix, in the case where all three particles are well separated ($j < k-1, k < l-1$), take the form

$$ \begin{aligned} E\psi_{\{\vec{j},k,l\}} &= \bar{\alpha}\psi_{\{\vec{j}-1,k,l\}} + \alpha\psi_{\{j,k-1,l\}}, \\ E\psi_{\{j,\vec{k},l\}} &= \bar{\alpha}\psi_{\{\vec{j}-1,k,l\}} + \alpha\psi_{\{j,k,l-1\}}, \\ E\psi_{\{j,k,\vec{l}\}} &= \bar{\alpha}\psi_{\{j,k,l-1\}} + \alpha\psi_{\{j-1,k,l\}}. \end{aligned} \quad (C1) $$

The Bethe ansatz for $\psi$ then reads

$$ \begin{aligned} \psi_{\{\vec{j},k,l\}} &= A_{\bullet\bullet\circ} z_1^j z_2^k z_3^l + B_{\bullet\bullet\circ} z_1^j z_2^l z_1^k + \cdots + F_{\bullet\bullet\circ} z_1^l z_2^k z_3^j, \\ \psi_{\{j,\vec{k},l\}} &= A_{\circ\bullet\circ} z_1^j z_2^k z_3^l + B_{\circ\bullet\circ} z_1^j z_2^l z_1^k + \cdots + F_{\circ\bullet\circ} z_1^l z_2^k z_3^j, \\ \psi_{\{j,k,\vec{l}\}} &= A_{\circ\circ\vec{l}} z_1^j z_2^k z_3^l + B_{\circ\circ\vec{l}} z_1^j z_2^l z_1^k + \cdots + F_{\circ\circ\vec{l}} z_1^l z_2^k z_3^j. \end{aligned} \quad (C2) $$

The three sectors $\{\bullet\circ\circ\}$, $\{\circ\bullet\circ\}$, and $\{\circ\circ\bullet\}$ correspond to the position of the active particle, while the six coefficients $A, B, \dots, F$ correspond to the permutations of the indices $i, j, k$. Inserting Eq. (C2) into Eq. (C1) and comparing coefficients yields

$$ \begin{aligned} \mathcal{E}_1 A_{\bullet\bullet\circ} &= \frac{\alpha}{z_2} A_{\circ\circ\circ}, & \mathcal{E}_1 B_{\bullet\bullet\circ} &= \frac{\alpha}{z_3} B_{\circ\circ\circ}, \dots, \\ \mathcal{E}_2 A_{\circ\circ\bullet} &= \frac{\alpha}{z_3} A_{\circ\circ\circ}, & \mathcal{E}_2 B_{\circ\circ\bullet} &= \frac{\alpha}{z_2} B_{\circ\circ\circ}, \dots, \\ \mathcal{E}_3 A_{\circ\circ\circ} &= \frac{\alpha}{z_1} A_{\bullet\bullet\circ}, & \mathcal{E}_2 B_{\circ\circ\circ} &= \frac{\alpha}{z_1} B_{\bullet\bullet\circ}, \dots, \end{aligned} \quad (C3) $$

where we define

$$ \mathcal{E}_a = E - \frac{\bar{\alpha}}{z_a}. \quad (C4) $$

Multiplying these equations for each of the $A, B, \dots, F$ coefficients yields

$$ \prod_{a=1}^{3} \left( E - \frac{\bar{\alpha}}{z_a} \right) = \frac{\alpha^3}{z_1 z_2 z_3}, \quad (C5) $$

which is identical for $N=3$ to the second expression in Eq. (46).

Next, we consider the eigenvalue equations in the case when two particles are on neighboring sites and the third particle is detached:

$$ \begin{aligned} E\psi_{\{\vec{k}-1, k, l\}} &= \bar{\alpha}\psi_{\{\vec{k}-2, k, l\}} + \alpha\psi_{\{\vec{k}-1, k, l\}}, \\ E\psi_{\{k-1, \vec{k}, l\}} &= \bar{\alpha}\psi_{\{\vec{k}-1, k, l\}} + \alpha\psi_{\{k-1, k, l-1\}}, \\ E\psi_{\{j, l-1, l\}} &= \bar{\alpha}\psi_{\{j, l-2, l\}} + \alpha\psi_{\{j, l-1, l\}}, \\ E\psi_{\{j, l-1, \vec{l}\}} &= \bar{\alpha}\psi_{\{j, l-1, l\}} + \alpha\psi_{\{j-1, l-1, l\}}. \end{aligned} \quad (C6) $$

Substituting Eq. (C2) into Eq. (C6) and comparing coefficients gives

$$ \begin{aligned} \frac{A_{\bullet\bullet\circ}}{C_{\bullet\bullet\circ}} &= S_{12}, & \frac{B_{\bullet\bullet\circ}}{D_{\bullet\bullet\circ}} &= S_{13}, & \frac{E_{\bullet\bullet\circ}}{F_{\bullet\bullet\circ}} &= S_{23}, \\ \frac{A_{\circ\circ\bullet}}{B_{\circ\circ\bullet}} &= S_{23}, & \frac{C_{\circ\circ\bullet}}{E_{\circ\circ\bullet}} &= S_{13}, & \frac{F_{\circ\circ\bullet}}{D_{\circ\circ\bullet}} &= S_{21}, \end{aligned} \quad (C7) $$

where

$$ S_{kl} = - \frac{z_k}{z_l} \frac{E - \alpha - \bar{\alpha}/z_l}{E - \alpha - \bar{\alpha}/z_k}. \quad (C8) $$

Proceeding analogously with the second and third lines in Eq. (C6) yields

$$ \begin{aligned} E \left( \frac{A_{\circ\circ\vec{l}}}{z_2} + \frac{B_{\circ\circ\vec{l}}}{z_3} \right) &= \bar{\alpha} \left( \frac{A_{\circ\circ\vec{l}}}{z_2} + \frac{B_{\circ\circ\vec{l}}}{z_1} \right) + \alpha \left( \frac{A_{\circ\circ\vec{l}}}{z_1 z_2} + \frac{B_{\circ\circ\vec{l}}}{z_3 z_1} \right), \\ E \left( \frac{A_{\circ\circ\vec{l}}}{z_1} + \frac{C_{\circ\circ\vec{l}}}{z_2} \right) &= \bar{\alpha} \left( \frac{A_{\circ\circ\vec{l}}}{z_1} + \frac{C_{\circ\circ\vec{l}}}{z_2} \right) + \alpha \left( \frac{A_{\circ\circ\vec{l}}}{z_3 z_1} + \frac{C_{\circ\circ\vec{l}}}{z_3 z_2} \right), \end{aligned} $$
and five additional equations that start with $B, ..., F$. The solution to these conditions is

$$ 
\begin{aligned}
B_{\cdot\cdot\cdot} &= \frac{z_2}{z_3} S_{32} A_{\cdot\cdot\cdot}, \\
C_{\cdot\cdot\cdot} &= S_{21} A_{\cdot\cdot\cdot}, \\
D_{\cdot\cdot\cdot} &= \frac{z_2}{z_3} S_{31} S_{32} A_{\cdot\cdot\cdot}, \\
E_{\cdot\cdot\cdot} &= \frac{z_1}{z_3} S_{21} S_{31} A_{\cdot\cdot\cdot}, \\
F_{\cdot\cdot\cdot} &= \frac{z_1}{z_3} S_{21} S_{31} S_{32} A_{\cdot\cdot\cdot},
\end{aligned}
\quad (C9) $$

which fixes all amplitudes $B, ..., F$ in terms of $A$.

Finally, the sector with three adjacent occupied sites yields the following eigenvalue equations:

$$ E\psi_{\{\vec{j},j+1,j+2\}} = \bar{\alpha}\psi_{\{j-1,j+1,j+2\}} + \alpha\psi_{\{\vec{j},j+1,j+2\}}, \quad (C10) $$

$$ E\psi_{\{j,j+1,j+2\}} = \bar{\alpha}\psi_{\{\vec{j},j+1,j+2\}} + \alpha\psi_{\{j,j+1,j+2\}}, \quad (C11) $$

$$ E\psi_{\{j,j+1,j+2\}} = \bar{\alpha}\psi_{\{j,j+1,j+2\}} + \alpha\psi_{\{j-1,j+1,j+2\}}, \quad (C12) $$

Of these equations, only Eq. (C11) leads to a new relation:

$$ 
\begin{aligned}
& \{[E - \alpha]A_{\cdot\cdot\cdot} - \bar{\alpha}A_{\cdot\cdot\cdot}\}_{\frac{z_3}{z_1}} + \{[E - \alpha]B_{\cdot\cdot\cdot} - \bar{\alpha}B_{\cdot\cdot\cdot}\}_{\frac{z_2}{z_1}} + \{[E - \alpha]C_{\cdot\cdot\cdot} - \bar{\alpha}C_{\cdot\cdot\cdot}\}_{\frac{z_3}{z_2}} + \\
& \{[E - \alpha]D_{\cdot\cdot\cdot} - \bar{\alpha}D_{\cdot\cdot\cdot}\}_{\frac{z_2}{z_1}} + \{[E - \alpha]E_{\cdot\cdot\cdot} - \bar{\alpha}E_{\cdot\cdot\cdot}\}_{\frac{z_1}{z_2}} + \{[E - \alpha]F_{\cdot\cdot\cdot} - \bar{\alpha}F_{\cdot\cdot\cdot}\}_{\frac{z_1}{z_3}} = 0.
\end{aligned}
\quad (C13) $$

Crucially, this relation is identically satisfied as a consequence of Eqs. (C3) and (C9).

The periodic boundary conditions

$$ \psi_{\{j,k,\vec{L}\}} = \psi_{\{\vec{0},j,k\}} \quad (C14) $$

give the following relations for the spectral parameters $\{z_1, z_2, z_3\}$:

$$ z_a^{L-1} = \left(\frac{E - \bar{\alpha}/z_a}{\alpha}\right) \prod_{b \neq a}^{3} T(z_a, z_b), \quad a = 1, 2, 3, \quad (C15) $$

where [cf. Eq. (38)]

$$ T(z_a, z_b) = - \frac{E - \alpha - \bar{\alpha}/z_b}{E - \alpha - \bar{\alpha}/z_a}. \quad (C16) $$

[1] D. Landau and K. Binder, *A Guide to Monte Carlo Simulations in Statistical Physics* (Cambridge University Press, Cambridge, England, 2013).

[2] W. Krauth, *Statistical Mechanics: Algorithms and Computations* (Oxford University, New York, 2006).

[3] E. L. Pollock and D. M. Ceperley, *Path-integral computation of superfluid densities*, Phys. Rev. B **36**, 8343 (1987).

[4] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, *Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions*, Rev. Mod. Phys. **68**, 13 (1996).

[5] J. B. Kogut, *An introduction to lattice gauge theory and spin systems*, Rev. Mod. Phys. **51**, 659 (1979).

[6] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, *Equation of state calculations by fast computing machines*, J. Chem. Phys. **21**, 1087 (1953).

[7] R. J. Glauber, *Time-dependent statistics of the Ising model*, J. Math. Phys. (N.Y.) **4**, 294 (1963).

[8] M. Creutz, *Monte Carlo study of quantized SU(2) gauge theory*, Phys. Rev. D **21**, 2308 (1980).

[9] S. Geman and D. Geman, *Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images*, IEEE Trans. Pattern Anal. Mach. Intell. **PAMI-66**, 721 (1984).

[10] E. P. Bernard and W. Krauth, *Two-step melting in two dimensions: First-order liquid-hexatic transition*, Phys. Rev. Lett. **107**, 155704 (2011).

[11] B. Li, Y. Nishikawa, P. Höllmer, L. Carillo, A. C. Maggs, and W. Krauth, *Hard-disk pressure computations—A historic perspective*, J. Chem. Phys. **157**, 234111 (2022).

[12] R. H. Swendsen and J.-S. Wang, *Nonuniversal critical dynamics in Monte Carlo simulations*, Phys. Rev. Lett. **58**, 86 (1987).

[13] U. Wolff, *Collective Monte Carlo updating for spin systems*, Phys. Rev. Lett. **62**, 361 (1989).

[14] C. Dress and W. Krauth, *Cluster algorithm for hard spheres and related systems*, J. Phys. A **28**, L597 (1995).

[15] J. Liu and E. Luijten, *Rejection-free geometric cluster algorithm for complex fluids*, Phys. Rev. Lett. **92**, 035504 (2004).

[16] J. A. Anderson, E. Jankowski, T. L. Grubb, M. Engel, and S. C. Glotzer, *Massively parallel Monte Carlo for many-particle simulations on GPUs*, J. Comput. Phys. **254**, 27 (2013).

[17] M. Engel, J. A. Anderson, S. C. Glotzer, M. Isobe, E. P. Bernard, and W. Krauth, *Hard-disk equation of state: First-order liquid-hexatic transition in two dimensions with three simulation methods*, Phys. Rev. E **87**, 042134 (2013).
[18] B. Li, S. Todo, A. Maggs, and W. Krauth, Multithreaded event-chain Monte Carlo with local times, Comput. Phys. Commun. **261**, 107702 (2021).

[19] Y. Nishikawa, W. Krauth, and A. C. Maggs, Liquid-hexatic transition for soft disks, Phys. Rev. E **108**, 024103 (2023).

[20] B. Derrida and M. R. Evans, Bethe ansatz solution for a defect particle in the asymmetric exclusion process, J. Phys. A **32**, 4833 (1999).

[21] G. M. Schütz, Exactly solvable models for many-body systems far from equilibrium, in Phase Transitions and Critical Phenomena (Elsevier, New York, 2001), Vol. 19, pp. 1–251.

[22] F. Chen, L. Lovász, and I. Pak, Lifting Markov chains to speed up mixing, in Proceedings of the 17th Annual ACM Symposium on Theory of Computing (Association for Computing Machinery, New York, 1999), p. 275.

[23] P. Diaconis, S. Holmes, and R. M. Neal, Analysis of a nonreversible Markov chain sampler, Ann. Appl. Probab. **10**, 726 (2000).

[24] E. P. Bernard, W. Krauth, and D. B. Wilson, Event-chain Monte Carlo algorithms for hard-sphere systems, Phys. Rev. E **80**, 056704 (2009).

[25] M. Michel, S. C. Kapfer, and W. Krauth, Generalized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps, J. Chem. Phys. **140**, 054116 (2014).

[26] J. Bierkens, P. Fearnhead, and G. Roberts, The zig-zag process and super-efficient sampling for Bayesian analysis of big data, Ann. Stat. **47**, 1288 (2019).

[27] A. C. Maggs and W. Krauth, Large-scale dynamics of event-chain Monte Carlo, Phys. Rev. E **105**, 015309 (2022).

[28] M. Klement and M. Engel, Efficient equilibration of hard spheres with Newtonian event chains, J. Chem. Phys. **150**, 174108 (2019).

[29] E. A. J. F. Peters and G. de With, Rejection-free Monte Carlo sampling for general potentials, Phys. Rev. E **85**, 026703 (2012).

[30] T. A. Kampmann, D. Müller, L. P. Weise, C. F. Vorsmann, and J. Kierfeld, Event-chain Monte-Carlo simulations of dense soft matter systems, Front. Phys. **9**, 96 (2021).

[31] M. F. Faulkner, L. Qin, A. C. Maggs, and W. Krauth, All-atom computations with irreversible Markov chains, J. Chem. Phys. **149**, 064113 (2018).

[32] P. Höllmer, L. Qin, M. F. Faulkner, A. C. Maggs, and W. Krauth, JELLYFYSH-VERSION1.0—A PYTHON application for all-atom event-chain Monte Carlo, Comput. Phys. Commun. **253**, 107168 (2020).

[33] F. Spitzer, Interaction of Markov processes, Adv. Math. **5**, 246 (1970).

[34] H. Lacoin, The cutoff profile for the simple exclusion process on the circle, Ann. Probab. **44**, 3399 (2016).

[35] H. Lacoin, The simple exclusion process on the circle has a diffusive cutoff window, Ann. Inst. H. Poincaré Probab. Statist. **53**, 1402 (2017).

[36] C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin, Kinetics of biopolymerization on nucleic acid templates, Biopolymers **6**, 1 (1968).

[37] C. T. MacDonald and J. H. Gibbs, Concerning the kinetics of polypeptide synthesis on polyribosomes, Biopolymers **7**, 707 (1969).

[38] D. Chowdhury, L. Santen, and A. Schadschneider, Statistical physics of vehicular traffic and some related systems, Phys. Rep. **329**, 199 (2000).

[39] L.-H. Gwa and H. Spohn, Six-vertex model, roughened surfaces, and an asymmetric spin Hamiltonian, Phys. Rev. Lett. **68**, 725 (1992).

[40] L.-H. Gwa and H. Spohn, Bethe solution for the dynamical scaling exponent of the noisy Burgers equation, Phys. Rev. A **46**, 844 (1992).

[41] D. Kim, Bethe ansatz solution for crossover scaling functions of the asymmetric XXZ chain and the Kardar-Parisi-Zhang-type growth model, Phys. Rev. E **52**, 3512 (1995).

[42] B. Derrida, An exactly soluble non-equilibrium system: The asymmetric simple exclusion process, Phys. Rep. **301**, 65 (1998).

[43] K. Mallick, Some exact results for the exclusion process, J. Stat. Mech. (2011) P01024.

[44] T. Chou, K. Mallick, and R. K. P. Zia, Non-equilibrium statistical mechanics: From a paradigmatic model to biological transport, Rep. Prog. Phys. **74**, 116601 (2011).

[45] D. Aldous and P. Diaconis, Shuffling cards and stopping times, Am. Math. Mon. **93**, 333 (1986).

[46] D. Dhar, An exactly solved model for interfacial growth, Phase Transitions **9**, 51 (1987).

[47] J. Baik and Z. Liu, TASEP on a ring in sub-relaxation time scale, J. Stat. Phys. **165**, 1051 (2016).

[48] M. Kardar, G. Parisi, and Y.-C. Zhang, Dynamic scaling of growing interfaces, Phys. Rev. Lett. **56**, 889 (1986).

[49] Z. Lei, W. Krauth, and A. C. Maggs, Event-chain Monte Carlo with factor fields, Phys. Rev. E **99**, 043301 (2019).

[50] W. Krauth, Event-chain Monte Carlo: Foundations, applications, and prospects, Front. Phys. **9**, 229 (2021).

[51] F. C. Alcaraz, M. Droz, M. Henkel, and V. Rittenberg, Reaction-diffusion processes, critical dynamics, and quantum chains, Ann. Phys. (N.Y.) **230**, 250 (1994).

[52] F. H. L. Essler, J. Gipouloux, and W. Krauth, The generalized lifted TASEP (to be published).

[53] S. C. Kapfer and W. Krauth, Irreversible local Markov chains with rapid convergence towards equilibrium, Phys. Rev. Lett. **119**, 240603 (2017).

[54] Z. Lei and W. Krauth, Mixing and perfect sampling in one-dimensional particle systems, Europhys. Lett. **124**, 20003 (2018).

[55] Z. Lei and W. Krauth, Irreversible Markov chains in spin models: Topological excitations, Europhys. Lett. **121**, 10008 (2018).

[56] A. Sokal, Monte Carlo methods in statistical mechanics: Foundations and new algorithms, in Functional Integration: Basics and Applications, edited by C. DeWitt-Morette, P. Cartier, and A. Folacci (Springer, New York, 1997), pp. 131–192.

[57] D. A. Levin and Y. Peres, Markov Chains and Mixing Times, second ed. (American Mathematical Society, Providence, 2017).

[58] X. Zhang and W. Krauth, Diameters of symmetric and lifted simple exclusion models, J. Stat. Phys. **191**, 103 (2024).
[59] O. Golinelli and K. Mallick, *Bethe ansatz calculation of the spectral gap of the asymmetric exclusion process*, J. Phys. A **37**, 3321 (2004).

[60] S. F. Edwards and D. R. Wilkinson, *The surface statistics of a granular aggregate*, Proc. R. Soc. A **381**, 17 (1982).

[61] P. Caputo, T. Liggett, and T. Richthammer, *Proof of Aldous' spectral gap conjecture*, J. Am. Math. Soc. **23**, 831 (2010).

[62] N. Madras and A. D. Sokal, *The pivot algorithm: A highly efficient Monte Carlo method for the self-avoiding walk*, J. Stat. Phys. **50**, 109 (1988).

[63] J. de Gier and F. H. L. Essler, *Bethe ansatz solution of the asymmetric exclusion process with open boundaries*, Phys. Rev. Lett. **95**, 240601 (2005).

[64] J. de Gier and F. H. L. Essler, *Exact spectral gaps of the asymmetric exclusion process with open boundaries*, J. Stat. Mech. (2006) P12011.

[65] J. de Gier and F. H. L. Essler, *Slowest relaxation mode of the partially asymmetric exclusion process with open boundaries*, J. Phys. A **41**, 485002 (2008).

[66] V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin, *Quantum Inverse Scattering Method and Correlation Functions* (Cambridge University Press, Cambridge, England, 1997), Vol. 3.

[67] B. Derrida, M. R. Evans, V. Hakim, and V. Pasquier, *Exact solution of a 1D asymmetric exclusion model using a matrix formulation*, J. Phys. A **26**, 1493 (1993).