# Dynamics of Fluctuations in Quantum Simple Exclusion Processes

Denis Bernard, Fabian H. L. Essler, Ludwig Hruza, Marko Medenjak

► To cite this version:

Denis Bernard, Fabian H. L. Essler, Ludwig Hruza, Marko Medenjak. Dynamics of Fluctuations in Quantum Simple Exclusion Processes. SciPost Physics, 2022, 12 (1), pp.42. 10.21468/SciPost-Phys.12.1.042. [hal-03365177](https://HAL.science/hal-03365177v1)

HAL Id: hal-03365177

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

Submitted on 23 Aug 2024

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
# Dynamics of fluctuations in quantum simple exclusion processes

Denis Bernard¹, Fabian H.L. Essler², Ludwig Hruza¹ and Marko Medenjak³

1 Laboratoire de Physique de l'École Normale Supérieure, CNRS, ENS & PSL University, Sorbonne Université, Université de Paris, 75005 Paris, France

2 The Rudolf Peierls Centre for Theoretical Physics, Oxford University, Parks Road, Oxford OX1 3PU, UK

3 Institut de Physique Théorique Philippe Meyer, École Normale Supérieure, PSL University, Sorbonne Universités, CNRS, 75005 Paris, France

## Abstract

We consider the dynamics of fluctuations in the quantum asymmetric simple exclusion process (Q-ASEP) with periodic boundary conditions. The Q-ASEP describes a chain of spinless fermions with random hoppings that are induced by a Markovian environment. We show that fluctuations of the fermionic degrees of freedom obey evolution equations of Lindblad type, and derive the corresponding Lindbladians. We identify the underlying algebraic structure by mapping them to non-Hermitian spin chains and demonstrate that the operator space fragments into exponentially many (in system size) sectors that are invariant under time evolution. At the level of quadratic fluctuations we consider the Lindbladian on the sectors that determine the late time dynamics for the particular case of the quantum symmetric simple exclusion process (Q-SSEP). We show that the corresponding blocks in some cases correspond to known Yang-Baxter integrable models and investigate the level-spacing statistics in others. We carry out a detailed analysis of the steady states and slow modes that govern the late time behaviour and show that the dynamics of fluctuations of observables is described in terms of closed sets of coupled linear differential-difference equations. The behaviour of the solutions to these equations is essentially diffusive but with relevant deviations, that at sufficiently late times and large distances can be described in terms of a continuum scaling limit which we construct. We numerically check the validity of this scaling limit over a significant range of time and space scales. These results are then applied to the study of operator spreading at large scales, focusing on out-of-time ordered correlators and operator entanglement.

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

Copyright D. Bernard *et al.*

This work is licensed under the Creative Commons

Attribution 4.0 International License.

Published by the SciPost Foundation.

Received 19-07-2021

Accepted 22-11-2021

Published 28-01-2022

doi:10.21468/SciPostPhys.12.1.042

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

Check for
updates

## Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>2</td><td>Q-ASEP & Q-SSEP dynamics as spin chain dynamics</td><td>4</td></tr><tr><td>2.1</td><td>Dynamics of fluctuations in the Q-ASEP & Q-SSEP</td><td>4</td></tr><tr><td>2.2</td><td>Q-ASEP & Q-SSEP as non-hermitian gl(2R) spin chains</td><td>6</td></tr><tr><td>2.2.1</td><td>Super-operator formalism</td><td>8</td></tr><tr><td>2.2.2</td><td>Hilbert space doubling formalism</td><td>10</td></tr></table>
<table><tbody><tr><td>3</td><td>U(1) symmetries, sectors and fragmentation</td><td>13</td></tr><tr><td>3.1</td><td>Symmetries</td><td>13</td></tr><tr><td>3.2</td><td>Sectors and fragmentation</td><td>13</td></tr><tr><td>4</td><td>Integrability and non-integrability</td><td>16</td></tr><tr><td>4.1</td><td>Spin chain identification</td><td>16</td></tr><tr><td>4.2</td><td>Numerics</td><td>17</td></tr><tr><td>5</td><td>Steady states and slow modes</td><td>20</td></tr><tr><td>5.1</td><td>Q-SSEP steady states</td><td>20</td></tr><tr><td>5.2</td><td>Q-SSEP low lying magnon states</td><td>21</td></tr><tr><td>5.3</td><td>Cat state preparation in the non trivial integrable sector</td><td>21</td></tr><tr><td>6</td><td>Q-SSEP dynamics on the lattice and in the scaling limit</td><td>23</td></tr><tr><td>6.1</td><td>Two point functions for two replicas</td><td>23</td></tr><tr><td>6.2</td><td>Numerical solution</td><td>25</td></tr><tr><td>6.2.1</td><td>CDW initial state</td><td>25</td></tr><tr><td>6.2.2</td><td>Domain wall initial state</td><td>26</td></tr><tr><td>6.2.3</td><td>Correlated domain wall initial state</td><td>27</td></tr><tr><td>6.3</td><td>Scaling limit</td><td>28</td></tr><tr><td>6.3.1</td><td>Late time regime Dt/Ω² ≫ 1</td><td>30</td></tr><tr><td>6.3.2</td><td>Intermediate time regime Dt/Ω² ≪ 1 ≪ Dt/a₀²</td><td>31</td></tr><tr><td>6.3.3</td><td>Numerical test of the scaling limit</td><td>31</td></tr><tr><td>6.3.4</td><td>"Diagonal" correlations in the scaling limit</td><td>33</td></tr><tr><td>6.4</td><td>Evolution of higher point correlation functions</td><td>34</td></tr><tr><td>6.5</td><td>Dictionary for 2-replica operators</td><td>35</td></tr><tr><td>7</td><td>Application: Operator spreading</td><td>36</td></tr><tr><td>7.1</td><td>Out-of-time ordered correlation functions</td><td>37</td></tr><tr><td>7.2</td><td>Operator entanglement</td><td>37</td></tr><tr><td>8</td><td>Summary and conclusions</td><td>40</td></tr><tr><td>A</td><td>Dictionary: Super-operator versus Hilbert space doubling approach</td><td>41</td></tr><tr><td>A.1</td><td>One replica</td><td>41</td></tr><tr><td>A.2</td><td>Two replicas</td><td>42</td></tr><tr><td>B</td><td>Steady states, low-lying states and gaps</td><td>42</td></tr><tr><td>C</td><td>GL(2R) and U(L) interplay on steady states</td><td>44</td></tr><tr><td>D</td><td>Unfolding and level-spacing for a Poisson process.</td><td>47</td></tr><tr><td>E</td><td>Two point functions in the Lindbladian formulation</td><td>48</td></tr><tr><td>F</td><td>Spectrum of the two-replica Lindbladian in the two-particle sector</td><td>48</td></tr><tr><td>G</td><td>Derivation of the time evolution of n-point bubbles</td><td>51</td></tr><tr><td>G.1</td><td>Three-point bubble</td><td>51</td></tr><tr><td>G.2</td><td>n-point bubbles</td><td>53</td></tr><tr><td></td><td>References</td><td>54</td></tr></tbody></table>
# 1 Introduction

Environmental effects on many-particle quantum dynamics continue to attract a great deal of interest for experimental [1,2] as well as fundamental theoretical reasons [3,4]. In many cases of interest the environment is approximately Markovian, which allows for a description of the average system dynamics by means of the Lindblad formalism [5-7]. In the many-particle context Lindblad equations are generally quite difficult to solve and a natural question that arises is whether there exist paradigmatic models for which exact results can be obtained. The simplest such class of Lindblad equations can be mapped to imaginary-time Schrödinger equations with non-Hermitian "Hamiltonians" that are quadratic in creation/annihilation operators [8-11]. More recently it has been shown that there exist Lindblad equations whose evolution operators are related to interacting integrable quantum spin chains [12-24]. This opens the door to obtaining exact results on the average dissipative dynamics of many-particle Lindblad equations by employing methods from quantum integrability.

While historically the focus of attention has mainly been on the average system dynamics, there are interesting questions pertaining to fluctuations of system degrees of freedom induced by their coupling to the environment. While classical stochastic dynamics has been a well studied domain for many decades, stochastic dynamical quantum extended systems remained largely unexplored, to the best of our knowledge, until the very recent developments in random quantum circuit models [25-29]. These advances proceeded in parallel, but independently, to work on quantum exclusion processes and more generally stochastic dynamics in quantum many-body systems such as spin chains coupled to Markovian "noise" [30]. These model systems provide quantum extensions of discrete versions of the classical fluctuating hydrodynamics whose impact on our understanding of transport in low dimensional systems cannot be underestimated [31-33]. Moreover, the average late time dynamics in such models typically reduces to classical stochastic processes, which have been intensively analyzed over the last several decades [34-37]. A key motivation for studying fluctuating quantum dynamics, as in [30], is the search for a quantum extension of the classical macroscopic fluctuation theory (MFT) [38, 39], see [40]. The latter is a general setting, whose formulation relies on stochastic hydrodynamic equations and gives rise to a precise description of the fluctuations in classical, diffusive, out-of-equilibrium systems. Quantum simple exclusion processes are paradigmatic models for fluctuating diffusive quantum dynamics. The closed quantum symmetric simple exclusion process (Q-SSEP), was introduced in [41] (under another name), its open version in [42], and its asymmetric analogue, the quantum asymmetric simple exclusion process (Q-ASEP) was proposed in [43]. The average dynamics of both the Q-ASEP and the Q-SSEP was recently shown to be integrable [22] and it reduces to their classical analogues (which have been studied in great detail, see e.g. [44-55]) at late times. However, fluctuations in the Q-ASEP and Q-SSEP are not classical and reveal patterns which are possibly generic for mesoscopic diffusive quantum systems. In a recent work, fluctuations in the steady state were analysed in the framework of a continuum limit of the open Q-SSEP [56] by constructing the associated invariant probability measure.

In this work we carry out a detailed analysis of the (finite time) dynamics of fluctuations in Q-ASEP and Q-SSEP.

In section 2 we introduce the Q-ASEP and Q-SSEP and formulate the dynamics of fluctuations in terms of Lindblad equations. The corresponding Lindbladians are shown to take the form of non-Hermitian $gl(2R)$ quantum spin chains. In section 3 we establish that the Lindbladians fragment [22, 23] into an exponential (in system size) number of sectors characterized by local conserved charges. The remaining part of the paper focuses on the case of the Q-SSEP. In section 4 we focus on the five sectors that describe the late time dynamics of quadratic fluctuations in the Q-SSEP and show that two of them are Yang-Baxter integrable
and two others are trivial. We then investigate the level-spacing statistics of the Lindbladian eigenvalues in the remaining sector (as well as the two integrable ones) in order to ascertain whether it may be integrable as well. In section 5 we identify the algebraic structure of all steady states and low-lying "magnon" excitations in terms of the underlying gl(2R) algebra. In section 6 we turn to a detailed analysis of the dynamics of correlation functions on the lattice, and show that there is a scaling regime in which correlators are described by hierarchies of partial differential equations that take the form of diffusion equations with source terms that couple the different levels of the hierarchy. We then apply these findings to the large scale dynamics of operator spreading in section 7. Section 8 contains our conclusions. A number of technical details are discussed in six appendices.

# 2 Q-ASEP & Q-SSEP dynamics as spin chain dynamics

## 2.1 Dynamics of fluctuations in the Q-ASEP & Q-SSEP

The quantum asymmetric (Q-ASEP) or symmetric (Q-SSEP) simple exclusion processes [41, 42] are models of quantum many-body stochastic dynamics describing fermions hopping stochastically along a one dimensional chain. The chain can be open, with injection/extraction of fermions at the two ends of the chain, or closed, without any external injection/extraction processes. We shall deal with the latter case, in which the state converges to a steady state dressed by sub-leading fluctuations at late time [42,56]. The purpose of the following analysis is to provide the description for the dynamics of those fluctuations in terms of Lindbladian time evolution.

In the closed setup, the Q-SSEP & Q-ASEP dynamics is unitary but noisy. This means that the system density matrix $\rho_t$ evolves unitarily, $\rho_{t+dt} = e^{-idH_t} \rho_t e^{idH_t}$, with the Hamiltonian increment between times $t$ and $t+dt$ $dH_t$, which depends on the external noise<fnref n="1" />. For Q-SSEP & Q-ASEP, on the chain of length $L$ and periodic boundary conditions, the Hamiltonian increment reads [41–43]

$$dH_t := \sum_{j=1}^{L} (c_{j+1}^{\dagger} c_j dW_t^j + c_j^{\dagger} c_{j+1} d\overline{W}_t^j), \quad (1)$$

where $c_j$ and $c_j^\dagger$ are canonical fermion operators at site $j$ of the chain,

$$\{c_j, c_k^\dagger\} = \delta_{j,k}, \qquad (2)$$

while $W_t^j$ and $\overline{W}_t^j$ are classical complex Brownian motions in the Q-SSEP case and complex quantum Brownian motions in the case of the Q-ASEP. Specifically, they correspond to Gaussian processes with zero mean $\mathbb{E}[dW_t^j] = 0$, $\mathbb{E}[d\overline{W}_t^j] = 0$ and covariances

$$\mathbb{E}[dW_t^j d\overline{W}_t^k] = pJ \delta^{j;k} dt, \quad \mathbb{E}[d\overline{W}_t^j dW_t^k] = qJ \delta^{j;k} dt, \quad (3)$$

and independent increments. Here $p, q$ are real positive dimensionless parameters, $J$ is a dimensionfull rate parameter, and $\mathbb{E}$ denotes the average over the Brownian motions. In the

<footnote n="1">¹A better notation would be to label the noise by a variable, say $\omega$, taking values in some adapted probability space, and in denoting the Hamiltonian increments as $dH_t(\omega)$ to indicate its noise dependence. As a consequence, the density matrix, which also noise dependent, should be denoted as $\rho_t(\omega)$. We lighten the notation making the $\omega$-dependence implicit.</footnote>
symmetric case (Q-SSEP), $p = q = 1$, the processes $W_t^j$ and $\overline{W}_t^k$ commute, while in the asymmetric case $p \neq q$ (Q-ASEP) the commutator is nontrivial,

$$[dW_t^j, d\overline{W}_t^k] = (p-q)J \delta^{j;k} dt. \quad (4)$$

Periodic boundary conditions are assumed throughout.

The mean density matrix $\bar{\rho}_t := \mathbb{E}[\rho_t]$, which is obtained by averaging over all possible realizations of the Brownian motion, satisfies a dissipative Lindblad dynamics [30]

$$\partial_t \bar{\rho}_t = \mathcal{L}(\bar{\rho}_t), \quad (5)$$

where the Lindbladian $\mathcal{L} = \sum_{j=1}^L \mathcal{L}_{j;j+1}$ is the sum of local terms $\mathcal{L}_{j;j+1}$, which can be decomposed as

$$\mathcal{L}_{j;j+1} = J(p \overrightarrow{\mathcal{L}}_{j;j+1} + q \overleftarrow{\mathcal{L}}_{j;j+1}), \quad (6)$$

where $\overrightarrow{\mathcal{L}}_{j;j+1}$ corresponds to the hopping from site $j$ to site $j+1$ and $\overleftarrow{\mathcal{L}}_{j;j+1}$ the other way around:

$$\overrightarrow{\mathcal{L}}_{j;j+1}(\rho) := \ell_j^- \rho \ell_j^+ - \frac{1}{2}(\ell_j^+ \ell_j^- \rho + \rho \ell_j^+ \ell_j^-), \quad (7a)$$

$$\overleftarrow{\mathcal{L}}_{j;j+1}(\rho) := \ell_j^+ \rho \ell_j^- - \frac{1}{2}(\ell_j^- \ell_j^+ \rho + \rho \ell_j^- \ell_j^+). \quad (7b)$$

Here the jump operators are given by $\ell_j^+ := c_{j+1}^\dagger c_j$ and $\ell_j^- := c_j^\dagger c_{j+1}$. As a consequence, the total Lindbladian can be decomposed as

$$\mathcal{L} = J(p \overrightarrow{\mathcal{L}} + q \overleftarrow{\mathcal{L}}), \quad (8)$$

where the Lindbladian $\overrightarrow{\mathcal{L}}$ (resp. $\overleftarrow{\mathcal{L}}$) generates fermionic hopping to the right (resp. left). The cases where $p = 0$ or $q = 0$ correspond to the totally asymmetric processes (Q-TASEP).

Let $O$ be an operator acting on the system degrees of freedom only. We denote the associated Heisenberg picture operator acting on the full Hilbert space, comprised of the system and noise, by $O_t$, and its average over the noise by $\bar{O}_t := \mathbb{E}[O_t]$. The dynamics of $\bar{O}_t$ follows from duality $\text{Tr}(\bar{\rho}_t O) = \text{Tr}(\rho_0 \bar{O}_t)$, and is of the Lindblad type as well,

$$\partial_t \bar{O}_t = \mathcal{L}^*(O)_t, \quad (9)$$

where the dual Lindbladian $\mathcal{L}^* = J(p \overrightarrow{\mathcal{L}}^* + q \overleftarrow{\mathcal{L}}^*)$, is also the sum of local contributions, $\overrightarrow{\mathcal{L}}^* = \sum_j \overrightarrow{\mathcal{L}}_{j;j+1}^*$ and $\overleftarrow{\mathcal{L}}^* = \sum_j \overleftarrow{\mathcal{L}}_{j;j+1}^*$. We have

$$\overrightarrow{\mathcal{L}}_{j;j+1}^*(O) := \ell_j^+ O \ell_j^- - \frac{1}{2}(\ell_j^+ \ell_j^- O + O \ell_j^+ \ell_j^-), \quad (10a)$$

$$\overleftarrow{\mathcal{L}}_{j;j+1}^*(O) := \ell_j^- O \ell_j^+ - \frac{1}{2}(\ell_j^- \ell_j^+ O + O \ell_j^- \ell_j^+). \quad (10b)$$

Although $\overrightarrow{\mathcal{L}}_{j;j+1}^*$ (resp. $\overleftarrow{\mathcal{L}}_{j;j+1}^*$) is by construction the dual of $\overrightarrow{\mathcal{L}}_{j;j+1}$ (resp. $\overleftarrow{\mathcal{L}}_{j;j+1}$), it is different from $\overleftarrow{\mathcal{L}}_{j;j+1}$ (resp. $\overrightarrow{\mathcal{L}}_{j;j+1}$) since we have $\overrightarrow{\mathcal{L}}_{j;j+1}^*(O) - \overleftarrow{\mathcal{L}}_{j;j+1}(O) = \frac{1}{2}\{O, n_j - n_{j+1}\}$. However, summing over all sites of the chain, with periodic boundary conditions, implies that

$$\overrightarrow{\mathcal{L}}^* = \overleftarrow{\mathcal{L}}, \quad \overleftarrow{\mathcal{L}}^* = \overrightarrow{\mathcal{L}}. \quad (11)$$

In particular, exchanging $p$ and $q$, by reversing the orientation of the chain, amounts to exchange the total Lindbladian with its dual.
The above Lindbladians describe the time evolution of averages $E[\mathrm{Tr}(\rho_t O)]$ of operators $O$ acting on the system degrees of freedom but not of their fluctuations. Quadratic fluctuations associated with two operators $O_1$ and $O_2$ are given by $E[\mathrm{Tr}(\rho_t O_1)\mathrm{Tr}(\rho_t O_2)]$, and higher moments are defined analogously. The dynamics of the quadratic fluctuations $\mathrm{Tr}(\rho_t O_1)\mathrm{Tr}(\rho_t O_2) = \mathrm{Tr}[(\rho_t \otimes \rho_t)(O_1 \otimes O_2)]$ can be obtained from a replicated system with density matrix $\rho_t^{(2)} := E[\rho_t \otimes \rho_t]$. Here the two replicas of the system are coupled to the same realisation of the noise, which is then averaged over. Physically this is equivalent to a three-leg ladder, where the outer legs corresponding to the fermion degrees of freedom are coupled in an identical fashion to the central leg, which describes the noise degrees of freedom.

As a consequence, the dynamics of the quadratic fluctuations is still of Lindblad type,

$$ \partial_t \rho_t^{(2)} = \mathcal{L}(\rho_t^{(2)}), \qquad (12) $$

where $\mathcal{L} = J(p\vec{\mathcal{L}} + q\overleftarrow{\mathcal{L}})$ and $\vec{\mathcal{L}} = \sum_j \overleftarrow{\mathcal{L}}_{j;j+1}$ and $\overleftarrow{\mathcal{L}} = \sum_j \overleftarrow{\mathcal{L}}_{j;j+1}$ as in eqn (7), but with altered jump operators<fnref n="2" />

$$ \ell_j^{\pm} \rightarrow \ell_j^{\pm} \otimes \mathbb{I} + \mathbb{I} \otimes \ell_j^{\pm} = \ell_{1;j}^{\pm} + \ell_{2;j}^{\pm}, \qquad (13) $$

which correspond to the sum of jumps on both replicas. Here the indices 1,2 refer to the labelling of the replicas. A similar structure applies for higher moments, and hence higher number of replicas

$$ \ell_j^{\pm} = \sum_{a=1}^{R} \ell_{a;j}^{\pm}, \qquad (14) $$

where the index $a$ labels the different replicas.

By duality, the dynamics of replicated observables is also of Lindblad type. Quadratic fluctuations are encoded in $E[O_1 \otimes O_2]$, where as before the averaging is over the different realizations of the Brownian motions. Then we have

$$ \partial_t E[O_1 \otimes O_2]_t = \mathcal{L}^*(E[O_1 \otimes O_2])_t, \qquad (15) $$

where the dual Lindbladian $\mathcal{L}^*$ has the same structure as in eqn (10) but with the replacements $\ell_j^{\pm} \rightarrow \ell_j^{\pm} \otimes \mathbb{I} + \mathbb{I} \otimes \ell_j^{\pm} = \ell_{1;j}^{\pm} + \ell_{2;j}^{\pm}$. Similar remarks apply to higher moments.

## 2.2 Q-ASEP & Q-SSEP as non-hermitian $gl(2R)$ spin chains

The aim of this section is to represent the Q-SSEP & Q-ASEP replicated dynamics as non-Hermitian quantum spin chains. Viewing the replicas as different internal "color" degrees of freedom, the dynamics will include intra- and inter- replica interactions. As a consequence, the Q-SSEP & Q-ASEP dual Lindbladian for two replicas, i.e. for quadratic fluctuations, can be decomposed as

$$ \vec{\mathcal{L}}^* = \vec{\mathcal{L}}_1^* + \vec{\mathcal{L}}_2^* + \vec{\mathcal{L}}_{12}^*, \quad \overleftarrow{\mathcal{L}}^* = \overleftarrow{\mathcal{L}}_1^* + \overleftarrow{\mathcal{L}}_2^* + \overleftarrow{\mathcal{L}}_{12}^*, $$

where $\vec{\mathcal{L}}_a^*(\overleftarrow{\mathcal{L}}_a^*)$, $a = 1,2$ and $\overleftarrow{\mathcal{L}}_{12}^*(\vec{\mathcal{L}}_{12}^*)$ represent respectively intra and inter-replica Lindbladians. Each of them is the sum of local Lindbladians acting on edges connecting two sites, say $\vec{\mathcal{L}}_1^* = \sum_j \overleftarrow{\mathcal{L}}_{1;j,j+1}^*$ and similarly for the other terms. For a higher number $R$ of replicas, the dual Lindbladians can be decomposed as

$$ \vec{\mathcal{L}}^* = \sum_a \vec{\mathcal{L}}_a^* + \sum_{a<b} \vec{\mathcal{L}}_{ab}^*, \quad \overleftarrow{\mathcal{L}}^* = \sum_a \overleftarrow{\mathcal{L}}_a^* + \sum_{a<b} \overleftarrow{\mathcal{L}}_{ab}^*, \qquad (16) $$

<footnote n="2">²To be precise, we should have added to the Lindbladians a label referring to the number of replicas, say $\mathcal{L}^{(2)}$. We made it implicit to lighten the writing.</footnote>
where $\vec{\mathcal{L}}_{ab}^*(\vec{\mathcal{L}}_{ab}^*)$ are the inter-replica Lindbladians and $\vec{\mathcal{L}}_a^*(\vec{\mathcal{L}}_a^*)$, $a = 1, 2$, the intra-replica ones. In the following we will focus on the two replica case but mention in passing what are the appropriate generalisations to an arbitrary number of replicas.

As we will now describe, these Linbladians can be written in terms of super-operators, i.e. linear maps acting on operators defined over the system Hilbert space, which fulfil the commutation relations of the gl(4) algebra, or the gl(2R) algebra in case of R replicas.

* First, there are local generators $C_{a; j}$, $j = 1, \dots, L$ and $a = 1, \dots, R$ of $u(1)$ charges counting the fermionic degrees of the operators at site $j$. Their action on an arbitrary operator $X$ is defined by

$$C_{a; j} X := [n_{a; j}, X], \qquad (17)$$

where $n_{a; j}$ are the local number operators in the a-th replica.

* Second, on each replica there are generators $J_{a; j}^z$, $J_{a; j}^{\pm}$, $j = 1, \dots, L$, and $a = 1, \dots, R$ forming an sl(2) algebra of super-operators

$$[J_{a; j}^z, J_{a; k}^{\pm}] = \pm 2 \delta_{j; k} J_{a; j}^{\pm}, \quad [J_{a; j}^{+}, J_{a; k}^{-}] = \delta_{j; k} J_{a; j}^{z}, \qquad (18)$$

and commuting with the $u(1)$ charges defined above, $[C_{a; j}, J_{a; k}^*] = 0$. These super-operators are defined below in eqn (31).

The intra-replica Lindbladians are expressed in terms of these super-operators as :

$$\begin{aligned} \vec{\mathcal{L}}_{a; j, j+1}^* &= (J_{a; j+1}^+ J_{a; j}^- + \frac{1}{4}(J_{a; j+1}^z + 1)(J_{a; j}^z - 1) + \frac{1}{4}C_{a; j+1} C_{a; j}), \\ \vec{\mathcal{L}}_{a; j, j+1}^* &= (J_{a; j+1}^- J_{a; j}^+ + \frac{1}{4}(J_{a; j+1}^z - 1)(J_{a; j}^z + 1) + \frac{1}{4}C_{a; j+1} C_{a; j}). \end{aligned} \qquad (19)$$

In particular, in the symmetric case the sum

$$\vec{\mathcal{L}}_a^* + \vec{\mathcal{L}}_a^* = \sum_{j=1}^{L} (J_{a; j+1}^+ J_{a; j}^- + J_{a; j+1}^- J_{a; j}^+ + \frac{1}{2} J_{a; j+1}^z J_{a; j}^z + \frac{1}{2} (C_{a; j+1} C_{a; j} - 1))$$

is the isotropic Heisenberg Hamiltonian. This formula was also observed in [57] in the one replica case.

* Finally, we define in eqn (30) a set of super-operators $G_j^{AB}$, $A, B = (1r), (1l), (2r), (2l)$, satisfying the gl(4) commutation relations,

$$[G_j^{AB}, G_k^{CD}] = \delta_{j,k}(\delta^{BC}G_j^{AD} - \delta^{DA}G_j^{CB}), \qquad (20)$$

such that (for two replicas)

$$\begin{aligned} \vec{\mathcal{L}}_{12; j, j+1}^* &= (G_{j+1}^{1l; 2r} G_j^{2r; 1l} + G_{j+1}^{2l; 1r} G_j^{1r; 2l} \\ &\quad + \frac{1}{2}(G_{j+1}^{1l; 2l} G_j^{2l; 1l} + G_{j+1}^{2l; 1l} G_j^{1l; 2l} + G_{j+1}^{1r; 2r} G_j^{2r; 1r} + G_{j+1}^{2r; 1r} G_j^{1r; 2r})), \end{aligned} \qquad (21a)$$

$$\begin{aligned} \vec{\mathcal{L}}_{12; j, j+1}^* &= (G_{j+1}^{2r; 1l} G_j^{1l; 2r} + G_{j+1}^{1r; 2l} G_j^{2l; 1r} \\ &\quad + \frac{1}{2}(G_{j+1}^{1l; 2l} G_j^{2l; 1l} + G_{j+1}^{2l; 1l} G_j^{1l; 2l} + G_{j+1}^{1r; 2r} G_j^{2r; 1r} + G_{j+1}^{2r; 1r} G_j^{1r; 2r})). \end{aligned} \qquad (21b)$$

The relation between the above gl(2) generators and the gl(4) generators reads $J_{a; j}^z = G_j^{al; al} - G_j^{ar; ar}$, $J_{a; j}^+ = (-)^{R_a^0} G_j^{al; ar}$, $J_{a; j}^- = G_j^{ar; al} (-)^{R_a^0}$, with $(-)^{R_a^0}$ some Wigner-Klein factors to be defined below, and $C_{a; j} = G_j^{al; al} + G_j^{ar; ar} - 1$.
The above remarks extend to the case of an arbitrary number R of replicas. In that case, the generators $G_j^{AB}$ form a representation of $gl(2R)$ with index A, B running over the pairs (ar) and (al) with a = 1, ..., R referring to the replica index and r/l to right/left action. The dual Lindbladians $\vec{\mathcal{L}}_a^*$ ($\vec{\mathcal{L}}_a^*$) and $\vec{\mathcal{L}}_{ab}^*$ ($\vec{\mathcal{L}}_{ab}^*$) are simply given by the previous expressions but replacing the replica indices 1, 2 by a, b.

In the symmetric Q-SSEP case, with R replicas, the total Lindbladian $\mathcal{L}_{\text{ssep}} = J(\vec{\mathcal{L}}^* + \vec{\mathcal{L}}^*)$ can be written as (for Q-SSEP $\mathcal{L}_{\text{ssep}}^* = \mathcal{L}_{\text{ssep}}$)

$$ \mathcal{L}_{\text{ssep}} = J \sum_{j} \left( \sum_{A,B} G_{j+1}^{AB} G_{j}^{BA} - \frac{1}{2} (C_{j+1} + C_{j}) - R \right), \qquad (22) $$

with

$$ C_j := \sum_a C_{a;j} = \sum_A G_j^{AA} - R, \qquad (23) $$

the total $u(1)$ charge over the different replicas at site j. This expression makes the $gl(2R)$ symmetry explicit. Alternatively, one can decompose the $gl(2R)$ algebra as the sum of a $sl(2R)$ algebra plus a $u(1)$ and introduce the $sl(2R)$ generators $J_j^{AB} := G_j^{AB} - \frac{\delta^{AB}}{2R}(C_j + R)$ with $\sum_A J_j^{AA} = 0$. Then

$$ \mathcal{L}_{\text{ssep}} = J \sum_{j} \left( \sum_{A,B} J_{j+1}^{AB} J_{j}^{BA} + \frac{1}{2R} C_{j+1} C_{j} - \frac{R}{2} \right). \qquad (24) $$

We shall present the derivation of these results using two complementary approaches: one defining the appropriate super-operators by their actions on operators via specific left/right multiplications, the other one introducing explicit basis on the space of operators and defining the super-operators via their actions on this basis. They have different advantages depending on the question raised.

Since starting from Section 4 this manuscript deals with Q-SSEP only we shall then drop the label “ssep” attached to the Lindbladian.

### 2.2.1 Super-operator formalism

When dealing with many replicas, we have as many fermion species as there are replicas. We denote them by $c_{a,j}$, $c_{a;j}^\dagger$, with a the replica index and j the site index. They satisfy:

$$ \begin{gathered} \{c_{a;j}^{\dagger}, c_{a;k}\} = \delta_{j;k}, \quad \{c_{a;j}, c_{a;k}\} = 0, \quad \{c_{a;j}^{\dagger}, c_{a;k}^{\dagger}\} = 0, \\ [c_{a;j}^{\dagger}, c_{b;k}] = 0, \quad [c_{a;j}, c_{b;k}] = 0, \quad [c_{a;j}^{\dagger}, c_{b;k}^{\dagger}] = 0, \quad b \neq a, \end{gathered} \qquad (25) $$

i.e. fermion operators anti-commute when in the same replica but commute when in different replicas.

We then define super-operators by left/right multiplications with $c_{a,j}$ or $c_{a;j}^\dagger$. That is, for any operator O on the physical system Hilbert space, we set:

$$ R_{a;j}^{+} O := O c_{a;j}^{\dagger}, \quad R_{a;j}^{-} O := O c_{a;j}, \qquad (26a) $$

$$ L_{a;j}^{+} O := c_{a;j}^{\dagger} O, \quad L_{a;j}^{-} O := c_{a;j} O. \qquad (26b) $$

For instance, the local $u(1)$ charge $C_{a;j}$, whose action is specified by $C_{a;j} O := [n_{a;j}, O]$, can be represented as

$$ C_{a;j} = L_{a;j}^{+} L_{a;j}^{-} - R_{a;j}^{-} R_{a;j}^{+} = L_{a;j}^{+} L_{a;j}^{-} + R_{a;j}^{+} R_{a;j}^{-} - 1. \qquad (27) $$
We recall that the total charge is defined as $C_j = \sum_a C_{a;j}$. The operators defined in this way satisfy mixed commutation/anti-commutation relations (for instance left/right multiplication commute). It is convenient to transform them into a set of fermionic super-operators, anti-commuting for all replica indices and all sites by introducing Klein factors, which we implement through a Wigner like transformation on the replica space (it is not intrinsic, we could have made different choices). We set (for two replicas, but the generalisation to any number of replicas is straightforward)

$$ \begin{aligned} \hat{R}_{1;j}^{\pm} &:= R_{1;j}^{\pm}, & \hat{L}_{1;j}^{\pm} &:= (-)^{R_1^0} L_{1;j}^{\pm}, \\ \hat{R}_{2;j}^{\pm} &:= (-)^{R_1^0+L_1^0} R_{2;j}^{\pm}, & \hat{L}_{2;j}^{\pm} &:= (-)^{R_1^0+L_1^0+R_2^0} L_{2;j}^{\pm}, \end{aligned} $$

where $R_a^0 := \sum_j R_{a;j}^+ R_{a;j}^-$ and $L_a^0 := \sum_j L_{a;j}^+ L_{a;j}^-$ are the (total) number $R/L$ operators. By construction we then have the anti-commutation relations

$$ \begin{aligned} \{\hat{R}_{a;j}^s, \hat{R}_{b;k}^{s'}\} &= \delta_{j,k} \delta_{a,b} \delta^{s+s';0}, \\ \{\hat{R}_{a;j}^s, \hat{L}_{b;k}^{s'}\} &= 0, \\ \{\hat{L}_{a;j}^s, \hat{L}_{b;k}^{s'}\} &= \delta_{j,k} \delta_{a,b} \delta^{s+s';0}. \end{aligned} \qquad (28) $$

We denote these fermions by $f_j^{A\dagger}$, $f_j^A$ with $A = (1r), (1l), (2r), (2l)$, where the labelling $r/l$ refers to right/left multiplication (the generalization to $R$ replicas is $A = (ar)$ or $A = (al)$ with $a = 1, \dots, R$), with

$$ \{f_j^{A\dagger}, f_k^{B\dagger}\} = 0, \quad \{f_j^{A\dagger}, f_k^B\} = \delta_{j,k} \delta^{AB}, \quad \{f_j^A, f_k^B\} = 0. \qquad (29) $$

There are $2R$ such effective fermions (or super-fermions as they are linear maps acting on operators).

The $gl(2R)$ generators are then defined by

$$ G_j^{AB} := f_j^{A\dagger} f_j^B. \qquad (30) $$

By construction, they satisfy the $gl(2R)$ commutation relations. In particular, they include as many $sl(2)$ sub-algebras as the number of replicas which are generated by the super-operators $J_{a;j}^{\pm}$ and $J_{a;j}^{z}$ defined by

$$ J_{a;j}^{+} O = c_{a;j}^{\dagger} O c_{a,j}, \quad J_{a;j}^{-} O = c_{a;j} O c_{a,j}^{\dagger}, \quad J_{a;j}^{z} O = n_{a} O + O n_{a}^{-} - O, \qquad (31) $$

for any operator $O$. We have $J_{a;j}^z = G_j^{al;al} - G_j^{ar;ar}$, $J_{a;j}^+ = (-)^{R_a^0} G_j^{al;ar}$, $J_{a;j}^- = G_j^{ar;al} (-)^{R_a^0}$.

The proof of eqns (19) and (21) is then just a matter of expanding the expressions (10) of the dual Lindbladian in terms of the jump operators $\ell_j^\pm$ and expressing them in terms of the effective fermions $f_j^{A\dagger}$ and $f_j^A$. We shall not reproduce all computations here, as they are slightly lengthy, but only reproduce two examples. Let us first look at one intra-replica contribution to $\mathcal{L}^*(O)$ for some operator $O$. One such contribution is $\ell_{a;j}^+ O \ell_{a;j}^-$, for instance. We have:

$$ \begin{aligned} \ell_{a;j}^{+} O \ell_{a;j}^{-} &= c_{a;j+1}^{\dagger} c_{a,j} O c_{a;j}^{\dagger} c_{a,j+1} \\ &= L_{a;j+1}^{+} R_{a;j+1}^{-} L_{a;j}^{-} R_{a;j}^{+} O = G_{j+1}^{al;ar} G_{j}^{ar;al} O. \end{aligned} $$

An inter-replica contribution is for instance $\ell_{a;j}^+ \ell_{b;j}^- O$. We have:

$$ \begin{aligned} \ell_{a;j}^{+} \ell_{b;j}^{-} O &= c_{a;j+1}^{\dagger} c_{a,j} c_{b;j}^{\dagger} c_{b,j+1} O \\ &= L_{a;j+1}^{+} L_{b;j+1}^{-} L_{b;j}^{+} L_{a;j}^{-} O = G_{j+1}^{al;bl} G_{j}^{bl;al} O. \end{aligned} $$

It is clear that all terms contributing to the (dual) Lindbladians are quadratic in the $gl(2R)$ generators. Gathering all terms proves eqns (19) and (21).
### 2.2.2 Hilbert space doubling formalism

A basis of quantum states over site j of our one-dimensional lattice is

$$|0\rangle_j, \quad |1\rangle_j = c_j^\dagger |0\rangle_j. \qquad (32)$$

A basis of the Hilbert space of states $\mathcal{H}$ is thus

$$|\sigma\rangle = \otimes_{j=1}^{L} |\sigma_j\rangle_j, \quad \sigma_j \in \{0, 1\}. \qquad (33)$$

In the one-replica case the noise-averaged density matrix can be expanded in this basis as

$$\bar{\rho} = \sum_{\alpha, \beta} \rho_{\alpha, \beta} |\alpha\rangle\langle\beta|. \qquad (34)$$

In the Hilbert-space doubling formalism this is considered as a vector on the linear vector space of operators on $\mathcal{H}$

$$\|\bar{\rho}\rangle = \sum_{\alpha, \beta} \rho_{\alpha, \beta} |\alpha\rangle \otimes |\beta\rangle. \qquad (35)$$

A basis of the 4-dimensional space of states associated with site j is given by

$$\begin{aligned} \|1\rangle_j &= |1\rangle_j \otimes |1\rangle_j, & \|2\rangle_j &= |0\rangle_j \otimes |1\rangle_j, \\ \|3\rangle_j &= |1\rangle_j \otimes |0\rangle_j, & \|4\rangle_j &= |0\rangle_j \otimes |0\rangle_j. \end{aligned} \qquad (36)$$

Introducing notations such that

$$\mathcal{O} = \mathcal{O} \otimes 1, \quad \tilde{\mathcal{O}} = 1 \otimes \mathcal{O}, \qquad (37)$$

and using that the jump operators $\ell_j^+ := c_{j+1}^\dagger c_j$ fulfil $(\ell_j^\pm)^T = \ell_j^\mp$, the Lindblad equation (5), (7) is then recast in the form

$$\frac{d}{dt} \|\bar{\rho}_t\rangle = \mathcal{L} \|\bar{\rho}_t\rangle, \qquad (38)$$

where

$$\mathcal{L} = J \sum_j p \left[ \ell_j^- \tilde{\ell}_j^- - \frac{1}{2} (\ell_j^+ \ell_j^- + \tilde{\ell}_j^+ \tilde{\ell}_j^-) \right] + q \left[ \ell_j^+ \tilde{\ell}_j^+ - \frac{1}{2} (\ell_j^- \ell_j^+ + \tilde{\ell}_j^- \tilde{\ell}_j^+) \right]. \qquad (39)$$

In order to represent $\mathcal{L}$ in the basis (36) we define Hubbard operators in the usual way by

$$E_j^{ab} \equiv \|a\rangle_j \langle b\|. \qquad (40)$$

Keeping track of minus signs arising from fermionic anticommutation relations we have

$$\begin{aligned} \ell_j^- &= (E_j^{12} + E_j^{34})(E_{j+1}^{21} + E_{j+1}^{43}), \\ \tilde{\ell}_j^- &= (E_j^{13} - E_j^{24})(E_{j+1}^{31} - E_{j+1}^{42}). \end{aligned} \qquad (41)$$

This results in

$$\mathcal{L} = J \sum_j p \left[ J_j^+ J_{j+1}^- + \frac{1}{4} (J_j^z J_{j+1}^z - 1) \right] + q \left[ J_j^- J_{j+1}^+ + \frac{1}{4} (J_j^z J_{j+1}^z - 1) \right] + \frac{p+q}{4} C_j C_{j+1}, \qquad (42)$$

where we have defined

$$J_j^+ = E_j^{14}, \quad J_j^- = E_j^{41}, \quad J_j^z = E_j^{11} - E_j^{44}, \quad C_j = E_j^{22} - E_j^{33}. \qquad (43)$$
Quantities of physical interest are expressed in this formalism as follows

$$ \mathrm{Tr}[\bar{\rho}_t \mathcal{O}] = \langle 1 \| \mathcal{O} \| \bar{\rho}_t \rangle, \qquad (44) $$

where

$$ \langle 1 \| = \otimes_{j=1}^{L} [j \langle 1 \| + j \langle 4 \|]. \qquad (45) $$

In the two-replica case we proceed analogously. Our starting point are two sets of fermion creation and annihilation operators corresponding to the two replicas that commute with one another, cf. (25). To proceed we define new fermions that fulfil pure anticommutation relations

$$ \bar{c}_{1;j} = c_{1;j}, \quad \bar{c}_{2;j} = c_{2;j}(-1)^{N_1}, \quad N_1 = \sum_{j=1}^{L} c_{1;j}. \qquad (46) $$

In the Hilbert space doubling formalism the Lindblad equation (12) for $\rho_t^{(2)}$ then takes the form

$$ \frac{d}{dt} \|\rho_t^{(2)}\rangle = \hat{\mathcal{L}}_2 \|\rho_t^{(2)}\rangle. \qquad (47) $$

Now the space associated with site j is 16-dimensional and a convenient basis is given by the graded tensor product of the states (36)

$$ |\alpha\rangle_j \otimes |\beta\rangle_j, \quad \alpha, \beta \in \{1, 2, 3, 4\}. \qquad (48) $$

The jump operators (13) in the Lindblad equation (12) for $\rho_t^{(2)}$ are

$$ \ell_{a;j}^{-} = (\ell_{a;j}^{+})^{\dagger} = \bar{c}_{a;j}^{\dagger} \bar{c}_{a;j+1}, \quad a = 1, 2. \qquad (49) $$

The Lindbladian is given by

$$ \hat{\mathcal{L}}_2 = J \sum_{j=1}^{L} \sum_{a=1}^{2} (\mathcal{L}_{a,j} + \mathcal{L}_{a,j}^{\mathrm{int}}), \qquad (50) $$

with

$$ \begin{aligned} \mathcal{L}_{a,j} &= p \left[ \ell_{a,j} \tilde{\ell}_{a,j} - \frac{1}{2} \left( \ell_{a,j}^{\dagger} \ell_{a,j} + \tilde{\ell}_{a,j}^{\dagger} \tilde{\ell}_{a,j} \right) \right] + q \left[ \ell_{a,j}^{\dagger} \tilde{\ell}_{a,j}^{\dagger} - \frac{1}{2} \left( \ell_{a,j} \ell_{a,j}^{\dagger} + \tilde{\ell}_{a,j} \tilde{\ell}_{a,j}^{\dagger} \right) \right], \\ \mathcal{L}_{a,j}^{\mathrm{int}} &= p \left[ \ell_{a,j} \tilde{\ell}_{a,j} - \frac{1}{2} \left( \ell_{a,j}^{\dagger} \ell_{a,j} + \tilde{\ell}_{a,j}^{\dagger} \ell_{a,j} \right) \right] + q \left[ \ell_{a,j}^{\dagger} \tilde{\ell}_{a,j}^{\dagger} - \frac{1}{2} \left( \ell_{a,j} \ell_{a,j}^{\dagger} + \tilde{\ell}_{a,j}^{\dagger} \ell_{a,j} \right) \right], \end{aligned} \qquad (51) $$

where we defined $\bar{1} = 2, \bar{2} = 1$. We now construct explicit matrix representations of $\hat{\mathcal{L}}_2$. We start by decomposing the 16 basis states on site j according to

$$ \begin{aligned} \|0\rangle_j &= |3\rangle_j \otimes |2\rangle_j, & \|1\rangle_j &= |1\rangle_j \otimes |1\rangle_j, & \|2\rangle_j &= |1\rangle_j \otimes |4\rangle_j, \\ \|3\rangle_j &= |4\rangle_j \otimes |1\rangle_j, & \|4\rangle_j &= |4\rangle_j \otimes |4\rangle_j, & \|5\rangle_j &= |2\rangle_j \otimes |3\rangle_j, \\ \|6\rangle_j &= |1\rangle_j \otimes |2\rangle_j, & \|7\rangle_j &= |2\rangle_j \otimes |1\rangle_j, & \|8\rangle_j &= |4\rangle_j \otimes |2\rangle_j, & \|9\rangle_j &= |2\rangle_j \otimes |4\rangle_j, \\ \|10\rangle_j &= |1\rangle_j \otimes |3\rangle_j, & \|11\rangle_j &= |3\rangle_j \otimes |1\rangle_j, & \|12\rangle_j &= |4\rangle_j \otimes |3\rangle_j, & \|13\rangle_j &= |3\rangle_j \otimes |4\rangle_j, \\ \|14\rangle_j &= |2\rangle_j \otimes |2\rangle_j, & \|15\rangle_j &= |3\rangle_j \otimes |3\rangle_j. \end{aligned} \qquad (53) $$

We then define Hubbard operators by

$$ \mathcal{E}_j^{a,b} = \|a\rangle_j \langle b\|, \qquad (54) $$
and operators $(G_j^{ab})^\dagger = G_j^{ba}$

$$
\begin{align*}
G_j^{12} &= \varepsilon_j^{1,3} + \varepsilon_j^{2,4} + \varepsilon_j^{6,8} + \varepsilon_j^{10,12}, \\
G_j^{13} &= \varepsilon_j^{0,4} + \varepsilon_j^{1,5} - \varepsilon_j^{6,9} - \varepsilon_j^{11,12}, \\
G_j^{14} &= \varepsilon_j^{0,3} - \varepsilon_j^{2,5} - \varepsilon_j^{6,7} + \varepsilon_j^{13,12}, \\
G_j^{23} &= -\varepsilon_j^{0,2} + \varepsilon_j^{3,5} - \varepsilon_j^{8,9} + \varepsilon_j^{11,10}, \\
G_j^{24} &= -\varepsilon_j^{0,1} - \varepsilon_j^{4,5} - \varepsilon_j^{8,7} - \varepsilon_j^{13,10}, \\
G_j^{34} &= \varepsilon_j^{2,1} + \varepsilon_j^{4,3} + \varepsilon_j^{9,7} + \varepsilon_j^{13,11}, \\
G_j^{11} &= \varepsilon_j^{0,0} + \varepsilon_j^{1,1} + \varepsilon_j^{2,2} + \varepsilon_j^{6,6} + \varepsilon_j^{10,10} + \varepsilon_j^{11,11} + \varepsilon_j^{13,13} + \varepsilon_j^{15,15}, \\
G_j^{22} &= \varepsilon_j^{0,0} + \varepsilon_j^{3,3} + \varepsilon_j^{4,4} + \varepsilon_j^{8,8} + \varepsilon_j^{11,11} + \varepsilon_j^{12,12} + \varepsilon_j^{13,13} + \varepsilon_j^{15,15}, \\
G_j^{33} &= \varepsilon_j^{2,2} + \varepsilon_j^{4,4} + \varepsilon_j^{5,5} + \varepsilon_j^{7,7} + \varepsilon_j^{10,10} + \varepsilon_j^{11,11} + \varepsilon_j^{12,12} + \varepsilon_j^{15,15}, \\
G_j^{44} &= \varepsilon_j^{1,1} + \varepsilon_j^{3,3} + \varepsilon_j^{5,5} + \varepsilon_j^{9,9} + \varepsilon_j^{10,10} + \varepsilon_j^{12,12} + \varepsilon_j^{13,13} + \varepsilon_j^{15,15}. \tag{55}
\end{align*}
$$

The operators $G_j^{ab}$ fulfil the gl(4) commutation relations

$$ [G_j^{ab}, G_j^{cd}] = \delta_{b,c} G_j^{ad} - \delta_{a,d} G_j^{ca}. \qquad (56) $$

In terms of (55) the two-replica Lindbladian becomes

$$ \hat{\mathcal{L}}_2 = J \sum_{j=1}^{L} p \overleftarrow{\mathcal{L}}_{j,j+1} + q \overrightarrow{\mathcal{L}}_{j,j+1} \equiv J p \overleftarrow{\mathcal{L}}_2 + J q \overrightarrow{\mathcal{L}}_2, \qquad (57) $$

where $\overleftarrow{\mathcal{L}}_2 = (\overrightarrow{\mathcal{L}}_2)^\dagger$ and

$$
\begin{align}
\overleftarrow{\mathcal{L}}_{j,j+1} = & G_j^{12} G_{j+1}^{21} + G_j^{43} G_{j+1}^{34} + G_j^{13} G_{j+1}^{31} + G_j^{42} G_{j+1}^{24} + \frac{1}{2} [G_j^{41} G_{j+1}^{14} + G_j^{32} G_{j+1}^{23} + \text{h.c.}] \nonumber \\
& + \frac{1}{2} \sum_{a=1}^{4} G_j^{aa} G_{j+1}^{aa} - \frac{1}{4} (C_j + C_{j+1} + 4). \tag{58}
\end{align}
$$

Here we have defined $C_j = \sum_{a=1}^4 G_j^{aa} - 2$, as above. This is in complete agreement with (21) once we account for the relation between the Lindbladian and its dual (11).

We note that the states (53) transform under gl(4) as follows:

* $\{\|a\rangle_j|a = 0, \dots, 5\}$ form the six-dimensional representation [6] of gl(4) and have $C_j = 0$;

* $\{\|a\rangle_j|a = 6, \dots, 9\}$ form the fundamental representation [4] of gl(4) and have $C_j = 1$;

* $\{\|a\rangle_j|a = 10, \dots, 13\}$ form the fundamental representation [4] of gl(4) and have $C_j = -1$;

* $\{\|14\rangle_j\}$ and $\{\|15\rangle_j\}$ form one-dimensional representations of gl(4) and have $C_j = \mp 2$.

The quantities of physical interest are expressed as follows

$$ \mathrm{Tr}[\bar{\rho}_t^{(2)} \mathcal{O}_1 \otimes \mathcal{O}_2] = \langle 1_2 \| \mathcal{O}_1 \otimes \mathcal{O}_2 \| \bar{\rho}_t^{(2)} \rangle, \qquad (59) $$

where

$$ \langle 1_2 \| = \otimes_{j=1}^{L} [j \langle 1 \| + j \langle 2 \| + j \langle 3 \| + j \langle 4 \|]. \qquad (60) $$
# 3 U(1) symmetries, sectors and fragmentation

## 3.1 Symmetries

The Q-ASEP & Q-SSEP models possess $u(1)$ symmetries, one for each site, which echo the fact that two pairs of complex Brownian motions ($W_t, \overline{W}_t$) and ($\tilde{W}_t, \overline{\tilde{W}}_t$) differing by a phase, $\tilde{W}_t = e^{i\theta} W_t$, have identical distributions.

For any realization of the Brownian motions $W_t^j$, let us consider the representation of $U(1)^L$ on the physical Hilbert space generated by the local particle numbers $n_j$. That is: elements of $U(1)^L$ are represented by $e^{i\sum_j \theta_j n_j}$ with arbitrary parameter $\theta_j$. The two Hamiltonian increments, $dH_t$ from eqn (1), and $d\tilde{H}_t$, obtained by conjugation with these $U(1)$'s, with

$$d\tilde{H}_t := e^{-i\sum_j \theta_j n_j} dH_t e^{i\sum_j \theta_j n_j},$$

have identical distributions. This holds because $d\tilde{H}_t$ coincides with eqn (1) but with Brownian increments $d\tilde{W}_t^j = e^{i(\theta_j - \theta_{j+1})} dW_t^j$ whose distribution is identical to that of $dW_t^j$. As a consequence, the two density matrices $\rho_t$ and $\tilde{\rho}_t = e^{-i\sum_j \theta_j n_j} \rho_t e^{i\sum_j \theta_j n_j}$, time evolved with the Hamiltonian increments $dH_t$ and $d\tilde{H}_t$, respectively, have identical distributions. This property holds for Q-ASEP and Q-SSEP as well.

This implies that, for any $R$ replicas, acting with the Lindbladian or conjugating with the total local particle numbers are commuting operations:

$$[\mathcal{L}, C_j] = 0, \quad [\mathcal{L}^*, C_j] = 0, \qquad (61)$$

for all $j = 1, \dots, L$ with $C_j = \sum_a C_{a; j}$. Recall that the $C_{a; j}$'s act on any operators as $C_{a; j} O = [n_{a; j}, O]$ with $n_{a; j} = c_{a; j}^\dagger c_{a; j}$ the local number operator in the $a$-th replica. They count the number of fermion insertions in the operators $O$. Since the latter are made of product of the local fermionic operators, the spectrum of $C_j$ runs from $-R$ to $R$:

$$C_j = -R, -R+1, \dots, 0, \dots, R-1, R.$$

For instance, the operator $\prod_a c_{a; j}^\dagger$ (resp. $\prod_a c_{a; j}$) is only the operator with $C_j$-charge $+R$ (resp. $-R$).

## 3.2 Sectors and fragmentation

Since the $C_j$'s commute with the Lindbladian, their eigenvalues are good quantum numbers preserved by the dynamics. Phrased differently, the Q-ASEP & Q-SSEP dynamics only involves mixing between operators with identical $C$-charges. Note that one $C$-charge is assigned per site. We denote the set of $C$-charges by $\{c_j | j = 1, \dots, L\}$ by $c$.

We call the subspace of the operators of given $C$-charges a sector. Each sector is labelled by the corresponding set $c$. The Q-ASEP & Q-SSEP dynamics takes place within given sector and does not induce transitions between different sectors. In particular, since the extremal $C$-charges $\pm R$ are non-degenerate, the dynamics inside those sectors is trivial in the sense that operators with extremal $C$ charge do not move. The splitting of the Q-ASEP & Q-SSEP dynamics into sectors has been termed “fragmentation” in Ref. [22].

Let us analyze the operator content of the different sectors depending on the number of replicas. For a single replica, the different values of the $C$-charge are $\pm 1$ and 0 with the following content, at each given site $j$ (with $n_j = c_j^\dagger c_j$):

$$C_j^{[R=1]} = 0 : \mathbb{I}_j, n_j, \rightsquigarrow [2], \\ C_j^{[R=1]} = \pm 1 : c_j^\dagger (c_j), \rightsquigarrow [1], ([\bar{1}]).$$
In order to describe time evolution in a given sector it is useful to define projection operators onto the different subspaces at a given site

$$P_j^{[\alpha]}, \quad \alpha = 0, \pm 1. \qquad (62)$$

In terms of the explicit representation (36) these read

$$P_j^{[0]} = |1\rangle\langle 1| + |4\rangle\langle 4|, \quad P_j^{[1]} = |2\rangle\langle 2|, \quad P_j^{[-1]} = |3\rangle\langle 3|. \qquad (63)$$

The Lindbladian in a given sector c is then given by

$$\mathcal{L}_c = \left( \prod_{j=1}^{L} P_j^{[c_j]} \right) \mathcal{L} \left( \prod_{k=1}^{L} P_j^{[c_k]} \right). \qquad (64)$$

The time evolution of the single-particle Green's function $\text{Tr}[\rho(t)c_j^\dagger c_k]$, with $k > j + 1$, occurs in sector

$$c_1 = (0, \dots, 0, \underbrace{1}_{j}, 0, \dots, 0, \underbrace{-1}_{k}, 0, \dots, 0). \qquad (65)$$

The corresponding Lindbladian is (with periodic boundary conditions)

$$\mathcal{L}_{c_1} = P_j^{[1]} + \mathcal{L}_{[j+1,k-1]} + P_k^{[-1]} + \mathcal{L}_{[k+1,j-1]}, \qquad (66)$$

where [22]

$$\mathcal{L}_{[1,l]} = \sum_{j=1}^{l-1} p \left[ E_j^{14} E_{j+1}^{41} - E_j^{44} E_{j+1}^{11} \right] + q \left[ E_j^{41} E_{j+1}^{14} - E_j^{11} E_{j+1}^{44} \right] \\ - \frac{q}{2} E_l^{11} - \frac{p}{2} E_l^{44} - \frac{p}{2} E_1^{11} - \frac{q}{2} E_1^{44}. \qquad (67)$$

The dynamics described by (66) is that of two stationary “defects” at sites j and k, which affect the dynamics in the spatial regions between them through the boundary terms in (67). The dynamics is integrable [22], but so far has only been solved for the case of the SSEP.

For two replicas, the possible C-charges are $\pm 2, \pm 1$ and 0 with contents (with the convention, $c^+ \equiv c^\dagger, c^- \equiv c, n = c^\dagger c$):

$$C_j^{[R=2]} = 0 : \mathbb{I}_j \otimes \mathbb{I}_j, \mathbb{I}_j \otimes n_j, n_j \otimes \mathbb{I}_j, n_j \otimes n_j, c_j^+ \otimes c_j^-, c_j^- \otimes c_j^+ \rightsquigarrow [6]$$

$$C_j^{[R=2]} = \pm 1 : \mathbb{I}_j \otimes c_j^\pm, c_j^\pm \otimes n_j, n_j \otimes c_j^\pm, c_j^\pm \otimes \mathbb{I}_j, \rightsquigarrow [4], ([\bar{4}])$$

$$C_j^{[R=2]} = \pm 2 : c_j^\pm \otimes c_j^\pm, \rightsquigarrow [1], ([\bar{1}]).$$

The operators in the $c_0 = (0, \dots, 0)$ sector form an algebra (because the product of two such operators still has C charge zero). There is an action of this algebra on the other sectors, or more precisely the other sectors form modules of the $c_0$ algebra under multiplication, because the C-charge is additive under products.

The above numbers [·] refer to the multiplicity of the corresponding sectors. They have a group theoretical interpretation. The space of local operators at site j is isomorphic to that of the Fock space of the $2R$ effective fermions $f_j^{A\dagger}$. Its dimension $2^{2R} = 2^R \cdot 2^R$ indeed coincides with that of the space of endomorphisms on the local physical Hilbert space (generated by the R physical fermions $c_{a:j}^\dagger$ at site j). As explained above, there is an action of $gl(2R)$ on the operators localised at site j whose generators are the super-operators $G_j^{AB} = f_j^{A\dagger} f_j^B$, see eqn (30). The algebra $gl(2R)$ has a central $u(1)$, generated by the C-charge $C_j = \sum_A f_j^{A\dagger} f_j^A - R$,
# 4 Integrability and non-integrability

## 4.1 Spin chain identification

Thanks to eqn (22) or eqn (24), $\mathcal{L} \propto \sum_j (\sum_{A,B} J_{j+1}^{AB} J_j^{BA} + \frac{1}{2R} C_{j+1} C_j - \frac{R}{2})$, the Q-SSEP Lindblad dynamics has been re-written as a spin chain dynamics. In the one replicas case $R = 1$, the spin chain is equivalent to the isotropic XXX Heisenberg model and thus integrable [58]. The purpose of this section is to determine whether the $R = 2$ replica dynamics is integrable. The structure of the Lindbladian density (69) is different from that of integrable spin chains obtained by varying the representations of the symmetry algebra acting on a given site, see. e.g. [59–62] in that the latter constructions generically lead to three-"spin" interactions, while in our case the Lindbladian density only involves two sites. This suggests that the dynamics in a general sector $c$ is not integrable in a standard Yang-Baxter fashion. We therefore focus on the $c_{\pm 1}$ and $c_0$ sectors.

*Integrability in the $c_{\pm 1}$ sector.*

The quantity $\sum_{A,B} J_{j+1}^{AB} J_j^{BA}$ is the tensor Casimir acting on the tensor product of the $sl(2R)$ representations located at site $j$ and $j + 1$. By $sl(2R)$ invariance, its diagonalisation is obtained by decomposing this tensor product into irreducible $sl(2R)$ sub-representations. As a consequence, $\sum_{A,B} J_{j+1}^{AB} J_j^{BA}$ can be written as a linear combination of the projectors on these irreducible sub-representations.

In the $c_{\pm 1}$ sector, for two replicas, the representation on each site is the $sl(4)$ vector representation $\square = [4]$. Its tensor product with itself, $\square \otimes \square$, decomposes into two irreducible sub-representations, the symmetric and anti-symmetric tensors. The associated projectors are $\frac{1}{2}(1 \pm P_{j;j+1})$ with $P_{j;j+1}$ the permutation operator. As a consequence, the Q-SSEP Lindbladian in the $c_0$ sector reads

$$ \mathcal{L}^{[C=\pm 1, R=2]} \propto \sum_j P_{j;j+1} + \text{const.} $$

It is the $sl(2R)$ version of the isotropic Heisenberg spin chain, and it is known to be integrable [63]. In other words, the dynamics of Q-SSEP quadratic fluctuations in the $c \pm 1$ sector is integrable.

*Integrability in the $c_0$ sector?*

The formula eqn (24) also holds in the $c_0$ sector. In this sector and for two replicas, the on-site representation is the [6] representation of $sl(4)$ or $so(6)$, that of rank 4 anti-symmetric tensors of $sl(4)$ or that of 6-dimensional vectors of $so(6)$. We have the decomposition rule, $[6] \otimes [6] = [1] + [15] + [20]$ as $sl(4)$ or $so(6)$ modules, with [15] the rank two antisymmetric $so(6)$ tensors and [20] the rank two traceless symmetric $so(6)$ tensors. The projectors on these representations are (with $d = 6$)

$$ P_A = \frac{1}{2}(1-P), \quad P_S = \frac{1}{2}(1+P) - \frac{1}{d}Q, \quad P_\bullet = \frac{1}{d}Q, $$

with $P$ the permutation operator and $Q$ the so-called trace operator. Thus $P$, $Q$ and the identity are the only $sl(4) \equiv so(6)$ intertwiners in $[6] \otimes [6]$. As a consequence, in the $c_0$ sector or for two replicas, $\sum_{A,B} J_{j+1}^{AB} J_j^{BA}$ is a linear combination of $P_{j;j+1}$, $Q_{j;j+1}$, up to an additive constant. Computing explicitly the tensor Casimir, see the end of Appendix B, we find

$$ \mathcal{L}^{[C=0, R=2]} \propto \sum_j (P_{j;j+1} - Q_{j;j+1}) + \text{const.} $$

The Hamiltonian of the known $sl(4) \equiv so(6)$ integrable model [64] is $H = \sum_j (2P_{j;j+1} - Q_{j;j+1})$. Note the differences in the coefficients.
Surprisingly, and despite this absence of direct connection to a known integrable models, we found numerical evidences that the dynamics in the $c_0$ sector could be integrable. It therefore remains a challenge to analytically prove that the $c_0$ dynamics is integrable and to decipher the underlying structures responsible for this property.

A similar analysis applies to higher number of replicas. In particular, because they associated to the $sl(2R)$ vector representations, the dynamics in the sectors with all $C$ equal to $\pm(R-1)$ are always integrable in the usual sense (they are mapped to the $sl(2R)$ analogues of the isotropic Heisenberg spin chain). Whether the other sectors are integrable is an open question.

<figure><img src="image_4.png" /><figcaption>(a) Spacings $s_n = |e_n - e_{n-1}|$ of adjacent eigenvalues</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Spacings $\hat{s}_n = |\hat{e}_n - \hat{e}_{n-1}|$ of unfolded eigenvalues $\hat{e}_n$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) Ratio $r_n = s_n/s_{n-1}$ of adjacent spacings</figcaption></figure>

**Figure 1:** Histograms for the spacings and ratios of the $N = 1077$ eigenvalues of the $c_1$ Lindbladian on $L = 11$ sites after the degeneracies have been removed. The chosen symmetry sector is defined by $k = 2\pi/11$ ($e^{ik}$ is the eigenvalue under translation by one site) and $(n_1, n_2, n_3, n_4) = (1, 2, 3, 5)$ where $n_i$ is the number of times the local state $i = 1, 2, 3, 4$ appears in the tensor-product-state (the $c_1$ sector has local dimension 4 on each site). This corresponds to the Cartan charges $(J_1^z, J_2^z, \frac{C_1-C_2}{2}) = (-2, -3, \frac{3}{2})$. The average ratio of consecutive spacings is $\langle \tilde{r} \rangle = 0.3826$.

## 4.2 Numerics

Here we investigate the integrability of the two-replica Lindbladian numerically in sectors $c_0$ and $c_1$ by considering the level-spacing statistics of its eigenvalues. Most studies that have been conducted on the level-spacing statistics of integrable spin chains deal with the case where the global symmetry algebra is $u(1)$ or $sl(2)$. In contrast to this, the symmetry algebra of our spin chain model is of higher rank, namely $sl(4)$. This leads to a new source of degeneracies of eigenvalues, which makes the level statistics differ from the usual $sl(2)$ case at first sight. To our knowledge, this is the first time this problem has been addressed in the literature.

Our analysis of the level-spacing statistics suggests that the $c_{\pm 1}$ sectors are indeed integrable, while for the $c_0$ sector it provides weak signs of integrability breaking, which sometimes lack an overall consistency. A key problem in the analysis is that the effects of integrability breaking perturbations become visible only for sufficiently large system sizes [65, 66]. We suspect that the maximal system size we could achieve in practice is simply too small to get a completely consistent picture. Moreover, the presence of a higher rank symmetry algebra may further obscure signs of integrability breaking in the level-spacing statistics.

Let us start by recalling some known results on the level-spacing statistics of integrable and non-integrable models. Integrable Hamiltonians possess the very particular property that their eigenvalues are i.i.d. random variables, as if the Hamiltonian was just a random diagonal matrix. This was first conjectured by Berry and Tabor [67] and has been confirmed in many explicit examples such as the XXX Heisenberg chain [68, 69]. Importantly, the spacing
<figure><img src="image_7.png" /><figcaption>Figure 2: Histograms of the spacings between consecutive raw eigenvalues (first row) and unfolded eigenvalues (second row) for a deformation of the Q-SSEP Linbladian for quadratic fluctuations on L = 10 sites in the c₀ sector. The deformations are parametrised by g and defined as $\mathcal{L}(g) \propto \sum_j (P_{j;j+1} - gQ_{j;j+1} - 1)$ such that g = 1 corresponds to the original Q-SSEP case. The chosen symmetry sector is defined by k = 2π/10 and Cartan charges $(J_1^z, J_2^z, \frac{C_1-C_2}{2}) = (5, 3, 3)$. The value for N provides the number of eigenvalues left after removing the degeneracies.</figcaption></figure>

$s_n = e_n - e_{n-1}$ between adjacent eigenvalues follows an exponential distribution

$$p(s) = e^{-s}. \qquad (71)$$

To be precise, this holds only for the so-called "unfolded spectrum" of the Hamiltonian, where one performs a local change of variable on the eigenvalues $e_n$ such that the density of the new variables $\hat{e}_n$ is uniform (see Appendix D). Instead, as showed in [70], one can also consider the ratio of consecutive spacings $r_n = s_n/s_{n-1}$ whose distribution is independent of the local density of eigenvalues and is given by

$$p(r) = \frac{1}{1 + r^2}. \qquad (72)$$

In contrast to integrable Hamiltonians the eigenvalues of a generic Hamiltonian - a random Matrix - tend to repel each other. The spacing between eigenvalues of a 2 × 2 random matrix in the GOE (Gaussian Orthogonal Ensemble) - which would be the appropriate ensemble to deal with since the Q-SSEP Linbladian is symmetric and real - has a probability distribution known as *Wigners surmise*

$$p(s) = \frac{\pi s}{2} e^{-\pi s^2/4}. \qquad (73)$$

This turns out to be good approximation also for the level-spacing of large GOE random matrices. In particular, there is a zero probability to find consecutive eigenvalues with spacing zero. The same is true if instead of the spacing one again considers the ratio of adjacent spacings r, which behaves as $p(r) \sim r^\beta$ for small r, where β is the Dyson index of the matrix ensemble [70].

The ratio of adjacent spacings has the nice property, that the average of

$$\tilde{r} = \min(r, 1/r) \in [0, 1] \qquad (74)$$
over the given ensemble takes is a constant, and can therefore be used to classify the en-
semble to which some numerical distribution might belongs. For Poisson random variables<fnref n="3" />,
$\langle \tilde{r} \rangle = 2 \ln 2 - 1 \approx 0.3863$, while for the eigenvalues of a matrix in the GOE, $\langle \tilde{r} \rangle \approx 0.5359$.

Before discussing the results, let us also comment on the reduction of the Lindbladian
to its remaining symmetry sectors. The eigenvalues in each symmetry sector are statistically
independent and therefore one should treat each sector independently. In practice, we bring
the Lindbladian to block diagonal form with respect to all its mutually commuting symmetries
$I_i, [\mathcal{L}, I_i] = 0, [I_i, I_j] = 0$. After fixing to $\mathbf{c}_{\pm 1}$ or $\mathbf{c}_0$, the maximally set of commuting symmetries
are translation $T$, the three Cartan elements ($J_1^z = \sum_j J_{1,j}^z, J_2^z = \frac{c_1-c_2}{2}$) of $gl(4)$ (which are the
analogues of the magnetization for $sl(2)$) and depending on the choice of the three Cartan
elements, a permutation $F$ of the states (which is the generalization of a spin flip in the $m=0$
sector for $sl(2)$). However, once all these charges have been fixed, there is still a degeneracy
in the eigenvalues of the Lindbladian left, which would lead to an artificial large peak at zero
in the level-spacing statistic. This is because the Cartan subalgebra for $gl(4)$ consists of more
than one element, hence there are more than one "lowering-operator" and therefore weight-
spaces in an irreducible $gl(4)$ representation can be more than one-dimensional. But the $gl(4)$
symmetry of the Lindbladian ensures that all states in an irreducible $gl(4)$ representation have
the same eigenvalue and hence, selecting a weight-space (i.e. fixing the Cartan charges) will
not lift all degeneracies (as it would do for the $sl(2)$ spin chain). We therefore manually
deleted all the degenerate copies of eigenvalues from the complete set of eigenvalues in a
given symmetry sector and analysed the level-spacing and ratio statistics for the remaining
eigenvalues. The procedure is not entirely correct, because there can also be degeneracies in
the spectrum solely due to integrability, which would be neglected in our procedure. But it
turns out that for the overall statistics this only plays a minor role.

Fig. 1 shows the results for the **c**₁ sector and both the shape of distribution and the value for ⟨**r**⟩ suggest that this model is integrable. The results for the **c**₀ sectors are less clear. In some of the sectors with fixed Cartan elements there is no sign of integrability breaking. In others, such as in Fig. 2 there are weak signs. Increasing the perturbation of the known integrable model with the trace operator **Q**, parametrised by the coupling strength **g**, one observes a gradual deviation from the exponential distribution, which is correctly reproduced for the integrable case (g = 1/2). However, a real Wigner-like distribution is not visible, even for higher values of **g**, and we suspect this to be due to the small system size we are able to achieve in practice. Also note, that the value for ⟨**r**⟩ seems to suggest integrability breaking for the Q-SSEP case. However, in our example, it does not consistently increase with **g** as one would expect: Value for **g** = 2 is lower than that for **g** = 1. Again we think this is due to the limited system size. Finally, it might be interesting to look at the value of *N*, which describes the number of eigenvalues left after removing the degeneracies. For the known integrable case (a), *N* = 1205 is lower than for the other two cases (b) and (c) where the number almost coincides, *N* = 1335 and *N* = 1336, respectively<fnref n="4" /> This hints, that in the integrable case (a), there were more degenerate eigenvalues than the gl(4) symmetry with its higher dimensional weight spaces could explain. We think that these additional degeneracies are probably a result of integrability - or reversely, their absence (as in (b) and (c)) is a sign of integrability breaking. To sum up, all these observations suggest that the quadratic fluctuations for Q-SSEP are not integrable in the **c**₀ sector.

<footnote n="3">³A short derivation of this fact: Since the probability of consecutive spacings $s_1$ followed by $s_2$ is equivalent the inverse order, i.e. $P(s_1, s_2) = P(s_2, s_1)$, it follows that $r = s_1/s_2$ and $1/r = s_2/s_1$ have the same distributions. Therefore $R(\tilde{r}) = 2p(r)\Theta(1-r)$, from which the average value $\langle \tilde{r} \rangle = 2\ln 2 - 1$ can be computed.</footnote>

<footnote n="4">⁴The slight difference could arise due to the numerical imprecision.</footnote>
# 5 Steady states and slow modes

We now look at the steady states of the Q-SSEP. The Q-SSEP invariant measure has been described in [41] using a mapping to random matrix theory. We here adopt an alternative approach which yields a better access to their structure in terms of the underlying $gl(2R)$ symmetry. The case of the Q-ASEP steady states will be discussed elsewhere [71].

## 5.1 Q-SSEP steady states

The steady states of the Q-SSEP dynamics are the zero modes of the Q-SSEP Lindbladian $\mathcal{L}$. Since $\mathcal{L}$ commutes the $gl(2R)$ action, these modes form $gl(2R)$ multiplets. Recall that the Q-SSEP Lindbladian is self-dual, $\mathcal{L} = \mathcal{L}^*$.

We first study the zero modes in the $c_0$ sector (i.e. with all charges $C_j$ along the chain set to zero). The identity belongs to this sector and is an eigenstate with eigenvalue zero since $\mathcal{L}^*\mathbb{I} = 0$ by definition and $\mathcal{L} = \mathcal{L}^*$. The identity is a factorised operator $\mathbb{I} = \mathbb{I}_1\mathbb{I}_2\cdots\mathbb{I}_L$ with $\mathbb{I}_j$ the identity on the $j$-th site. As a consequence of the $sl(2R)$ symmetry, all factorised operators $x_1x_2\cdots x_L$ with $x_j$ localised on site $j$ and with $x$ in the $sl(2R)$ orbit of the identity are also zero modes. This orbit is isomorphic to the middle $sl(2R)$ fundamental representation $[\mu_R]$ and the $sl(2R)$ action on this representation is faithful. As a consequence, all factorised operators of the form

$$ \|\Omega_x^0\rangle := x_1x_2\cdots x_L, \quad x \in [\mu_R], \qquad (75) $$

are zero eigenstates of the Q-SSEP Lindbladian in the $c_0$ sector: $\mathcal{L}\|\Omega_x\rangle = 0$. Of course, linear combinations of such states are also steady states. They form an irreducible representation of $sl(2R)$ with highest weight $L\mu_R$ (with a slight abuse of notation consisting in denoting by $\mu_k$ the highest weight of the representation $[\mu_k]$). Indeed, the representation formed by the above operators $\|\Omega_x^0\rangle$, $x \in [\mu_R]$, is irreducible with the highest weight vector $\|\Omega_{\mu_R}^0\rangle$ equal to the tensor product of $L$ identical highest weight vectors, each of weight $\mu_R$. The Young tableau of the representation of highest weight $L\mu_R$ is a rectangular tableau with $R$ lines and $L$ columns. Its dimension can be computed using the so-called hook rule (see below for examples).

As detailed in the Appendix B, this remark holds true in any sector. If $|\mu_{R+k}\rangle$ is the highest weight vector of the representation $[\mu_{R+k}]$ corresponding to local operator with $C$ charge equals to $k$, then $|\mu_{R+k}\rangle \otimes \cdots \otimes |\mu_{R+k}\rangle$ is a zero mode of the Q-SSEP Lindbladian. By $sl(2R)$ invariance, this implies that all states,

$$ \|\Omega_x^k\rangle := x_1x_2\cdots x_L, \quad x \in [\mu_{R+k}], \qquad (76) $$

are zero eigenstates of the Q-SSEP Lindbladian in the $c_k$ sector. Linear combinations of such states are also zero modes. They form the $sl(2R)$ representation whose Young tableau is rectangular with $R+k$ lines and $L$ columns.

Let us give examples for one or two replicas. For $R=1$ replica, the structure of the zero modes as $gl(2)$ modules is: $c_0 \rightsquigarrow [L+1]$ (spin $L/2$) and $c_{\pm 1} \rightsquigarrow [1]$ (spin 0).

For $R=2$ replicas, the structure of the zero modes as $gl(4)$ modules is: $c_0 \rightsquigarrow [L\mu_2]$ (Young tableau with 2 lines of length $L$), $c_{\pm 1} \rightsquigarrow [L\mu_1]$ (Young tableau with 1 line of length $L$) and $c_{\pm 2} \rightsquigarrow [L\mu_0] = [1]$ the trivial representation. The number of zero modes are:

$$ c_0 \rightsquigarrow \frac{(L+1)(L+2)^2(L+3)}{12}, \quad c_{\pm 1} \rightsquigarrow \frac{(L+1)(L+2)(L+3)}{6}, \quad c_{\pm 2} \rightsquigarrow 1. $$

We have checked that these zero modes form all possible Q-SSEP steady states for $L=1$ and $L=2$, and we conjecture that this is true for all $L$. It is known from [41] that the Q-SSEP steady distribution is $SL(L)$ invariant, so that there is a peculiar (nice) interplay between the $SL(L)$ and $GL(2R)$ perspectives on the zero modes. See Appendix C for details.
All Q-SSEP steady operators have homogeneous C-charge along the chain. In contrast, a domain wall operator connecting sectors with different C charges $c_1$ and $c_2$, say of the form $|\mu_{R+c_1}\rangle^{\otimes M} \otimes |\mu_{R+c_2}\rangle^{\otimes (L-M)}$ and their $gl(2R)$ orbits, are eigenstates of the Q-SSEP Lindbladian with eigenvalues $-\frac{J}{2}|c_1 - c_2|$. Hence, domain walls decay exponentially fast in time, cf. [22]. More generally, for an inhomogeneous sector with domain walls at a series of edges $w_1, \dots, w_d$ connecting sectors with charges $c_\alpha^-, c_\alpha^+(c_\alpha^- \neq c_\alpha^+)$, $\alpha = 1, \dots, d$, the gap is the sum of the gaps associated to each of these domain walls, i.e. $-\frac{J}{2}\sum_\alpha |c_\alpha^- - c_\alpha^+|$. A proof of this statement is given in Appendix B.

## 5.2 Q-SSEP low lying magnon states

Let us first consider the one replica case ($R = 1$). In the $c_0$ sector the zero modes form a spin $L/2$ representation of $sl(2R)$. The other zero modes are $sl(2R)$ singlets. With respect to the $sl(2R)$ basis, $J^\pm = \sum_j J_j^\pm$, $J^z = \sum_j J_j^z$, the highest weight vector of this spin $L/2$ representation is the operator $\|\Omega_n^1\rangle := n_1 n_2 \cdots n_L$ with $n_j$ the number operator at site $j$. Since the one replica Lindbladian has been identified with the XXX Heisenberg spin chain Hamiltonian, one may construct magnon excitations by reversing the polarization of this fully polarized state in a coherent way. Algebraically, this is done by acting with the super-operator $\sum_j e^{ip_j} J_j^-$ depending on a momentum $p$. It yields the eigenstates

$$ \sum_j e^{ip_j} n_1 n_2 \cdots (J_j^- n_j) \cdots n_L, \qquad (77) $$

with $J_j^- n_j = I_j - n_j$. Their eigenvalues are $\omega_p := J(\cos p - 1)$. At small momenta, the dispersion relation is diffusive, $\omega_p \approx -J p^2$, as expected.

This observation generalizes to any number of replicas. Since the zero modes $\|\Omega_x^k\rangle = x_1 x_2 \cdots x_L$ with $x \in [\mu_{R+k}]$, given in eqn (76) with R the number of replicas and k the C charge, are also ferromagnetic like states, one can construct one-magnon eigenstates. The latter are obtained by considering the fully polarized state $\|\Omega_{\mu_{R+k}}\rangle$ with x chosen to be the highest weight vector $\mu_{R+k}$ of the $sl(2R)$ representation $[\mu_{R+k}]$ and by acting on it with the appropriate $gl(2R)$ generator. Namely, one considers the operators

$$ \|p; \gamma; \mu_{R+k}\rangle := \sum_j e^{ip_j} G_j^{-\gamma} \|\Omega_{\mu_{R+k}}\rangle, \qquad (78) $$

with $\gamma$ a positive root of $sl(2R)$ and $G^{-\gamma}$ the associated root generator. As shown in Appendix B, these states are eigenstates of the R replica Lindbladian iff $\gamma$ is one of the simple roots of $sl(2R)$. There are $2R - 1$ such simple roots but only one of them (dual to the weight $\mu_{R+k}$) yields a non vanishing state so that the labelling $\gamma$ is actually redundant since it is linked to $\mu_{R+k}$. The dispersion relation is also diffusive:

$$ \omega_p = 2J(\cos p - 1) \approx -J p^2. \qquad (79) $$

By $sl(2R)$ invariance, any state obtained by acting with $SL(2R)$ on $\|p; \gamma; \mu_{R+k}\rangle$ is an eigenstate. That is: all states $\otimes_j g_j \|p; \alpha; \mu_{R+k}\rangle$, with $g \in SL(2R)$ and $\alpha$ a simple root, are eigenstates. We have $\otimes_j g_j \|p; \alpha; \mu_n\rangle = \|p; g \alpha g^{-1}; x\rangle$ with $x = g \mu_{R+k}$, they form an irreducible $sl(2R)$ representation isomorphic to $[\mu_{R+k}]$.

## 5.3 Cat state preparation in the non trivial integrable sector

In order to demystify the meaning of different sectors and demonstrate how integrability could be used for addressing interesting physical questions, we discuss the dynamics of the cat states, which partially belong to the integrable $c_{\pm 1}$ subsectors for two replicas.
Imagine initially preparing the system in the pure state $|\psi\rangle$. Its initial $R = 2$ replicated density matrix is $\rho_0^{(2)} = \rho_0 \otimes \rho_0$, with $\rho_0 = |\psi\rangle\langle\psi|$.

In order for $\rho^{(2)}$ to have a non zero component on the $c_{\pm 1}$ sectors, one has to consider $|\psi\rangle$ to be a cat state, i.e. a particular linear combination of states with macroscopically different numbers of particles. The simplest state $|\psi\rangle$ to consider would be the sum $|\psi_0\rangle + |\psi_M\rangle$ with $|\psi_0\rangle$ (resp. $|\psi_M\rangle$) a state with 0 (resp. M) particles. However, this is actually not enough for producing quadratic fluctuations in the $c_{\pm 1}$ sector. One has to consider three position cat states $|\psi\rangle$ which are the sum of three states with different particle numbers of the form

$$|\psi\rangle = |\psi_0\rangle + |\psi_M\rangle + |\psi_{L-M}\rangle \qquad (80)$$

with $|\psi_0\rangle$, $|\psi_M\rangle$ and $|\psi_{L-M}\rangle$ having 0, M and $L-M$ particles respectively. We assume $0 < M < L/2$. (We could also consider more general cat states but this is the simplest with the required properties).

When considering quadratic fluctuations, we have to consider the initial tensor product $|\psi\rangle \otimes |\psi\rangle$. To ensure that the quadratic fluctuation dynamics survives the large time limit, we have to decipher on which $SU(L)$ representations this product is decomposed onto in order to verify that it has a non trivial component on $c_{\pm 1}$ steady states. See Appendix C for details.

The state $|\psi_0\rangle$ belongs to the trivial $SU(L)$ representation $\Lambda_0$, while the states $|\psi_M\rangle$ (resp. $|\psi_{L-M}\rangle$) belong the fundamental representation $\Lambda_M$ (resp. $\Lambda_{L-M}$), with Young tableau with one column and M (resp. $L-M$) boxes. The tensor product $|\psi\rangle \otimes |\psi\rangle$ thus contains $|\psi_0\rangle \otimes |\psi_0\rangle$ which belongs to the scalar representation $\Lambda_0 \otimes \Lambda_0 \equiv \Lambda_0$. It also contains the products $|\psi_M\rangle \otimes |\psi_{L-M}\rangle$ (resp. $|\psi_{L-M}\rangle \otimes |\psi_M\rangle$) which are elements of $\Lambda_M \otimes \Lambda_{L-M}$ (resp. $\Lambda_{L-M} \otimes \Lambda_M$). The decomposition of the latter tensor product of $SU(L)$ representations contains the fully antisymmetric representation $\Lambda_0^+$ which, as $SU(L)$ module, is isomorphic to $\Lambda^0$. However, $\Lambda_0^+$ and $\Lambda_0$ are isomorphic as $SU(L)$ modules but not as $U(L)$ modules since they differ by their C charges.

As a consequence, the initial quadratic density matrix $\rho^{(2)} := |\psi\rangle\langle\psi| \otimes |\psi\rangle\langle\psi|$ contains blocks intertwining $\Lambda_0 \to \Lambda_0^+$ and $\Lambda_0^+ \to \Lambda_0$. These blocks are in the $c_{\pm 1}$ sectors. They are sub-blocks in the matrices $|\psi_0\rangle\langle\psi_M| \otimes |\psi_0\rangle\langle\psi_{L-M}|$ or $|\psi_0\rangle\langle\psi_{L-M}| \otimes |\psi_0\rangle\langle\psi_M|$ and their Hermitian conjugates.

To check that they have a non trivial projection on the $c_{\pm 1}$ steady states, let us compute their overlaps with the $c_{\pm 1}$ steady operator (that is: operators in the kernel of the dual Lindbladian in the sector $c_{\pm 1}$). These operators, say with $c_{-1}$, have exactly one fermionic annihilation operator per site in either of the two replicas. The simplest such operator is of the form

$$\mathfrak{C} := \prod_{m \in \mathfrak{M}} c_{1;m} \cdot \prod_{n \in \mathfrak{M}^c} c_{2;n},$$

with $\mathfrak{M}$ a sub-set of $\{1, \dots, L\}$ and $\mathfrak{M}^c$ its complement. Any other operator in the $c_{\pm 1}$ sector is obtained by acting on the latter with $c_0$ operators (which amounts to multiply $\mathfrak{C}$ by the number operators $n_{a;j}$ at different sites or replicas). Let us choose $\text{Card}\,\mathfrak{M} = M$ and hence $\text{Card}\,\mathfrak{M} = L - M$. Then,

$$\text{Tr}(\rho^{(2)} \mathfrak{C}) = \langle \psi_0 | \prod_{m \in \mathfrak{M}} c_m | \psi_M \rangle \langle \psi_0 | \prod_{n \in \mathfrak{M}^c} c_n | \psi_{L-M} \rangle \neq 0.$$

This is non zero if the state $|\psi_M\rangle$ (resp. $|\psi_{L-M}\rangle$) has particles inserted at the site selected by $\mathfrak{M}$ (resp. $\mathfrak{M}^c$). The same holds true if we multiply $\mathfrak{C}$ by an operator in the $c_0$ sector.

Hence, the three position cat states, $|\psi\rangle = |\psi_0\rangle + |\psi_M\rangle + |\psi_{L-M}\rangle$, provide examples of physical states such that the dynamics of their quadratic fluctuations has a non trivial component on the integrable sector $c_{\pm 1}$. Of course, part of the dynamics of those fluctuations is also in the non integrable sector $c_0$.
# 6 Q-SSEP dynamics on the lattice and in the scaling limit

In this section we consider the dynamics of (quadratic) fluctuations in the Q-SSEP in some detail. We show that the dynamics is essentially diffusive, but at sufficiently late time and large distances the deviations from diffusive behaviour are captured by a continuum description. More precisely, we show that there is a scaling regime in which correlation functions are described by hierarchies of partial differential equations that take the form of diffusion equations with source terms that couple the different levels of the hierarchy. By numerically solving the equations of motion for some examples we show that the scaling limit gives an accurate description of an intermediate time regime of the lattice dynamics.

## 6.1 Two point functions for two replicas

In this section we investigate the equations of motion for averages involving two replicas, i.e. quantities of the form $E[O_1 \otimes O_2]_t$. We first show that in the Q-SSEP the mean expectations evolve diffusively, cf. eqn (84), while the quadratic fluctuations satisfy a set of linear equations encoding diffusion plus interactions, cf. eqn (87).

Averaging the Heisenberg equations of motion for system plus noise over the latter gives

$$E[O_1 \otimes O_2]_{t+\Delta t} = E[O_1 \otimes O_2]_t - \frac{1}{2} E\left[(dH_t^{(2)}, [dH_t^{(2)}, O_1 \otimes O_2])\right]_t + \dots, \quad (81)$$

where the two-replica Hamiltonian increment reads

$$dH_t^{(2)} = \sum_{j=1}^{L} \sum_{a=1}^{2} c_{a,j+1}^{\dagger} c_{a,j} dW_t^j + c_{a,j}^{\dagger} c_{a,j+1} d\overline{W}_t^j. \quad (82)$$

Here $dW_t^j, d\overline{W}_t^j$ are complex Brownian noises, cf. eqn (3), and $c_{a,j}$ and $c_{a,j}^{\dagger}$ are the canonical fermion annihilation and creation operators defined in (25) (we recall that the fermion operators of different species commute). As we will see, for the Q-SSEP the operators of interest are

$$n_{a,j} = c_{a,j}^{\dagger} c_{a,j}, \quad S_j^+ = c_{1,j}^{\dagger} c_{2,j} = S_j^{-\dagger}, \quad S_j^- = c_{1,j} c_{2,j}^{\dagger} = (S_j^+)^{\dagger}. \quad (83)$$

Working out the relevant commutators and noise-averages in (81) we obtain

$$\begin{aligned} \frac{d}{dt} E[n_{a,j}]_t &= J \Delta_{j,k} E[n_{a,k}]_t, \\ \frac{d}{dt} E[S_j^{\pm}]_t &= J \Delta_{j,k} E[S_k^{\pm}]_t, \end{aligned} \quad (84)$$

where $\Delta_{j,k} = \delta_{j,k+1} + \delta_{j,k-1} - 2\delta_{j,k}$ is the lattice Laplacian and repeated indices are summed over. This establishes that one-point functions exhibit purely diffusive dynamics. Assuming that the initial density matrix is invariant under the exchange of the two replicas, the average fermion density must be independent of the replica index $a$. For later convenience we define it as

$$\rho_j(t) = \mathrm{Tr}[\rho^{(2)}(0)E[n_{a,j}]_t]. \quad (85)$$

The fermion density fulfils a simple diffusion equation on the lattice

$$\frac{d}{dt} \rho_j(t) = J \Delta_{j,k} \rho_k(t). \quad (86)$$
The analogous calculation for operators involving two operators $n_{a,j}$, $S_j^\pm$ gives

$$ \frac{d}{dt} \mathbb{E}[n_{1,j} n_{2,k}]_t = J(\Delta_{j,n} + \Delta_{k,m}) \mathbb{E}[n_{1,n} n_{2,m}]_t - J \delta_{j,k} (\mathbb{E}[S_{j+1}^+ S_j^-]_t + \mathbb{E}[S_{j-1}^+ S_j^-]_t + \text{h.c.}) \\ + J(\delta_{j,k-1} + \delta_{j,k+1}) (\mathbb{E}[S_j^- S_k^+]_t + \mathbb{E}[S_j^+ S_k^-]_t), \quad (87) $$

$$ \frac{d}{dt} \mathbb{E}[S_j^- S_k^+]_t = J(\Delta_{j,n} + \Delta_{k,m}) \mathbb{E}[S_n^- S_m^+]_t \\ + J \delta_{j,k-1} \mathbb{E}[\bar{n}_{1,j} n_{2,j+1} + n_{2,j} \bar{n}_{1,j+1}]_t + J \delta_{j,k+1} \mathbb{E}[\bar{n}_{1,j} n_{2,j-1} + n_{2,j} \bar{n}_{1,j-1}]_t \\ - J \delta_{j,k} \mathbb{E}[\bar{n}_{1,j} (n_{2,j-1} + n_{2,j+1}) + n_{2,j} (\bar{n}_{1,j-1} + \bar{n}_{1,j+1})]_t, \quad (88) $$

where we have defined $\bar{n}_{a,j} = n_{a,j} - 1$.

We now take the trace with an initial density matrix that is invariant under exchange of the replicas and define the following averages of the (replicated) system degrees of freedom

$$ \begin{aligned} g_+(j, k; t) &:= \operatorname{Tr}[\rho^{(2)}(0)\mathbb{E}[n_{1,j} n_{2,k}]_t], \\ g_-(j, k; t) &:= \begin{cases} -\operatorname{Tr}[\rho^{(2)}(0)\mathbb{E}[S_j^- S_k^+]_t] & \text{if } j \neq k, \\ \operatorname{Tr}[\rho^{(2)}(0)\mathbb{E}[n_{1,j} n_{2,j}]_t] & \text{if } j = k. \end{cases} \end{aligned} \quad (89) $$

Alternatively, we can express these as

$$ g_+(j, k; t) = \mathbb{E}[G_{jj} G_{kk}], \quad g_-(j, k; t) = \mathbb{E}[G_{jk} G_{kj}], \quad G_{ij} = \operatorname{Tr}(\rho_t c_i^\dagger c_j). \quad (90) $$

This shows that $g_+(j, k; t)$ and $g_-(j, k; t)$ encode respectively density and coherence correlation fluctuations.

Treating the case $j=k$ carefully, we obtain

$$ \frac{d}{dt} g_{\sigma}(j, k; t) = J(\Delta_{j,n} g_{\sigma}(n, k; t) + \Delta_{k,m} g_{\sigma}(j, m; t)). \quad (91) $$

In the following we will also consider the connected correlation functions,

$$ \begin{aligned} g_{+,c}(j, k; t) &\equiv g_+(j, k; t) - \rho_j(t)\rho_k(t), \\ g_{-,c}(j, k; t) &\equiv g_-(j, k; t) - \delta_{j,k}(\rho_j(t))^2, \end{aligned} \quad (92) $$

where $\rho_j(t)$ is the average fermion density in replica $a$ defined above in eqn (85) (it is independent of $a$ as a result of our choice of initial conditions).

It is easy to see that (91) have time-independent sum rules

$$ A_{\sigma} = \sum_{j,k=1}^{L} g_{\sigma}(j,k;t), \quad \frac{dA_{\sigma}}{dt} = 0. \quad (93) $$

These conservation laws are consequences of the fact that the dynamics is unitary for any given realization of the noise, and hence it preserves the spectrum of the system density matrix. For the initial states we consider here the constants $A_\sigma$ have regular expansions in inverse powers of the system size

$$ A_{\sigma} = a_{\sigma}^{(0)} + \frac{a_{\sigma}^{(1)}}{L} + \frac{a_{\sigma}^{(2)}}{L^2} + \dots \quad (94) $$

The calculation of correlation functions can of course also be formulated in the Hilbert space doubling approach. For example, we have

$$ g_+(j, k; t) = \langle 1_2 \| n_{1,j} n_{2,k} \| \rho^{(2)}(t) \rangle, \quad (95) $$
where $\langle 1_2 |$ has been defined in (60). The equation of motion thus reads

$$ \frac{d}{dt} g_+(j, k; t) = \langle 1_2 \| n_{1,j} n_{2,k} \hat{\mathcal{L}}_2 \| \rho^{(2)}(t) \rangle. \qquad (96) $$

This, and the analogous equation for $g_-(j, k; t)$, can be cast in the form of an imaginary time two-particle Schrödinger equation for a non-Hermitian "Hamiltonian" obtained from $\hat{\mathcal{L}}_2$ by exploiting the existence of the two $gl(2)$ sub-algebras. This is sketched in Appendix E, and a discussion of the spectrum of the Lindbladian in the two-particle sector is given in Appendix F.

## 6.2 Numerical solution

The coupled equations (91) can be straightforwardly solved numerically. We take the initial two-replica density matrix to encode no correlations between the replicas

$$ \mathbb{E}[\rho^{(2)}(t=0)] = \rho^{(1)} \otimes \rho^{(1)}. \qquad (97) $$

### 6.2.1 CDW initial state

As our first example we choose an initial product state that is invariant under translations by three lattice sites

$$ \rho^{(1)} = |\text{CDW}\rangle\langle\text{CDW}|, \quad |\text{CDW}\rangle = \prod_{j=1}^{L/3} c_{3j-2}^{\dagger}|0\rangle. \qquad (98) $$

Here $L$ is taken to be a multiple of 3 so that the charge-density order is commensurate with the system size. This gives the initial conditions

$$ \begin{aligned} g_+(j, k, 0) &= D_j D_k, & D_j &= \sum_{i=1}^{L/3} \delta_{j,3i-2}, \\ g_-(j, k, 0) &= D_j^2 \delta_{j,k}. & & \end{aligned} \qquad (99) $$

The fermion density in the steady state is $\rho_j(\infty) = 1/3 \equiv \rho$ and connected density and coherence fluctuations behave as

$$ \begin{aligned} g_{+,c}(j, k; \infty) &= \frac{\rho(1-\rho)}{L} \left[ \delta_{j,k} - \frac{1}{L} \right] + \mathcal{O}(L^{-3}), \\ g_{-,c}(j, k; \infty) &= \frac{\rho(1-\rho)}{L} \left[ 1 - \frac{\delta_{j,k}}{L} \right] + \mathcal{O}(L^{-3}). \end{aligned} \qquad (100) $$

The results of numerical integration of the equations of motion is shown in Figs 3 and 4. The fermion density is seen to relax quite quickly to its steady state value 1/3 as shown in Fig. 3. The qualitative behaviour of the connected two-point functions is as follows:

* The connected correlations initially vanish, then build up over a time scale $Jt \sim 1$, and subsequently decay towards their steady state values (100);

* The decay of $g_+(j, k; t)$ towards its steady-state values is well-characterised by a power law $t^{-1/2}$ on times scales $Jt \ll L$;

* The region in which connected correlations are non-negligible spreads diffusively, i.e. there is a parabolic “envelope”.
<figure><img src="image_8.png" /><figcaption>Figure 3: Density $\rho_n(t)$ for the CDW initial state on the sublattice $n = 3k-2$ (blue) and the other two sublattices $n = 3k-1, n = 3k$ (yellow) as functions of time. We observe that the density relaxes quickly to its steady state value $1/3$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 4: Connected correlators $g_{\sigma,c}(1, n)$ as functions of time for $n = 1$ (solid blue), $n = 2$ (dashed yellow), $n = 3$ (dotted green) on a ring of $L = 36$ sites. The correlations are very short-ranged.</figcaption></figure>

### 6.2.2 Domain wall initial state

The next example we consider is that of an uncorrelated domain wall state at half-filling, where all fermions are initially located on sites $[1, L/2]$, i.e.

$$ \rho^{(1)} = |DW\rangle\langle DW|, \qquad |DW\rangle = \prod_{j=1}^{L/2} c_j^\dagger |0\rangle. \qquad (101) $$

The initial conditions are thus

$$ \begin{aligned} g_+(j, k, 0) &= \begin{cases} 1 & \text{if } 1 \le j, k \le \frac{L}{2}, \\ 0 & \text{else.} \end{cases} \\ g_-(j, k, 0) &= \begin{cases} 1 & \text{if } 1 \le j = k \le \frac{L}{2}, \\ 0 & \text{else.} \end{cases} \end{aligned} \qquad (102) $$

In Fig. 5 we show the ensuing time evolution of the density and connected two-point function $g_{-c}(\frac{L}{2}+1, n, t)$. The initial state is unentangled, but connected correlations are seen to build up in time and are negligible outside a causal region that spreads diffusively in time, i.e. there is a parabolic "envelope".
<figure><img src="image_10.png" /><figcaption>Figure 5: (a) Density $\rho_j(t)$ for the uncorrelated domain-wall initial state (101) as a function of $j$ and time on a ring of $L = 96$ sites. (b) Connected correlator $g_{-,c}(\frac{L}{2} + 1, n, t)$ for the same parameters.</figcaption></figure>

### 6.2.3 Correlated domain wall initial state

In the two previous examples the initial states were unentangled. We now consider the half-filled ground state of a tight-binding model on half the system, i.e. sites $1 \le j \le L/2$ (we assume $L$ to be even). Defining

$$ \tilde{c}(k_n) = \sqrt{\frac{2}{L/2+1}} \sum_{j=1}^{L/2} c_j \sin(k_n j), \quad k_n = \frac{\pi n}{L/2+1}, \quad n=1, \dots, \frac{L}{2}, \quad (103) $$

our initial state is

$$ |\Psi(0)\rangle = \prod_{n=1}^{L/4} \tilde{c}^{\dagger}(k_n) |0\rangle. \quad (104) $$

The initial two-point functions are then easily calculated using

$$ c_j = \sqrt{\frac{2}{L/2+1}} \sum_{k_n} \tilde{c}(k_n) \sin(k_n j). \quad (105) $$

Defining $2k_F = \frac{\pi}{2}$ we have for $j, k \in [1, L/2]$

$$ g_{j,k} = \langle \Psi(0) | c_j^\dagger c_k | \Psi(0) \rangle = f(j-k) - f(j+k), \quad j \neq k, \\ f(n) = \begin{cases} \frac{\sin(2k_F n)}{(L+2)\sin(\frac{\pi n}{L+2})} & \text{if } n \neq L+2, \\ \frac{1}{2} & \text{if } n = L+2. \end{cases} \quad (106) $$

The corresponding initial values for our one and two-point functions are then

$$ \rho_j(0) = \begin{cases} \frac{1}{2} & \text{if } 1 \le j \le \frac{L}{2}, \\ 0 & \text{if } \frac{L}{2} \le j \le L, \end{cases} \quad g_+(j,k,0) = \rho_j(0)\rho_k(0), \\ g_-(j,k,0) = \begin{cases} g_{j,k} g_{k,j} & \text{if } 1 \le j \neq k \le \frac{L}{2}, \\ \frac{1}{4} & \text{if } 1 \le j = k \le \frac{L}{2}. \\ 0 & \text{else.} \end{cases} \quad (107) $$
The stationary behaviour at late times can again be determined explicitly

$$ 
\begin{aligned}
\lim_{t \to \infty} g_+(j, k, t) &= \frac{1}{16} \left[ 1 - \frac{3}{L^2} + \frac{3}{L} \delta_{j,k} \right], \\
\lim_{t \to \infty} g_-(j, k, t) &= \frac{1}{16} \left[ \frac{3}{L} + \delta_{j,k} \left( 1 - \frac{3}{L^2} \right) \right].
\end{aligned} 
\qquad (108) $$

In Fig. 6 we show the dynamics of the average density $\rho_j(t)$ as a function of time.

<figure><img src="image_11.png" /><figcaption>Figure 6: (a) Density $\rho_j(t)$ for the correlated domain-wall initial state (104) as a function of $j$ and time on a ring of $L = 96$ sites. (b) Connected correlator $g_{-,c}(\frac{L}{2} + 1, n, t)$ for the same parameters.</figcaption></figure>

The domain wall is seen to spread diffusively in a way that is qualitatively quite similar to the case of an uncorrelated domain wall initial state considered above. The connected correlations are non-vanishing in the initial state and as a result the region in which they are sizeable spreads in an asymmetric way.

## 6.3 Scaling limit

The aim of this section is to show that there exists a scaling limit in which a continuum description applies and to derive a set of partial differential equations describing the dynamics of quadratic fluctuations in this limit. We then establish the regime of applicability of results obtained in the scaling limit for the domain wall initial conditions considered above by comparing numerical result for the continuum and lattice descriptions. The scaling limit is defined as

$$ a_0 \to 0, \quad J \to \infty, \quad \mathcal{D} = J a_0^2 \text{ fixed,} \qquad (109) $$

where $a_0$ is the lattice spacing. The limit (109) has proved to provide very useful insights into the stationary behaviour [56], and we now show that this extends to the full dynamics. By introducing co-ordinates

$$ x = ja_0, \quad y = ka_0, \qquad (110) $$

we can turn difference equations into differential equations for the functions $G_\sigma(x, y, t) = g_\sigma(j, k; t)$ and $\rho(x, t) = \rho_j(t)$. The scaling limit of (86) is simply

$$ \left[ \frac{d}{dt} - \mathcal{D} \frac{d^2}{dx^2} \right] \rho(x, t) = 0. \qquad (111) $$
To work out the scaling limits of (91) it is useful to rewrite the lattice equations in the form

$$ \frac{d}{dt} g_{\sigma}(j, k; t) = J \sum_{n} \Delta_{j,n} g_{\sigma}(n, k; t) + \Delta_{k,n} g_{\sigma}(j, n; t) \\ + J \delta_{j,k} \sum_{n} \Delta_{j,n} g_{-\sigma}(n, k; t) + \Delta_{k,n} g_{-\sigma}(j, n; t) - 2J D_{j,k} g_{-\sigma}(j, k; t), \quad (112) $$

where $D_{j,k} = \delta_{j+1,k} + \delta_{j-1,k} - 2\delta_{j,k}$. Taking $h(n)$ to be a discretization of a test function we have

$$ \sum_{k} \delta_{j,k} h(k) = h(j) \Rightarrow \delta_{j,k} \rightarrow a_0 \delta(x-y), \\ \sum_{k} D_{j,k} h(k) = \Delta_{j,n} h(n) \Rightarrow D_{j,k} \rightarrow a_0^3 \delta''(x-y). \quad (113) $$

This then gives the following result for the scaling limit

$$ \frac{d}{dt} G_{\sigma}(x, y; t) = \mathcal{D} \nabla^2 G_{\sigma}(x, y; t) + \mathcal{D} a_0 [\delta(x-y)\nabla^2 - 2\delta''(x-y)] G_{-\sigma}(x, y; t) \\ + \mathcal{O}(J a_0^4). \quad (114) $$

We observe that in the limit $a_0 \to 0$ we obtain a simple diffusion equation. This indicates that for large separations we do not have connected correlations, in agreement with the numerical solution of the lattice model. The sub-leading term in (114) should be understood in the sense that we keep the overall physical length $\mathcal{L} = L a_0$ of our system fixed while taking the lattice spacing $a_0$ to zero and the number of sites $L$ to infinity. Hence we have

$$ \mathcal{L} = a_0 L. \quad (115) $$

The corrections to the diffusion equations are indeed proportional to $1/L$ as expected, say from the known properties of the steady state [41, 42]. Eqns (114) are linear and therefore can be solved by expanding

$$ G_{\sigma}(x, y; t) = G_{\sigma}^{(0)}(x, y; t) + \frac{a_0}{\mathcal{L}} G_{\sigma}^{(1)}(x, y; t) + \dots \quad (116) $$

The functions $G_{\sigma}^{(n)}(x, y; t)$ then fulfil a hierarchy of partial differential equations

$$ \begin{aligned} \partial_t G_{\sigma}^{(0)}(x, y; t) &= \mathcal{D} \nabla^2 G_{\sigma}^{(0)}(x, y; t), \\ \partial_t G_{\sigma}^{(1)}(x, y; t) &= \mathcal{D} \nabla^2 G_{\sigma}^{(1)}(x, y; t) + h_{\sigma}(x, y; t), \end{aligned} \quad (117) $$

where the inhomogeneities are given by

$$ h_{\sigma}(x, y; t) = \mathcal{D} \mathcal{L} [\delta(x-y)\nabla^2 - 2\delta''(x-y)] G_{-\sigma}^{(0)}(x, y; t). \quad (118) $$

Importantly, for generic initial conditions (assuming the absence of long-range order) the connected Green's functions on the lattice will decay to zero with distance, which implies that

$$ G_{+}^{(0)}(x, y; 0) = \rho(x, 0)\rho(y, 0), \quad G_{-}^{(0)}(x, y; 0) = 0. \quad (119) $$

This results in

$$ G_{+}^{(0)}(x, y; t) = \rho(x, t) \rho(y, t), \quad G_{-}^{(0)}(x, y; t) = 0. \quad (120) $$

The sub-leading contributions fulfil (driven) diffusion equations

$$ \left[ \frac{d}{dt} - \mathcal{D} \nabla^2 \right] G_{\sigma}^{(1)}(x, y; t) = h_{\sigma}(x, y; t). \quad (121) $$

These are heat equations with a time-dependent external source acting along the line $x = y$. Let us recapitulate the main results obtained in this section.
• The average fermion density $\rho(x, t)$ evolves diffusively

$$ \partial_t \rho(x, t) = \mathcal{D} \nabla_x^2 \rho(x, t). \qquad (122) $$

• The density correlations $G_+(x, y, t)$ and the coherence fluctuations $G_-(x, y; t)$ admit expansions in powers of $a_0/\mathcal{L}$, cf. eqn (116).

• As a consequence of the absence of long range order in the initial state, the off-diagonal coherences vanish at leading order, i.e.

$$ G_{-}^{(0)}(x, y; t) = 0. \qquad (123) $$

• Since the initial density correlations are factorized $G_+^{(0)}(x, y, 0) = \rho(x, 0)\rho(y, 0)$ and then evolve diffusively, the density correlations remain uncorrelated at leading order in the expansion for all times $t$, i.e.

$$ G_{+}^{(0)}(x, y; t) = \rho(x, t)\rho(y, t). \qquad (124) $$

• As a consequence of the triangular structure of the equations of motion, the coherence fluctuations, which are of order $\mathcal{O}(L^{-1})$, evolve diffusively but are subject to a source term depending on the spatial variation of the fermion density

$$ \partial_t G_{-}^{(1)}(x, y; t) = \mathcal{D}(\nabla_x^2 + \nabla_y^2)G_{-}^{(1)}(x, y; t) + h_{-}(x, y; t), \qquad (125) $$

where $h_{-}(x, y, t)$ is given by (118) or equivalency by

$$ h_{-}(x, y, t) = 2 \mathcal{D} \mathcal{L} \nabla_x \nabla_y (\delta(x - y) \rho(x, t) \rho(y, t)). $$

• As the densities in the two replicas are by construction initially uncorrelated they remain uncorrelated at order $a_0/\mathcal{L}$ at any time

$$ G_{+}^{(1)}(x, y; t) = 0, \qquad (126) $$

where $h_{+}(x, y, t)$ is given by (118) or equivalency by

$$ h_{+}(x, y, t) = 2 \mathcal{D} \mathcal{L} \nabla_x \nabla_y (\delta(x - y) \rho(x, t) \rho(y, t)). $$

• As the densities in the two replicas are by construction initially uncorrelated they remain uncorrelated at order $a_0/\mathcal{L}$ at any time

$$ G_{+}^{(1)}(x, y, t) = 0. \qquad (127) $$

This means that the connected density correlations are of order $\mathcal{O}((a_0/\mathcal{L})^2)$, in agreement with known properties of the steady state [56].

### 6.3.1 Late time regime $\mathcal{D}t/\mathcal{L}^2 \gg 1$

In this regime we have

$$ \rho(x, t) = \rho(x, \infty) + \mathcal{O}(e^{-\gamma t}), \quad \gamma \ge \frac{4\pi^2 \mathcal{D}}{\mathcal{L}^2}. \qquad (128) $$

Here we have used the fact that there are low-lying magnon excitations with dispersion (79) and the smallest non-zero momentum is $2\pi/L$. At asymptotically late times the relaxation to the steady state is exponential in time, and only the lowest “excited” mode contributes to the dynamics.
### 6.3.2 Intermediate time regime $\mathcal{D}t/\mathcal{L}^2 \ll 1 \ll \mathcal{D}t/a_0^2$

In this regime the continuum approximation applies as $\mathcal{D}t \gg a_0^2$, but correlations have not yet spread throughout the entire volume and all low-lying excited modes (with momenta close to zero) contribute to the dynamics. In order to work out the solution of our system of equations in this regime we require the Green's function

$$ \left[ \frac{d}{dt} - \mathcal{D}\nabla^2 \right] G_2(\mathbf{x}, t; \mathbf{y}, t') = \delta^{(2)}\left(\frac{\mathbf{x} - \mathbf{y}}{\mathcal{L}}\right) \delta(t - t'). \quad (129) $$

With $k_{n,m} = \frac{2\pi}{\mathcal{L}}(n, m)^T$, we have

$$ \begin{aligned} G_2(\mathbf{x}, t; \mathbf{y}, t') &= \theta(t-t') \sum_{n,m \in \mathbb{Z}} e^{-ik_{n,m} \cdot (\mathbf{x}-\mathbf{y}) - \mathcal{D}(t-t')k_{n,m}^2} \\ &= \frac{\theta(\tau)e^{-\frac{\xi^2}{4\tau}}}{4\pi\tau} \vartheta_3(i\frac{\xi_1}{4\tau}, e^{-\frac{1}{4\tau}}) \vartheta_3(i\frac{\xi_2}{4\tau}, e^{-\frac{1}{4\tau}}), \end{aligned} \quad (130) $$

where $\vartheta_3$ is an elliptic theta function whose arguments were transformed by a modular transformation and we have defined

$$ \xi = \frac{\mathbf{x} - \mathbf{y}}{\mathcal{L}}, \quad \tau = \frac{\mathcal{D}(t - t')}{\mathcal{L}^2}. \quad (131) $$

The correlation functions of interest can be decomposed as e.g.

$$ G^{(1)}(x, y; t) = \mathcal{G}(x, y; t) + \mathcal{H}(x, y; t), \quad (132) $$

where $\mathcal{G}(x, y, t)$ is the solution of the homogeneous equation with the appropriate initial conditions and $\mathcal{H}(x, y, t)$ is the solution of the inhomogeneous equation with vanishing initial conditions. With $\mathbf{x} = (x, y)$, the latter can be written as

$$ \begin{aligned} \mathcal{H}(x, y, t) &= \int_0^\infty dt' \int_0^\mathcal{L} \int_0^\mathcal{L} \frac{d^2\mathbf{z}}{\mathcal{L}^2} G_2(\mathbf{x}, t; \mathbf{z}, t') h_-(z_1, z_2, t') \\ &= \int_0^t dt' \int_0^1 dz \, K(\zeta_1 - \mathbf{z}, \zeta_2 - \mathbf{z}, \tau) \rho^2(z\mathcal{L}, t'), \end{aligned} \quad (133) $$

where

$$ \begin{aligned} K(\zeta_1, \zeta_2, \tau) &= \frac{2\mathcal{D}}{\pi\mathcal{L}^2(4\tau)^3} e^{-\frac{\zeta_1^2+\zeta_2^2}{4\tau}} g_1(\zeta_1, \tau)g_1(\zeta_2, \tau), \\ g_1(\zeta, \tau) &= 2\zeta\vartheta_3\left(i\frac{\zeta}{4\tau}, e^{-\frac{1}{4\tau}}\right) - i\vartheta'_3\left(i\frac{\zeta}{4\tau}, e^{-\frac{1}{4\tau}}\right). \end{aligned} \quad (134) $$

Here we have defined

$$ \tau(t') = \frac{\mathcal{D}(t-t')}{\mathcal{L}^2}, \quad \zeta_1 = \frac{x}{\mathcal{L}}, \quad \zeta_2 = \frac{y}{\mathcal{L}}. \quad (135) $$

The two equivalent expressions in eqn (133) are related by integration by parts.

### 6.3.3 Numerical test of the scaling limit

We now test the accuracy of the continuum approximation for the domain-wall initial state (6.2.2). In the continuum we replace the initial conditions (102) by

$$ \begin{aligned} G_+^{(0)}(x, y; 0) &= \theta(\frac{\mathcal{L}}{2} - x)\theta(\frac{\mathcal{L}}{2} - y), \\ G_-^{(1)}(x, y; 0) &= \delta(\frac{x-y}{\mathcal{L}})\theta(\frac{\mathcal{L}}{2} - x). \end{aligned} \quad (136) $$
Importantly these correctly reproduce the sum rules (93). Strictly speaking we should smear out the delta-function over some microscopic length scale associated with the high-energy cutoff of the continuum theory. We have verified that this essentially only modifies the short time behaviour. Using the initial conditions (136) we find

$$ \rho(x, t) = \frac{1}{2} + \frac{2}{\pi} \sum_{n=0}^{\infty} \frac{\sin(k_{2n+1}x)}{2n+1} e^{-\mathcal{D}t k_{2n+1}^2}, \quad (137) $$

$$ \begin{aligned} \mathcal{G}(x, y; t) = & \frac{\rho(\frac{x+y}{2}, t/2)e^{-\frac{(\xi_1 - \xi_2)^2}{8\tau}}}{\sqrt{8\pi\tau}} \vartheta_3\left(-\frac{\pi}{2} - i\frac{\xi_1 - \xi_2}{8\tau}e^{-\frac{1}{8\tau}}\right) \\ & + \frac{e^{-\frac{(\xi_1 - \xi_2)^2}{8\tau}}}{4\sqrt{2\pi\tau}} \vartheta_3\left(i\frac{\xi_1 - \xi_2}{8\tau}e^{-\frac{1}{8\tau}}\right) - \frac{e^{\frac{(\xi_1 - \xi_2)^2}{8\tau}}}{4\sqrt{2\pi\tau}} \vartheta_3\left(-\frac{\pi}{2} - i\frac{\xi_1 - \xi_2}{8\tau}e^{-\frac{1}{8\tau}}\right). \end{aligned} \quad (138) $$

The contribution $\mathcal{H}(x, y; t)$ is then obtained by numerically integrating (133) with the density (137). The result can then be directly compared to the full lattice model. Recalling that $x = ja_0$, $y = ka_0$ and the lattice spacing is $a_0 = \mathcal{L}/L$ we expect that for sufficiently large $L$ we have

$$ g_{-}(L\frac{x}{\mathcal{L}}, L\frac{y}{\mathcal{L}}; t) \approx \frac{1}{L} G_{-}^{(1)}(x, y; t) + \mathcal{O}(L^{-2}). \quad (139) $$

In Fig. 7 we show a comparison of the lattice and continuum results for $(x/\mathcal{L}, y/\mathcal{L}) = (\frac{1}{2}, \frac{3}{8})$ and $L = 48, 96, 192$.

<figure><img src="image_12.png" /><figcaption>Figure 7: Comparison between $g_{-}(L/2, 3L/8; t)$ as a function of $Jt/L^2$ on a ring of $L = 48$ (blue dashed), $L = 96$ (yellow dotted), $L = 192$ (gree dot-dashed) sites with domain wall initial conditions and $G_{-}^{(1)}(x, y; t)$ (red solid). The continuum limit result is approached throughout the depicted time range as $L$ increases.</figcaption></figure>

We observe that as $L$ is increased, the lattice result approaches the continuum one, as expected. The agreement between lattice and continuum is excellent throughout the intermediate time range considered. In Fig. 8 we show the analogous comparison for $(x/\mathcal{L}, y/\mathcal{L}) = (\frac{1}{2}, \frac{23}{48})$ and $L = 96$. The lattice and continuum results match perfectly, and the $L$-dependence is negligible in the sense that the results for $L = 48, 192$ are indistinguishable on the scale of the plot. We also show the result for $\mathcal{G}(x, y; t)$ only, which is clearly different. This shows that the “interaction term” $h_{-}(x, y; t)$ in the continuum equations (117) is important.
<figure><img src="image_13.png" /><figcaption>Figure 8: Comparison between $g_{-}(L/2, 23L/48; t)$ as a function of $Jt/L^2$ on a ring of $L = 96$ (blue dashed) sites with domain wall initial conditions and $G_{-}^{(1)}(x, y; t)$ (red solid). The dotted brown line is the result for $G(x, y; t)$ only.</figcaption></figure>

### 6.3.4 "Diagonal" correlations in the scaling limit

So far we have considered the two-point functions $G_{\sigma}(x, y; t)$ for non-vanishing separations $x \neq y$. As we will see later on, the "diagonal" correlations $G_{\sigma}(\xi\sqrt{t}, \xi\sqrt{t}; t)$ are of interest in the study of operator spreading. We therefore turn to their analysis next.

In the following we present a numerical check that as a function $\xi = \frac{x}{\sqrt{t}}$ the diagonal correlations $G_{\sigma}(\xi\sqrt{t}, \xi\sqrt{t}; t)$ factorize in the leading order in $\frac{1}{L}$

$$ \mathbb{E}[G_{\sigma}(\xi\sqrt{t}, \xi\sqrt{t}; t)^{\alpha}] = \mathbb{E}[G_{\sigma}(\xi\sqrt{t}, \xi\sqrt{t}; t)]^{\alpha}, \quad (140) $$

in the limit $L \to \infty$, $t \to \infty$. In order to avoid finite size effects, in the study of operator spreading, we also take the limit $\mathcal{L} \to \infty$. In Fig. 9 we plot $\mathbb{E}[G(\xi)^2] = G_{+}^{(0)}(\xi \times \sqrt{t}, \xi \times \sqrt{t}, t)$ at different times and compare it to the asymptotic distribution, which is approached as $\frac{1}{\sqrt{t}}$, as presented in Fig. 10.

<figure><img src="image_14.png" /><figcaption>Figure 9: On the left we plot the comparison of $\mathbb{E}[G(\xi)^2]$ at time $t = 150$ (the blue line), and the square of expression (173) (the black dashed line). On the right we plot the difference between the asymptotic profile, and finite time results. The system size is fixed to $L = 120$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 10: In this plot we show the scaling of the maximal difference, with respect to the coordinate $\xi$, between the average of square of $\mathbb{E}[G(\xi)^2]$ for times between 50 and 150 and the square of expression (173) (blue dots). We can see that the difference decays as $1/\sqrt{t}$. Indeed, the linear extrapolation (black line) yields the value at the origin $8 \times 10^{-4}$.</figcaption></figure>

## 6.4 Evolution of higher point correlation functions

We now turn to multi-point correlation functions in the Q-SSEP. To ease notations we set $\mathcal{L} = 1 = \mathcal{D}$ (which implies $J = L^2$) in the remainder of this section. It is straightforward to reintroduce $\mathcal{L}, \mathcal{D}$ in the various equations by dimensional analysis. The main point we make here is to show that higher point correlation functions satisfy equations of motion with a triangular structure similar to eqn (114) in the scaling limit.

We consider n-point correlation functions of the form

$$g_n(l_1, \dots, l_n) = \mathbb{E}[G_{l_1 l_2} G_{l_2 l_3} \dots G_{l_n l_1}] =: \left[ \begin{smallmatrix} n \\ l_1 \\ \vdots \\ l_n \\ n \end{smallmatrix} \right] =: \left[ \begin{smallmatrix} n \\ n \end{smallmatrix} \right], \quad (141)$$

which we call n-point bubbles due to its the diagrammatic representation. The discrete indices $l_1, \dots, l_n$ take values in $\{1, \dots, L\}$, $L$ is the number of sites of the system, $G_{ij} = \text{Tr}(\rho c_j^\dagger c_i)$ and $\mathbb{E}[\cdot]$ is the expectation value w.r.t. the noise. In the diagrammatic representation each index $l_a$ is a node and $G_{l_a l_b}$ is a directed edge from node $l_a$ to $l_b$.

In the scaling limit ($J = L^2 \to \infty$) the n-bubbles become functions of the continuous variables $x_a := l_a/L \in [0, 1]$ and they have an expansion in inverse powers of $L$,

$$g_n = g_n^{(0)} + L^{-1} g_n^{(1)} + L^{-2} g_n^{(2)} + \dots \quad (142)$$

As shown in Appendix G, the first $n-1$ terms in this expansion are actually zero and the leading order of an n-bubble scales with $L^{-(n-1)}$. We furthermore show that any expectation value of products of bubbles will factorize at leading order into products of expectation value. Therefore, if interested in the leading order of any correlation function it is enough to study bubbles$^5$.

<footnote>$^5$Open diagrams such as $\mathbb{E}[G_{xy}G_{yz}]$ decay exponentially in time and are of no interest. Pinched diagrams such as $\mathbb{E}[G_{xy}G_{yx}G_{xx}]$ can be obtained by continuity from the three-bubble $\mathbb{E}[G_{xy}G_{yz}G_{zx}]$ as the limit where $z \to x$ (to be precise, this is true for the connected part only). Therefore, the leading order of any correlation function is encoded into bubbles.</footnote>
Denoting the leading contribution of any term by the superscript #, we find the time evolution of n-point bubbles to be

$$ \frac{d}{dt} \left[ \begin{smallmatrix} \left( \begin{smallmatrix} n \\ \bullet \end{smallmatrix} \right) \end{smallmatrix} \right]^\# = \sum_{a=1}^{L} \Delta_{x_a} \left[ \begin{smallmatrix} \left( \begin{smallmatrix} n \\ \bullet \end{smallmatrix} \right) \end{smallmatrix} \right]^\# \\ + \sum_{\substack{a,b=1 \\ a<b}}^{L} \partial_{x_a} \partial_{x_b} \left( \left[ \begin{smallmatrix} \left( \begin{smallmatrix} n-b+a \\ \bullet \end{smallmatrix} \right) & \bullet \\ \bullet & \bullet \end{smallmatrix} \right]^\# \left[ \begin{smallmatrix} \left( \begin{smallmatrix} n-b+a \\ \bullet \end{smallmatrix} \right) & \bullet \\ \bullet & \bullet \end{smallmatrix} \right]^\# \right] + (l_a \leftrightarrow l_b). \quad (143) $$

The equations suggests that the leading contribution of a n-point bubble diffuses and is sourced by the product of the leading contribution of smaller bubbles - hence the "triangular" structure. The source term arises if to indices $x_a$ and $x_b$ are in contact and its construction can be visualized as follows: The original n-bubbles is squeezed together such that the nodes $x_a$ and $x_b$ touch each other and the diagram forms an eight. Then the eight is split apart into two disconnected bubbles and we sum over the two possible ways how one can attribute $x_a$ and $x_b$ to the two nodes at the splitting junction. To illustrate this equation, we give its explicit form in the case of 3-point bubbles.

$$ \frac{d}{dt} \left[ \begin{smallmatrix} \bullet & \bullet \\ \bullet & \bullet \end{smallmatrix} \right]^\# = (\Delta_x + \Delta_y + \Delta_z) \left[ \begin{smallmatrix} \bullet & \bullet \\ \bullet & \bullet \end{smallmatrix} \right]^\# \\ + \partial_x \partial_y \left( \delta(x-y) \left[ \begin{smallmatrix} \bullet & \bullet \\ \bullet & \bullet \end{smallmatrix} \right]^\# \left[ \begin{smallmatrix} \bullet & \bullet \\ \bullet & \bullet \end{smallmatrix} \right]^\# + (x \leftrightarrow y) \right) \\ + \partial_y \partial_z \left( \delta(y-z) \left[ \begin{smallmatrix} \bullet & \bullet \\ \bullet & \bullet \end{smallmatrix} \right]^\# \left[ \begin{smallmatrix} \bullet & \bullet \\ \bullet & \bullet \end{smallmatrix} \right]^\# + (y \leftrightarrow z) \right) \\ + \partial_z \partial_x \left( \delta(z-x) \left[ \begin{smallmatrix} \bullet & \bullet \\ \bullet & \bullet \end{smallmatrix} \right]^\# \left[ \begin{smallmatrix} \bullet & \bullet \\ \bullet & \bullet \end{smallmatrix} \right]^\# + (z \leftrightarrow x) \right). \quad (144) $$

## 6.5 Dictionary for 2-replica operators

Since initially the aim was to find equations for the evolution of any correlation function on two replicas, we will show here how Wick's theorem can be used to find a correspondence between any connected 2-replica correlation function and the leading diagrams contributing to these correlation functions - which are bubbles by the argument in the last section.

Any 2-replica operator can be written in terms of $n_{a,i} = c_{a,i}^\dagger c_{a,i}$ (where $a = 1,2$) and pairs of $S_j^+ = c_{1,j}^\dagger c_{2,j}$ and $S_k^- = c_{1,k} c_{2,k}^\dagger$ (fermions on different replica commute in our convention). Note that we need the same number of $S^+$ and $S^-$, otherwise the quantum expectation number is zero for a Gaussian state. Wick's theorem applies if the state of the system is described by a Gaussian density matrix $\rho = e^{c^\dagger M_c} / \text{Tr}(e^{c^\dagger M_c})$, which contains all possible quadratic terms in fermionic creation and annihilation operators.

As an example, the four-point correlation function $\mathbb{E}[\langle S_x^+ S_y^- n_{1,a} n_{2,b} \rangle]$, where $\langle \cdot \rangle$ is the quantum expectation value, decomposes into

$$ - \left[ \begin{smallmatrix} \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \end{smallmatrix} \right] + \left[ \begin{smallmatrix} \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \end{smallmatrix} \right] + \left[ \begin{smallmatrix} \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \end{smallmatrix} \right] + \left[ \begin{smallmatrix} \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \end{smallmatrix} \right]. \quad (145) $$

We learned in the last section that products of bubbles factorize at leading order. Therefore, the correspondence two-replica correlation functions and bubbles (which we will in the following denote by ~) is on the level of the connected part. To leading order in $L$ we have

$$ \mathbb{E}[\langle S_x^+ S_y^- n_{1,a} n_{2,b} \rangle]^c \sim \left[ \begin{smallmatrix} \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \\ \bullet & \bullet & \bullet & \bullet \end{smallmatrix} \right]. \quad (146) $$
This correspondence can be carried further to any combination of operators. First consider,

$$ \mathbb{E} [ \langle S_{x_1}^+ S_{y_1}^- \dots S_{x_n}^+ S_{y_n}^- \rangle ]^c \sim - \sum_{\sigma \in S_n, \tau \in S_{n-1}} \left[ \begin{smallmatrix} 2n & & & \\ & \ddots & x_{\tau(1)} & \\ & & y_{\sigma(1)} & \\ & & & \ddots & \end{smallmatrix} \right], \quad (147) $$

where $S_n$ is the permutation group of $n$ elements.

Next, any correlation function of the form $\mathbb{E}[\langle n_{1,a_1} \cdots n_{1,a_l} n_{2,b_1} \cdots n_{2,b_k} \rangle]^c$ is in fact a correlation function of the product of 1-replica operators, $\mathbb{E}[\langle n_{a_1} \cdots n_{a_l} \rangle \langle n_{b_1} \cdots n_{b_k} \rangle]^c$. Therefore, the leading contribution of its connected part is not given by a single $l+k$ bubble, but rather by the connected part of the product of two bubbles with $l$ and $k$ nodes. These diagrams scale with even higher (negative) power of $L$ than a $k+l$-bubble, and their time evolution has not been considered in the last section. Still we can give the correspondence for 1-replica correlation functions of the form

$$ \mathbb{E} [ \langle n_{a_1} \dots n_{a_l} \rangle ]^c \sim (-1)^{l+1} \sum_{\sigma \in S_{l-1}} \left[ \begin{smallmatrix} l \\ a_{\sigma(l-1)} \\ \vdots \\ a_l \\ a_{\sigma(1)} \end{smallmatrix} \right] . \quad (148) $$

Finally, the most general correlation function $\mathbb{E}[\langle S_{x_1}^+ S_{y_1}^- \dots S_{x_n}^+ S_{y_n}^- n_{1,a_1} \dots n_{1,a_l} n_{2,b_1} \dots n_{2,b_k} \rangle]^c$ has as its leading contribution bubbles with $2n+l+k$ nodes that must be arranged on a circle according to the following rules:

* every $x$ is followed by either $b$ or $y$;

* every $y$ is followed by either $a$ or $x$;

* every $a$ is followed by either $a$ or $x$;

* every $b$ is followed by either $b$ or $y$.

One can convince oneself of these rules by looking at the example given above.

# 7 Application: Operator spreading

In this section, we will consider the large scale dynamics of operator spreading, which attracted a significant amount of attention in recent years as an indicator of quantum chaos [72, 73]. Research, thus far, has mostly focused on operator spreading in isolated quantum systems, or in random unitary circuits [28, 29, 74, 75, 76, 76–78]. It is, however, important to understand operator spreading in continuous (in time) systems, which are coupled to an environment. The Q-SSEP represents a perfect test-bed to address these questions.

We will focus on two aspects of operator spreading, the first one being the out-of-time ordered correlators (OTOC), which can in some cases be related to Lyapunov exponents [72, 73], and secondly on the hydrodynamics of the operator entanglement spreading [79–94], which has been conjectured to be able to distinguish between integrable and non-integrable isolated systems [95, 96], even in cases when OTOCs fail to distinguish between different types of dynamics [96]. Operator entanglement is also interesting due to its relation to the computational complexity in the matrix product ansatz based algorithms.

In what follows we will consider the spreading of the single site fermion creation and annihilation operators $c_i, c_i^\dagger$. Without loss of generality, we will focus on the time evolution of operator $c_0$. For every realization of the noise, the time evolution of these operators is free, meaning that

$$ c_0(t) = \sum_i A_i(t) c_i . \quad (149) $$
Nevertheless, we have to keep in mind that coefficients

$$A_i(t) = \mathrm{Tr}(c_i^\dagger c_0(t)) \qquad (150)$$

are fluctuating objects.

## 7.1 Out-of-time ordered correlation functions

In general, the weight of an operator $O(t)$ that has spread to site $x$ at time $t$ can be quantified by summing up the contributions of the commutators between the operator $O(t)$ and the basis elements of the local operator algebra $e_x^{(r)} \in \{c_x^\dagger, c_x, n_x, \mathbb{I} - n_x\}$

$$C(x, t) = \sum_r \frac{1}{\mathrm{Tr}(\mathbb{1})} \mathrm{Tr}([O(t), e_x^{(r)}]^\dagger [O(t), e_x^{(r)}]). \qquad (151)$$

For the operator $O(0) = c_0$, the OTOC $C(x, t)$ can be related to the amplitudes $A_x(t)$ defined in (150)

$$C(x, t) = 4 - |A_x(t)|^2. \qquad (152)$$

Note the constant 4, which arises due to the non-commutativity of the basis elements associated with different sites of the chain.

The OTOC amplitude (152) can also be accessed by considering the dynamics of the two point function, if the system is initially prepared in the state with a single fermion at position $x$. Namely, the two point function

$$G_{0,0}^{[x]}(t) = \mathrm{Tr}(c_0^\dagger c_0 \rho^{[x]}(t)) = \mathrm{Tr}(c_0^\dagger(t)c_0(t)\rho^{[x]}) = A_i^\dagger(t)A_j(t)\mathrm{Tr}(c_i^\dagger c_j \rho^{[x]}) \qquad (153)$$

gives us

$$C(x, t) = 4 - G_{0,0}^{[x]}(t), \qquad (154)$$

if we choose $\rho^{[x]}$ such that $\mathrm{Tr}(c_{i'}^\dagger c_i \rho^{[x]}) = \delta_{i,x} \delta_{i',x}$.

It is clear that at the level of averages the spreading of OTOCs is purely diffusive due to translation invariance

$$\mathbb{E}[C(x, t)] = 4 - \mathbb{E}[G_{-x,-x}^{[0]}(t)]. \qquad (155)$$

In the scaling limit $L \to \infty$ and for large times $t$, we have

$$\mathbb{E}[C(x, t)] = 4 - \frac{1}{\sqrt{4\pi t}} \exp(-x^2/4t) + o(t^{-1/2}). \qquad (156)$$

## 7.2 Operator entanglement

A complete operator basis of the system can be obtained by multiplying the basis operators associated with the local algebra $e_x^{(r)}$. Considering a bi-partition into the subsystem A and the subsystem B any operator can be represented as

$$O/\sqrt{\mathrm{Tr}(O^\dagger O)} = \sum_{i,j} M_{i,j} O_{A,i} O_{B,j}, \qquad (157)$$

where $O_{A,i}$ and $O_{B,j}$ are orthonormalized basis elements of the operator spaces associated with the subspaces A and B respectively

$$\mathrm{Tr}_{A/B}(O_{A/B,i}^\dagger O_{A/B,j}) = \delta_{ij}. \qquad (158)$$
Operator entanglement of the local operator O can then be obtained by performing the singular value decomposition of the matrix $[M]_{i,j} = M_{i,j}$, which provides the Schmidt decomposition of the operator O

$$O/\sqrt{\mathrm{Tr}(O^{\dagger}O)} = \sum_{i} \sqrt{\lambda_{i}} O_{A,i} O_{B,i}. \qquad (159)$$

Schmidt coefficients $\lambda_i$ are normalized, $\sum_i \lambda_i = 1$, and $\alpha$-th Rényi operator entanglement entropy reads

$$S_{\alpha} = \frac{1}{1-\alpha} \log \left( \sum_{i} \lambda_{i}^{\alpha} \right). \qquad (160)$$

We will focus on the operator entanglement $S_\alpha(x, t)$ of the operator (149), for two connected bipartitions $A = [1, x]$ and $B = (x, L]$, which can be deduced from the decomposition

$$O/\sqrt{\mathrm{Tr}(O^{\dagger}O)} = \sum_{i,j} M_{(r_1, r_2, \dots, r_x), (r_{x+1}, \dots, r_L)} \left( \prod_{k=1}^{x} e_k^{(r_k)} \right) \left( \prod_{h=x+1}^{L} e_h^{(r_h)} \right). \qquad (161)$$

Note that $S_\alpha(x, t)$ is defined for a fixed realization of the noise, and in the following we will be concerned with the averaged quantity. The terms comprising the operator can be grouped in two parts

$$\begin{aligned} c_0(t) &= \sum_{i \in A} A_i(t) c_i + \sum_{i \in B} A_i(t) c_i \qquad (162) \\ &= \left( \sum_{i \in A} A_i(t) c_i \right) \left( \prod_{i \in B} 1_i \right) + \left( \prod_{i \in A} 1_i \right) \left( \sum_{i \in B} A_i(t) c_i \right). \end{aligned}$$

Clearly the two parts satisfy the orthogonality condition, and the Schmidt coefficients can be obtained simply by computing the norms of associated operators:

$$\lambda_1(x, t) = \sum_{i=-L/2+1}^{x} |A_i(t)|^2; \quad \lambda_2(x, t) = \sum_{i=x+1}^{L/2} |A_i(t)|^2. \qquad (163)$$

Operator entanglement is upper bounded by the logarithm of the number of non-zero Schmidt coefficients $S_\alpha(x, t) \le \log(2)$, which can be understood as a direct consequence of the free dynamics of fermions for any realization of the noise. Similar upper bounds can be obtained for any composite operator [85].

Also in this case it proves useful to relate the dynamics of operator entanglement to the two point function $G_{i,j}^{[x^{-}]}(t) = \mathrm{Tr}(c_{i}^{\dagger} c_{j} \rho_{x}^{-}(t))$ for the domain wall initial conditions

$$\mathrm{Tr}(c_{l'}^{\dagger} c_{k'} \rho_{x}^{-}(0)) = \delta_{l'l'} \theta(x + 1/2 - l'), \quad x, l', k' \in [1, L], \qquad (164)$$

as

$$\lambda_1(x, t) = G_{L/2, L/2}^{[x^{-}]}(t). \qquad (165)$$

The second Schmidt coefficients can be related to the quench with a complementary domain wall initial conditions

$$\mathrm{Tr}(c_{l'}^{\dagger} c_{k'} \rho_{x}^{+}(0)) = \delta_{l'l'} \vartheta(l' - x - 1/2) \quad x, l', k' \in [1, L]. \qquad (166)$$

as

$$\lambda_2(x, t) = G_{L/2, L/2}^{[x^{+}]}(t). \qquad (167)$$
The operator entanglement therefore takes the following suggestive form

$$S_{\alpha}(x, t) = \frac{1}{1-\alpha} \log(G_{L/2,L/2}^{[x^{-}]}(t)^{\alpha} + G_{L/2,L/2}^{[x^{+}]}(t)^{\alpha}). \quad (168)$$

In what follows, we will be concerned with the large system size limit $L \to \infty$, while, for the moment, keeping $x$ and $t$ finite. Firstly, this allows us to relate operator entanglement to the spatio-temporal profile of two point functions, starting from the single initial domain wall condition, with the cusp at the origin

$$S_{\alpha}(x, t) = \frac{1}{1-\alpha} \log(G_{-x,-x}(t)^{\alpha} + G_{x+1,x+1}(t)^{\alpha}). \quad (169)$$

Here we used a simplified notation $G_{x,x}(t) \equiv G_{x,x}^{[\rho]}(t)$, for the domain wall at the origin

$$\mathrm{Tr}(c_{l'}^{\dagger} c_{k'} \rho_0(0)) = \delta_{l'k'} \theta(1/2 + l'). \quad (170)$$

The intuition of the large $L$ limit, which we built in the preceding section, lies at the heart of our conjecture that in the scaling limit $t, x \to \infty$, while $\xi = \frac{x}{\sqrt{t}}$ is kept fixed, the identification (140) holds.

Equation (140), which asserts that products of Green functions are self-averaging at leading order in $1/L$, implies that the averaged operator entanglement

$$S_{\alpha}(\xi) := \lim_{t \to \infty} S_{\alpha}(\xi \sqrt{t}, t) \quad (171)$$

can be, in the infinite size limit limit$^6$, reduced to the entanglement of the averaged behaviour

$$\mathbb{E}[S_{\alpha}(\xi)] = \frac{1}{1-\alpha} \log(\mathbb{E}[G(\xi)]^{\alpha} + \mathbb{E}[G(-\xi)]^{\alpha}). \quad (172)$$

In the continuum limit $\mathbb{E}[G(\xi)]$ satisfies the diffusion equation, which implies that

$$\mathbb{E}[G(\xi)] = \frac{1}{2}(\mathrm{Erf}(\xi/2) + 1). \quad (173)$$

The hydrodynamics of the operator entanglement is therefore described by

$$\mathbb{E}[S_{\alpha}(\xi)] = \frac{1}{1-\alpha} \log \left[ \left( \frac{\mathrm{Erf}(\xi/2) + 1}{2} \right)^{\alpha} + \left( \frac{\mathrm{Erf}(-\xi/2) + 1}{2} \right)^{\alpha} \right]. \quad (174)$$

Clearly the operator entanglement at the origin is maximal, $\mathbb{E}[S_{\alpha}(0)] = \log(2)$, for any value of $\alpha$. There are two particularly interesting limits, which we might wish to consider. The first one is the Von Neumann entanglement entropy, which is obtained in the limit $\alpha \to 1$, and the $\infty$ entanglement entropy $\alpha \to \infty$, which corresponds to the logarithm of the largest Schmidt coefficient

$$\mathbb{E}[S_1(\xi)] = \mathrm{erfc}\left(\frac{\xi}{\sqrt{2}}\right) \tanh^{-1}\left[\mathrm{erf}\left(\frac{\xi}{\sqrt{2}}\right)\right] - \log\left[\mathrm{erf}\left(\frac{\xi}{\sqrt{2}}\right) + 1\right] + \log 2, \quad (175)$$

$$\mathbb{E}[S_{\infty}(\xi)] = \log(2) - \log(1 + \mathrm{erf}(|\xi|/\sqrt{2})). \quad (176)$$

In Fig. 11 we plot the entanglement entropy for different Rényi indexes $\alpha$.

In this section we have demonstrated that spreading of operators in Q-SSEP model falls in a different "dynamical universality class" from isolated quantum systems or random unitary circuits. While in the latter two cases operators spread ballistically with diffusive corrections,

<footnote>$^6$To leading order in the system size.</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 11: Asymptotic profile of operator entanglement for different values of Rényi indexes α. The dashed curves correspond to the extremal values of parameters.</figcaption></figure>

we observed that in Q-SSEP operator spreading is purely diffusive. This observation is im-
portant as it begs the question whether in open quantum systems, where the bath degrees
of freedom are much faster than the dynamics in the quantum system under consideration,
interactions with bath can in general reduce operator spreading. While this is a fundamental
question in its own right, it is also related to the possibility of using weak dissipation in order to
facilitate classical simulations of quantum dynamics. Similar behaviour was already exploited
in a number of numerical studies of quantum transport employing the dissipative boundary
driving, e.g. [97–99].

**8 Summary and conclusions**

In this work we have conducted a detailed analysis of the dynamics of fluctuations in the
quantum asymmetric simple exclusion process (Q-ASEP) with periodic boundary conditions.
We have demonstrated that fluctuations of the fermionic degrees of freedom obey evolution
equations of Lindblad type, and shown that the corresponding Lindbladians can be repre-
sented as non-Hermitian quantum spin chains that can be expressed in terms of the generators
of a gl(2R) algebra. In case of the Q-SSEP this algebra is a symmetry of the Lindbladian. The
operator space in our model fragments into exponentially many (in system size) sectors that
are invariant under time evolution. This extends recent findings for the average time evolu-
tion to fluctuations. Focusing on the five sectors that describe the late time dynamics of the
Q-SSEP we showed that two of them (c_±1) correspond to a Yang-Baxter integrable model.
Numerical checks of the c_0 block revealed signatures consistent with integrability or weak in-
tegrability breaking in the level-spacing statistics. However, they are not conclusive given the
limitations on the system sizes we were able to simulate. In the particular case of Q-SSEP we
determined the algebraic structure underlying the steady states and slow modes that govern
the late time behaviour. We showed that the former can be understood in terms of "ferromag-
netic" gl(2R) multiplets, while the latter can be viewed as diffusive magnon-like excitations.
We then showed that the dynamics of fluctuations of observables in the Q-SSEP is described by
a closed sets of coupled linear differential-difference equations. The behaviour of the solutions
to these equations is essentially diffusive but with relevant deviations, that at sufficiently late
times and large distances can be described in terms of a continuum scaling limit which we
constructed. We established the applicability of this scaling limit over a significant range of
time and space scales by comparing it to numerical solutions of the corresponding equations
of motion for the lattice model. We finally applied this continuum description to the study of operator spreading at large scales, focusing on out-of-time ordered correlators and operator entanglement. In contrast to operator spreading in random unitary circuits and isolated many-particle Hamiltonian systems, where operators spread ballistically with diffusive corrections, we observe purely diffusive spreading in the Q-SSEP.

Our work raises a number of interesting questions that warrant further enquiry. First, it would be interesting to extend the results reported here for the Q-SSEP to the Q-ASEP. The determination of the steady state manifold will be addressed in a forthcoming publication, but the nature of excited states, dynamics of correlations and dynamics of operator spreading are significantly harder to address. Second, there should be a field theory that gives rise to the continuum description of the equations of motion for correlation functions and it is an open problem to construct it. Third, it would be interesting to extend our analysis to the case of open boundaries with particle injection and extraction. Fourth, one ought to go beyond purely dissipative dynamics and investigate the effects of a Hamiltonian part of the Lindbladian, at least in some limiting cases. Finally, in order to analyze the level-spacing statistics we introduced a conjecture of how to treat additional degeneracies arising from the presence of a higher-rank symmetry. It would be interesting to further test the validity of this conjecture by considering larger system sizes.

**Acknowledgements:** This work was initiated during the thematic trimester program Systems out of equilibrium at the Institut Henri Poincaré. We are grateful to the IHP for hospitality and to Sorbonne University for support. D.B. acknowledges Michel Bauer and Jean-Bernard Zuber for regular discussions. M.M. thanks Tony Jin for initial collaboration on the project and, in particular, for his contributions to analysing the spectrum in the two-particle sector. This work was in part supported by CNRS, by the ENS, by the ANR project "ESQuisses", contract number ANR-20-CE47-0014-01 and by the EPSRC under grant EP/S020527/1.

# A Dictionary: Super-operator versus Hilbert space doubling approach

## A.1 One replica

The basis of states (36) in the Hilbert space doubling approach correspond to physical operators as follows

$$ \|1\rangle_j = n_j, \quad \|2\rangle_j = c_j, \quad \|3\rangle_j = c_j^\dagger, \quad \|4\rangle_j = 1-n_j, \qquad (177) $$

where $n_j = c_j^\dagger c_j$. In terms of the Hubbard operators (40) we have

$$ \begin{aligned} c_j^\dagger &= c_j^\dagger \otimes 1 = E_j^{12} + E_j^{34}, \\ \tilde{c}_j^\dagger &= 1 \otimes c_j^\dagger = -E_j^{13} + E_j^{24}. \end{aligned} \qquad (178) $$

Here the "doubled" operator $\tilde{c}_j^\dagger$ describes the right action of the annihilation operator $c_j$, so that the correspondence with the fermions $f_j^{l/r}$ defined in eqn (29) in the super-operator formalism is

$$ c_j^\dagger \leftrightarrow (f_j^l)^\dagger, \quad \tilde{c}_j^\dagger \leftrightarrow f_j^r. \qquad (179) $$

The two states $\|1\rangle_j$ and $\|4\rangle_j$ form an $su(2)$ doublet with generators $J_j^+ = E^{14}$, $J_j^- = E^{41}$ and $J_j^z = E^{11} - E^{44}$. The identity operator $\|1\rangle_j + \|4\rangle_j$ on this subspace is not an eigenstate of $J_j^z$ but of $J_j^x = J_j^+ + J_j^- = E^{14} + E^{41}$.
## A.2 Two replicas

In the two-replica case the correspondence between the sixteen basis states (53) and physical operators is as follows

<table><thead><tr><th>||0⟩</th><th>||1⟩</th><th>||2⟩</th><th>||3⟩</th><th>||4⟩</th><th>||5⟩</th><th>||6⟩</th><th>||7⟩</th></tr></thead><tbody><tr><td>c<sub>1</sub><sup>†</sup>c<sub>2</sub></td><td>n<sub>1</sub>n<sub>2</sub></td><td>n<sub>1</sub>(1 − n<sub>2</sub>)</td><td>(1 − n<sub>1</sub>)n<sub>2</sub></td><td>(1 − n<sub>1</sub>)(1 − n<sub>2</sub>)</td><td>c<sub>1</sub>c<sub>2</sub><sup>†</sup></td><td>n<sub>1</sub>c<sub>2</sub></td><td>c<sub>1</sub>n<sub>2</sub></td></tr></tbody></table>

<table><thead><tr><th>||8⟩</th><th>||9⟩</th><th>||10⟩</th><th>||11⟩</th><th>||12⟩</th><th>||13⟩</th><th>||14⟩</th><th>||15⟩</th></tr></thead><tbody><tr><td>(1 − n<sub>1</sub>)c<sub>2</sub></td><td>c<sub>1</sub>(1 − n<sub>2</sub>)</td><td>n<sub>1</sub>c<sub>2</sub><sup>†</sup></td><td>c<sub>1</sub><sup>†</sup>n<sub>2</sub></td><td>(1 − n<sub>1</sub>)c<sub>2</sub><sup>†</sup></td><td>c<sub>1</sub><sup>†</sup>(1 − n<sub>2</sub>)</td><td>c<sub>1</sub>c<sub>2</sub></td><td>c<sub>1</sub><sup>†</sup>c<sub>2</sub><sup>†</sup></td></tr></tbody></table>

In terms of the Hubbard operators (54) we have

$$ 
\begin{aligned}
c_{1,j}^{\dagger} &= \varepsilon_j^{0,8} + \varepsilon_j^{1,7} + \varepsilon_j^{2,9} + \varepsilon_j^{6,14} + \varepsilon_j^{10,5} + \varepsilon_j^{11,3} + \varepsilon_j^{13,4} + \varepsilon_j^{15,12}, \\
c_{2,j}^{\dagger} &= \varepsilon_j^{2,6} + \varepsilon_j^{3,8} - \varepsilon_j^{5,9} - \varepsilon_j^{7,14} + \varepsilon_j^{10,2} - \varepsilon_j^{11,0} + \varepsilon_j^{12,4} - \varepsilon_j^{15,13}, \\
\tilde{c}_{1,j}^{\dagger} &= -\varepsilon_j^{1,11} - \varepsilon_j^{2,13} + \varepsilon_j^{5,12} - \varepsilon_j^{6,0} + \varepsilon_j^{7,3} + \varepsilon_j^{9,4} - \varepsilon_j^{10,15} + \varepsilon_j^{14,8}, \\
\tilde{c}_{2,j}^{\dagger} &= -\varepsilon_j^{0,13} - \varepsilon_j^{1,10} - \varepsilon_j^{3,12} + \varepsilon_j^{6,2} + \varepsilon_j^{7,5} + \varepsilon_j^{8,4} + \varepsilon_j^{11,15} - \varepsilon_j^{14,9}.
\end{aligned} 
\quad (180) $$

The correspondence with the fermions $f_j^{a;l/r}$ defined in eqn (29) in the super-operator formalism is

$$ c_{a,j}^{\dagger} \leftrightarrow (f_j^{a;l})^{\dagger}, \quad \tilde{c}_{a,j}^{\dagger} \leftrightarrow f_j^{a;r}. \quad (181) $$

## B Steady states, low-lying states and gaps

To simplify notation we set $J = 1$ in this Appendix.

Via eqn (22) or eqn (24), the (local) Q-SSEP Lindbladian is expressed in terms of the so-called tensor Casimir $C$ for the $gl(2R)$ or $sl(2R)$ algebras, up to simple terms depending on the $C$ charges. Namely, eqn (22) can be written as (we set $D_0 = 1$ to simplify the notation)

$$ \mathcal{L}_{j+1; j}^* = C_{j; j+1} - \frac{1}{2}(C_{j+1} + C_j) - R, \quad (182) $$

with $C_{j; j+1} := \sum_{AB} G_{j+1}^{AB} G_j^{BA}$ the $gl(2R)$ tensor Casimir.

Imagine considering the tensor product of two representations $W_1 \otimes W_2$ with highest weights $\mu_1$ and $\mu_2$ of some Lie algebra (and in particular $gl(2R)$ or $sl(2R)$). The tensor Casimir $C_{12}$ is defined by

$$ C_{12} := \frac{1}{2}(C_{W_1 \otimes W_2} - C_{W_1} \otimes I - I \otimes C_{W_2}), \quad (183) $$

with $C_W$ the Casimir in the representation $W$. Now, if $|\mu_{1,2}\rangle$ are the highest weight vectors in $W_{1,2}$ then $|\mu_1\rangle \otimes |\mu_2\rangle$ is a highest weight vector in $W_1 \otimes W_2$ with highest weights $\mu_1 + \mu_2$. Thus, it is also an eigenvector of the tensor Casimir $C_{12}$ $|\mu_1\rangle \otimes |\mu_2\rangle = \kappa_{\mu_1; \mu_2} |\mu_1\rangle \otimes |\mu_2\rangle$, with $\kappa_{\mu_1; \mu_2} = \frac{1}{2}(\mathrm{Cas}_{\mu_1+\mu_2} - \mathrm{Cas}_{\mu_1} - \mathrm{Cas}_{\mu_2})$, with $\mathrm{Cas}_\mu$ the value of the Casimir in the representation of weight $\mu$. Up to a normalisation factor (to be determined according the chosen convention), $\mathrm{Cas}_\mu = (\mu, \mu+2\rho)$ with $\rho$ half of the sum over the positive roots. Here $(\cdot, \cdot)$ denotes the scalar product on the weight space induced by the Killing form. Hence, $\kappa_{\mu_1; \mu_2} = (\mu_1, \mu_2)$ and

$$ C_{12} |\mu_1\rangle \otimes |\mu_2\rangle = (\mu_1, \mu_2) |\mu_1\rangle \otimes |\mu_2\rangle. \quad (184) $$

To apply this to our problem, we have to check the normalisation of the Casimir in the fermionic representation of $gl(2R)$ we are using. We have $G^{AB} := f^{A\dagger}f^B$ with $\hat{C} := \sum_A G^{AA}$
the total number of f fermions. Recall that $\hat{C}$ is related the C charge via $\hat{C} = C+R$. Using the fermionic anti-commutation rules, we compute:

$$ \sum_{A,B} G^{AB} G^{BA} = \sum_{A,B} f^{A\dagger} f^{B} f^{B\dagger} f^{A} = \hat{C}(2R+1-\hat{C}). $$

This coincides with $(\mu_p, \mu_p + 2\rho)$ for $\mu_p = \epsilon_1 + \cdots + \epsilon_p$, with $\epsilon_i$ ortho-normalized, the highest weight of the p antisymmetric tensor and $\hat{C} = p$ or equivalently $C = p - R$. In particular $(\mu_p, \mu_q) = \min(p, q)$ and this yields the eigien-value of the tensor Casimir in the corresponding tensor product of $gl(2R)$ representations.

Consider now the fully polarized state in the sector with C charge $p-R$ with $0 \le p \le 2R$. This state is an eigenstate of the Q-SSEP Lindbladian. Since $\mathcal{L}_{j+1; j}^* = \mathbb{C}_{j; j+1} - \frac{1}{2}(C_{j+1} + C_j) - R$, its eigenvalue is $p - (p-R) - R = 0$. Thus, as claimed in the main text, the state $\|\Omega_{\mu_p}\rangle := |\mu_p\rangle \otimes \cdots \otimes |\mu_p\rangle$ are zero modes of the Lindbladian $\mathcal{L}^* \|\Omega_{\mu_p}\rangle = 0$. By action of $GL(2R)$ this generates $\|\Omega_x\rangle = x_1 x_2 \cdots x_L$ with $x \in [\mu_p]$ and $\mathcal{L}^* \|\Omega_x\rangle = 0$.

Furthermore, for the tensor product states $|\mu_p\rangle \otimes |\mu_q\rangle$, with $p \ne q$, the eigen-value of the tensor Casimir is $\min(p, q)$. As a consequence of domain wall states $|\mu_p\rangle^{\otimes M} \otimes |\mu_q\rangle^{\otimes (L-M)}$ between sectors with C charges $p-R$ and $q-R$ are eigenstates of the Q-SSEP Lindbladian with eigenvalues $-\frac{1}{2}|p-q|$ (because the eigenvalue is $(\mu_p, \mu_q) - \frac{1}{2}(p+q)$ and $(\mu_p, \mu_q) = \min(p, q)$).

We can use this observation to prove that the spectrum of the Lindbladian is gapped in any "non-homogeneous" sector (i.e. any c such that $\exists j, k$ s.t. $c_j \ne c_k$). Any such sector is made of a series of homogeneous portions separated by 'domain wall edges' with different C-charges on its two ends. Let us denote these domain wall edges by $w_\alpha$, $\alpha = 1, \cdots, d$, where $d$ is the number of homogeneous segments, and denote by $c_\alpha^\pm$ the two distinct C-charges on either side of the corresponding edge. Since the total Q-SSEP Lindbladian is the sum of Lindbladian attached to each of the edges, and since any Lindbladian is a non-positive operator, we have the operator inequality

$$ \mathcal{L}_{\text{ssep}}^* \le \sum_{\alpha} \mathcal{L}_{j_{\alpha}, j_{\alpha}+1}^*, \qquad (185) $$

where $(j_\alpha, j_\alpha + 1)$ are the two vertices connected by the edge $w_\alpha$. The two representations $[\mu_{R+c_\alpha^+}]$ and $[\mu_{R+c_\alpha^-}]$ are attached to these vertices, so that $\mathcal{L}_{j_\alpha, j_\alpha+1}^*$ is acting on their tensor product $[\mu_{R+c_\alpha^+}] \otimes [\mu_{R+c_\alpha^-}]$, as in eqn (182). As a consequence, its maximum eigenvalue is that of the tensor Casimir on $[\mu_{R+c_\alpha^+}] \otimes [\mu_{R+c_\alpha^-}]$. The latter is diagonalized by decomposing the tensor product $[\mu_{R+c_\alpha^+}] \otimes [\mu_{R+c_\alpha^-}]$ into $gl(2R)$ irreducible representations, $[\mu_{R+c_\alpha^+}] \otimes [\mu_{R+c_\alpha^-}] = \oplus_\lambda [\mu_\lambda]$. From eqn (183), its maximal eigenvalue is given by the representation $\mu_\lambda$ maximising its Casimir. Since $\text{Cas}_\mu = (\mu, \mu + 2\rho) = |\mu + \rho|^2 - |\rho|^2$ the maximum of the Casimir among the representation $[\mu_\lambda]$ occurring in the decomposition of the tensor product $[\mu_{R+c_\alpha^+}] \otimes [\mu_{R+c_\alpha^-}]$ is for $\mu_\lambda = \mu_{R+c_\alpha^+} + \mu_{R+c_\alpha^-}$. All other weights $\mu_\lambda + \rho$ in this decomposition have smaller norm. Hence, from the computation of the previous paragraph, we have

$$ \mathcal{L}_{j_{\alpha}, j_{\alpha}+1}^{*} \leq -\frac{1}{2}|c_{\alpha}^{+} - c_{\alpha}^{-}|. \qquad (186) $$

This proves that all inhomogeneous sectors are gapped with a gap given by the sum of the gaps associated with each domain wall edge.

Let us present here the proof that the one-magnons states $\|p; \gamma, \mu_n\rangle$ defined in eqn (78) are indeed eigenstates of the Q-SSEP Lindbladian $\mathcal{L} = \mathcal{L}^*$. Since the local charges $G_j^{-\gamma}$ is one of the conserved currents, commuting it with the Lindbladian yields,

$$ [\mathcal{L}^*, G_j^{-\gamma}] = V_j^{-\gamma} - V_{j-1}^{-\gamma}, \qquad (187) $$

with

$$ V_j^{-\gamma} := -(\vec{\gamma} \cdot \vec{H}_{j+1}) G_j^{-\gamma} + G_{j+1}^{-\gamma} (\vec{\gamma} \cdot \vec{H}_j) + 2 \sum_{\alpha: \alpha-\gamma \in \text{roots}} G_{j+1}^{-\alpha} G_j^{\alpha-\gamma}. $$
For $||p; \gamma, \mu_n||$ to be an eigenstate of $\mathcal{L}^*$, the last sum in the above expression has to annihilate it. This amounts to say that either $\alpha - \gamma$ or $-\alpha$ have to be a positive root for all $\alpha$ such that $\alpha - \gamma$ is a root. We now show that this is indeed true iff $\gamma$ is a simple root. For a specific choice of the Weyl chamber, the positive (resp. negative) roots of $sl(2R)$ can be written as $\epsilon_k - \epsilon_l$ with $k < l$ (resp. $k > l$) with $\epsilon_k$ orthonormal basis. Choose $\gamma = \epsilon_i - \epsilon_j$, $i < j$. Then, either $\alpha = \epsilon_k - \epsilon_j$ (and $\alpha - \gamma = \epsilon_k - \epsilon_i$) or $\alpha = \epsilon_i - \epsilon_l$ (and $\alpha - \gamma = \epsilon_j - \epsilon_l$). It is then easy to check (by analysing each four cases) that we fulfil the condition that either $\alpha - \gamma$ or $-\alpha$ is positive, for all $\alpha$, if and only if $j = i + 1$, i.e. iff $\gamma = \epsilon_i - \epsilon_{i+1}$ is a simple root. Thus, for $\gamma = \alpha_i$ a simple root, we have

$$V_j^{-\alpha_i} \|\Omega_\mu\rangle = (\alpha_i, \mu)(G_{j+1}^{-\alpha_i} - G_j^{-\alpha_i}) \|\Omega_\mu\rangle,$$

for any highest weight vector $\mu$. Finally, for $\mu$ one of the fundamental weight of $sl(2R)$, the scalar product $(\alpha_i, \mu)$ is equal either to 0 or to 1.

Finally, we describe here the computations needed to compare the two-replica Lindbladian with the so(6) integrable spin chain Hamiltonian in the vector representation. The so(6) generators in the vector representation can be written as

$$J_{ij} = E_{ij} - E_{j'i'},$$

with $E_j = |i\rangle\langle j|$ and $i, j = 1, \dots, 6$. The permutation $P$ and trace $Q$ operators are

$$P = \sum_{ij} E_{ij} \otimes E_{ji}, \quad Q = \sum_{ij} E_{ij} \otimes E_{i'j'}.$$

The tensor Casimir is $C = \sum_{ij} J_{ij} \otimes J_{ji}$. Computing, we get

$$\begin{aligned}
C &= \sum_{ij} (E_{ij} - E_{j'i'}) \otimes (E_{ji} - E_{i'j'}) \\
&= E_{ij} \otimes E_{ji} - E_{ij} \otimes E_{i'j'} - E_{j'i'} \otimes E_{ji} + E_{j'i'} \otimes E_{i'j'} \\
&= 2(P - Q).
\end{aligned}$$

Hence the Lindbladian is proportional to $\sum_j (P_{j;j+1} - Q_{j;j+1})$, up to a constant, which differs from the known Yang-Baxter so(6) integrable Hamiltonian $H \propto \sum_j (2P_{j;j+1} - Q_{j;j+1})$.

## C $GL(2R)$ and $U(L)$ interplay on steady states

The analysis of the R replica Q-SSEP zero modes presented in the main text tells us how they decompose into representations of $gl(2R)$: steady states in the sector with C charge equals to $p-R$ form an irreducible $gl(2R)$ isomorphic to that associated with a rectangular Young tableau with $p$ boxes vertically and $L$ boxes horizontally. The purpose of this Appendix is to understand the nature of these steady states in terms of the physical degrees of freedom.

The system Hilbert space, generated by the physical fermions $c_j^\dagger$, decomposes into the direct sum of sub-spaces with fixed particle numbers:

$$\mathcal{H}_{\text{sys}} = \Lambda_1 \oplus \Lambda_2 \oplus \cdots \oplus \Lambda_L,$$

where $L$ is the number of sites. Each such subspace $\Lambda_m$, with $m$ particles, forms an irreducible $SU(L)$ representation isomorphic to arank $m$ antisymmetric tensor in dimension $L$. The two representations with $m = 0$ and $m = L$ are isomorphic, and isomorphic to the scalar representation, i.e. $\Lambda_0 \cong \Lambda_L \cong C$. The others are not isomorphic. Thus, as a $SU(L)$ module, the system Hilbert space decomposes as

$$\mathcal{H}_{\text{sys}} \equiv \Lambda_0 \otimes \mathbb{C}^2 \oplus \Lambda_1 \oplus \Lambda_2 \cdots \oplus \Lambda_{L-1},$$
where the first factor encodes for the multiplicity of the trivial representation, i.e. $\Lambda_0 \otimes \mathbb{C}^2 \equiv \Lambda_0 \oplus \Lambda_L$. We shall denote by $\mathbb{P}_m$ the projector on the subspace $\Lambda_m$ with $m$ particles.

Let us first recall a few properties of the Q-SSEP dynamics and invariant measure:

1. For each realisation of the noise, the dynamics preserves the total number of particles $\hat{M} = \sum_j c_j^\dagger c_j$ and the spectrum of the density matrix (in particular it preserves pure states).

2. The Q-SSEP invariant measure is $SU(L)$ invariant. On each sector with fixed number of particle, it is that induced by the $SU(L)$ Haar measure on that orbit, see Ref. [41].

As a consequence the R replica steady states have to be covariant under the $SU(L)$ action and this is the route that we shall follow to understand them. We shall only discuss in details the cases of $R = 1$ or $R = 2$ replicas.

* For one replica ($R = 1$).

Here the symmetry algebra is $gl(2)$. In the $c_0$ sector, the zero modes form a multiplet of dimension $L + 1$ (spin $L/2$). The $c_{\pm 1}$ sectors have each a single zero mode, scalar under the $gl(2)$ symmetry. The $SU(L)$ invariance of the steady measure implies that the steady mean density matrix $\bar{\rho} := E_{\text{inv}}[\rho]$ is invariant under the $SU(L)$ adjoint action. As a consequence, by Schur's lemma, it has to read:

$$ \bar{\rho} := E_{\text{inv}}[\rho] = q_0 \mathbb{P}_0 \otimes \bar{\sigma}_{2 \times 2} + q_1 \mathbb{P}_1 + \cdots + q_{L-1} \mathbb{P}_{L-1}, $$

with $\bar{\sigma}_{2 \times 2}$ a $2 \times 2$ matrix acting on the $\mathbb{C}^2$ factor of the component $\Lambda_0 \otimes \mathbb{C}^2$. The diagonal entries corresponding to the $L+1$ projectors $\mathbb{P}_m$, $m = 0, \cdots, L$ (with $\mathbb{P}_0$ and $\mathbb{P}_L$ the two diagonal entries in $\bar{\sigma}_{2 \times 2}$) are zero modes belonging to the $c_0$ sector. The off-diagonal ones in $\bar{\sigma}_{2 \times 2}$ belong in the $c_{\pm 1}$ sector and couple the empty ($m = 0$) and full ($m = L$) states. These off-diagonal terms are non-zero if and only if there were already present in the initial mean density matrix (because these terms are frozen under the Q-SSEP dynamics, cf. eqn (69)).

The $gl(2)$ symmetry acts on the diagonal parts (the off-diagonal terms are scalar under $gl(2)$), i.e. it acts on the projectors $\mathbb{P}_m$, $m = 0, \cdots, L$. It is easy to verify that,

$$ J^+ \cdot \mathbb{P}_m = (m+1)\mathbb{P}_{m+1}, \quad J^- \cdot \mathbb{P}_m = (L+1-m)\mathbb{P}_{m-1}, \quad J^z \cdot \mathbb{P}_m = (2m-L)\mathbb{P}_m, $$

for $m = 0, \cdots, L$. Thus, these projectors form a spin $L/2$ representation of $gl(2)$ as they should. Phrased differently, the $gl(2)$ action on one replica steady state induces moves between steady states with different particle numbers.

* For two replicas ($R = 2$).

By $SU(L)$ invariance, the second moment $\bar{\rho}^{(2)} := E_{\text{inv}}[\rho \otimes \rho]$ is also $SU(L)$ invariant under the diagonal adjoint action. Hence, by Schur's lemma, it is block diagonal on irreducible $sl(2R)$ sub-representations of $\mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}}$. That is: it decomposes as the sum over the projectors of the irreducible $SU(L)$ components in $\mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}}$.

To decompose $\mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}}$ as a sum of $SU(L)$ irreducible representations is a (nice) exercise in group theory$^7$. The output is the following decomposition:

$$ \mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}} \equiv (L+3)\Lambda_0 + \bigoplus_{k=1}^{L-1} (L+2)\Lambda_k + \bigoplus_{p=1}^{L-1} \bigoplus_{q=0}^{L-p-1} (q+1)\Lambda_{p+q;p}, \quad (188) $$

where $\Lambda_r$ is the r-th fundamental representation of $SU(L)$, with Young tableau made of one column with r boxes, and $\Lambda_{r;s;\cdots}$ the representation with Young tableau with r boxes on the first

<footnote>$^7$We thank Jean-Bernard Zuber for his help in clarifying this problem.</footnote>
column, s boxes on the second column, etc. The proof of the above formula (188) follows from using rules for decomposing of the tensor product of two fundamental SU(L) representations: for $r \ge s$, $\Lambda_r \otimes \Lambda_s = \bigoplus_{k=0}^{\min(s,L-r)} \Lambda_{r+k;s-k}$.

Let us give a few examples (at least to check the general formula):

$$- \text{For } L = 2, \mathcal{H}_{\text{sys}} = \Lambda_0 \oplus \Lambda_1 \oplus \Lambda_2, \text{ with } \Lambda_2 \equiv \Lambda_0, \text{ and}$$

$$\mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}} = 5\Lambda_0 \oplus 4\Lambda_1 \oplus \Lambda_{1;1}.$$

$$- \text{For } L = 3, \mathcal{H}_{\text{sys}} = \Lambda_0 \oplus \Lambda_1 \oplus \Lambda_2 \oplus \Lambda_3, \text{ with } \Lambda_3 \equiv \Lambda_0, \text{ and}$$

$$\mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}} = 6\Lambda_0 \oplus 5\Lambda_1 \oplus 5\Lambda_2 \oplus \Lambda_{1;1} \oplus 2\Lambda_{2;1} \oplus \Lambda_{2;2}.$$

$$- \text{For } L = 4, \mathcal{H}_{\text{sys}} = \Lambda_0 \oplus \Lambda_1 \oplus \Lambda_2 \oplus \Lambda_3 \oplus \Lambda_4, \text{ with } \Lambda_4 \equiv \Lambda_0, \text{ and}$$

$$\begin{aligned} \mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}} = & 7\Lambda_0 \oplus 6\Lambda_1 \oplus 6\Lambda_2 \oplus 6\Lambda_3 \\ & \oplus \Lambda_{1;1} \oplus 2\Lambda_{2;1} \oplus 3\Lambda_{3;1} \\ & \oplus \Lambda_{2;2} \oplus 2\Lambda_{3;2} \\ & \oplus \Lambda_{3;3}. \end{aligned}$$

We see a simple (triangular) pattern emerging. Formula (188) is a compact way to encode it.

As a consequence, the quadratic fluctuations $\mathbb{E}_{\text{inv}}[\rho \otimes \rho]$ in the steady measure decompose as:

$$\mathbb{E}_{\text{inv}}[\rho \otimes \rho] = \mathbb{P}_0 \otimes \bar{\sigma}^{(0)}_{(L+3) \times (L+3)} + \sum_{m=1}^{L-1} \mathbb{P}_m \otimes \bar{\sigma}^{(m)}_{(L+2) \times (L+2)} + \sum_{p=1}^{L-1} \sum_{q=0}^{L-p-1} \mathbb{P}_{p+q; q} \otimes \bar{\sigma}^{(p+q; q)}_{(q+1) \times (q+1)}, \quad (189)$$

where $\mathbb{P}_{r;s;...}$ are the projectors on $\Lambda_{r;s;...}$. The square matrices $\bar{\sigma}^{(*)}_{n \times n}$ are $n \times n$ matrices parametrizing the two replica steady states. The total number of parameters (without taking the trace normalization into account) is

$$\begin{aligned} \# \text{ parameters} &= (L+3)^2 + (L+2)^2(L-1) + \sum_{k=1}^{L-1} k^2(L-k) \\ &= \frac{(L+1)(L+2)^2(L+3)}{12} + 2\frac{(L+1)(L+2)(L+3)}{6} + 2. \end{aligned}$$

It matches the total dimension of the zero mode spaces that we identify in the main text using $gl(4)$ arguments. This shows that we didn't miss any zero modes using these arguments and that, given their C charge, these zero modes form irreducible $gl(4)$ multiplets.

The above construction and counting do not distinguish between the different C charge sectors. To do that we have to refine the decomposition (188) of $\mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}}$ by taking into account the $u(1)$ charges in $U(L)$ and not only the $SU(L)$ content. Indeed, some of the Young tableaux appearing in the decomposition of $\mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}}$ are of the form $\Lambda_{L;s}$ with a first column full with L boxes. As $SU(L)$ representation $\Lambda_{L;s}$ and $\Lambda_s$ are isomorphic but not as $U(L)$ representations. Going from $\Lambda_s$ to $\Lambda_{L;s}$ amounts to add +1 to all occupation numbers $n_j$ (recall that the $U(L)$ Cartan generators are the occupation numbers, $H_j = n_j$, whereas those of $SU(L)$ are $h_j = H_j - H_{j-1}$). We denote $\Lambda_s^+ := \Lambda_{L;s}$ (for instance $\Lambda_0^+ = \Lambda_L$) and $\Lambda_s^{++} := \Lambda_{L;L;s}$.

Using $\Lambda_r \otimes \Lambda_s = \bigoplus_{k=0}^{\min(s,L-r)} \Lambda_{r+k;s-k}$, for $r \ge s$, the detailed analysis of the decomposition of $\mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}}$ is:
- $\Lambda_0$ comes either (a) from the tensor products $\Lambda_0 \otimes \Lambda_0 \rightsquigarrow \Lambda_0$, $\Lambda_0 \otimes \Lambda_L \rightsquigarrow \Lambda_0^+$, $\Lambda_L \otimes \Lambda_0 \rightsquigarrow \Lambda_0^+$ and $\Lambda_L \otimes \Lambda_L \rightsquigarrow \Lambda_0^{++}$, or (b) from $\Lambda_k \otimes \Lambda_{L-k} \rightsquigarrow \Lambda_0^+$, with $k = 1, \dots, L-1$;
- $\Lambda_k$, for any $k = 1, \dots, L-1$, comes either (a) from the tensor products $\Lambda_r \otimes \Lambda_s \rightsquigarrow \Lambda_{r+s=k}$ with
$r = 0, \dots, [\frac{k}{2}]$, $s = [\frac{k+1}{2}]$, $\dots, k$ or with the role of $r$ and $s$ exchanged, and there are $(k+1)$ such possibilities, or (b) from $\Lambda_r \otimes \Lambda_s \rightsquigarrow \Lambda_{r+s=L+k}$ with $r, s = k, \dots, L$, and there are $(L+1-k)$ such possibilities.

$-\Lambda_{p+q;q}$ comes from $\Lambda_r \otimes \Lambda_s \rightsquigarrow \Lambda_{p+q;q}$, with $r+s = 2p+q$ and $r, s = p, \dots, p+q$. They are all of the form $\Lambda_{p+q;q}$ and not $\Lambda_{p+q;q}^+$ because their Young tableaux have two columns (so that a full column could not have been produced from the tensor product $\mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}}$).

Thus, the $U(L)$ decomposition of $\mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}}$ is

$$ \mathcal{H}_{\text{sys}} \otimes \mathcal{H}_{\text{sys}} \equiv (\Lambda_0 + (L+1)\Lambda_0^+ + \Lambda_0^{++}) + \bigoplus_{m=1}^{L-1} ((m+1)\Lambda_m + (L+1-m)\Lambda_m^+) \\ + \bigoplus_{p=1}^{L-1} \bigoplus_{q=0}^{L-p-1} (q+1)\Lambda_{p+q;p}. \qquad (190) $$

In particular, $\Lambda_0$ is the empty state and $\Lambda_0^{++}$ the totally (on both replicas) full state.

The zero modes in the $c_{+1}$ sector are the intertwiners from $\Lambda_r$ to $\Lambda_r^+$, or from $\Lambda_r^+$ to $\Lambda_r^{++}$ (i.e. they are operators intertwining representations whose Young tableaux differ by a full column). Explicitly, in the $c_{+1}$ sector, we have the intertwiners:

$$ \begin{aligned} \Lambda_0 &\rightarrow \Lambda_0^+ [(L+1)\text{--times}], \quad \Lambda_0^+ \rightarrow \Lambda_0^{++} [(L+1)\text{--times}] \\ \Lambda_m &\rightarrow \Lambda_m^+ [(m+1)(L+1-m)\text{--times}], \quad m = 1, \dots, L-1. \end{aligned} $$

The total number of $C=1$ intertwiners is:

$$ \# \text{off.diag} = 2(L+1) + \sum_{m=1}^{L-1} (m+1)(L+1-m) = \frac{(L+1)(L+2)(L+3)}{6}. $$

It matches the dimension of the $c_{+1}$ (or $c_{-1}$) zero modes, as it should. The $c_{-1}$ intertwiners are the transposed. The $c_{\pm 2}$ intertwiners are the empty-to-full interwiners $\Lambda_0^{++} \rightarrow \Lambda_0$ and $\Lambda_0 \rightarrow \Lambda_0^{++}$, respectively.

As a consequence, the decomposition (189) of $\mathbb{E}_{\text{inv}}[\rho \otimes \rho]$ can be refined. The matrices $\bar{\sigma}^{(0)}_{(L+3)\times(L+3)}$ can be decomposed into blocks of sizes 1, $L+1$, 1, according to the decomposition $\Lambda_0 + (L+1)\Lambda_0^+ + \Lambda_0^{++}$ in eqn (190). The matrices $\bar{\sigma}^{(m)}_{(L+1)\times(L+1)}$ can be decomposed into blocks of sizes $m+1$, $L+1-m$, according to the decomposition $(m+1)\Lambda_m + (L+1-m)\Lambda_m^+$ in eqn (190). The off-diagonal blocks correspond to the $c_{\pm 1}$ sector (and $c_{\pm 2}$). The $gl(4)$ algebra act (faithfully on each sector of given charge $C = 0, \pm 1, \pm 2$) on these intertwiners.

## D Unfolding and level-spacing for a Poisson process.

Given $N$ i.i.d random variables $X_i$ (the eigenvalues of an integrable Hamiltonian) taking values in $\mathbb{R}$ with density $p_X(x)$, we perform a local change of variables such that the new variable $\hat{x}(x)$ describes the average number of old variables below $x$,

$$ x \rightarrow \hat{x}(x) = N \int_{-\infty}^{x} p_X(x') dx'. \qquad (191) $$

This procedure is called "unfolding the spectrum" and it ensures that the density of the new variables $\hat{p}_X(\hat{x})d\hat{x} = p_X(x)dx$ is indeed uniform, $\hat{p}_X(\hat{x}) = 1/N$.

The probability to find $N_s = k$ of the new random variables in the interval $[0, s]$ is now is given by the Poisson distribution

$$ \mathbb{P}[N_s = k] = \frac{(\lambda s)^k}{k!} e^{-\lambda s}, \qquad (192) $$
with average "even rate" $\lambda = 1$. We can define the probability $p_S(s)$ to observe a spacing $s$ between adjacent eigenvalues by

$$ \mathbb{P}[N_s = 0] = \int_s^\infty p_S(s')ds'. \qquad (193) $$

Deriving w.r.t $s$ provides us with the exponential distribution $p_S(s) = e^{-s}$.

## E Two point functions in the Lindbladian formulation

As the state $\|1_2\rangle$ belongs to the sector $c_0$ we can project the equations of motion (96) to this sector. This means that at each site we need to consider only the states $\|0\rangle_j, \|1\rangle_j, \|2\rangle_j, \|3\rangle_j, \|4\rangle_j, \|5\rangle_j$, which form a six dimensional gl(4) multiplet with C charge zero. A representation of the gl(4) generators on the corresponding subspace is obtained by retaining only the Hubbard operators $\epsilon_j^{k,l}$, with $k, l = 1, \dots, 6$ in eqn (55). The two gl(2) sub-algebras associated with each of the replicas are spanned by generators leaving operators on the opposite replica invariant. The action of $J_1^\alpha$ (resp. $J_2^\alpha$) is as in the 1-replica case, i.e.

$$ \begin{aligned} J_{1,j}^+ &= E_j^{14} \otimes 1 = G_j^{12}, \quad J_{2,j}^+ = 1 \otimes E_j^{14} = G_j^{43}, \\ J_{1,j}^z &= (E_j^{11} - E_j^{44}) \otimes 1 = G_j^{11} - G_j^{22}, \\ J_{2,j}^z &= 1 \otimes (E_j^{11} - E_j^{44}) = G_j^{44} - G_j^{33}. \end{aligned} \qquad (194) $$

On this subspace we have $2n_{a,j} - 1 = J_{a,j}^z$, which allows us to rewrite (96) as

$$ \frac{d}{dt} \langle 1_2 \|(J_{1,j}^z + 1)(J_{2,k}^z + 1)\| \rho^{(2)}(t) \rangle = \langle 1_2 \|(J_{1,j}^z + 1)(J_{2,k}^z + 1) \hat{\mathcal{L}}_{2,\epsilon_0} \| \rho^{(2)}(t) \rangle. \qquad (195) $$

We now observe that while the identity $\|1_2\rangle$ is not an eigenstate of $J_{a,j}^z$, it is in fact a simultaneous eigenstate of $J_{a,j}^x = J_{a,j}^+ + J_{a,j}^-$, $a = 1, 2$. To proceed we therefore "rotate the quantization axis" and work in a basis of eigenstates of $J_{a,j}^x$, rather than $J_{a,j}^z$. This basis is given by

$$ \begin{aligned} \|1\rangle_j &= \frac{1}{2}, \quad \|2\rangle_j = n_{1,j} - \frac{1}{2}, \quad \|3\rangle_j = n_{2,j} - \frac{1}{2}, \quad \|4\rangle_j = 2(n_{1,j} - \frac{1}{2})(n_{2,j} - \frac{1}{2}), \\ \|0\rangle_j &= c_{1,j}^\dagger c_{2,j}, \quad \|5\rangle_j = c_{1,j} c_{2,j}^\dagger. \end{aligned} $$

The eigenvalues of $J_{1,j}^x, J_{2,j}^x$ on these states are $(1, 1), (-1, 1), (1, -1), (-1, -1), (0, 0)$ and $(0, 0)$ respectively. We note that $\|1\rangle = 2^L \otimes_{j=1}^L |1\rangle_j$. It is now apparent that $J_{a,j}^z$ as well as $S_j^\pm$ defined in (83) create single-particle excitations over the state $\|1\rangle$. As the Lindbladian $\hat{\mathcal{L}}_{2,\epsilon_0}$ commutes with $J_{1,2}^x$, the corresponding particle numbers are good quantum numbers and (195) as well as its analog for $g_{-}(j, k; t)$ can therefore be viewed as Schrödinger equations in the two-particle sector.

## F Spectrum of the two-replica Lindbladian in the two-particle sector

As a consequence of the free dynamics governing the behaviour of Q-SSEP model for every realization of the noise, the equations for the n-point functions wrt the noise close, as discussed in the main text. In this appendix we will focus on the spectral properties in the case of two
point functions. Furthermore, we will specialize our discussion to the closed loop diagrams, on which, once again, the dynamics closes. Correlation functions associated with disconnected diagrams decay exponentially fast and are, as such, of lesser importance when considering the hydrodynamics of the model.

On the level of two-point functions it proves advantageous to define two sets of observables, comprised of diagonal and off-diagonal correlations

$$A_{jk}^{\pm} := \mathbb{E}[G_{jj}G_{kk} \pm G_{kj}G_{jk}]. \qquad (196)$$

Interestingly, the dynamics for correlations (196) decouples and is governed by equations

$$\frac{d}{dt} A_{ij}^{\pm} = (\Delta_i + \Delta_j) A_{ij}^{\pm} \pm 2 \delta_{i,j} (A_{i,i+1}^{\pm} + A_{i-1,i}^{\pm}) \mp 2 (\delta_{j,i+1} + \delta_{i,j+1}) A_{ij}^{\pm}. \qquad (197)$$

We can consider equations (197) as a two particle problem, where the positions of two particles is associated with indices i and j, which can be solved through the scattering wave ansatz

$$A_{j,k}^{\pm}(t) = (e^{i(p_j+q_k)} + S^{\pm}(p,q)e^{i(q_j+p_k)})e^{\omega t}, \quad \text{for } j < k,$$

$$A_{j,j}^{\pm}(t) = C^{\pm}(p,q)e^{i(p+q)_j+\omega t},$$

$$A_{j,k}^{\pm}(t) = (e^{i(p_k+q_j)} + S(p,q)e^{i(q_k+p_j)})e^{\omega t}, \quad \text{for } j > k.$$

The decay rate $\omega$ is governed by the non-interacting, i.e. purely diffusive, part of equations (197) ($\Delta_i + \Delta_j$), which yields

$$\omega(p,q) = 2(\cos p + \cos q - 2). \qquad (198)$$

Coefficients $C^{\pm}(p,q)$ and $S^{\pm}(p,q)$ can be obtained by considering the equations for $A_{i,i+1}^{\pm}(p,q)$ and $A_{i,i}^{\pm}(p,q)$, yielding

$$S^{+}(p,q) = -e^{i(q-p)} \frac{(\cos p + \cos q)(e^{ip} + e^{-iq} + 2) - 4(1 + \cos(p+q))}{(\cos p + \cos q)(e^{-ip} + e^{iq} + 2) - 4(1 + \cos(p+q))}, \qquad (199)$$

$$C^{+}(p,q) = \frac{4}{(\omega + 4)}(e^{-ip} + e^{iq} + S^{+}(p,q)(e^{-iq} + e^{ip})), \qquad (200)$$

and

$$S^{-}(p,q) = -\frac{e^{i(p+q)} - 2e^{iq} + 1}{e^{i(p+q)} - 2e^{ip} + 1}, \qquad (201)$$

$$C^{-}(p,q) = 0. \qquad (202)$$

Note that the coefficients $S^{-}$ and $C^{-}$, as well as the equations for $A_{ij}^{-}(t)$ match perfectly the Bethe ansatz coefficients and equations for the two particle sector of the Heisenberg XXX chain. Finally, periodicity implies the quantization conditions

$$S^{\pm}(p,q) e^{ipL} = S^{\pm}(p,q) e^{-iqL} = 1. \qquad (203)$$

Due to the translational invariance the total momentum takes the values $k = p+q = \frac{2\pi n}{L}$, $n \in \mathbb{Z}$, which means that we can reduce the eigenvalue problem at the fixed total momentum k to a single transcendental equation. The quantization condition can be represented in terms of a new variable $\xi$

$$z := e^{iu/2}, \quad \xi := z + z^{-1} = 2 \cos(u/2), \qquad (204)$$
which parametrizes the difference of particle momenta $u = p-q$, and gives a compact expression for the eigenvalue

$$ \omega_n(z) = 2(z + z^{-1})c_n - 4 = 2(c_n \xi - 2), \qquad (205) $$

where we introduced $c_n = \cos(\frac{2\pi n}{L})$, and eigenfunctions<fnref n="8" />

$$ A_{ij}(t) = e^{ik(i+j)/2} [z^{i-j+N} + (-)^n z^{-(i-j+N)}] e^{\omega t}, \quad \text{for } i < j, $$

$$ A_{jj}(t) = e^{ikj} 4 \left[ \frac{z^{N-1} + (-)^n z^{1-N}}{z^{-1} + z} \right] e^{\omega t}. $$

For $k = \frac{2\pi n}{L} \neq \pi$ (i.e. for $\cos(k/2) \neq 0$), the S-matrix, for a total moment $k = \frac{2\pi n}{L}$, reads<fnref n="9" />

$$ S_n^+(z) = \frac{z^{-1} R_n(z)}{z R_n(z^{-1})}; \quad R_n^+(z) := 4c_n - (z+z^{-1})(1+z c_n), \qquad (206) $$

$$ S_n^-(z) = \frac{c_n - 1/z}{c_n - z} \qquad (207) $$

and since $e^{ikL/2} = (-1)^n$, the quantization condition reduces to

$$ S_n^{\pm}(z) z^L = (-1)^n. \qquad (208) $$

The condition can be neatly represented in terms of Chebyshev-like polynomials

$$ P_p^+(\xi) := z^{2p} + z^{-2p}, \quad Q_p^+(\xi)\xi := z^{2p-1} + z^{1-2p} \qquad (209) $$

$$ Q_p^-(\xi)(z-z^{-1})\xi := z^{2p} - z^{-2p}, \quad P_p^-(\xi)(z-z^{-1}) := z^{2p+1} - z^{-1-2p}. \qquad (210) $$

Namely, for chains of even length $L = 2N$, we have to consider the following cases for $A^+$:

* $N = 2M$ even, $n$ even

$$ (4c_n - \xi) \xi Q_M^+(\xi) = c_n \xi P_M^+(\xi), \qquad (211) $$

* $N = 2M + 1$ odd, $n$ even

$$ (4c_n - \xi) P_M^+(\xi) = c_n \xi^2 Q_{M+1}^+(\xi), \qquad (212) $$

* $N = 2M + 2$ even, $n$ odd

$$ (4c_n - \xi) P_M^-(\xi) = c_n \xi^2 Q_{M+1}^-(\xi), \qquad (213) $$

* $N = 2M + 1$ odd, $n$ odd

$$ (4c_n - \xi) \xi Q_M^-(\xi) = c_n \xi P_M^-(\xi). \qquad (214) $$

And similarly for $A^-$:

* $N = 2M$ even

$$ c_n P_M^+(\xi) = Q_M^+(\xi), \qquad (215) $$

* $N = 2M + 1$ odd

$$ c_n Q_{M-1}^+(\xi) = P_M^+(\xi). \qquad (216) $$

A numerical analysis of these equations reveals that solutions $z$ always lie on the unit circle $|z|=1$, except for two reciprocal values, see Fig. 12.

<footnote n="8">⁸up to the multiplicative constant</footnote>

<footnote n="9">⁹Note that for $A^+$ the special case with $k = \pi$ has to be analysed separately.</footnote>
<figure><img src="image_17.png" /><figcaption>Figure 12: Solutions of the equation (211) for $N = 12$ and $n = 4, 6, 8$. All of the solutions lie on the unit circle (red dashed line) except for two solutions at any $n$, if $n \neq N/2$.</figcaption></figure>

# G Derivation of the time evolution of n-point bubbles

Here we show how to arrive at (143) in the main text that describes time evolution of the n-point bubble $g_n(l_1, ..., l_n) = [G_{l_1 l_2} G_{l_2 l_3} ... G_{l_n l_1}]$, where [...] is shorthand for $\mathbb{E}[...]$. The time evolution of the two-point functions $G_{ij} = \text{Tr}(\rho c_j^\dagger c_i)$ is given by the stochastic differential equation [41]

$$dG_{ij} = (\delta_{ij}(G_{i+1,i+1} + G_{i-1,i-1}) - 2G_{i,j})Jdt + i(G_{i,j-1}d\bar{W}^{j-1} + G_{i,j+1}dW^j - G_{i-1,j}dW^{i-1} - G_{i+1,j}d\bar{W}^i), \quad (217)$$

where we used the convention for the Brownian motion $dW$ from (3) with $p = q$ (Q-SSEP), which makes $dW$ dimensionless by introducing the rate parameter $J$. Using Ito convention we have

$$\frac{[dG_{ij}dG_{kl}]}{Jdt} = \delta_{jk}[G_{i,j\pm1}G_{k\pm1,l}] + \delta_{il}[G_{i\pm1,j}G_{k,l\pm1}] - (\delta_{i,k\pm1} + \delta_{j,l\pm1})[G_{il}G_{kj}], \quad (218)$$

where we introduced the shorthand notation $G_{i,j\pm1}G_{k\pm1,l} := G_{i,j\pm1}G_{k\pm1,l} + G_{i,j-1}G_{k-1,l}$ and similarly $\delta_{i,k\pm1} := \delta_{i,k\pm1} + \delta_{i,k-1}$.

## G.1 Three-point bubble

To illustrate the derivation, we will first explicitly derive the evolution of the three-point bubble,

$$dg_3(i, j, k) = [dG_{ij}G_{jk}G_{ki}] + [G_{ij}dG_{jk}G_{ki}] + [G_{ij}dG_{jk}G_{ki}] + [dG_{ij}dG_{jk}G_{ki}] + [G_{ij}dG_{jk}dG_{ki}] + [dG_{ij}G_{jk}dG_{ki}]. \quad (219)$$

In particular we need

$$\frac{[dG_{ij}]}{Jdt} = \delta_{ij}[G_{i\pm1,i\pm1}] - 2[G_{ij}], \quad (220)$$

$$\frac{[dG_{ij}dG_{jk}G_{ki}]}{Jdt} = [G_{i,j\pm1}G_{j\pm1,k}G_{ki}] + \delta_{ik}[G_{i\pm1,j}G_{j,i\pm1}G_{ii}] - (\delta_{i,j\pm1} + \delta_{j,k\pm1})[G_{ik}G_{kk}G_{ki}]. \quad (221)$$
The terms in (219) without Kronecker deltas yield discrete Laplacians ($\Delta_i + \Delta_j + \Delta_k$)g$_3(i, j, k)$. The contact terms (terms with Kronecker delta) are built by products of two bubbles which appear in a symmetric manner. We therefore introduce the notation

$$g_3(i|j, k) := [G_{ii}G_{jk}G_{kj}] = [\begin{smallmatrix} \circled{2} \\ \bullet_i & \bullet_j & \bullet_k \end{smallmatrix}], \qquad (222)$$

$$f_3(i|j) := g_3(i|j, k) + g_3(j|i, k) = [\begin{smallmatrix} \circled{2} \\ \bullet_i & \bullet_j & \bullet_k \end{smallmatrix}] + [\begin{smallmatrix} \circled{2} \\ \bullet_j & \bullet_i & \bullet_k \end{smallmatrix}]. \qquad (223)$$

The idea behind this notation is that the vertical bar separates bubbles into disconnected pieces. In writing $f_n(\text{var.})$ we mean the symmetrization of $g_n(\text{var.} + \text{other var.})$ in the displayed variables. Of course, $f_n(\text{var.})$ also depends on the other variables (in total there must be $n$ variables), but for simplicity they are suppressed in this notation. Denoting $\mathcal{D}^{ij}$ the operator that produces the contact terms between $i$ and $j$, we find

$$\begin{aligned} \mathcal{D}^{ij} g_3(i, j, k) &= \delta_{ij} (f_3(i+1|j) + f_3(i-1|j)) - (\delta_{i+1,j} + \delta_{i-1,j}) f_3(i|j) \\ &= \delta_{ij} (\Delta_i f_3(i|j)) - (\Delta_i \delta_{ij}) f_3(i|j). \end{aligned} \qquad (224)$$

To summarize the discrete evolution, we have found that

$$\frac{dg_3(i, j, k)}{Jdt} = \left( \sum_{a \in \{i, j, k\}} \Delta_a + \sum_{\substack{a, b \in \{i, j, k\} \\ \text{unordered pairs}}} \mathcal{D}^{a,b} \right) g_3(i, j, k). \qquad (225)$$

This sum over nodes and unordered pairs of nodes can be easily generalized to $n$-point bubbles with an appropriate generalization of $f_3(i|j)$, see (233) to understand the spirit of this equation.

The scaling limit of (225) is taken by introducing continuous variables

$$x = i/L, y = j/L, z = k/L,$$

while sending $J = L^2$ to infinity (note that we set $\mathcal{L} = \mathcal{D} = 1$ for simplicity). For simplicity, will denote the 3-point bubbles in the scaling limit by the same name $^{10}$. Expanding in inverse powers of $L$, the Laplacian becomes (the factor $L^2$ on the rhs. is there to cancel with $1/J$)

$$L^2 \Delta_i g_3(i, j, k) \rightarrow \Delta_x g_3(x, y, z) + \mathcal{O}(L^{-2} g_3), \qquad (226)$$

while a Kronecker delta scales as

$$\delta_{ij} \rightarrow L^{-1} \delta(x-y) + \mathcal{O}(L^{-2}). \qquad (227)$$

Therefore the contact term becomes

$$\begin{aligned} L^2 \mathcal{D}^{ij} g_3(i, j, k) &\rightarrow L^{-1} (\delta(x-y) \Delta_x f_3(x|y) - \delta''(x-y) f_3(x|y)) + \mathcal{O}(L^{-2} g_3) \\ &= L^{-1} \partial_x \partial_y (\delta(x-y) f_3(x|y)) + \mathcal{O}(L^{-2} g_3). \end{aligned} \qquad (228)$$

$$\left| \sum_{i=1}^{L} \frac{1}{L} f(x_i) - \int_{0}^{1} f(x) dx \right| \le \frac{1}{24L^2} \max_{x \in [0,1]} f''(x),$$

<footnote>$^{10}$There is a small subtlety in taking the scaling limit. Formally we should contract the discrete equation with a test function (here done for a one-point function), $(h, g_1) = \sum_{i=1}^L h(i)g_1(i)$, and then approximate the sum as an integral, $(h, g_1) \rightarrow \int_0^1 dx \tilde{h}(x) \tilde{g}_1(x)$ (we denote the continuous function by a tilde in this footnote). However, the approximation leads to an error that also scales with some power of $L$. It can be quantified to scale as $\mathcal{O}(L^{-2})$ using the "mid-point rule",</footnote>

<footnote>where $x_i = (i - 1/2)/L$. To have an error of at most $\mathcal{O}(L^{-2}\tilde{g}_3)$ in (226) and (228), one should therefore define the continuous 3-point function as $\tilde{g}_3((i - 1/2)/L, \dots) = g_3(i, \dots)$. When $L \rightarrow \infty$ this becomes the same as the definition in the main text.</footnote>
The last step is obtained by partial integration w.r.t $x$, which is allowed, since all of these equations should be thought of evaluated against a test function, i.e.

$$ \int dx dy \varphi(x, y) \partial_x \partial_y (\delta(x-y) f_3(x|y)). \qquad (229) $$

Then, we can summarize the continuous evolution equation as

$$ \frac{dg_3(x, y, z)}{dt} = \sum_{a \in \{x, y, z\}} \Delta_a g_3(x, y, z) + L^{-1} \sum_{\substack{a, b \in \{x, y, z\} \\ \text{unordered pairs}}} \partial_x \partial_y (\delta(a-b) f_3(a|b)). \qquad (230) $$

We now discuss the scaling of the 3-point bubble in inverse powers of $L$ that, a priori, could have an expansion according to

$$ g_3(x, y, z) = g_3^{(0)} + L^{-1} g_3^{(1)} + L^{-2} g_3^{(2)} + \dots $$

If we assume that the initial density matrix $\rho_0$ does not involve long range correlations, then all n-point bubbles are originally zero, except for the average fermion density $g_1(x)$. At a later time, $g_3(x, y, z)$ can only be non-zero if the source term $L^{-1}f_3(x|y)$ is non-zero. But $f_3(x|y)$ also satisfies a diffusion equation with new source terms<fnref n="11" /> of which one is $L^{-1}g_3(x|y|z)$. At 0th order in $1/L$, $g_3(x|y|z)$ satisfies a pure diffusion equation with non-zero initial condition, because at time zero (and therefore at all times) it factorizes, $g_3(x|y|z) = g_1(x)g_1(y)g_1(z)$, and $g_1(x)$ is non-zero as argued above. In the original equation for $g_3(x, y, z)$, this gives rise to a non-zero source term only at order $L^{-2}$ (consecutively "descending" the source terms we encountered two factors of $L^{-1}$) and the lower orders are actually zero, $g_3^{(0)} = g_3^{(1)} = 0$. More general, the leading order of an n-bubble is $n-1$ and we denote it by the superscript # in the following. This argument also shows that at leading order the expectation value of products of bubbles factorizes. Rewriting (230) at leading order in terms of diagrams gives rise to (144) in the main text.

## G.2 n-point bubbles

We denote the discrete indices by $l_1, ..., l_n \in \{1, ..., L\}$ and the continuous variables by $x_1, ..., x_n \in [0, 1]$. The generalization of (225) is

$$ \frac{dg_n(l_1, ..., l_n)}{Jdt} = \left( \sum_{a=1}^{L} \Delta_{l_a} + \sum_{\substack{a, b=1 \\ a < b}}^{L} \mathcal{D}^{a,b} \right) g_n(l_1, ..., l_n), \qquad (231) $$

where

$$ \mathcal{D}^{a,b} = \delta_{l_a, l_b} (\Delta_{l_a} f_n(l_a|l_b)) - (\Delta_{l_a} \delta_{l_a, l_b}) f_n(l_a|l_b), \qquad (232) $$

$$ \frac{d}{dt} g_3(x|y,z) = (\Delta_x + \Delta_y + \Delta_z)g_3(x|y,z) + L^{-1}2\partial_y \partial_z (\delta(y-z)g_3(x|y|z)) \\ + L^{-1}(\partial_x \partial_y (\delta(x-y)\bullet) + \partial_x \partial_z (\delta(x-z)\bullet))(f_3(x,y,z)). $$

Here we made full use of the notation introduced in (222): $g_3(x|y|z)$ is a product of three one-bubbles and $f_3(x, y, z) = g_3(x, y, z) + g_3(x, z, y)$ is the symmetrization of a single three-bubble in all variables.

<footnote n="11">¹¹Explicitly, one finds</footnote>
and for $a < b$ we generalize

$$f_n(l_a|l_b) = g_n(l_a, l_{a+1}, \dots, l_{b-1}|l_b, l_{b+1}, \dots, l_{a-1}) + g_n(l_b, l_{a+1}, \dots, l_{b-1}|l_a, l_{b+1}, \dots, l_{a-1}) \quad (233)$$

$$= \left[ \begin{tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}] \node at (0,0) {$l_{b-1}$}; \node at (0,1) {$l_a$}; \node at (1,0) {$l_{a+1}$}; \node at (1,1) {$l_b$}; \node at (2,0) {$l_{a-1}$}; \node at (2,1) {$l_{b+1}$}; \draw[->] (0,0) -- (1,0); \draw[->] (0,1) -- (1,1); \draw[->] (0,0) -- (2,0); \draw[->] (0,1) -- (2,1); \node at (0.5,0) {$b-a$}; \node at (0.5,1) {$n-b+a$}; \node at (1.5,0) {$l_{a-1}$}; \node at (1.5,1) {$l_{b+1}$}; \end{tikzpicture} \right] + \left[ \begin{tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}] \node at (0,0) {$l_{b-1}$}; \node at (0,1) {$l_a$}; \node at (1,0) {$l_{a+1}$}; \node at (1,1) {$l_b$}; \node at (2,0) {$l_{a-1}$}; \node at (2,1) {$l_{b+1}$}; \draw[->] (0,0) -- (1,0); \draw[->] (0,1) -- (1,1); \draw[->] (0,0) -- (2,0); \draw[->] (0,1) -- (2,1); \node at (0.5,0) {$b-a$}; \node at (0.5,1) {$n-b+a$}; \node at (1.5,0) {$l_{a-1}$}; \node at (1.5,1) {$l_{b+1}$}; \end{tikzpicture} \right]. \quad (234)$$

In terms of diagrams the contact term $f_n(l_a|l_b)$ has an easy visualisation: The original n-bubbles is squeezed together such that the nodes $l_a$ and $l_b$ touch each other and the diagram forms an eight. Then the eight is split apart into two disconnected bubbles and we sum over the two possible ways how one can attribute $l_a$ and $l_b$ to the two nodes at the splitting junction.

Then it is easy to see, that the generalization of (230) is just

$$\frac{dg_n(x_1, \dots, x_n)}{dt} = \sum_{a=1}^{L} \Delta_{x_a} g_n(x_1, \dots, x_n) + L^{-1} \sum_{\substack{a, b=1 \\ a < b}}^{L} \partial_{x_a} \partial_{x_b} (\delta(x_a - x_b) f_n(x_a|x_b)). \quad (235)$$

In particular, at leading order in inverse powers of $L$, we recover (143) in the main text.

## References

[1] A. Johnson, S. S. Szigeti, M. Schemmer and I. Bouchoule, Long-lived nonthermal states realized by atom losses in one-dimensional quasicondensates, Phys. Rev. A 96, 013623 (2017), doi:10.1103/PhysRevA.96.013623.

[2] Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan and B. L. Lev, Thermalization near integrability in a dipolar quantum Newton's cradle, Phys. Rev. X 8, 021030 (2018), doi:10.1103/PhysRevX.8.021030.

[3] D. Bernard, T. Jin and O. Shpielberg, Transport in quantum chains under strong monitoring, Europhys. Lett. 121, 60006 (2018), doi:10.1209/0295-5075/121/60006.

[4] B. Skinner, J. Ruhman and A. Nahum, Measurement-induced phase transitions in the dynamics of entanglement, Phys. Rev. X 9, 031009 (2019), doi:10.1103/PhysRevX.9.031009.

[5] V. Gorini, A. Kossakowski and E. C. G. Sudarshan, Completely positive dynamical semi-groups of N-level systems, J. Math. Phys. 17, 821 (1976), doi:10.1063/1.522979.

[6] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys. 48, 119 (1976), doi:10.1007/BF01608499.

[7] H.-P. Breuer and F. Petruccione, The theory of open quantum systems, Oxford University Press, Oxford, UK, ISBN 9780199213900 (2006), doi:10.1093/acprof:oso/9780199213900.001.0001.

[8] T. Prosen, Third quantization: A general method to solve master equations for quadratic open Fermi systems, New J. Phys. 10, 043026 (2008), doi:10.1088/1367-2630/10/4/043026.

[9] B. Horstmann, J. I. Cirac and G. Giedke, Noise-driven dynamics and phase transitions in fermionic systems, Phys. Rev. A 87, 012108 (2013), doi:10.1103/PhysRevA.87.012108.
[10] E. Vernier, *Mixing times and cutoffs in open quadratic fermionic systems*, SciPost Phys. **9**, 049 (2020), doi:10.21468/SciPostPhys.9.4.049.

[11] V. Alba and F. Carollo, *Spreading of correlations in Markovian open quantum systems*, Phys. Rev. B **103**, L020302 (2021), doi:10.1103/PhysRevB.103.L020302.

[12] M. V. Medvedyeva, F. H. L. Essler and T. Prosen, *Exact Bethe ansatz spectrum of a tight-binding chain with dephasing noise*, Phys. Rev. Lett. **117**, 137202 (2016), doi:10.1103/PhysRevLett.117.137202.

[13] D. A. Rowlands and A. Lamacraft, *Noisy coupled qubits: Operator spread-ing and the Fredrickson-Andersen model*, Phys. Rev. B **98**, 195125 (2018), doi:10.1103/PhysRevB.98.195125.

[14] N. Shibata and H. Katsura, *Dissipative quantum Ising chain as a non-Hermitian Ashkin-Teller model*, Phys. Rev. B **99**, 224432 (2019), doi:10.1103/PhysRevB.99.224432.

[15] N. Shibata and H. Katsura, *Dissipative spin chain as a non-Hermitian Kitaev ladder*, Phys. Rev. B **99**, 174303 (2019), doi:10.1103/PhysRevB.99.174303.

[16] A. A. Ziolkowska and F. H. Essler, *Yang-Baxter integrable Lindblad equations*, SciPost Phys. **8**, 044 (2020), doi:10.21468/SciPostPhys.8.3.044.

[17] M. Nakagawa, N. Kawakami and M. Ueda, *Exact Liouvillian spectrum of a one-dimensional dissipative Hubbard model*, Phys. Rev. Lett. **126**, 110404 (2021), doi:10.1103/PhysRevLett.126.110404.

[18] B. Buca, C. Booker, M. Medenjak and D. Jaksch, *Dissipative Bethe ansatz: Exact solutions of quantum many-body dynamics under loss*, arXiv:2004.05955.

[19] P. Ribeiro and T. Prosen, *Integrable quantum dynamics of open collective spin models*, Phys. Rev. Lett. **122**, 010401 (2019), doi:10.1103/PhysRevLett.122.010401.

[20] S. Lerma-Hernández, A. Rubio-García and J. Dukelsky, *Trigonometric SU(N) Richardson-Gaudin models and dissipative multi-level atomic systems*, J. Phys. A: Math. Theor. **53**, 395302 (2020), doi:10.1088/1751-8121/abab54.

[21] D. Yuan, H.-R. Wang, Z. Wang and D.-L. Deng, *Solving the Liouvillian gap with artificial neural networks*, Phys. Rev. Lett. **126**, 160401 (2021), doi:10.1103/PhysRevLett.126.160401.

[22] F. H. L. Essler and L. Piroli, *Integrability of one-dimensional Lindbladi-ans from operator-space fragmentation*, Phys. Rev. E **102**, 062210 (2020), doi:10.1103/physreve.102.062210.

[23] J. Robertson and F. H. L. Essler, *Exact solution of a quantum asymmetric exclusion process with particle creation and annihilation*, J. Stat. Mech. 20102 (2020), doi:10.1088/1742-5468/ac22f8.

[24] M. de Leeuw, C. Paletta and B. Pozsgay, *Constructing integrable Lindblad superoperators*, Phys. Rev. Lett. **126**, 240403 (2021), doi:10.1103/physrevlett.126.240403.

[25] W. Brown and O. Fawzi, *Decoupling with random quantum circuits*, Commun. Math. Phys. **340**, 867-900 (2015), doi:10.1007/s00220-015-2470-1.
[26] F. G. Brandao, A. W. Harrow and M. Horodecki, *Local random quantum circuits are approximate polynomial-designs*, Commun. Math. Phys. **346**, 397-434 (2016), doi:10.1007/s00220-016-2706-8.

[27] A. Nahum, J. Ruhman, S. Vijay and J. Haah, *Quantum entanglement growth under random unitary dynamics*, Phys. Rev. X **7**, 031016 (2017), doi:10.1103/physrevx.7.031016.

[28] A. Nahum, S. Vijay and J. Haah, *Operator spreading in random unitary circuits*, Phys. Rev. X **8**, 021014 (2018), doi:10.1103/physrevx.8.021014.

[29] A. Chan, A. De Luca and J. T. Chalker, *Solution of a minimal model for many-body quantum chaos*, Phys. Rev. X **8**, 041019 (2018), doi:10.1103/physrevx.8.041019.

[30] M. Bauer, D. Bernard and T. Jin, *Stochastic dissipative quantum spin chains (I): Quantum fluctuating discrete hydrodynamics*, SciPost Phys. **3**, 033 (2017), doi:10.21468/sciPostPhys.3.5.033.

[31] H. Spohn, *Nonlinear fluctuating hydrodynamics for anharmonic chains*, J. Stat. Phys. **154**, 1191 (2014), doi:10.1007/s10955-014-0933-y.

[32] H. Spohn, *Large scale dynamics of interacting particles*, Springer Berlin Heidelberg, Berlin, Heidelberg, ISBN 9783642843716 (1991), doi:10.1007/978-3-642-84371-6_1.

[33] C. Kipnis, S. Olla and S. R. S. Varadhan, *Hydrodynamics and large deviation for simple exclusion processes*, Comm. Pure Appl. Math. **42**, 115 (1989), doi:10.1002/cpa.3160420202.

[34] F. Spitzer, *Interaction of Markov processes*, Adv. Math. **5**, 246 (1970), doi:10.1016/0001-8708(70)90034-4.

[35] T. M. Liggett, *Interacting particle systems*, Springer, New York, ISBN 9783540269625 (1985), doi:10.1007/b138374.

[36] B. Derrida, *An exactly soluble non-equilibrium system: The asymmetric simple exclusion process*, Phys. Rep. **301**, 65 (1998), doi:10.1016/S0370-1573(98)00006-4.

[37] G. M. Schütz, *Exactly solvable models for many-body systems far from equilibrium*, Phase transitions and critical phenomena **19**, 1 (2001), doi:10.1016/S1062-7901(01)80015-X.

[38] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, *Current fluctuations in stochastic lattice gases*, Phys. Rev. Lett. **94**, 030601 (2005), doi:10.1103/PhysRevLett.94.030601.

[39] L. Bertini, A. D. Sole, D. Gabrielli, G. Jona-Lasinio and C. Landim, *Macroscopic fluctuation theory*, Rev. Mod. Phys. **87**, 593 (2015), doi:10.1103/RevModPhys.87.593.

[40] D. Bernard, *Can the macroscopic fluctuation theory be quantized?*, To appear.

[41] M. Bauer, D. Bernard and T. Jin, *Equilibrium fluctuations in maximally noisy extended quantum systems*, SciPost Phys. **6**, 045 (2019), doi:10.21468/sciPostPhys.6.4.045.

[42] D. Bernard and T. Jin, *Open quantum symmetric simple exclusion process*, Phys. Rev. Lett. **123**, 080601 (2019), doi:10.1103/PhysRevLett.123.080601.

[43] T. Jin, A. Krajenbrink and D. Bernard, *From stochastic spin chains to quantum Kardar-Parisi-Zhang dynamics*, Phys. Rev. Lett. **125**, 040603 (2020), doi:10.1103/PhysRevLett.125.040603.
[44] L.-H. Gwa and H. Spohn, Six-vertex model, roughened surfaces, and an asymmetric spin Hamiltonian, Phys. Rev. Lett. 68, 725 (1992), doi:10.1103/PhysRevLett.68.725.

[45] L.-H. Gwa and H. Spohn, Bethe solution for the dynamical-scaling exponent of the noisy Burgers equation, Phys. Rev. A 46, 844 (1992), doi:10.1103/PhysRevA.46.844.

[46] B. Derrida, M. R. Evans, V. Hakim and V. Pasquier, Exact solution of a 1d asymmetric exclusion model using a matrix formulation, JOJ. Phys. A: Math. Gen. 26, 1493 (1993), doi:10.1088/0305-4470/26/7/011.

[47] G. Schütz and E. Domany, Phase transitions in an exactly soluble one-dimensional exclusion process, J. Stat. Phys. 72, 277 (1993), doi:10.1007/BF01048050.

[48] 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), doi:10.1103/PhysRevE.52.3512.

[49] T. Sasamoto, One-dimensional partially asymmetric simple exclusion process with open boundaries: Orthogonal polynomials approach, J. Phys. A: Math. Gen. 32, 7109 (1999), doi:10.1088/0305-4470/32/41/306.

[50] 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), doi:10.1103/PhysRevLett.95.240601.

[51] R. A. Blythe and M. R. Evans, Nonequilibrium steady states of matrix-product form: A solver's guide, J. Phys. A: Math. Theor. 40, R333 (2007), doi:10.1088/1751-8113/40/46/r01.

[52] J. de Gier and F. H. L. Essler, Large deviation function for the current in the open asymmetric simple exclusion process, Phys. Rev. Lett. 107, 010602 (2011), doi:10.1103/PhysRevLett.107.010602.

[53] A. Lazarescu and K. Mallick, An exact formula for the statistics of the current in the TASEP with open boundaries, J. Phys. A: Math. Theor. 44, 315001 (2011), doi:10.1088/1751-8113/44/31/315001.

[54] K. Mallick, Some exact results for the exclusion process, J. Stat. Mech. P01024 (2011), doi:10.1088/1742-5468/2011/01/P01024.

[55] N. Crampé, E. Ragoucy and D. Simon, Matrix coordinate Bethe ansatz: Applications to XXZ and ASEP models, J. Phys. A: Math. Theor. 44, 405003 (2011), doi:10.1088/1751-8113/44/40/405003.

[56] D. Bernard and T. Jin, Solution to the quantum symmetric simple exclusion process: The continuous case, Commun. Math. Phys. 384, 1141 (2021), doi:10.1007/s00220-021-04087-x.

[57] R. Frassek, C. Giardinà and J. Kurchan, Duality in quantum transport models, SciPost Phys. 10, 135 (2021), doi:10.21468/SciPostPhys.10.6.135.

[58] V. E. Korepin, N. M. Bogoliubov and A. G. Izergin, Quantum inverse scattering method and correlation functions, Cambridge University Press, Cambridge, UK, ISBN 9780521373203 (1993), doi:10.1017/CBO9780511628832.
[59] N. Andrei and H. Johannesson, *Heisenberg chain with impurities (an integrable model)*, Phys. Lett. A **100**, 108 (1984), doi:10.1016/0375-9601(84)90675-3.

[60] K.-J.-B. Lee and P. Schlottmann, *Integrable spin-1 Heisenberg chain with impurity*, Phys. Rev. B **37**, 379 (1988), doi:10.1103/PhysRevB.37.379.

[61] H. J. de Vega and F. Woynarovich, *New integrable quantum chains combining different kinds of spins*, J. Phys. A: Math. Gen. **25**, 4499 (1992), doi:10.1088/0305-4470/25/17/012.

[62] G. Bedürftig, F. H. L. Eßler and H. Frahm, *Integrable impurity in the supersymmetric $t-J$ model*, Phys. Rev. Lett. **77**, 5098 (1996), doi:10.1103/PhysRevLett.77.5098.

[63] B. Sutherland, *Model for a multicomponent quantum system*, Phys. Rev. B **12**, 3795 (1975), doi:10.1103/PhysRevB.12.3795.

[64] N. Y. Reshetikhin, *Integrable models of quantum one-dimensional magnets with $O(n)$ and $Sp(2k)$ symmetry*, Theor. Math. Phys. **63**, 555 (1985), doi:10.1007/BF01017501.

[65] D. Szász-Schagrin, B. Pozsgay and G. Takacs, *Weak integrability breaking and level spacing distribution*, SciPost Phys. **11**, 037 (2021), doi:10.21468/SciPostPhys.11.2.037.

[66] R. Modak, S. Mukerjee and S. Ramaswamy, *Universal power law in crossover from integrability to quantum chaos*, Phys. Rev. B **90**, 075152 (2014), doi:10.1103/physrevb.90.075152.

[67] M. V. Berry and M. Tabor, *Level clustering in the regular spectrum*, Proc. R. Soc. Lond. A, **356**, 375 (1977), doi:10.1098/rspa.1977.0140.

[68] D. Poilblanc, T. Ziman, J. Bellissard, F. Mila and G. Montambaux, *Poisson vs. GOE statistics in integrable and non-integrable quantum Hamiltonians*, Europhys. Lett. **22**, 537 (1993), doi:10.1209/0295-5075/22/7/010.

[69] K. Kudo and T. Deguchi, *Level statistics of XXZ spin chains with discrete symmetries: Analysis through finite-size effects*, J. Phys. Soc. Jpn. **74**, 1992 (2005), doi:10.1143/jpsj.74.1992.

[70] Y. Y. Atas, E. Bogomolny, O. Giraud and G. Roux, *Distribution of the ratio of consecutive level spacings in random matrix ensembles*, Phys. Rev. Lett. **110**, 084101 (2013), doi:10.1103/physrevlett.110.084101.

[71] D. Bernard, F. H. L. Essler, L. Hruza and M. Medenjak, *Steady fluctuations in $q-assep$*, In Preparation.

[72] A. I. Larkin and Y. N. Ovchinnikov, *Quasiclassical method in the theory of superconductivity*, Sov. J. Exp. Theor. Phys. **28**, 1200 (1969).

[73] J. Maldacena, S. H. Shenker and D. Stanford, *A bound on chaos*, J. High Energy Phys. **08**, 106 (2016), doi:10.1007/jhep08(2016)106.

[74] V. Khemani, D. A. Huse and A. Nahum, *Velocity-dependent Lyapunov exponents in many-body quantum, semiclassical, and classical chaos*, Phys. Rev. B **98**, 144304 (2018), doi:10.1103/physrevb.98.144304.

[75] S. Gopalakrishnan, *Operator growth and eigenstate entanglement in an interacting integrable Floquet system*, Phys. Rev. B **98**, 060302 (2018), doi:10.1103/physrevb.98.060302.
[76] S. Gopalakrishnan, D. A. Huse, V. Khemani and R. Vasseur, *Hydrodynamics of operator spreading and quasiparticle diffusion in interacting integrable systems*, Phys. Rev. B **98**, 220303 (2018), doi:10.1103/physrevb.98.220303.

[77] C. W. von Keyserlingk, T. Rakovszky, F. Pollmann and S. L. Sondhi, *Operator hydrodynamics*, *OTOCs, and entanglement growth in systems without conservation laws*, Phys. Rev. X **8**, 021013 (2018), doi:10.1103/physrevx.8.021013.

[78] J. De Nardis, D. Bernard and B. Doyon, *Hydrodynamic diffusion in integrable systems*, Phys. Rev. Lett. **121**, 160603 (2018), doi:10.1103/physrevlett.121.160603.

[79] P. Zanardi, *Entanglement of quantum evolutions*, Phys. Rev. A **63**, 040304 (2001), doi:10.1103/PhysRevA.63.040304.

[80] T. Prosen and I. Pižorn, *Operator space entanglement entropy in a transverse Ising chain*, Phys. Rev. A **76**, 032316 (2007), doi:10.1103/PhysRevA.76.032316.

[81] M. Žnidarič, T. Prosen and I. Pižorn, *Complexity of thermal states in quantum spin chains*, Phys. Rev. A **78**, 022103 (2008), doi:10.1103/PhysRevA.78.022103.

[82] I. Pižorn and T. Prosen, *Operator space entanglement entropy in XY spin chains*, Phys. Rev. B **79**, 184416 (2009), doi:10.1103/physrevb.79.184416.

[83] M. J. Hartmann, J. Prior, S. R. Clark and M. B. Plenio, *Density matrix renormalization group in the Heisenberg picture*, Phys. Rev. Lett. **102**, 057202 (2009), doi:10.1103/physrevlett.102.057202.

[84] D. Muth, R. G. Unanyan and M. Fleischhauer, *Dynamical simulation of integrable and nonintegrable models in the Heisenberg picture*, Phys. Rev. Lett. **106**, 077202 (2011), doi:10.1103/physrevlett.106.077202.

[85] J. Dubail, *Entanglement scaling of operators: A conformal field theory approach, with a glimpse of simulability of long-time dynamics in 1+1d*, J. Phys. A: Math. Theor. **50**, 234001 (2017), doi:10.1088/1751-8121/aa6f38.

[86] T. Zhou and D. J. Luitz, *Operator entanglement entropy of the time evolution operator in chaotic systems*, Phys. Rev. B **95**, 094206 (2017), doi:10.1103/physrevb.95.094206.

[87] C. Jonay, D. A. Huse and A. Nahum, *Coarse-grained dynamics of operator and state entanglement*, arXiv:1803.00089.

[88] S. Xu and B. Swingle, *Accessing scrambling using matrix product operators*, Nat. Phys. **16**, 199 (2019), doi:10.1038/s41567-019-0712-4.

[89] E. van Nieuwenburg and O. Zilberberg, *Entanglement spectrum of mixed states*, Phys. Rev. A **98**, 012327 (2018), doi:10.1103/physreva.98.012327.

[90] R. Pal and A. Lakshminarayan, *Entangling power of time-evolution operators in integrable and nonintegrable many-body systems*, Phys. Rev. B **98**, 174304 (2018), doi:10.1103/physrevb.98.174304.

[91] T. Takayanagi, *Holographic spacetimes as quantum circuits of path-integrations*, J. High Energy Phys. **12**, 048 (2018), doi:10.1007/jhep12(2018)048.

[92] L. Nie, M. Nozaki, S. Ryu and M. Tian Tan, *Signature of quantum chaos in operator entanglement in 2d CFTs*, J. Stat. Mech. 093107 (2019), doi:10.1088/1742-5468/ab3a29.
[93] M. Žnidarič, *Entanglement growth in diffusive systems*, Commun. Phys. **3**, 100 (2020), doi:10.1038/s42005-020-0366-7.

[94] V. Alba, *Diffusion and operator entanglement spreading*, Phys. Rev. B **104**, 094410 (2021), doi:10.1103/PhysRevB.104.094410.

[95] T. Prosen and M. Žnidarič, *Is the efficiency of classical simulations of quantum dynamics related to integrability?*, Phys. Rev. E **75**, 015202 (2007), doi:10.1103/PhysRevE.75.015202.

[96] V. Alba, J. Dubail and M. Medenjak, *Operator entanglement in interacting integrable quantum systems: The case of the rule 54 chain*, Phys. Rev. Lett. **122**, 250603 (2019), doi:10.1103/physrevlett.122.250603.

[97] M. Žnidarič, *Spin transport in a one-dimensional anisotropic Heisenberg model*, Phys. Rev. Lett. **106**, 220601 (2011), doi:10.1103/physrevlett.106.220601.

[98] M. Žnidarič, *Transport in a one-dimensional isotropic Heisenberg model at high temperature*, J. Stat. Mech. P12008 (2011), doi:10.1088/1742-5468/2011/12/p12008.

[99] M. Žnidarič, A. Scardicchio and V. Kerala Varma, *Diffusive and subdiffusive spin transport in the ergodic phase of a many-body localizable system*, Phys. Rev. Lett. **117**, 040601 (2016), doi:10.1103/physrevlett.117.040601.