Integrability of two-species partially asymmetric exclusion processes

Ivan Lobaskin^1, Martin R Evans^1 and Kirone Mallick^2

$^1$ School of Physics and Astronomy, University of Edinburgh, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, United Kingdom
$^2$ Institut de Physique Théorique, Université Paris-Saclay, CEA and CNRS, 91191 Gif-sur-Yvette, France
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
August 27, 2024

We work towards the classification of all one-dimensional exclusion processes with two species of particles that can be solved by a nested coordinate Bethe Ansatz.
Using the Yang-Baxter equations, we obtain conditions on the model parameters that ensure that the underlying system is integrable.
Three classes of integrable models are thus found.
Of these, two classes are well known in literature, but the third has not been studied until recently, and never in the context of the Bethe ansatz.
The Bethe equations are derived for the latter model as well as for the associated dynamics encoding the large deviation of the currents.

§ INTRODUCTION

In statistical physics far from equilibrium, the one-dimensional asymmetric simple exclusion process (ASEP) is a paradigmatic model, playing a similar role to the Ising model for equilibrium systems [1, 2, 3].
Its simplicity means that it can be used as a generic model in many fields, including biophysics [4, 5, 6], traffic modelling [7] and microfluidics [8].
At the same time, a multitude of exact results have been derived, which is rare for interacting $N$-body problems.
The two main approaches for obtaining exact results have been the matrix product ansatz [9] and the Bethe ansatz (BA) [10, 11, 12].
For general reviews of these methods in the ASEP context, see [2] for the matrix product approach (which will not be used  in the present work) and [13] for the BA.

The BA was first developed to diagonalize one-dimensional quantum spin chain Hamiltonians [14], but has since been extended to many other fields, such as 2D vertex models [15] and ASEPs.
It was first shown that it can be applied to the  exclusion process on a ring [10, 11, 12].
Since then, it has also been used on systems with various boundary conditions [16, 17, 18, 19], partially asymmetric (PASEP) models [20, 21, 22, 23, 24, 25] and multi-species cases [26, 21, 22, 27].
Using a deformation of the Markov operator, it was first shown in [28] that the BA can be used very effectively to derive higher order particle displacement statistics, such as the diffusion coefficient.
For some particularly simple cases, it has even been possible to calculate the full statistics to all orders [28, 25].

The question of which systems can be solved by the BA, i.e. are integrable,  has a long history in the quantum and classical equilibrium contexts [29].
An answer to this problem is through the Yang-Baxter equations (YBE), which can be interpreted as consistency conditions for the BA; indeed, by checking whether the YBE hold on a general class of models, it is possible to select values of the  parameters  for which  exact solutions exist.
Though this approach was  used for quantum [30] and vertex model  [31], it has been applied less systematically to  study classical interacting particles systems, see for example  [32, 27, 33, 34].

In this paper, we  conduct such a systematic study for two-species exclusion processes on a ring.

By postulating a very general coordinate BA and deriving the YBE, we are able to assess the integrability of a wide range of models, which suggests a classification of integrable two-species systems.

We find that there are three families of integrable models: two of them are well-known in the literature, but there exists another, special, model that has only come to attention recently [35] and has not yet been studied from the viewpoint of integrability.
There, its steady state was solved using the matrix product approach, and the phase diagram, currents and density profiles were calculated.
This model has the peculiar feature that it is only integrable when there is only one particle of one of the species, which can be seen clearly in the YBE.
We conjecture that our study entirely exhausts all integrable two-species exclusion processes.

The remainder of the paper is structured as follows.
In BAderivation, we show how the time evolution problem of a two-species PASEP can be solved by a general form of the coordinate BA. In addition, we show that  a very similar ansatz can be used to solve the conditioned (or deformed) time evolution problem, which allows direct calculation of the large deviations of the current.
In YBE, we derive the YBE  and present  the three classes of solutions.
Two of these are well-known in literature [36, 37, 28, 21, 27], but the third has only come to attention recently [35] and has not yet been studied in the context of the BA.
In newmodelBA, we focus on this new model and derive the Bethe equations for the deformed time evolution problem.
In <ref>, we explain the procedure that we used to find the solutions of the YBE.
Finally, in <ref>, we fill in some algebraic steps in the calculation of newmodelBA.

§ BETHE ANSATZ SOLUTION OF TWO-SPECIES PASEP PROBLEM

§.§ Problem statement

We consider an $L$ site ring with particles of two-species hopping stochastically in continuous time.
Let there be $M_1$ particles of species $1$ and $M_2$ particles of species $2$, with $M_1+M_2=M$.
For the most general two-species partially asymmetric problem with overtaking, the following processes take place with the given rates:

\begin{eqnarray*}
\underset{q_1}{\overset{p_1}{\rightleftarrows}}
01 \quad ; \quad
\underset{q_2}{\overset{p_2}{\rightleftarrows}}
02 \quad ; \quad
\underset{b _{21}}{\overset{b _{12}}{\rightleftarrows}}
21 \quad ;
\end{eqnarray*}

where $0,1,2$ represent empty sites, and particles of species $1$ and $2$, respectively.

Fixing the reference frame, we denote by $Q_j$ and $x_j$ respectively the species and position of particle $j$, counting from the left, where $Q_j\in \{1,2\}$ and $x_j \in {\mathbb Z} _L$.
As a shorthand, we denote the configuration of the system as ${\cal C} = \{ Q_1 \dots Q_M ; x_1 \dots x_M \}$.
The time evolution of the probability of a configuration ${\cal C}$ at time $t$ is given by the master equation

\begin{equation}
\frac{\partial P_t ({\cal C})}{\partial t} =
\sum\limits_{{\cal C}'} {\cal M}({\cal C},{\cal C}') P_t({\cal C}') \;,
\label{master}
\end{equation}

where $\cal M$ is the continuous-time  Markov operator containing the transition rates from configuration ${\cal C}'$ to ${\cal C}$.
We wish to solve the eigenvalue-eigenvector problem for the Markov operator.
It can be written in the form

\begin{eqnarray}
\fl
\lambda \psi ^{Q_1 \dots Q_M} (x_1 \dots x_M) =
-[M_1 (p_1+q_1) + M_2(p_2+q_2)]\psi ^{Q_1 \dots Q_M} (x_1\dots x_M) \nonumber
\\
\quad +\sum\limits_{j=1}^{M} [p_{Q_j} \psi ^{Q_1 \dots Q_M} (\dots  x_j -1\dots  )
+q_{Q_j} \psi ^{Q_1 \dots Q_M} (\dots  x_j +1\dots )] \;,
\label{evproblem}
\end{eqnarray}

where $\lambda$ is an eigenvalue and $\psi$ is the corresponding eigenvector.

This equation holds in configurations in which no two particles are on adjacent sites.
It can be assumed to hold in all configurations, provided an additional set of equations, called collision equations, also hold.
Considering the $j$-th and $(j+1)$-th particles, with $x_{j+1}=x_{j}+1$, these are given by the following, where $Q_{j,j+1}=1,2$,

\begin{eqnarray}
\fl -(p_{Q_j} + q_{Q_{j+1}})\psi^{\dots Q_jQ_{j+1}\dots}(\dots x_j,x_j+1 \dots)
+ p_{Q_{j+1}} \psi^{\dots Q_jQ_{j+1}\dots}(\dots x_j,x_j \dots) \nonumber \\
\fl \qquad + q_{Q_j} \psi^{\dots Q_jQ_{j+1}\dots}(\dots x_j+1,x_j+1 \dots)
= -b_{Q_j Q_{j+1}}\psi^{\dots Q_jQ_{j+1}\dots}(\dots x_j,x_j+1 \dots)\nonumber \\
\fl \qquad + b_{Q_{j+1}Q_j}\psi^{\dots Q_{j+1}Q_{j}\dots}(\dots x_j,x_j+1 \dots)\;,
\label{collpsi}
\end{eqnarray}

where $b_{22}=b_{11}=0$.

Finally, the periodic boundary conditions are imposed via the requirement

\begin{equation}
\psi ^{Q_1 Q_2 \dots Q_M} (x_1, x_2 \dots x_M) =
\psi ^{Q_2\dots Q_M Q_1} (x_2\dots x_M, x_1+L) \;. \label{pbc}
\end{equation}

Then  equations evproblem, collpsi, pbc define the eigenvalue-eigenvector problem to be solved.

§.§ Bethe ansatz solution

We use the following modified nested coordinate BA,

\begin{equation}
\psi ^{Q_1 \dots Q_M} (x_1 \dots x_M) =
r_{Q_1} ^{x_1} \dots r_{Q_M} ^{x_M} \sum\limits _\sigma A^{Q_1 \dots Q_M}_\sigma
z_{\sigma (1)}^{x_1} \dots  z_{\sigma (M)}^{x_M} \;,
\label{ansatz}
\end{equation}

where $\sigma $ indicates the permutations of the integers $1, \dots ,M$; the $A$'s and $z$'s are undetermined complex amplitudes and Bethe roots respectively; and the $r$'s are species dependent pre-factors, to be determined.
This is the most general known coordinate ansatz and is a straightforward nested generalization of the ansatz introduced in [26].
There, it was used for a totally asymmetric problem with one defect particle.

The ansatz ansatz is found to solve evproblem provided the following condition holds:

\begin{equation}
p_1q_1=p_2q_2\;. \label{pqpq}
\end{equation}

In this case, we can set

\begin{equation}
p_1 = \alpha p_2\;,\qquad q_1=q_2/\alpha \;, \label{alpha}
\end{equation}

for some constant $\alpha $.
Then the ansatz ansatz works, with

\begin{equation}
r_1 = \alpha \;,\qquad r_2 =1\;.\label{r1r2}
\end{equation}

We assume $\alpha\neq0$, as otherwise the ansatz vanishes.
Notably, this excludes the model studied by Arndt [38, 39], which can be defined with our notation as $p_1=0,q_2=0$.
However, as this is believed to be non-integrable [40], this is not of major concern.

Putting ansatz with alpha into evproblem gives the following expression for the eigenvalue in terms of the Bethe roots,

\begin{equation}
\lambda = -[M_1 (\alpha p_2+q_2/\alpha ) + M_2(p_2+q_2)] + \sum\limits_{j=1}^{M}( p_2z_j ^{-1} + q_2z_j )\;. \label{lambda}
\end{equation}

We remark that condition pqpq comes from the form of the ansatz ansatz.
However, a similar condition was derived from the the algebraic Bethe ansatz for the $15$-vertex model [31].
These two cases are linked, as with certain choices of parameters, the $15$-vertex model can be mapped to a two-species PASEP with overtaking.
This suggests that condition pqpq is generic in the context of exactly solvable models and not a special feature of the ansatz ansatz.

In the subsequent calculations, it is convenient to use the following dimensionless parameters

\begin{eqnarray}
x=q_2/p_2 \;,\\
\beta_{12}=b_{12}/p_2 \;,\qquad \beta_{21}=b_{21}/p_2 \;.
\end{eqnarray}

Putting the ansatz ansatz into the collision equations collpsi gives

\begin{eqnarray}
[ {\cal D}_{ij} -(\alpha + x/\alpha )  z_j ] A^{11}_{ij}
+ (i\leftrightarrow j)= 0 \;, \label{coll_nested_a} \\ [5pt]
{[ {\cal D}_{ij}-(1 + x )  z_j ]} A^{22}_{ij}
+ (i\leftrightarrow j) = 0 \;, \label{coll_nested_b} \\ [5pt]
{[ ({\cal D}_{ij}-\kappa_2  z_j ) A^{12}_{ij}
- \alpha \beta_{21} z_j A^{21}_{ij}]}
+ (i\leftrightarrow j) = 0 \;, \label{coll_nested_c} \\ [5pt]
{[ ({\cal D}_{ij}-\kappa_1  z_j ) A^{21}_{ij}
- \alpha^{-1}\beta_{12} z_j A^{12}_{ij}]}
+ (i\leftrightarrow j) = 0 \;,\label{coll_nested_d}
\end{eqnarray}

\begin{eqnarray}
{\cal D}_{ij} = 1+xz_iz_j \;, \\
\kappa _1 = 1+x/\alpha -\beta_{21} \;,\qquad \kappa _2 = \alpha +x-\beta_{12} \;,
\end{eqnarray}

and the notation $(i\leftrightarrow j)$ indicates a term similar to the explicitly written one but with the indices $i,j$ swapped.

The system of  equations coll_nested_a–coll_nested_b can be written compactly in vector form

\begin{eqnarray}
[{\cal D}_{ij}1_{4\times 4}- z_j H]{\vec A}_{ij} + [{\cal D}_{ij}1_{4\times 4}- z_i H]{\vec A}_{ji} = 0 \;,
\end{eqnarray}

\begin{eqnarray}
H = \left(
\begin{array}{cccc}
\alpha + x/\alpha & 0 & 0 & 0  \\
0 & \kappa _2 & \alpha \beta_{21} & 0 \\
0 & \beta_{12}/\alpha  & \kappa _1 & 0 \\
0 & 0 & 0 & 1 + x
\end{array}
\right) \;, \\[5pt]
{\vec A}_{ij} = (A^{11}_{ij},A^{12}_{ij},A^{21}_{ij},A^{22}_{ij})^{\rm T} \;.
\end{eqnarray}

This can be formally rewritten as

\begin{eqnarray}
{\vec A}_{ji} =
-\frac{{\cal D}_{ij}1_{4\times 4}- z_j H}{{\cal D}_{ij}1_{4\times 4}- z_i H}
{\vec A}_{ij} \;, \label{collnestedfra}
\end{eqnarray}

where division is to be understood as multiplication by the matrix inverse from the left.
Performing this calculation explicitly, we get relations of the form

\begin{equation}
A^{Q_1Q_2}_{ji} =
\sum\limits_{Q_1',Q_2'=1,2}S^{Q_1Q_2}_{Q_1'Q_2'}(z_i , z_j)
A^{Q_1'Q_2'}_{ij} \;, \label{formfactors}
\end{equation}

in which $S^{Q_1Q_2}_{Q_1'Q_2'}(z_i , z_j)$ are called form factors.
Of these, only the following are non-zero

\begin{eqnarray}
S^{11}_{11}(z_i,z_j) = -\frac{{\cal D}_{ij}-(\alpha +x/\alpha)z_j}{{\cal D}_{ij}-(\alpha +x/\alpha)z_i} \;, \label{s1111}
\\
S^{22}_{22}(z_i,z_j) = -\frac{{\cal D}_{ij}-(1+x)z_j}{{\cal D}_{ij}-(1+x)z_i} \;, \label{s2222}\\
S^{12}_{12}(z_i,z_j) = -\Phi_{ij}[{\cal D}_{ij}^2-{\cal D}_{ij}(\kappa_1 z_j + \kappa _2 z_i)+(\kappa _1 \kappa_2 - \beta_{12} \beta_{21})z_iz_j] \;, \label{s1212}\\
S^{12}_{21}(z_i,z_j) = -\Phi_{ij}\alpha \beta_{21} {\cal D}_{ij}(z_i-z_j) \;, \label{s1221}\\
S^{21}_{12}(z_i,z_j) = -\Phi_{ij}\alpha^{-1}\beta_{12} {\cal D}_{ij}(z_i-z_j) \;, \label{s2112}\\
S^{21}_{21}(z_i,z_j) = -\Phi_{ij}[{\cal D}_{ij}^2-{\cal D}_{ij}(\kappa_2 z_j + \kappa _1 z_i)+(\kappa _1 \kappa_2 - \beta_{12} \beta_{21})z_iz_j] \label{s2121}\;,
\end{eqnarray}

\begin{eqnarray}
\Phi_{ij} =
[({\cal D}_{ij}-\kappa _1z_i)({\cal D}_{ij}-\kappa _2z_i)
\end{eqnarray}

§.§ Bethe ansatz for deformed Markov operator

Large deviations  in interacting particle systems can be derived by conditioning the dynamics  on a given, atypical value, of the observable under consideration.
The resulting evolution is governed by a deformation of the original Markov operator.
This method goes back to Donsker and Varadhan (see [41, 42] and references therein).
In the context of ASEP, this idea has been has been used systematically in conjunction with the BA, leading to exact results for the large deviation function of the current [28, 26, 24, 43, 44].
We now outline this calculation for the two-species case, merely for completeness, since much of the detail is similar to the undeformed case, analyzed in the previous section.

We define observables $Y^{1}_t,Y^{2}_t,Y^{12}_t$, which count the number of processes $10\to 01\;,20\to 02\;,12\to 21$ minus the inverse processes respectively, up to time $t$.
Thus, for example, $Y^1_t+Y^{12}_t$ is the net total displacement of particles of species $1$ to the right.
We can split the Markov operator ${\cal M}$ into

\begin{eqnarray}
{\cal M} = {\cal M}_0 + \sum\limits_{i=1,2,12} ({\cal M}_{i,+}+{\cal M}_{i,-}) \;,
\end{eqnarray}

where ${\cal M}_0$ contains the diagonal entries and ${\cal M}_{i,+(-)}$ contains the rates of the processes that increase (decrease) $Y^i_t$ by 1.
Let $P_t({\cal C},Y^1,Y^2,Y^{12})$ be the joint probability of configuration ${\cal C}$ with $Y^i_t=Y^i$ at time $t$.
Its generating function is defined as

\begin{equation}
F_t ({\cal C},\gamma _1,\gamma _2, \gamma _{12}) =
\sum\limits _{Y^1,Y^2,Y^{12} = -\infty} ^{\infty}
\rme ^{\gamma _1 Y^1 +\gamma _2 Y^2 +\gamma _{12} Y^{12}} P_t({\cal C},Y^1,Y^2,Y^{12}) \;.
\end{equation}

Note that by marginalizing over configurations, we obtain the moment generating function of $Y^{1}_t,Y^2_t,Y^{12}_t$,

\begin{equation}
\sum\limits_{\cal C} F_t ({\cal C},\gamma _1,\gamma _2, \gamma _{12}) =
\langle \rme ^{\gamma _1 Y^1_t+\gamma _2 Y^2_t +\gamma _{12} Y^{12}_t}\rangle \;.
\end{equation}

A key property is that this satisfies a large deviation principle as $t\to \infty$,

\begin{equation}
\langle \rme ^{\gamma _1 Y^1_t+\gamma _2 Y^2_t +\gamma _{12} Y^{12}_t}\rangle
\sim \rme ^{\mu (\gamma _1,\gamma _2,\gamma _{12})t}\;, \label{ldp}
\end{equation}

with some rate function $\mu (\gamma _1,\gamma _2,\gamma _{12})$.

From master, we get the evolution equation

\begin{equation}
\partial _t F_t ({\cal C},\gamma _1,\gamma _2, \gamma _{12}) =
\sum\limits_{{\cal C}'} {\cal M}_{\gamma }({\cal C},{\cal C}')F_t ({\cal C'},\gamma _1,\gamma _2, \gamma _{12}) \;,
\end{equation}

where ${\cal M}_\gamma $ is the deformed Markov operator

\begin{equation}
{\cal M} _\gamma = {\cal M}_0 + \sum\limits_{i=1,2,12}
(\rme ^{\gamma _i}{\cal M}_{i,+}+\rme ^{-\gamma _i}{\cal M}_{i,-}) \;.
\end{equation}

The eigenvalue-eigenvector problem for ${\cal M}_\gamma $ then becomes

\begin{eqnarray}
\fl
\lambda (\gamma _1,\gamma_2,\gamma_{12})\psi ^{Q_1 \dots Q_M} (x_1 \dots x_M) =
-[M_1 (p_1+q_1) + M_2(p_2+q_2)]\psi ^{Q_1 \dots Q_M} (x_1\dots x_M) \nonumber \\
\fl \qquad +\sum\limits_{j=1}^{M}
[p_{Q_j} \rme^{\gamma _{Q_j}}\psi ^{Q_1 \dots Q_M} (\dots  x_j -1\dots  )
+q_{Q_{j}} \rme^{-\gamma _{Q_j}} \psi ^{Q_1 \dots Q_M} (\dots  x_j +1\dots )] \;.
\end{eqnarray}

As with the undeformed case, this equation is strictly valid only for configurations with no two particles being on adjacent sites, which is why $\gamma_{12}$ does not appear here.
The problem has to be supplemented with collision equations, which are given below.

A particularly important eigenvalue, is the one with the largest real part, which can be identified with the rate function $\mu (\gamma _1,\gamma _2,\gamma _{12})$ in ldp.
From the Perron-Frobenius theorem it follows that this eigenvalue is unique.
In the limit $\gamma _1\to 0,\;\gamma_2\to 0,\;\gamma_{12}\to 0$, it converges to the unique $0$ eigenvalue of the (undeformed) Markov operator ${\cal M}$, whose eigenvector is the stationary measure.

The collision equations read

\begin{eqnarray}
\fl -(p_{Q_j} + q_{Q_{j+1}})\psi^{\dots Q_jQ_{j+1}\dots}(\dots x_j,x_j+1 \dots)
+p_{Q_{j+1}}  \rme^{\gamma_{Q_{j+1}}} \psi^{\dots Q_jQ_{j+1}\dots}(\dots x_j,x_j \dots)
\nonumber \\
\fl \qquad + q_{Q_j} \rme^{-\gamma_{Q_j}}\psi^{\dots Q_jQ_{j+1}\dots}(\dots x_j+1,x_j+1 \dots)
= -b_{Q_j Q_{j+1}}\psi^{\dots Q_jQ_{j+1}\dots}(\dots x_j,x_j+1 \dots)\nonumber \\
\fl \qquad + b_{Q_{j+1}Q_j}\rme^{\gamma_{Q_{j+1}Q_j}}\psi^{\dots Q_{j+1}Q_{j}\dots}(\dots x_j,x_j+1 \dots)\;,
\label{collpsidef}
\end{eqnarray}

where $\gamma_{21}=-\gamma_{12}$.
This problem is solved by the ansatz

\begin{eqnarray}
\psi ^{Q_1 \dots Q_M} (x_1 \dots x_M) =
r_{Q_1} ^{x_1} \dots r_{Q_M} ^{x_M} \sum\limits _\sigma A^{Q_1 \dots Q_M}_\sigma
z_{\sigma (1)}^{x_1} \dots  z_{\sigma (M)}^{x_M} \;, \label{ansatzdefa}\\
r_1 = \alpha \rme ^{\gamma _1} \;,\;
r_2 = \rme ^{\gamma _2}
\;. \label{ansatzdefb}
\end{eqnarray}

Plugging the ansatz ansatzdefa-ansatzdefb into the collision equations collpsidef gives after some simplification an equation similar to collnestedfra, with the only difference being that now

\begin{eqnarray}
H = \left(
\begin{array}{cccc}
\alpha + x/\alpha & 0 & 0 & 0  \\
0 & \kappa _2 & \alpha \beta_{21}\rme^{\bar \gamma} & 0 \\
0 & \alpha ^{-1}\beta_{12} \rme^{-{\bar \gamma}} & \kappa _1 & 0 \\
0 & 0 & 0 & 1 + x
\end{array}
\right) \;,
\end{eqnarray}

\begin{equation}
{\bar \gamma} = \gamma_1-\gamma_2-\gamma_{12}\;.
\end{equation}

Then repeating the calculation of the form factors, we get that most of them are the same as in the undeformed case, with only the following modifications

\begin{eqnarray}
S^{12}_{21}(z_i,z_j) = -\Phi_{ij}\alpha \beta_{21} \rme^{{\bar\gamma}} {\cal D}_{ij}(z_i-z_j) \;,\\
S^{21}_{12}(z_i,z_j) = -\Phi_{ij}\alpha^{-1}\beta_{12}\rme^{-{\bar\gamma}} {\cal D}_{ij}(z_i-z_j) \;.
\end{eqnarray}

§ YANG-BAXTER EQUATIONS

§.§ Derivation of YBE

The Yang-Baxter relations are obtained by analyzing  collisions that involve 3 particles.
In this case, two sequences of binary collisions  can  result in the same permutation: $ijk\to ikj \to kij \to kji$ and $ijk \to jik \to jki \to kji$.
In order for the nested coordinate BA to be consistent, the products of form factors obtained from these two sequences must be the same.
This requirement yields the YBE (see [21] for a similar analysis).
We shall follow this procedure both for the regular and deformed cases, with the only difference being the explicit forms of the form factors.
Somewhat surprisingly, for the systems at hand, the results are insensitive to the deformation parameters.

The consistency condition gives some relations between the amplitudes $A^{Q_1Q_2Q_3}_{ijk}$ and $A^{Q_1'Q_2'Q_3'}_{kji}$, which for the two-species case is in principle an equation for $8\times 8$ matrices, called $R$-matrices.
The YBE can be formulated as an algebraic relation for the $R$-matrices.
In the context of driven diffusive systems this was done for instance in [27].
However, due to particle conservation, the relations can be simplified by dividing them into disjoint blocks with a fixed number of particles of each species.
It is not hard to verify that the $111$ and $222$ blocks are trivially satisfied for any choice of parameters.
For the $112$ block, we define

\begin{eqnarray}
\fl
\Xi (z_i,z_j)=
\left(
\begin{array}{ccc}
S^{11}_{11} & 0 & 0 \\
0 & S^{12}_{12} & S^{12}_{21} \\
0 & S^{21}_{12} & S^{21}_{21}
\end{array}
\right) _{(z_i,z_j)} \;, \;
{\tilde\Xi} (z_i,z_j)=
\left(
\begin{array}{ccc}
S^{12}_{12} & S^{12}_{21} & 0 \\
S^{21}_{12} & S^{21}_{21} & 0 \\
0 & 0 & S^{11}_{11}
\end{array}
\right) _{(z_i,z_j)} \;.
\end{eqnarray}

Then the YBE are given by

\begin{equation}
\Xi  (z_j,z_k) {\tilde\Xi} (z_i,z_k) \Xi  (z_i,z_j)-
{\tilde\Xi} (z_i,z_j) \Xi (z_i,z_k) {\tilde\Xi} (z_j,z_k)=0 \;.
\label{yb}
\end{equation}

The YBE for the $122$ block can be derived similarly.

We remark that all the entries of the matrices $\Xi$ are rational functions of the Bethe roots.
Thus, when we combine the two matrix products in yb, the result is a $3\times 3$ matrix whose entries are rational functions of the Bethe roots.
Their explicit forms are evidently very complicated but they can be analyzed with the help of a symbolic programming language, such as Mathematica.

§.§ Solutions of YBE

We wish to find conditions on the model parameters $(\alpha ,x, \beta _{12},\beta_{21} )$ that ensure that yb are satisfied identically, for all values of the Bethe roots.

Firstly, we note that the numerators of the entries are polynomials in $z_i,z_j,z_k$, which must vanish identically.
We may therefore examine individual coefficients of $z_i,z_j,z_k$ in the numerators and require them to vanish.
This can be done in many ways, as there are many possible coefficients to look at.
Fortunately, taking only a few is sufficient to obtain conditions that are strong enough to satisfy the full set of equations.
One possible method is outlined in <ref>.

We obtain three general solutions.
For the sake of conciseness, we do not list one-species models (which trivially form subsets of the models listed below and are known to be integrable), and models which are equivalent to the given ones up to permutation of species labels (including empty sites as a “zeroth species").

We remark that all the following models satisfy the relation

\begin{equation}
b_{12}b_{21}=p_2q_2 \;,
\end{equation}

which was not imposed a priori.
This shows that, the condition pqpq cannot be circumvented by a relabelling of the species, which further supports the claim that it might be a necessary condition.

§.§.§ Solution 1: TASEP case.

The first solution is defined by the conditions

\begin{equation}
x=0 \;, \qquad \beta_{21}=0 \;,
\end{equation}

which corresponds to the process

\begin{eqnarray*}
\overset{\alpha p_2}{\to}
01 \quad;\quad
\overset{p_2}{\to}
02 \quad;\quad
\overset{b_{12}}{\to}
21 \;.
\end{eqnarray*}

This is a totally asymmetric process, with particles of species  1 hopping at rate $\alpha p_2$ and overtaking  particles  of species 2 at rate $b_{12}$.
The matrix product solution for the steady state was first given in [36].
The case of a single particle of species 2 (referred to as a defect) has been much studied in the literature [37, 26, 27].

§.§.§ Solution 2: PASEP with first- and second-class particles.

The second solution is defined by the conditions

\begin{equation}
\alpha = 1 \;, \qquad \beta_{12} = 1 \;, \qquad \beta_{21} = x \; ,
\end{equation}

which corresponds to the process

\begin{eqnarray*}
\underset{q_2}{\overset{p_2}{\rightleftarrows}}
01 \quad;\quad
\underset{q_2}{\overset{p_2}{\rightleftarrows}}
02 \quad ; \quad
\underset{q_2}{\overset{p_2}{\rightleftarrows}}
21 \; .
\end{eqnarray*}

This is a two-species partially asymmetric process, with both species hopping at the same rates and the species 1 particles overtaking species 2 particles with their usual hopping rates.
Since species 1 does not distinguish between species 2 and a vacancy, species 1 can be seen as having priority in the dynamics and is therefore referred to in the literature as first-class and species 2 as second-class.
This case has been studied in the literature as well, with the matrix product solution for the steady state given in [36] and  BA equations derived in [21, 22].
We remark that although the Bethe equations have been derived, the statistics of the current fluctuations have not yet been calculated.

§.§.§ Solution 3: PASEP with a single first-class defect.

The third solution is defined by the conditions

\begin{equation}
\beta_{12} = \alpha \;, \qquad \beta_{21} = x/\alpha \;, \qquad M_1 = 1 \; ,
\end{equation}

which corresponds to the process

\begin{eqnarray*}
\underset{q_2/\alpha }{\overset{\alpha p_2}{\rightleftarrows}}
01 \quad ; \quad
\underset{q_2}{\overset{p_2}{\rightleftarrows}}
02 \quad ; \quad
\underset{q_2/\alpha }{\overset{\alpha p_2}{\rightleftarrows}}
21  \;.
\end{eqnarray*}

This is a partially asymmetric process, with a single first-class defect particle with hopping rates $\alpha p_2,q_2/\alpha $ that overtakes the second-class particles with its usual hopping rates.
This model has not been analyzed in the literature until recently, when its steady state was solved using a matrix product approach [35].

We remark that the constraint on particle number $M_1=1$ is unusual in the context of YBE.
It emerges here as the given conditions on the rates are sufficient to satisfy yb for the $122$ block but not the $112$ block.
This means that the ansatz works for this general case only if there is exactly one particle of species 1 in the system.
To satisfy yb for the $112$ block as well, it is necessary to impose either $\alpha =1$ or $\alpha=x$ or $x=0$.
In all these cases, this solution reduces down to solution 1 or 2.

§ DERIVATION OF BETHE EQUATIONS FOR SOLUTION 3

As solution 3 has not appeared much in the literature, it presents a case of particular interest.
As mentioned previously, it was studied using a matrix product approach in [35].
The phase diagram was shown to consist of localized phases, in which the defect particle has no macroscopic effect on the system, and a shock phase, in which the defect creates two extended regions with different bulk densities.
The density profiles and currents were also calculated.
The BA approach presented here complements these results, as it should allow one to calculate the fluctuations of the current as well.

To that end, we now derive the Bethe equations for the deformed eigenvalue problem from BAdeformed.
As we have only one particle of species $1$, it is not necessary to use the full nested ansatz.
Instead, we can fix the particle indexing so that the species $1$ particle is always the leftmost one.
We can always do this using the boundary condition pbc, which gives the following relation for amplitudes

\begin{eqnarray}
A^{Q_1Q_2\dots }_{ij\dots} = \rme^{\gamma _{Q_1} L}z_i^L A^{Q_2\dots Q_1}_{j\dots i}\;.
\end{eqnarray}

Now that the order of species is always fixed to be $12\dots 2$, we can drop the upper indices without ambiguity.
Then we get for 22, 12 and 21 collisions respectively,

\begin{eqnarray}
A_{\dots ji\dots} =
-\frac{{\cal D}_{ij}-(1+x)z_j}{{\cal D}_{ij}-(1+x)z_i} A_{\dots ij\dots}
\;, \label{coll_a} \\ [5pt]
{[ ({\cal D}_{ij} -xz_j )A_{ij\dots}
- x \rme^{\bar \gamma} z_i (\rme ^{\gamma _2}z_j)^L A _{i\dots j} ]}
+ (i\leftrightarrow j) =0\;, \label{coll_b}\\ [5pt]
{[ ({\cal D}_{ij}-z_i )(\rme ^{\gamma _2}z_j) ^L A_{i\dots j}
- \rme^{-{\bar\gamma}} z_j A_{ij\dots } ]} + (i\leftrightarrow j) =0\;. \label{coll_c}
\end{eqnarray}

We will denote the form factor in coll_a as $S^{22}_{22}(z_i,z_j)$ for conciseness.
Applying coll_a $M_2-1$ times to coll_b and coll_c gives

\begin{eqnarray}
\left[ {\cal D}_{ij}-xz_j
+ x z_i S^{22}_{22}(z_i,z_j)E_j \right]
A_{ij\dots} + (i\leftrightarrow j) = 0 \;, \label{beq_0a}\\ [5pt]
\left[ ({\cal D}_{ij}-z_i) S^{22}_{22}(z_i,z_j)E_j + z_j \right] A_{ij\dots}
+ (i\leftrightarrow j)=0 \label{beq_0b}\;,
\end{eqnarray}

\begin{equation}
E_j = \rme^{\bar\gamma}(\rme^{\gamma _2}z_j) ^L
\prod\limits _{k=1}^{M}S^{22}_{22}(z_j,z_k) \;.
\end{equation}

In order for the eigenvector not to vanish, the system of equations beq_0a-beq_0b for $A_{ij\dots},A_{ji\dots}$ must be degenerate.
The equation resulting from this requirement can be simplified, although this procedure is algebraically somewhat cumbersome and not very illuminating.
We give the details in <ref> for completeness.
In this calculation, it becomes convenient to consider the following transformation of the Bethe roots [20, 24],

\begin{equation}
z_i = \frac{1-y_i}{1-xy_i} \;. \label{ztoy}
\end{equation}

Ultimately, one obtains a set of equations, called the Bethe equations, for the transformed Bethe roots $y_i$,

\begin{equation}
\rme ^{{\bar\gamma}+L\gamma_2}
\left( \frac{1-y_i}{1-xy_i}\right) ^L
\prod\limits_{k=1}^{M}\frac{xy_i-y_k}{y_i-xy_k}
= -\frac{B-y_i}{Bx-y_i} \;,
\label{beq_1}
\end{equation}

where $B$ is a constant.
The constant $B$ can be fixed by multiplying the Bethe equations for all $i$,

\begin{equation}
\prod\limits_{i=1}^{M}\frac{B-y_i}{Bx-y_i}
= \alpha ^{-L}\rme ^{-(L-M){\bar\gamma}-L\gamma_{12}} \;,
\label{bfix}
\end{equation}

where we have implemented the periodic boundary condition pbc as

\begin{equation}
\alpha ^L \rme ^{L(\gamma_1+M_2\gamma_2)}z_1 ^L \dots z_M ^L = 1 \;. \label{periodicity}
\end{equation}

It is interesting to compare beq_1 to the corresponding equations for the well-known cases of a one-species PASEP [24] and the TASEP with a defect [26, 27].
The one-species PASEP can also be simplified using a change of variable similar to ztoy.
The resulting equation involves a product of fractions of linear differences of the transformed Bethe roots $y_i$, like beq_1.

Meanwhile, the TASEP with a defect also features an additional constant, like $B$, that has to be fixed as part of solving the Bethe equations.
As will be shown in a forthcoming publication, this leads to one of the Bethe roots converging to a different point in the complex plane to the other roots in the $\gamma_1\to 0,\gamma_2\to0,\gamma_{12}\to0$ limit.

We now briefly comment on how one would proceed to solve the Bethe equations.
beq_1 is a system of polynomial equations of degree $L+M+1$ for each Bethe root $y_i$, in which the coefficients in turn depend on all the roots.
One would first need to solve beq_1 consistently with bfix, then choose $M$ roots to satisfy periodicity for some $L$-th root of unity.
This would then give an eigenvalue through lambda.
In practice this is very challenging both analytically and numerically.
However, special functional methods have been developed that allow to calculate perturbative expansions of the largest eigenvalue in terms of $\gamma_1,\gamma_2,\gamma_{12} $ [28, 26, 24, 25].
In a forthcoming publication, we will show how these methods can be applied to the problem at hand.

§ CONCLUSION

We have presented a nested coordinate BA that solves the two-species PASEP problem with overtaking, in both the deformed and simple (undeformed) cases, provided the condition pqpq is satisfied.

The YBE were derived, which allowed us to assess the BA integrability of all models in this class.
Although the derivation presented here is not entirely general, as a particular form of the nested coordinate Bethe wave function ansatz was used, it is conjectured that the list of integrable models presented is exhaustive.

It was shown that the solutions fall into three cases,
namely a TASEP with two species that move at different rates, with one species overtaking the other with an arbitrary rate; a PASEP with two species moving at the same rate and one species treating the other as holes; and a PASEP with a single defect particle that moves with a different rate and does not distinguish between particles and holes.
These results are summarized in tab:model_list.
Two of these cases are well-known in literature.
Namely, the TASEP with defect particles, and the second-class PASEP.

Solution 3 is a PASEP with a single first-class defect particle and overtaking.
This model has only been investigated recently, when its steady state was solved using the matrix product approach [35].
It is peculiar among the solutions listed here in that it is only integrable for a single defect particle.
This follows from the fact that the YBE are satisfied for three-body collisions in which the particle species are $122$ and $222$ but not $112$.

Model     Parameter restrictions     Free parameters

TASEP with defect species     $q_1=q_2=b_{21}=0$     $p_1,\; b_{12}\;, M_1,\; M_2$

Second-class PASEP     $p_1=b_{12}=p_2,\; p_2=b_{21}=q_2$     $q_2,\; M_1,\; M_2$

PASEP with one defect     $p_1 = b_{12} = \alpha p_2 ,\; q_1=b_{21}=q_2/\alpha,\; M_1=1$
$q_2,\; \alpha,\; M_2$

All Bethe ansatz integrable sub-classes of the two-species partially asymmetric exclusion process with overtaking on a ring.
The rate $p_2$ is not included as a free parameter as it can be seen as just setting the timescale.
The particle numbers $M_1,M_2$ are included as free parameters to highlight that $M_1$ is fixed in the third case.

It is interesting to note that that among the models examined here, all the models that are integrable in the normal (undeformed) case remain integrable after a deformation of the Markov operator.
This is a priori not obvious and it would be of interest to investigate whether integrability is preserved under deformation in other contexts as well.

Finally, we derived the Bethe equations for the deformed time evolution problem of solution 3 of the YBE.
This should allow one to calculate the long-time current fluctuations in this model.
This would complement the result for the mean current, which has already been calculated using a matrix product ansatz [35].
The Bethe equations have similar features to the well-known one-species PASEP and two-species TASEP, which makes it amenable to the techniques developed for those cases.

Another possible direction for future work would be to consider some extensions of the traditional YBE, which have been proposed in the literature.
These include the dynamical YBE [45] and the braided YBE [46, 47].
It would be interesting to investigate whether those generalizations can provide further insight in the context of asymmetric exclusion processes.

§ ACKNOWLEDGEMENTS

IL acknowledges studentship funding from EPSRC under Grant No. EP/R513209/1.
The work of KM has been supported by the project RETENU ANR-20-CE40-0005-01 of the French National Research Agency (ANR).
For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

§ DERIVATION OF SOLUTIONS OF YANG-BAXTER EQUATIONS

We wish to obtain conditions on the model parameters $\alpha ,x,\beta _{12},\beta_{21}$ such that yb are satisfied for the 112 and 122 blocks identically for any values of $z_i,z_j,z_k$.
As explained previously, we can do this by looking at individual coefficients of $z_i^{n_i}z_j^{n_j}z_k^{n_k}$ for any powers $n_i,n_j,n_k$ and requiring them to vanish.
We perform the following calculations with the help of Mathematica to evaluate the needed expressions symbolically.
We show the expressions for the deformed case, as the undeformed limit can be obtained simply by setting $\gamma_1=\gamma_2=\gamma_{12}=0$.
The procedure presented here is shown schematically in fig:tree for clarity.

First, we note that all entries have a factor of $(z_i-z_j)(z_i-z_k)(z_j-z_k)$, so we divide by this for simplicity.
This will be assumed to have been done in all subsequent expressions.
Then, looking at the constant term (i.e. taking $z_i\to0,z_j\to0,z_k\to0$), in the (2,1) entry in the 112 block and the (1,2) entry in the 122 block, we get respectively,

\begin{equation}
\rme^{\gamma_2+\gamma_{12}}\alpha^4\beta_{12}(\beta_{12}\beta_{21}-x)=0\;,
\qquad \rme^{\gamma_1}\alpha^4\beta_{21}(\beta_{12}\beta_{21}-x)=0 \;.
\end{equation}

We cannot set $\alpha=0$, as, due to r1r2, this would make the ansatz vanish.
Hence we can either set $\beta_{21}=x/\beta_{12}$ or $\beta_{12}=\beta_{21}=0$.

§.§ Case 1: $\beta_{21}=x/\beta_{12}$

With this choice, we now look at the constant term of the numerator in the (2,1) entry in the 122 block.
This gives us the additional condition

\begin{equation}
\rme^{\gamma_2+\gamma_{12}}\alpha x(\alpha-1)(\alpha-\beta_{12})(1+\beta_{12})=0\;.
\end{equation}

As before, we cannot set $\alpha=0$.
This leaves us with three options: $x=0$ or $\alpha=1$ or $\beta_{12}=\alpha$.
The choice $x=0$ is sufficient to satisfy the YBE fully for the 112 and 122 blocks, giving solution 1.

With the choice $\alpha=1$, we now look at the coefficient of $z_i^0z_j^1z_k^1$ in the $(1,1)$ entry of the $122$ block.
This gives the additional condition

\begin{eqnarray}
\end{eqnarray}

Setting $x=0$ gives solution 1.
Setting $\beta_{12}=1$ or $x$ is sufficient to satisfy the YBE for the $112$ and $122$ block.
Both options give solution 2 up to a relabelling of species.

The choice $\beta_{12}=\alpha$ is sufficient to satisfy the YBE for the 122 block.
The 112 block is not satisfied.
We can either limit ourselves to have only one particle of species 1, which gives solution 3, or try to satisfy the 112 block as well.
In the 112 block, looking at the constant term in the numerator of the (1,2) entry gives the condition

\begin{equation}
\rme^{\gamma_1}\alpha^3x(\alpha-1)(\alpha^2+x)(\alpha-x)=0\;.
\end{equation}

If we wish to satisfy this, we must set either $x=0$ or $\alpha=1$ or $\alpha=x$.
All of those choices lead either to solution 1 or 2.

[style = draw,->,thick]
[draw](2,1) in block 112, (1,2) in block 122[sibling distance = 9cm,level distance = 2.5cm]
child node [draw](2,1) in block 122[sibling distance = 3.5cm]
child node [draw,ultra thick]S 1 edge from parent node [left,xshift=-0.3cm]$x=0$
child node [draw](1,1) in block 112 [sibling distance = 2cm]
child node [draw, ultra thick]S2 edge from parent node [left]$\beta_{12}=1$ or $x$
child node [draw, ultra thick]S1 edge from parent node [right,xshift=0.1cm,yshift=-0.2cm]$x=0$
edge from parent node[left,xshift=0.1cm]$\alpha=1$
child node [draw]block 122 solved[sibling distance = 2.5cm]
child node [draw,ultra thick]S 3
edge from parent node [left,xshift=0.1cm,yshift=0.2cm]$M_1=1$
child node [draw](1,2) in block 112
child node [draw,ultra thick]S 1 edge from parent node [left,xshift=-0.3cm]$x=0$
child node [draw,ultra thick]S 2 edge from parent node [right,xshift=0.3cm]$\alpha=1$ or $x$
edge from parent node [right,xshift=0.3cm]solve block 112
edge from parent node [right,xshift=0.3cm]$\beta_{12}=\alpha$
edge from parent node [left,xshift=-0.5cm]$\beta_{21}=x/\beta_{12}$
child node [draw](1,1) in block 122[sibling distance = 4cm]
child node [draw,ultra thick]1-species PASEP edge from parent node [left,xshift=-0.3cm]$\alpha=1$
child node [draw,ultra thick]S 1 (with $\beta_{12}=0$) edge from parent node [right,xshift=0.3cm]$x=0$
edge from parent node [right,xshift=0.5cm]$\beta_{12}=\beta_{21}=0$;

Representation of the procedure to obtain solutions to the Yang-Baxter equations described in <ref>.
Solutions are outlined with thick boundaries.
S stands for solution.

§.§ Case 2: $\beta_{12}=\beta_{21}=0$

With this choice, we now look at the constant term in the numerator of the (1,1) entry of the 122 block.
This gives the further condition

\begin{equation}
\alpha^3x(\alpha-1)=0\;.
\end{equation}

Setting $\alpha=1$ is enough to satisfy all YBE but it gives a essentially a one-species PASEP, which is not of interest here, as we are looking for two-species solutions.

Choosing $x=0$ gives a special case of solution 1 with $\beta_{12}=0$.

§ ALGEBRAIC STEPS IN DERIVATION OF BETHE EQUATIONS FOR SOLUTION 3

Requiring the system of equations beq_0a-beq_0b to be degenerate, we get

\begin{eqnarray}
\fl \left[ {\cal D}_{ij}-xz_j + x z_i S^{22}_{22}(z_i,z_j)E_j \right]
\left[ ({\cal D}_{ij}-z_j) S^{22}_{22}(z_j,z_i)E_i + z_i \right]
- (i\leftrightarrow j)=0 \;.
\end{eqnarray}

Expanding and simplifying this, we get

\begin{eqnarray}
\fl ({\cal D}_{ij}-(1+x)z_j) S^{22}_{22}(z_j,z_i)E_i
- ({\cal D}_{ij}-(1+x)z_i) S^{22}_{22}(z_i,z_j)E_j \nonumber \\
+ (z_i-z_j)(1+xE_iE_j) = 0 \;.
\end{eqnarray}

Now plugging in the explicit forms of $S^{22}_{22}$, this becomes

\begin{equation}
\fl -({\cal D}_{ij}-(1+x)z_i) E_i + ({\cal D}_{ij}-(1+x)z_j) E_j
+ (z_i-z_j)(1+xE_iE_j) = 0 \;.
\end{equation}

With the change of variable ztoy, we get, after some simplification,

\begin{equation}
(y_i-xy_j) E_i -(y_j-xy_i) E_j
+ (y_i-y_j)(1+xE_iE_j) = 0 \;,
\end{equation}

where now

\begin{equation}
E_j = \rme^{{\bar\gamma} +L\gamma_2}
\left(\frac{1-y_j}{1-xy_j}\right) ^L
\prod\limits _{k=1}^{M} \frac{xy_j-y_k}{y_j-xy_k} \;.
\end{equation}

The terms with index $i$ and $j$ can now be separated to give the following equation,

\begin{equation}
y_i \frac{1+E_i}{1+xE_i} = y_j \frac{1+E_j}{1+xE_j} \;.
\end{equation}

As this procedure can be carried out for any two indices $i,j$, this means that this expression must equal some constant $B$ for all $i,j$.
Then equating the left-hand side with $B$ and rearranging to isolate $E_i$, we finally obtain beq_1.

§ REFERENCES

[1]
Derrida B 2007 Journal of Statistical Mechanics: Theory and Experiment
2007 P07023

[2]
Blythe R A and Evans M R 2007 Journal of Physics A: Mathematical and
Theoretical 40 R333

[3]
Chou T, Mallick K and Zia R K 2011 Reports on progress in physics 74 116601

[4]
MacDonald C T, Gibbs J H and Pipkin A C 1968 Biopolymers: Original
Research on Biomolecules 6 1–25

[5]
Szavits-Nossan J, Ciandrini L and Romano M C 2018 Phys. Rev. Lett. 120(12) 128101

[6]
Scott S and Szavits-Nossan J 2019 Physical biology 17 015004

[7]
Wolf D E, Schreckenberg M and Bachem A 1996 Traffic and granular flow
(World Scientific)

[8]
Cividini J, Mukamel D and Posch H A 2017 Physical Review E 95

[9]
Derrida B, Evans M R, Hakim V and Pasquier V 1993 Journal of Physics A:
Mathematical and General 26 1493

[10]
Alexander S and Holstein T 1978 Phys. Rev. B 18(1) 301–302

[11]
Dhar D 1987 Phase transitions 1 51

[12]
Gwa L H and Spohn H 1992 Physical Review A 46 844

[13]
Golinelli O and Mallick K 2006 Journal of Physics A: Mathematical and
General 39 12679

[14]
Bethe H 1931 Zeitschrift für Physik 71 205–226

[15]
Baxter R J 2016 Exactly solved models in statistical mechanics

[16]
Schütz G M 1997 Journal of statistical physics 88 427–445

[17]
Sasamoto T and Wadati M 1998 Physical Review E 58 4181

[18]
De Gier J and Essler F H 2005 Physical review letters 95 240601

[19]
Simon D 2009 Journal of Statistical Mechanics: Theory and Experiment
2009 P07017

[20]
Kim D 1995 Phys. Rev. E 52(4) 3512–3524

[21]
Alcaraz F and Bariev R 1999 Physical Review E 60 79

[22]
Alcaraz F and Bariev R 2000 Brazilian Journal of Physics 30

[23]
de Gier J and Essler F H 2008 Journal of Physics A: Mathematical and
Theoretical 41 485002

[24]
Prolhac S and Mallick K 2008 Journal of Physics A: Mathematical and
Theoretical 41 175002

[25]
Prolhac S 2010 Journal of Physics A: Mathematical and Theoretical 43 105002

[26]
Derrida B and Evans M R 1999 Journal of Physics A: Mathematical and
General 32 4833

[27]
Cantini L 2008 Journal of Physics A: Mathematical and Theoretical 41 095001

[28]
Derrida B and Lebowitz J L 1998 Physical review letters 80 209

[29]
Jimbo M 1990 Yang-Baxter equation in integrable systems vol 10 (World

[30]
Kulish P and Sklyanin E 1982 Journal of Soviet Mathematics 19

[31]
Vieira R S and Lima-Santos A 2021 Journal of Statistical Mechanics: Theory
and Experiment 2021 053103

[32]
Schütz G M 2001 Exactly Solvable Models for Many-Body Systems Far from
Equilibrium, in Phase Transitions and Critical Phenomena, C. Domb and J. L.
Lebowitz (Eds) vol 19 (Academic Press)

[33]
Crampé N, Mallick K, Ragoucy E and Vanicat M 2015 J. Phys. A:
Mathematical and Theoretical 48 175002

[34]
Crampé N, Ragoucy E and Vanicat M 2014 Journal of Statistical
Mechanics: Theory and Experiment 2014 P11032

[35]
Lobaskin I, Evans M R and Mallick K 2022 Journal of Physics A:
Mathematical and Theoretical 55 205002

[36]
Derrida B, Janowsky S A, Lebowitz J L and Speer E R 1993 Journal of
statistical physics 73 813–842

[37]
Mallick K 1996 Journal of Physics A: Mathematical and General 29

[38]
Arndt P F, Heinzel T and Rittenberg V 1998 Journal of Physics A:
Mathematical and General 31 L45

[39]
Rajewsky N, Sasamoto T and Speer E 2000 Physica A: Statistical Mechanics
and its Applications 279 123–142

[40]
Cantini L and Sasamoto T 2022 Personal communication

[41]
Touchette H 2009 Physics Reports 478 1–69 ISSN 0370-1573

[42]
Chetrite R and Touchette H 2015 Ann. Henri Poincare 16 2005

[43]
Simon D 2009 J. Stat. Mech: Theor and Exp. 2009 P07017

[44]
Popkov V, Schütz G M and Simon D 2010 J. Stat. Mech: Theor and Exp.
2010 P10007

[45]
Avan J, Babelon O and Billey E 1996 Communications in mathematical
physics 178 281–299

[46]
Fioravanti D and Rossi M 2001 Journal of Physics A: Mathematical and
General 34 L567

[47]
Fioravanti D and Rossi M 2002 Journal of Physics A: Mathematical and
General 35 3647


#  Universal scaling near band-tuned metal-insulator phase transitions

Simone Fratini  Université Grenoble Alpes, CNRS, Grenoble INP, Institut Néel,
38000 Grenoble, France  Sergio Ciuchi  Dipartimento di Scienze Fisiche e
Chimiche, Università dell’Aquila, 67100 Coppito (AQ), Italy  Istituto dei
Sistemi Complessi, CNR, P.le Aldo Moro I-00185 Roma, Italy  Vladimir
Dobrosavljević  Department of Physics and National High Magnetic Field
Laboratory, Florida State University, Tallahassee, FL, USA  Louk Rademaker
Department of Quantum Matter Physics, University of Geneva, 1211 Geneva,
Switzerland

###### Abstract

We present a theory for band-tuned metal-insulator transitions based on the
Kubo formalism. Such a transition exhibits scaling of the resistivity curves,
in the regime where $T\tau>1$ or $\mu\tau>1$, where $\tau$ is the scattering
time and $\mu$ the chemical potential. At the critical value of the chemical
potential, the resistivity diverges as a power law, $R_{c}\sim 1/T$.
Consequently, on the metallic side there is a regime with negative $dR/dT$,
which is often misinterpreted as insulating. We show that scaling and this
‘fake insulator’ regime is observed in a wide range of experimental systems.
In particular, we show that Mooij correlations in high-temperature metals with
negative $dR/dT$ can be quantitatively understood with our scaling theory in
the presence of $T$-linear scattering.

Thanks to the advent of highly tunable ‘twisted’ Van der Waals
heterostructures,[1, 2, 3] the field of quantum matter physics is in a
position to study continuous zero-temperature phase transitions with an
unprecedented accuracy. Detailed (and smooth!) experimental results allow a
systematic comparison between different theoretical predictions, which is
particularly true for continuous metal-to-insulator transitions (MITs).

Interaction-induced MITs, such as the Mott transitions, display quantum
critical behavior, including scaling of the resistivity.[4, 5] A full
theoretical understanding of Mott criticality, which would include a precise
calculation of the scaling exponents, is still lacking.[6] One of the main
challenges lies in the fact that an MIT is, in general, not a transition
described by symmetry breaking, which makes it challenging to identify the
source of scaling.

Recently, scaling has been observed in a simple band-tuned MIT in a MoTe2/WSe2
bilayer at full filling of the first valence flat band.[7] By tuning the
displacement field, one can open a band gap to the second valence band. The
scaling behavior there has been analysed using a model with disorder and a
bosonic field,[8] inspired by earlier work on ‘Mooij’ correlations.[9, 10]
However, the observed scaling can also be interpreted in a much simpler
perspective.

From a theoretical viewpoint, calculating the conductivity is notoriously
difficult. An exception is the classical Drude formula,
$\sigma=\frac{ne^{2}\tau}{m}$, which can also be derived with fully quantum-
mechanical advanced methods such as the Kubo formula,[11, 12]. A natural
question is whether the observed scaling at a metal-insulator transition can
be explained with the same set of assumptions that is used to derive Drude
theory.

Indeed, in this Letter we show that only a small number of very natural
assumptions leads to scaling behavior near a band-tuned MIT. The only
assumptions are that the scattering time $\tau$ is large, parametrized by
$T\tau>1$ or $\mu\tau>1$ respectively on the insulating and metallic sides of
the transition (with $\mu$ the chemical potential measured from the band
edges), and that the electron self-energy is local and proportional to the
electron density of states. These conditions naturally arise in weakly
correlated, weakly disordered metals. With this, the critical resistivity at
the MIT is diverging as $R_{c}(T)\sim 1/T$, in contrast to oft-cited picture
that the critical resistivity curve is independent of temperature. We derive
an explicit scaling form, showing that in the scaling regime the resistivity
is given by a universal $R(T,\mu)=R_{c}(T)f(\mu/T)$. Contrary to the physics
of universality at continuous phase transitions, the scaling of the
resistivity breaks down very close to the MIT.

Figure 1: In a band-tuned metal-insulator transition (MIT), the system changes
from having overlapping valence (blue) and conduction (red) bands in the
metallic side (right) to having a gap on the insulating side (left). The
tuning parameter is the chemical potential $\mu$. When either $T\tau>1$ or
$\mu\tau>1$, the resistivity (in shades on the background) can be described by
a scaling form, as shown in Fig. 2. This scaling relation breaks down very
close to the transition, where localization and interaction effects will
change the picture. Figure 2: Theoretical resistance curves close to a band-
tuned metal-insulator transition. a. Resistance calculated using Eq. (3), for
a constant scattering time $\tau=25$ eV-1, and chemical potential $\mu$
ranging from $-0.8$ to $+0.8$ eV. The resistance at the critical point $\mu=0$
diverges as $R_{c}(T)\sim 1/T$. On the metallic side, the resistance decreases
as a function of temperature (a ‘fake’ insulator), whereas on the insulating
side the resistance is activated. b. Resistance for a temperature-dependent
scattering rate $\tau^{-1}=\tau_{0}^{-1}+bT$ with $\tau_{0}=25$ eV-1 and
$b=0.1$. On the metallic side, there is a transition from a positive $dR/dT$
to a negative $dR/dT$. At high temperatures, this gives rise to Mooij
correlations (see Fig. 4). c. When $T\tau>1$ or $\mu\tau>1$, the resistance
curves follow a simple scaling law $R(T,\mu)=R_{c}(T)f(\mu/T)$. This can be
verified by plotting $R/R_{c}$ versus $T/|\mu|$. All data points collapse onto
one of the two curves, associated with either metallic or insulating behavior.

Band-tuned MIT –  Consider a weakly interacting electron system described by a
band-structure. The system is metallic if there is a nonzero density of charge
carriers, characterized by a nonzero chemical potential $\mu$. The system is
an insulator if there is a gap towards exciting charge carriers. By
continuously changing the bandstructure we can induce a band-tuned MIT. This
can be achieved with pressure, displacement field, or even due to spontaneous
symmetry breaking such as ferromagnetic polarization. Without loss of
generality, the dispersion at a band edge is parabolic, with the dispersion
set by $\xi_{\bf k}=\frac{k^{2}}{2m}-\mu$ where $m$ is the effective mass.
With this notation, $\mu>0$ corresponds to the metal, $\mu<0$ to an insulator,
and $\mu=0$ is the critical point. The chemical potential $\mu$ is thus the
tuning parameter of the MIT, as shown in Fig. 1.

In general, the conductivity is determined by disorder, electron-electron
interactions and electron-phonon coupling. Nonzero resistivity from electron-
electron interactions requires Umklapp scattering, which becomes
asymptotically irrelevant at low carrier densities (though there might be
nontrivial vertex corrections)[13]. Similarly, at zero temperature there is no
thermal occupation of phonons, and therefore no electron-phonon contribution
to the resistivity. The zero-temperature behavior of a band-tuned MIT is
therefore completely dominated by disorder. In principle strong disorder might
push the system into Anderson insulation. However, in $d=2,3$ it is considered
that the combination of weak disorder and weak interactions generally
precludes true localization [14, 15, 16, 17]. Moreover, even in the absence of
interactions, quantum corrections to the conductivity are not relevant in the
regimes $\mu\tau>1$ and $T\tau>1$ considered here, and will therefore be
neglected throughout this work.

Figure 3:  Scaling near the band-tuned MIT is observed in a range of
materials. Here, we apply our scaling analysis to three material systems[18]:
(a) the moiré heterobilayer MoTe2/WSe2[5], (b) the heterostructure
WSe2/bilayer graphene/WSe2[19], and (c) GST amorphous phase change materials
[20]. The measured resistivities are shown in insets. In panels a and c we see
a genuine MIT, with data collapse on both an insulating and conducting branch.
The theoretical scaling curve of Eq. (4) is shown as a dashed black line, and
shows remarkable agreement with the experimental results.

Conductivity –  With these natural assumptions, the conductivity close to the
MIT is calculated using the Kubo formula for local self-energies [11, 21],
which reads for a single band in $d$ dimensions, per spin species,

$\sigma_{xx}=\pi\int d\xi\Phi_{x}(\xi)\int
dzA^{2}(\xi,z)\left(-f^{\prime}(z)\right),$ (1)

where $A(\xi,z)$ is the one-particle spectral function and $f$ is the Fermi
function (see [18], Sec. A). The entire momentum-dependence is included in a
transport function $\Phi_{x}(\xi)=\int\frac{d^{d}p}{(2\pi)^{d}}j^{2}_{x}({\bf
p})\,\delta(\xi-\xi_{\bf p})$. The transport function itself displays
universal behavior in the vicinity of a band-tuned MIT: given the parabolic
band dispersion, the current operator equals ${\bf j}({\bf p})=\frac{e}{m}{\bf
p}$. Consequently the transport function reads

$\Phi_{x}(\xi)=\frac{2e^{2}}{dm}\left(\xi+\mu\right)N(\xi),$ (2)

where $N(\xi)$ is the non-interacting density of states,
$N(\xi)=\int\frac{d^{d}p}{(2\pi)^{d}}\delta(\xi-\xi_{\bf p})$. Assuming a
constant, energy-independent scattering rate $\tau$ in $d=2$, the imaginary
part of the self-energy is $\mathrm{Im}\Sigma(z)=-\Theta(z+\mu)(2\tau)^{-1}$.
This scattering time is typically of the order $\tau\sim 10^{-12}-10^{-14}$ s
$\sim 10-10^{3}$ eV-1. When $\mu\tau>1$ or $T\tau>1$, the Kubo formula
radically simplifies, and we find

$\sigma(T,\mu)=\frac{e^{2}}{h}\tau T\log\left[1+e^{\mu/T}\right].$ (3)

This is our central result for the conductivity close to the band-tuned MIT.
Surprisingly, it contrasts a few commonly held convictions on metal-insulator
transitions. First, at the critical point, the conductivity is linear in
temperature, $\sigma_{c}(T)\equiv\sigma(0,T)=\frac{e^{2}}{h}\tau T\log 2$,
rather than temperature-independent. Furthermore, on the metallic side of the
transition $(\mu>0)$, the temperature derivative of the resistivity can be
negative: a ‘fake insulator’ regime that is commonly misinterpreted as
insulating. Furthermore, Eq. (3) satisfies a universal scaling form

$\displaystyle\sigma(\mu,T)$ $\displaystyle=\sigma_{c}(T)F(\mu/T),$ (4)

which allows the collapse of many resistivity curves onto a simple scaling
function $F(x)=\log_{2}\left[1+e^{x}\right]$. The theoretical resistance
curves near the band-tuned MIT, including the scaling properties, are shown in
Fig. 2.

Hidden in plain view is the fact that Eq. (3) is, at zero temperature on the
metallic side, equivalent to Drude theory. Explicitly, in $d=2$,
$\sigma(\mu)=\frac{e^{2}}{h}\tau\mu=\frac{ne^{2}\tau}{m}$. In fact, the $T=0$
limit of Eq. (7) yields $\sigma=\Phi(E_{F})\tau$ with $\Phi(E_{F})=ne^{2}/m$
in any dimension $d$.

At finite temperature the scaling regime persists, even with a temperature-
dependent scattering time $\tau(T)$, provided that $\tau^{-1}$ is still
proportional to the density of states. When $\tau$ is temperature-independent,
in fact, all resistivity curves on the metallic side are ‘fake insulators’
with $d\rho/dT\leq 0$ (cf. Fig. 2). Only when the scattering rate increases
with temperature, for example from electron-phonon interactions shown in Fig.
2b or from Umklapp scattering, we find traditional metallic behavior with
$d\rho/dT>0$. In this case, inside the metallic regime there exists a point
where the temperature-derivative of the resistivity $d\rho/dT$ changes sign.
We will discuss universal properties around this point later in the context of
Mooij correlations.[9, 22]

It is important to emphasize that the scaling form of Eq. (4) is limited to
regions not too close to the transition. This limitation is similar to the one
proposed by Mott-Ioffe-Regel (MIR) [23]. A common formulation of the MIR limit
in metals is $k_{F}\ell\sim 1$ where $\ell$ is the mean-free path. This can be
rewritten as $\mu\tau\sim 1$; we therefore find that, upon approaching the
transition from the metallic side, the scaling hypothesis breaks down
precisely at the MIR boundary. What happens close to the transition is non-
universal, and depending on model parameters one can find various different
violations of scaling (see [18], Sec. B).

Band-tuned MIT in moiré bilayers –  We are now in a position to verify our
universal scaling result of Eq. (3) in experimental results on real physical
systems. Inspired by the recent developments in moiré materials, let us first
focus on the MIT in MoTe2/WSe2 at full filling of the first valence flat band
($f=2$).[7] By tuning the perpendicular displacement field, a gap is opened
up, yielding a band-tuned MIT. In Fig. 3 we fit the observed resistance curves
as a function of displacement field using our theory. Indeed, the critical
resistance diverges as $R_{c}\sim 1/T$, and the resistance curves obey
scaling. As shown in Fig. 3a, the scaling curve itself quantitatively matches
the analytical form derived in Eq. (3). A similar scaling plot for these data
has been reported in Ref. [8], inspired by earlier work in Ref. [10], which
describes disorder-induced polaron formation when the chemical potential is
far from the band edges. When the chemical potential approaches the band
edges, the theory of Ref. [8] reduces to the simpler theory presented here,
where polaronic effects are irrelevant.

There are many claims of MITs in graphene-based Moiré materials, that upon
closer inspection seem to exhibit ”fake insulator” behavior. Consider, for
example, the WSe2/bilayer graphene (BLG)/WSe2 heterostructure measured in Ref.
[19]. At filling $\nu=0$, the resistivity turns up at low temperatures
reminiscent of an insulating gap. However, at around $T=20K$, the resistivity
seems to saturate, to a displacement-field dependent value. The absence of a
true diverging resistance at low temperature suggests that these systems
retain a nonzero density of charge carriers, either from a band overlap or
induced by potential inhomogeneities that are common in graphene systems.
Indeed, when performing the scaling analysis, we can collapse all the curves
of this system to the metallic branch of our scaling form, as shown in Fig.
3b.

Disordered metallic alloys –  While Eq. (3) was derived for $d=2$ and weak
disorder scattering, it is in fact far more universal. Often a momentum-
independent self-energy arises through the equation $\Sigma(z)=s^{2}G(z)$, for
example in iterated schemes such as the self-consistent Born approximation for
disorder scattering or electron-phonon scattering in the adiabatic limit. Here
$s$ is a (possibly temperature dependent) parameter quantifying the scattering
process. Under this scheme, the inverse scattering time is in weak coupling
proportional to the density of states ${\rm Im}\Sigma(z)\propto
N(z)\sim\mathrm{Im}G(z)$. This leads to a conductivity of the form
$\sigma(T)=\frac{e^{2}}{dms^{2}}T\log\left[1+e^{\mu/T}\right]$ in general
dimensions $d$, consistent with Eq. (3).

The universal scaling is indeed also observed in three-dimensional compounds
away from the weak disorder limit. In particular, we look at GST[20], a phase-
change compound where the annealing history affects the effective number of
charge carriers.[24] Here, at high temperatures, a smooth evolution from
positive $dR/dT$ to negative $dR/dT$ is observed depending on the precise
composition and history of the sample. Since the main effect of these
compositional changes is in fact a shift of the chemical potential, we show in
Fig. 3c that the experimental data on GST can be accurately described by our
scaling theory.

Figure 4: The dimensionless temperature coefficient of the resistance versus
the dimensionless resistance, for a variety of materials[9, 5, 25, 26, 27, 28,
29, 30, 20] (for details see [18]), compared to our theoretical result of Eq.
(6) (solid black line). For the experimental data the only fitting parameter
is $R^{\infty}$, the limit of the critical resistance at high temperature. We
find an excellent agreement of the experimental Mooij correlations and our
theory.

Mooij correlations –  Universal scaling implies the existence of a ‘fake
insulator’ regime: a metal characterized by a (dimensionless) negative
temperature coefficient of the resistance $\alpha=(T/R)dR/dT<0$. Historically,
the observation of a negative $\alpha$ various disordered metals, including
binary alloys (NixCr1-x, TixAl1-x, FexSi1-x, etc.)[9, 10] was considered a
‘high temperature anomaly’.[22] In a seminal paper, Mooij[9] discovered a
correlation between the temperature coefficient $\alpha$ and the resistivity
$\rho$ itself. There is currently no consensus on the origin of these Mooij
correlations, though they have been interpreted in terms of quantum
localization corrections to the conductivity [31, 22] or the disorder-driven
formation of polarons[10].

Interestingly, the scaling theory proposed in this Letter allows to
quantitatively describe Mooij correlations. To do so, we assume that at high
temperature the scattering time $\tau$ is linear in $T$:

$\tau^{-1}=\tau_{0}^{-1}+bT.$ (5)

This form occurs in many metals, where $b$ is either proportional to the
electron-phonon coupling strength, or a more complex, ”Planckian” quantum
scattering.[32] With this assumption, the critical curve becomes flat at high
temperature, $R_{c}(T)\rightarrow R^{\infty}\propto b$. This allows us to
introduce a dimensionless resistivity $R/R^{\infty}$. By taking the derivative
of the scaling relation Eq. (3), and inverting it with respect to the tuning
parameter $\mu$ at a fixed temperature $T$, we find that the temperature
coefficient $\alpha$ only depends on $R/R^{\infty}$,

$\alpha(R)=\frac{R}{R^{\infty}}\left(1-2^{R^{\infty}/R}\right)\log_{2}\left[2^{R^{\infty}/R}-1\right].$
(6)

In Fig. 4 we compare our analytical result with the original data presented by
Mooij[9] and those collected in [10], finding a good agreement between the
experimental results on binary alloys and Eq. (6). The recent data on Moiré
bilayers shows an even more striking quantitative equivalence between the
resistivity data in the high temperature range $T=26-60$ K: without a fitting
parameter the experimental results of Ref. [5] match Eq. (6).

Outlook – In this Letter we have shown that a simple theory of conductivity
predicts universal scaling near band-tuned MITs consistent with experimental
results in a wide range of materials, from recent Moiré materials to decades
old data on binary alloys.

The predicted scaling regime does not extend arbitrarily close to the MIT:
when $T\tau<1$ and $\mu\tau<1$ deviations from or a full breakdown of scaling
can appear. Note that the difference between scaling close and further away
from the transition has been discussed in Ref. [6]. The scaling described in
this Letter is thus not due to the divergence of a length scale, and is not
related to Landau order parameters, the renormalization group, or any other
theory of universality in symmetry-breaking (quantum) phase transitions. The
universal behavior of resistivity scaling near the MIT throughout many
materials is just the consequence of a generic weakly interacting electrons
with weak disorder, in spirit similar to the stability of the Fermi liquid.
The properties of Anderson and weak localization as well as Wigner
crystallization and the Mott MIT[6] are phenomena that, on the other hand, are
outside the scaling regime discussed here. It is an interesting open question
whether the scaling described in this Letter can extend, under certain
conditions, arbitrarily close to the MIT, thus connecting to the standard
theoretical framework of continuous phase transitions.[33]

Acknowledgements –  We thank Jie Shan, Kin Fai Mak, Péter Makk and Bálint
Szentpéteri for sharing their experimental data. We thank Yuting Tan,
Christophe Berthod, and Giacomo Morpurgo for fruitful discussions. LR is
funded by the Swiss National Science Foundation by Starting Grant
TMSGI2_211296. SC is funded by the European Union - NextGenerationEU under the
Italian Ministry of University and Research (MUR) National Innovation
Ecosystem grant ECS00000041 - VITALITY - CUP E13C22001060006.

## References

  * [1] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young,  _Superconductivity and strong correlations in moiré flat bands_, Nature Physics 16, 725 (2020).
  * [2] D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. N. Pasupathy, and A. Rubio,  _Moiré heterostructures as a condensed-matter quantum simulator_, Nature Physics 17, 155 (2021).
  * [3] K. F. Mak and J. Shan,  _Semiconductor moiré materials_, Nature Nanotechnology 17, 686 (2022).
  * [4] A. Ghiotto, E.-M. Shih, G. S. S. G. Pereira, D. A. Rhodes, B. Kim, J. Zang, A. J. Millis, K. Watanabe, T. Taniguchi, J. C. Hone, L. Wang, C. R. Dean, and A. N. Pasupathy,  _Quantum Criticality in Twisted Transition Metal Dichalcogenides_, Nature 597, 345 (2021).
  * [5] T. Li, S. Jiang, L. Li, Y. Zhang, K. Kang, J. Zhu, K. Watanabe, T. Taniguchi, D. Chowdhury, L. Fu, J. Shan, and K. F. Mak,  _Continuous Mott transition in semiconductor moiré superlattices_, Nature 597, 350 (2021).
  * [6] Y. Tan, V. Dobrosavljevic, and L. Rademaker,  _How to Recognize the Universal Aspects of Mott Criticality?_, Crystals 12, 932 (2022).
  * [7] T. Li, S. Jiang, B. Shen, Y. Zhang, L. Li, Z. Tao, T. Devakul, K. Watanabe, T. Taniguchi, L. Fu, J. Shan, and K. F. Mak,  _Quantum anomalous Hall effect from intertwined moiré bands_, Nature 600, 641 (2021).
  * [8] Y. Tan, P. K. H. Tsang, and V. Dobrosavljević,  _Disorder-dominated quantum criticality in moiré bilayers_, Nature Communications 13, 7469 (2022).
  * [9] J. H. Mooij,  _Electrical conduction in concentrated disordered transition metal alloys_, physica status solidi (a) 17, 521 (1973).
  * [10] S. Ciuchi, D. D. Sante, V. Dobrosavljević, and S. Fratini,  _The origin of Mooij correlations in disordered metals_, npj Quantum Materials 3, 1 (2018).
  * [11] G. D. Mahan,  _Many-Particle Physics_ , 3rd ed., Kluwer Academic, New York (Kluwer Academic, New York, 2000).
  * [12] P. Coleman,  _Introduction to Many-Body Physics_ (Cambridge University Press, 2015) ISBN 9780521864886\.
  * [13] A. Mu, Z. Sun, and A. J. Millis,  _Optical conductivity of the two-dimensional Hubbard model: Vertex corrections, emergent Galilean invariance, and the accuracy of the single-site dynamical mean field approximation_, Physical Review B 106, 085142 (2022).
  * [14] F. Evers and A. D. Mirlin,  _Anderson transitions_, Reviews of Modern Physics 80, 1355 (2008).
  * [15] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn,  _Colloquium: Many-body localization, thermalization, and entanglement_, Reviews of Modern Physics 91, 021001 (2019).
  * [16] A. Punnoose and A. M. Finkel’stein,  _Metal-Insulator Transition in Disordered Two-Dimensional Electron Systems_, Science 310, 289 (2005).
  * [17] C. Castellani, C. D. Castro, P. A. Lee, and M. Ma,  _Interaction-driven metal-insulator transitions in disordered fermion systems_, Physical Review B 30, 527 (1984).
  * [18] _See online Supplementary Materials_.
  * [19] M. Kedves, B. Szentpéteri, A. Márffy, E. Tóvári, N. Papadopoulos, P. K. Rout, K. Watanabe, T. Taniguchi, S. Goswami, S. Csonka, and P. Makk,  _Stabilizing the inverted phase of a WSe$_2$/BLG/WSe$_2$ heterostructure via hydrostatic pressure_, arXiv (2023).
  * [20] T. Siegrist, P. Jost, H. Volker, M. Woda, P. Merkelbach, C. Schlockermann, and M. Wuttig,  _Disorder-induced localization in crystalline phase-change materials_, Nature Materials 10, 202 (2011).
  * [21] 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_, Reviews of Modern Physics 68, 13 (1996).
  * [22] P. A. Lee and T. V. Ramakrishnan,  _Disordered electronic systems_, Reviews of Modern Physics 57, 287 (1985).
  * [23] N. E. Hussey, K. Takenaka, and H. Takagi,  _Universality of the Mott–Ioffe–Regel limit in metals_, Philosophical Magazine 84, 2847 (2004).
  * [24] W. Zhang, A. Thiess, P. Zalden, R. Zeller, P. H. Dederichs, J.-Y. Raty, M. Wuttig, S. Blügel, and R. Mazzarello,  _Role of vacancies in metal–insulator transitions of crystalline phase-change materials_, Nature Materials 11, 952 (2012).
  * [25] M. Aftab, G. H. Jaffari, S. K. Hasanain, T. A. Abbas, and S. I. Shah,  _Magnetic and transport properties of Co2Mn1-xCrxSi Heusler alloy thin films_, Journal of Applied Physics 114, 103903 (2013).
  * [26] R. Dynes, J. Rowell, and P. Schmidt, Ternary Superconductors,
  * [27] Y. Feng, J. Y. Rhee, T. A. Wiener, D. W. Lynch, B. E. Hubbard, A. J. Sievers, D. L. Schlagel, T. A. Lograsso, and L. L. Miller,  _Physical properties of Heusler-like Fe2VAl_, Physical Review B 63, 165109 (2001).
  * [28] M. Fukuhara, C. Gangli, K. Matsubayashi, and Y. Uwatoko,  _Pressure-induced positive electrical resistivity coefficient in Ni-Nb-Zr-H glassy alloy_, Applied Physics Letters 100, 253114 (2012).
  * [29] K. Jin, B. C. Sales, G. M. Stocks, G. D. Samolyuk, M. Daene, W. J. Weber, Y. Zhang, and H. Bei,  _Tailoring the physical properties of Ni-based single-phase equiatomic alloys by modifying the chemical complexity_, Scientific Reports 6, 20159 (2016).
  * [30] K. Sumiyama, M. Yamazaki, T. Yoneyama, K. Suzuki, K. Takemura, Y. Kurokawa, and T. Hihara,  _Electric and Magnetic Evolution in Sputter-Deposited FexSi1-x Alloy Films_, MATERIALS TRANSACTIONS 57, 907 (2016).
  * [31] C. C. Tsuei,  _Nonuniversality of the Mooij Correlation—the Temperature Coefficient of Electrical Resistivity of Disordered Metals_, Physical Review Letters 57, 1943 (1986).
  * [32] J. A. N. Bruin, H. Sakai, R. S. Perry, and A. P. MacKenzie,  _Similarity of scattering rates in metals showing T-linear resistivity_, Science 339, 804 (2013).
  * [33] S. Mahmoudian and V. Dobrosavljević,  _Landau-Ginzburg Theory for Anderson localization_, arXiv.org 1503, arXiv:1503.00420 (2015).

Supplementary information

##  Appendix A Derivation of conductivity using Kubo formula

In the absence of vertex corrections the Kubo formula for conductivity reads

$\sigma_{xx}=\pi\int\frac{d^{d}p}{(2\pi)^{d}}\int dz\;A({\bf p},z)j_{x}({\bf
p})A({\bf p},z)j_{x}({\bf p})\left(-f^{\prime}(z)\right).$ (7)

where ${\bf j}({\bf p})=\frac{e{\bf p}}{m}$ is the current operator,
$f(z)=[e^{z/T}+1]^{-1}$ is the Fermi function, and $A({\bf p},z)$ is the
spectral function

$A({\bf p},z)=\frac{1}{\pi}\frac{-\mathrm{Im}\Sigma({\bf p},z)}{(z-\xi_{\bf
p}-\mathrm{Re}\Sigma({\bf p},z))^{2}+(\mathrm{Im}\Sigma({\bf p},z))^{2}}.$ (8)

Here $\Sigma({\bf p},z)$ is the electron self-energy. When the self-energy is
independent of momentum, we can replace the momentum integral by an integral
over dispersion $\xi$ and introduce the transport function as discussed in the
main text,

$\displaystyle\Phi_{x}(\xi)$ $\displaystyle=$
$\displaystyle\int\frac{d^{d}p}{(2\pi)^{d}}j^{2}_{x}({\bf
p})\,\delta(\xi-\xi_{\bf p})$ (9) $\displaystyle=$
$\displaystyle\frac{2e^{2}}{dm}\left(\xi+\mu\right)N(\xi).$ (10)

so that the Kubo formula reads

$\sigma_{xx}=\pi\int d\xi\Phi_{x}(\xi)\int
dzA^{2}(\xi,z)\left(-f^{\prime}(z)\right)$ (11)

Now we include the real part of the self-energy in a renormalization of the
bare dispersion; and assume the imaginary part of the self-energy is ${\rm
Im}\Sigma(z)=-\frac{1}{2\tau}\Theta(z+\mu)$. In this case, relevant for the
novel moiré systems, we can exactly integrate over $\xi$ in $d=2$ dimensions,

$\displaystyle\sigma(T,\mu)$ $\displaystyle=$ $\displaystyle\frac{e^{2}}{2\pi
h}\int_{-\mu}dz\left(-f^{\prime}(z)\right)\left[1+2\tau(z+\mu)\left(\frac{\pi}{2}+\mathrm{atan}\left(2\tau(z+\mu)\right)\right)\right],$
(12)

The leading order term in the limit where $(\mu+z)\tau$ is large, becomes

$\displaystyle\sigma(T,\mu)$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{h}\tau\int_{-\mu}dz\left(-f^{\prime}(z)\right)(z+\mu)=\frac{e^{2}}{h}\tau
T\log\left[1+e^{\mu/T}\right]$ (13)

which is the central result of the main text.

##  Appendix B Nonuniversal behavior close to the transition

The scaling ansatz is in general not valid arbitrarily close to the
transition. Explicit breakdown of scaling can be seen in theoretical models
even within the weak-coupling Kubo formula without vertex corrections
(meaning: even if we ignore Anderson localization, Mott localization, Wigner
crystallization, and percolation).

The scaling ansatz $\sigma(\mu,T)=\sigma_{c}(T)F(\mu/T)$ in the limit of
$T\rightarrow 0$ implies, with $\sigma_{c}(T)\sim T$, that
$\sigma(T=0)\propto\mu$. This is trivially true for the Drude formula in
$d=2$, whereas in $d=3$ it requires ${\rm Im}\Sigma\propto\sqrt{\mu}$. The
breakdown of scaling close to the transition can thus be inferred from having
nonlinear behavior of $\sigma$ as a function of $\mu$ at zero temperature.

In 2nd order perturbation theory with disorder we obtain ${\rm
Im}\Sigma=-\frac{1}{2\tau_{0}}$. The zero-temperature conductivity is thus,
for $\mu>0$,

$\displaystyle\sigma(\mu,T=0)=\frac{e^{2}}{2\pi
h}\left[1+2\mu\tau_{0}\left(\frac{\pi}{2}+\arctan 2\mu\tau_{0}\right)\right]$
(14)

and $\sigma=0$ on the insulating side. This result features a jump at $\mu=0$
of magnitude $\Delta\sigma=\frac{e^{2}}{2\pi}$, breaking the scaling ansatz.

This jump persists even when the density of states do not have a discontinuity
at the band edge. In the self-consistent Born approximation (SCBA) the self-
energy is given by $\Sigma(z)=\gamma G(z)$ where $G(z)$ is the fully dressed
Green’s function and $\gamma$ the strength of the disorder. Note that within
this scheme, the density of states $A(z)=-\frac{1}{\gamma\pi}{\rm
Im}\Sigma(z)$, which is now continuous at the band edge. Nevertheless, there
is still a (nonuniversal) jump in the conductivity.

In $d=3$, integrating over momenta in the Kubo formula yields the zero-
temperature conductivity for $\mu>0$,

$\sigma(\mu,T=0)=\frac{\sqrt{m}\left(\text{Im{\hbox{\Sigma(\mu)}}}^{2}+\mu^{2}\right)^{3/4}\left(\sqrt{\frac{|\text{Im{\hbox{\Sigma(\mu)}}}|}{\sqrt{\text{Im{\hbox{\Sigma(\mu)}}}^{2}+\mu^{2}}}+1}+\sqrt{1-\frac{|\text{Im{\hbox{\Sigma(\mu)}}}|}{\sqrt{\text{Im{\hbox{\Sigma(\mu)}}}^{2}+\mu^{2}}}}\right)^{3}}{24\sqrt{2}\pi^{2}\text{Im{\hbox{\Sigma(\mu)}}}}$
(15)

Note that indeed when $\mu>{\rm Im}\Sigma$ this reduces to the Drude formula
$\sigma=\frac{e^{2}\sqrt{2m}\mu^{3/2}\tau}{3\pi^{2}}=\frac{e^{2}n\tau}{m}$. In
the main text we assumed for scaling to hold $\tau^{-1}$ to be proportional to
the density of states, which is proportional to $\sqrt{\mu}$ near the band
edge. This results in a conductivity proportional to $1/\sqrt{\tau}$ close to
the transition, violating the scaling ansatz

All three cases are visualized in Fig. 5.

Figure 5:  The theoretical zero-temperature conductivity displays a breakdown
of scaling (in black dashed lines) near the MIT, in $d=2$ for fixed $\tau=1$
(a.), the self-consistent Born approximation (b.), and in $d=3$ for
$\tau^{-1}=\sqrt{\mu}$ (c.). The red shaded area indicates the deviation from
scaling. The inset in b. shows the imaginary part of the self-energy for the
SCBA with $\gamma=1$ in the same energy range as the main figure.

##  Appendix C Analysis of the experimental data: Scaling

Figure 6:  In the experimental scaling plots in the main text, we divide for
each curve the temperature by a scaling parameter $T_{0}$. Here we show the
value of $T_{0}$ as a function of tuning parameter for the three considered
experimental systems. a. The $\nu=2$ metal-insulator transition observed in
MoTe2/WSe2[5], the tuning parameter is the vertical displacement field $D$
measured in V/nm. We identify the critical $D_{c}$ as where the extrapolated
$\sigma(T=0)$ vanishes. b. In WSe2/BLG/WSe2, the tuning parameter is the
vertical displacement field. We could not identify a critical value of
$D_{c}$. Instead, we found $R_{c}(T)$ by optimizing the data collapse of all
the experimental curves. Note that we restricted our dataset to a regime of
$D$ where the change in resistance is monotonic. The value $T_{0}$ can be
identified as a chemical potential for conducting charge carriers. These are
possibly the result of potential fluctuations present in bilayer graphene. c.
The different samples of GST each have a different annealing history and
composition of Ge, Sb and Te, as outlined in Ref. [20]. Here we ordered the
data sets based on the scaling parameter $T_{0}$ ranging from insulating
($T_{0}<0$) to metal ($T_{0}>0$).

In order to perform scaling, we used publicly available resistance/resistivity
data from Refs. [5, 19, 20]. For each data set, we first identified the
critical curve of the form

$R_{c}(T)=\frac{1}{\tau_{0}T}+b$ (16)

where $\tau_{0}$ is the disorder-induced scattering time and $b$ is some
linear-$T$ scattering rate that can come from phonons or general ‘Planckian’
dissipation. Written like Eq. (16), these two parameters have units that
depend on the dimension ($d=2,3$). For the three experimental systems, we used
the following parameters for the critical curve:

System | $\tau_{0}$ | $b$
---|---|---
MoTe2/WSe2[5] | 192 k$\Omega$/K | 2.72 k$\Omega$
WSe2/BLG/WSe2[19] | 84.1 k$\Omega$/K | –
GeSbTe[20] | 1.35 $\Omega$cm/K | 0.003 $\Omega$cm

We rescaled the temperature with a parameter $T_{0}$ (which can be interpreted
as the chemical potential $\mu$) as a function of tuning parameter to achieve
a collapse of all resistance/resistivity curves. The relevant values of
$T_{0}$ are shown in Fig. 6.

For MoTe2/WSe2 and WSe2/BLG/WSe2 (BLG stands for bilayer graphene), the tuning
parameter is a vertical displacement field. The data reported are at filling
$\nu=-2$ and $\nu=0$, respectively, relative to charge neutrality.

The label GST refers to a collection of different compounds GeSb4Te7,
GeSb2Te4, Ge2Sb2Te5, and Ge3Sb2Te6. Different resistivity curves correspond to
those four materials each annealed at a different temperature, as outlined in
the Supplementary Information of Ref. [20].

##  Appendix D Breakdown of scaling

Note that in the case of MoTe2/WSe2, as expected, the scaling breaks down at
low temperatures. Details of this are shown in Fig. 7. Note that the breakdown
of scaling in the epxerimental system is consistent with the theoretical
picture of Fig. 1.

Figure 7: At low temperatures ($T\tau\leq 1$) we expect a breakdown of the
scaling regime. This is indeed observed in MoTe2/WSe2[5]. a. Here we show only
the low-temperature data ($T\leq 3$ or $7$ K, respectively on the insulator or
metal side of the MIT), scaled with the same $T_{0}$ and $R_{c}$ as is used
for the figure in the main text. The experimental data points are clearly not
on the theoretical scaling curve (dashed line). b. We can quantify the
deviation from scaling by calculating the distance-squared between the
experimental data points and the theoretical curve. Here we show in colorscale
(blue equals good fit, red is bad fit) the quality of the scaling ansatz as a
function of displacement field and temperature. Using the estimate of
$\tau_{0}$, we show with a dashed line where one would expect scaling to break
down.

##  Appendix E Analysis of the experimental data: Mooij correlations

The data presented in Fig. 4 of the main manuscript is collected from a
variety of sources. The analysis to arrive at the dimensionless temperature
coefficient of the resistivity $\alpha$ is for most materials[25, 26, 27, 28,
29, 30, 20] based on our earlier analysis presented in Ref. [10].

The TiAl data presented is from Fig. 6 of Ref. [9]. The three different sets
of data presented there are shown in Fig. 4 of the main manuscript with filled
circles of different colors.

The Mooij correlations for MoTe2/WSe2 are based on the same data analyzed for
the scaling of Fig. 3[5], limited to the temperature range $T=26$ – 63 K, and
the range of displacement fields $D=0.476$ – 0.408 V/nm. The temperature
derivative of the resistance is calculated with a two-point forward finite
difference. The only fitting parameter is $R^{\infty}=1.44$ k$\Omega$ is
obtained by collapsing the data for different displacement fields onto the
same curve. In Fig. 4 of the main manuscript, the data for MoTe2 is presented
with an empty upward triangle where the different colors represent the
different displacement fields.