# Coherent Fluctuations in Noisy Mesoscopic Systems, the Open Quantum SSEP, and Free Probability

Ludwig Hruza* and Denis Bernard●,†

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

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

(Received 5 May 2022; revised 5 December 2022; accepted 17 January 2023; published 24 March 2023)

Quantum coherences characterize the ability of particles to quantum mechanically interfere within some given distances. In the context of noisy many-body quantum systems, these coherences can fluctuate. A simple toy model to study such fluctuations in an out-of-equilibrium setting is the open quantum symmetric simple exclusion process (Q-SSEP), which describes spinless fermions in one dimension hopping to neighboring sites with random amplitudes coupled between two reservoirs. Here, we show that the dynamics of fluctuations of coherences in Q-SSEP have a natural interpretation as free cumulants, a concept from free probability theory. Based on this insight, we provide heuristic arguments as to why we expect free probability theory to be an appropriate framework to describe coherent fluctuations in generic mesoscopic systems where the noise emerges from a coarse-grained description. In the case of Q-SSEP, we show how the link to free probability theory can be used to derive the time evolution of connected fluctuations of coherences as well as a simple steady-state solution.

DOI: 10.1103/PhysRevX.13.011045

Subject Areas: Condensed Matter Physics,
Quantum Physics, Statistical Physics

## I. INTRODUCTION

From the point of view of statistical physics, a system is completely characterized if the probability of all possible microscopic configurations is known. From this knowledge, one can obtain the probability distribution of macroscopic (or thermodynamic) variables such as the density or the current profile, which are, in principle, experimentally observable. In an equilibrium situation, where the probability distribution on microscopic configurations has been well known since Boltzmann, macroscopic variables satisfy a large deviation principle: They are strongly peaked at their mean value, around which their probability distribution decays exponentially with the number of particles in the system, proportional to the appropriate thermodynamic potential (e.g., the entropy or the free energy); see, e.g., Ref. [1].

In contrast to this, out-of-equilibrium situations might depend on a large variety of system-dependent details, and there are no general formulas for the probability distribution on microscopic configurations. Nevertheless, over the course of the last 30 years, a lot of progress has been made

to understand the statistics of macroscopic variables in out-of-equilibrium systems, which show rich and new features. For example, density correlations in nonequilibrium steady states extend to macroscopic distances, as has been experimentally observed by Dorfmann [2]. Furthermore, fluctuation relations that go beyond the linear response regime have been shown to be generically applicable to out-of-equilibrium systems [3,4]. In addition, the large deviation principle introduced above in the equilibrium context has been realized to hold also in the nonequilibrium setting—with the need to formulate appropriate out-of-equilibrium thermodynamic potentials; see, e.g., Ref. [5].

For classical systems, these efforts have culminated in the formulation of the so-called macroscopic fluctuation theory (MFT) [6,7], which applies to systems with diffusive transport. Usually the system is maintained out of equilibrium by coupling the boundaries to reservoirs at different chemical potentials. MFT allows one to specify the probability distribution of the density and current profile in such systems, and remarkably, this relies on only two system-dependent quantities—the diffusion constant and the mobility. The development of MFT has been strongly inspired by the study of microscopic toy models, more precisely, stochastic lattice gases, which revealed universal properties in the density and current fluctuations in the sense that they did not depend on the precise underlying microscopic dynamics. A very important role in this context has been played by the symmetric and asymmetric simple exclusion processes (SSEP and ASEP) since these toy

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>* ludwig.hruza@ens.fr</footnote>

<footnote>† denis.bernard@ens.fr</footnote>
models are exactly solvable [8–11] (for a review, see Refs. [5,12]).

The major question we are concerned with is whether MFT can be extended to a quantum setting [13]. Such a theory, which could be called the quantum mesoscopic fluctuation theory, should not only describe the statistical properties of the diffusive transport (i.e., density and current profiles) in out-of-equilibrium quantum systems but also those of quantum coherent effects such as interference or entanglement. These effects are inscribed into the coherences $G_{ij} = \text{Tr}(\rho c_i^\dagger c_j)$, the off-diagonal elements of the density matrix $\rho$ when expressed in the particle number basis on each site $i, j = 1, \dots, L$. A big difference from classical noisy systems (such as SSEP) is that the density matrix in a noisy quantum system—usually interpreted as the probability distribution of quantum states —is itself a dynamical variable that fluctuates. A quantity sensible to these fluctuations is the entanglement entropy since it is composed of partial traces of powers of $\rho$. Fluctuations of the entanglement entropy have already been studied in the context of unitary random circuits [14,15] (see Refs. [16,17] for a review), as well as certain toy models of driven diffusive systems [18]. In seeking a mesoscopic fluctuation theory, we hope to establish a general framework to deal with coherent effects in diffusive systems, based on coherences $G_{ij}$.

The picture (see Fig. 1) we have in mind is that such a quantum extension of MFT should apply to mesoscopic systems, i.e., to systems in which the system size $L$ is of the order of the coherence length $L_\phi$ of the dynamical degrees of freedom (e.g., for electrons in a disordered metal of size $L_\phi \approx 1$ µm [19]). Another important length scale is the mean free path $\ell$, below which transport is ballistic and the microscopic degrees of freedom change very rapidly with time. After an average over the microscopic degrees of freedom inside ballistic cells of the size of the mean free paths $\ell$, we expect that the system has a coarse-grained description in terms of a stochastic and unitary process. These ballistic cells should be assumed to be thermodynamically large but small compared to the system size, $1 \ll \ell \ll L$ (we identify the physical length with the number of lattice sites by setting the lattice constant to 1, $a_{\text{uv}} = 1$). At lengths above $\ell$, we expect diffusive transport, and at

<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic representation of a mesoscopic system coupled to two reservoirs of different chemical potentials $\mu_a$ and $\mu_b$. Here, $\ell$ denotes the ballistic length, above which the transport is diffusive, and $L_\phi$ is the coherence length, below which interference effects can be observed.</figcaption></figure>

lengths below $L_\phi$, we observe quantum mechanical interference. This is the mesoscopic regime we are interested in.

The decomposition of a many-body quantum system into thermodynamically large cells finds many examples in the literature. Notably, it appears as the main feature of generalized hydrodynamics (GHD), which has allowed one to study out-of-equilibrium dynamics in integrable systems [20,21] (see Refs. [22,23] for a review). Since this theory treats the system as a composition of many thermodynamically large fluid cells that are locally in equilibrium with respect to a generalized Gibbs ensemble, it loses the statistical properties of quantum correlations between the fluid cells. Note that an approach to restore these quantum correlations in GHD has recently been proposed in Ref. [24]. Furthermore, GHD is mostly concerned with ballistic transport due to the integrable nature of the systems to which it applies. Therefore, it does not seem to be a good candidate for a quantum version of MFT.

A toy model that allows one to study fluctuations in the mesoscopic regime is the quantum symmetric simple exclusion process (Q-SSEP) [25–30]. It describes the noisy and coherent hopping of spinless fermions on a one-dimensional and discrete lattice chain, which is maintained out of equilibrium by two particle reservoirs at different chemical potentials. While a simple average over the noise completely destroys coherent effects and reduces the model to the classical SSEP, the fluctuations of coherent effects survive in the nonequilibrium steady state [27], although they are subleading in the system size. In this sense, the Q-SSEP describes mesoscopic transport. Furthermore, fluctuations of coherences satisfy a large deviation principle, in analogy to macroscopic variables in MFT. This can be seen as the first evidence that the fluctuations in Q-SSEP show universal properties that might apply to a larger class of mesoscopic systems. From this point of view, a single lattice site in Q-SSEP would correspond to a ballistic cell of the coarse-grained mesoscopic system.

In this paper, we show that the mathematical structure of Q-SSEP can be described within the framework of free probability theory, an extension of classical probability theory to noncommutative random variables pioneered by Voiculescu in the 1990s [31]. More precisely, connected fluctuations of coherences in Q-SSEP correspond to free cumulants, which generalize the notion of cumulants from classical probability. Free cumulants have a combinatorial nature and can be obtained from the moments of a random variable as a sum over noncrossing partitions. This is a structure of which we make repeated use in this paper: First, it allows us to derive the time evolution of connected fluctuations in Q-SSEP in a diffusive scaling limit (for the closed Q-SSEP, i.e., without boundaries, this scaling limit has been discussed in Ref. [30]). Secondly, it also allows us to find a simple derivation for the steady-state solution of the fluctuations of coherences. The connection between free probability and the steady-state solution of Q-SSEP
has already been observed by the Biane [32]. Here, we generalize his result to all times $0 < t < \infty$ and give a more physical explanation for why free probability arises in our context. In particular, we are able to extend this argument to generic mesoscopic systems, which implies that free probability could be the natural mathematical framework to characterize fluctuations in such systems.

This is the main point of the paper: Tools from free probability might play a significant role in understanding fluctuating many-body quantum systems in the mesoscopic regime, both in and out of equilibrium.

Indeed, free probability theory finds an explicit realization of its concepts through large random matrices [33–35]. For a generic mesoscopic system, the coarse-grained view of Fig. 1 in terms of ballistic cells, which undergo a very rapid (unitary) evolution—much faster than the system’s evolution as a whole—can be modeled by large random matrices whose size scales with $\ell$ and whose distribution is invariant under the rapid evolution of the ballistic cell, i.e., by the same unitary group (for the moment, we leave open the question of which degrees of freedom inside the ballistic cells are modeled precisely). Insofar, it is not surprising that fluctuations on mesoscopic scales, i.e., between ballistic cells, show signs of free probability. We believe that their origin simply lies in a coarse-grained description of microscopic spatial scales and unitary invariance at these scales. Further evidence for this interpretation is provided by the recent observation of free probability in the context of ETH, where instead of spatial scales the coarse graining is over small energy scales [36] (submitted simultaneously with our work).

The structure of the paper is as follows. Section II provides heuristic arguments that show why the fluctuations of spatial coherences in generic mesoscopic systems could be appropriately described within free probability theory. In Sec. III, we give a small introduction to the subject of free probability and, in particular, explain the relevance of crossing and noncrossing partitions. Section IV deals with the toy model Q-SSEP. In Sec. IV B, we recall the main properties of Q-SSEP. Sections IV C and IV D are concerned with the formulation of the dynamical equations of the correlation functions of the open Q-SSEP in the diffusive scaling limit. Finally, Sec. IV E explains how to identify free cumulants within Q-SSEP, and Sec. IV F explains how the relation to free cumulants can be used to find a simple steady-state solution. Details and proofs are given in a few appendixes. In particular, Appendixes E 1 and E 2 constitute the major derivations of the free probability structure in Q-SSEP.

## II. GENERAL PICTURE

We expect that there is a very general relationship between coherent fluctuations in diffusive one-dimensional

fermionic systems in the mesoscopic regime and the mathematical framework of free probability. The aim of this section is to convey an intuition for this claim, not to give complete proofs. We stress that the actual chronological order of our work is the opposite. First, it is observed that the mathematical structure of the toy model Q-SSEP has a relation to free probability. Then, we realize that one can reduce this relation to three simple conditions that apply to generic mesoscopic systems.

Before stating these conditions, some explanation is necessary. In the spirit of Fig. 1, we need to assume that the system exhibits a separation of timescales responsible for fast ballistic transport on scales of the mean free path $\ell$ (inside ballistic cells) versus slow diffusive transport on scales of the system size $L$. In a coarse-grained description of the slow degrees of freedom, the fast ones will act as a source of noise. Therefore, mesoscopic observables become random variables with a (yet to be defined) expectation value $\mathbb{E}_t$. If such observables do not depend on long-range coherences, such as the local particle number $\hat{n}_i := c_i^\dagger c_i$ ($c_i^\dagger$ is a fermionic creation operator on site $i = 1, \dots, L$) or the associated particle current, then we assume that they are well described by classical MFT, with the replacement $\mathbb{E}_t[\text{Tr}(\rho^*)] = \langle \cdot \rangle_M^{\text{MFT}}$. In contrast, we are interested in the statistical properties of coherences,

$$G_{ij}(t) = \text{Tr}(\rho_t c_i^\dagger c_j),$$

purely quantum mechanical properties without a classical analogue—and, therefore, outside the scope of classical MFT. In Appendix A, we outline a procedure for how one could, in principle, experimentally measure $G_{ij}$, following an idea in Ref. [18].

The three conditions sufficient for the link with free probability concern the statistical properties of coherences with respect to the noise expectation value $\mathbb{E}_t$ in the limit $1 \ll \ell \ll L$. Here, we adopt the view that $G_{ij}$ is a random variable with a time-dependent probability distribution, hence the subscript $t$.

(1) Local $U(1)$ invariance: The expectation value is invariant under $G_{ij} \to e^{-i\theta_i} G_{ij} e^{i\theta_j}$, which is a multiplication with local phases. In other words, unless $\{i_1, \dots, i_n\}$ is a permutation of $\{j_n, \dots, i_n\}$, we have

$$\mathbb{E}_t[G_{i_1 j_1} \dots G_{i_n j_n}] = 0. \qquad (1)$$

(2) Expectation values of “loops” with distinct indices, $i_k \neq i_l$, $\forall k \neq l \in \{1, \dots, n\}$, scale as

$$\mathbb{E}_t[G_{i_1 i_2} G_{i_2 i_3} \dots G_{i_n i_1}] \sim L^{-n+1}. \qquad (2)$$
(3) Expectation values of products of loops factorize as

$$ \mathbb{E}_t[G_{i_1 i_2} \dots G_{i_m i_1} G_{j_1 j_2} \dots G_{j_n j_1}] \\ = \mathbb{E}_t[G_{i_1 i_2} \dots G_{i_m i_1}] \mathbb{E}_t[G_{j_1 j_2} \dots G_{j_n j_1}] \quad (3) $$

at leading order. This holds, in particular, if $i_1 = j_1$.

We speak about loops because we can associate diagrams with the expectation values of the various products of $G_{ij}$ by connecting nodes $i$ and $j$ by a directed edge. The diagram associated with Eq. (2) is indeed a loop,

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

In fact, the local $U(1)$ invariance implies that a diagram built in this way is nonzero only if at each node the number of incoming edges is equal to the number of outgoing edges.

For later use, we note that the factorization in Eq. (3) also implies that connected expectation values scale as

$$ \mathbb{E}_t[G_{i_1 i_2} G_{i_2 i_3} \dots G_{i_n i_1}]^c = \mathcal{O}(L^{-n+1}) \quad (4) $$

even if indices become equal. If, for example, $i_1 = i_k$ and no other indices coincide, then $\mathbb{E}_t[G_{i_1 i_2} \dots G_{i_n i_1}]$ factorizes into $\mathbb{E}_t[G_{i_1 i_2} \dots G_{i_{k-1} i_1}] \mathbb{E}_t[G_{i_{k+1} i_1} \dots G_{i_n i_1}]$ at leading order. This product of loop expectation values scales as $L^{-(k-1)+1} L^{-(n-k+1)-1} = L^{-n+2}$. The connected expectation value is constructed precisely by subtracting this product of expectation values from $\mathbb{E}_t[G_{i_1 i_2} \dots G_{i_n i_1}]$. What remains must then scale at least 1 order of magnitude lower in $L$. This explains Eq. (4). The same argument also shows why the connected expectation value of an arbitrary product of loops can never become more dominant in $L$ than the connected expectation value of single loops of the same length.

As an example, consider loops of two points for which Eq. (2) implies $\mathbb{E}[G_{ij}G_{ji}] \sim 1/L$ if $i \neq j$ while Eq. (3) implies $\mathbb{E}[G_{ii}^2] \sim \mathbb{E}[G_{ii}]^2$ where $\mathbb{E}[G_{ii}] = \mathcal{O}(1)$. Therefore, $\mathbb{E}[G_{ij}G_{ji}]^c := \mathbb{E}[G_{ij}G_{ji}] - \delta_{ij}\mathbb{E}[G_{ii}]^2$ scales at most with $1/L$ for any choice of $i$ and $j$, which nicely illustrates that a loop's expectation values jump by an order of $L$ if indices become equal while its connected part does not.

The following subsection gives a definition of the noise average $\mathbb{E}_t$ and shows how to derive the three conditions from this definition. Section II B shows how the three conditions imply that fluctuations of coherences are in relation with free cumulants from the theory of free probability.

## A. Fluctuations in mesoscopic systems

We consider a one-dimensional system of spinless interacting fermions without noise on $L$ discrete lattice sites described by a density matrix $\rho$, see Fig. 2. In order to talk about transport, the system is required to have a locally conserved charge that can be transported. In the continuous description of the system, this is the density of particles $n(x, t)$; it satisfies a local conservation law, the continuity equation $\partial_t n + \partial_x j = 0$. For the following discussion, it is not important if the system is open or closed.

The separation of timescales between fast and slow degrees of freedom ensures that ballistic cells do not exchange particles with each other during times smaller than the typical timescale $t_\ell$ of ballistic cells. In other words, for $t < t_\ell$, each ballistic cell evolves as a “closed” system, and its dynamics is described by unitary and particle-preserving transformations, $\rho \rightarrow U^{(i)}\rho U^{(i)\dagger}$. Here, $U^{(i)}$ acts only on the Hilbert space $\mathcal{H}^{(i)} \equiv \mathbb{C}^{(2\ell)}$ of the cell around $i$. Particle preservation means that

$$ [U^{(i)}, N^{(i)}] = 0, \quad (5) $$

where $N^{(i)} = \sum_{i' \sim i} \hat{n}_{i'}$ is the total number operator in the cell around site $i$. The sum carries over all sites $i'$ within distance $\ell/2$ from $i$. As a consequence, $U^{(i)}$ decomposes into subspaces $\Lambda_p$ ($p$th exterior product of $\mathbb{C}^\ell$) of dimension $(\ell_p)$ where the number of fermions is fixed to $p$. In other words, $U^{(i)}$ consists of blocks $U_p^{(i)}$ on the diagonal.

To construct the average over ballistic cells, we make the ergodic hypothesis that within a time interval $t_\ell$, the ballistic cell has undergone all possible unitary and particle-preserving transformations. Then, time averages over intervals $t_\ell$ can be replaced by a uniform Haar average $[\dots]_U$ over such unitaries $U$ that only mix degrees of freedom within a cell but not between cells.

The idea that fluctuations of a chaotic quantum-many-body system can be characterized through a definition of ergodicity by unitary invariance has already been put forward [37]. The consequences of restricting such a global unitary invariance to local sectors of fixed energy has been explored in Refs. [38,39] (the latter in the context of ETH, though both discuss very similar topics). Here, we use similar ideas, but instead of local in energy, we restrict the unitary invariance to be local in space, i.e., to unitary invariance within ballistic cells.

If the separation of two sites $i$ and $j$ is larger than $\ell$, then the average over the corresponding ballistic cells is independent, i.e., $U = U^{(i)}U^{(j)}$, where $U^{(i)}$ and $U^{(j)}$ act only on the cells around $i$ and $j$, respectively. Otherwise, if $i$ and $j$ are in the same cell (but not necessarily equal), we assume that they interact via the same Haar-distributed unitary $U^{(i)} = U^{(j)}$. With this convention, we define the expectation value of coherences to be
$$ \mathbb{E}_t[G_{ij}] := \frac{1}{t_\ell} \int_t^{t+t_\ell} G_{ij}(t')dt' = \text{Tr}(\rho_t[c_i^\dagger c_j]_U), \quad (6) $$

where

$$ [c_i^\dagger c_j]_U := \int d\mu(U) U^\dagger c_i^\dagger c_j U \quad (7) $$

denotes the Haar average. Note that by the cyclic property of the trace, we could also choose to average $\rho_t$ instead of $c_i^\dagger c_j$, which makes it clear that the average carries over all possible evolutions of the cell.

The time average over $t_\ell$ promotes $G_{ij}(t)$ to a random variable that is only sensitive to the long (with respect to $t_\ell$) time behavior of the system. Therefore, we also expect that the right-hand side will not depend on the time at which $\rho'_t$ is evaluated, as long as $t' \in [t, t + t_\ell]$. Here, we choose to evaluate at time $t$.

For quadratic fluctuations, this process generalizes to

$$ \begin{aligned} \mathbb{E}_t[G_{ij}G_{kl}] &:= \frac{1}{t_\ell} \int_t^{t+t_\ell} G_{ij}(t')G_{kl}(t')dt' \\ &= \text{Tr}(\rho_t \otimes \rho_t \cdot [c_i^\dagger c_j \otimes c_k^\dagger c_l]_{U\otimes U}) \end{aligned} \quad (8) $$

and similarly for higher-order fluctuations.

### 1. U(1) invariance

The averages appearing in Eqs. (6) and (8) are invariant under unitary transformations of ballistic cells by construction. For example, we are interested in an average such as $[c_i^\dagger \otimes c_i]_{U\otimes U}$, and unitary invariance means

$$ [c_i^\dagger \otimes c_i]_{U \otimes U} = V^{(i)} \otimes V^{(i)} [c_i^\dagger \otimes c_i]_{U \otimes U} V^{(i)\dagger} \otimes V^{(i)\dagger}, \quad (9) $$

where $V^{(i)}$ can be any particle-preserving unitary on the cell around $i$. In particular, we can choose this $V^{(i)}$ to belong to the subgroup of local $U(1)$ transformations on each site, $V = \exp(\sum_i \theta_i \hat{n}_i)$, which are generated by the particle number operator $\hat{n}_i = c_i^\dagger c_i$. The creation and annihilation operators $c_i^\dagger$ and $c_i$ carry the charges +1 and -1 with respect to this $U(1)$ transformation. In other words, $[\hat{n}_i, c_i^\dagger] = c_i^\dagger$ and $[\hat{n}_i, c_i] = -c_i$. Then, it is easy to see that any combination of operators inside an average $[\cdots]_{U\otimes U}$ must be composed of an equal number of $c_i^\dagger$'s and $c_i$'s in order for the total charge to add up to zero—and, as a consequence, to be invariant under the subgroup of local $U(1)$ transformations. But this is exactly the same as what we claim in Eq. (1). For example, Eq. (8) is nonzero only if $i = l$ and $j = k$.

Evaluating averages such as Eq. (8) can be done using Schur's Lemma, but this does not add a new physical insight to the picture. If we restrict ourselves to the one-particle sector, however, one can learn that the local unitary average is equal to a homogeneous sum of $G_{ij}$'s over all

<figure><img src="image_4.png" /><figcaption>FIG. 2. Cells around sites $i$ and $j$ in the mesoscopic system corresponding to single sites in the coarse-grained description. The toy model Q-SSEP seems to be a good candidate for such a coarse-grained description. Noise emerges by averaging the mesoscopic system over all unitary transformations that only act on the individual cells, and it conserves the number of particles inside each cell.</figcaption></figure>

sites in the cells. This case is shown explicitly in the remainder of this subsection. The discussion can be skipped without having an impact on the comprehension of the following.

Unitaries in the one-particle sector are of the form

$$ U^{(i)} = \exp(\mathbf{c}^{(i)\dagger} M^{(i)} \mathbf{c}^{(i)}), \quad (10) $$

where $\mathbf{c}^{(i)} = (c_{i-\ell/2}, \dots, c_{i+\ell/2})$ comprises all fermionic annihilation operators inside the cell around $i$, and $M^{(i)}$ is an anti-Hermitian $\ell \times \ell$ matrix. Using the identity $e^{-c^{\dagger}Mc}c_i e^{c^{\dagger}Mc} = \sum_j (e^M)_{ij}c_j$, where $\mathbf{c} = (c_1, \dots, c_L)$, and the definition $U = e^{c^{\dagger}Mc}$, where the $L \times L$ matrix $M$ is composed of the blocks $M^{(i)}$ and $M^{(j)}$ on the diagonal at positions $i$ and $j$, and zeros otherwise, we obtain

$$ \text{Tr}(\rho_t U^\dagger c_i^\dagger c_j U) = (uG(t)u^\dagger)_{ij}. \quad (11) $$

Here, $u := e^{M^*} = \text{diag}(1, u^{(i)}, 1, u^{(j)}, 1)$ is a block-diagonal matrix with two unitary blocks $u^{(i)} = e^{M^{(i)*}}$ and $u^{(j)} = e^{M^{(j)*}}$ at positions $i$ and $j$ over which the Haar average can be taken [40]. For example, evaluating Eq. (6) using Eq. (11),

$$ \begin{aligned} \mathbb{E}_t[G_{ij}] &= \sum_{\substack{i' \sim i \\ j' \sim j}} [u_{ii'}^{(i)} G_{i'j'} u_{jj'}^{(j)*}] \\ &= \begin{cases} \frac{1}{\ell} \sum_{i'} G_{i'i} & \text{if } j = i \\ 0 & \text{otherwise,} \end{cases} \end{aligned} \quad (12) $$

where the sums carry over all $i'$ and $j'$ in the cells around $i$ and $j$. Note that although expressed through $G_{i'i'}$, the right-hand side becomes nonrandom as a result of the law of large numbers. Evaluating Eq. (8), using appropriate Haar averages, for any choice of $i$ and $j$, up to corrections in $1/\ell$, one finds

$$ \mathbb{E}_t[G_{ij}G_{ji}] = \frac{1}{\ell^2} \sum_{\substack{i' \sim i \\ j' \sim j}} G_{i'j'} G_{j'i'} + \delta_{ij} \left( \frac{1}{\ell} \sum_{\substack{i' \\ i \sim i}} G_{i'i'} \right)^2. \quad (13) $$

### 2. Scaling with system size

The scaling of loop expectation values in Eq. (2) with system size $L$ can be derived from the assumption that
$$ \mathbb{E}_t[G_{ii}^2] = \left(\frac{1}{\ell}\sum_{i' \sim i} G_{i'i'}\right)^2 = \mathbb{E}_t[G_{ii}]^2. \quad (19) $$

In the last equality, we used Eq. (12). This calculation shows that for the given example, the expectation value of products of loops (here, $G_{ii}^2$) is indeed proportional to the product of the expectation values of loops (here, single loops $G_{ii}$) at leading order. Using Schur's Lemma, one can show that this statement remains true at all particle sectors and for higher-order fluctuations, thereby confirming Eq. (3).

## B. Link with free probability

The link with free probability occurs when expressing loop expectation values $\mathbb{E}[G_{i_1 i_2} \dots G_{i_n i_1}]$ in terms of their connected parts. We find that this expansion involves a sum over noncrossing partitions only—see Eq. (31). Such noncrossing partitions are at the heart of a mathematical theory of noncommuting random variables that is called free probability theory; it is introduced in the next section. By means of noncrossing partitions, one can define so-called free cumulants, which will arise naturally in our setting.

As an intermediate step and as usual in random matrix theory, we consider the measure of $G$, the matrix of coherences $G_{ij}$, defined by the expectation value of its traces,

$$ \phi(G^n) := \frac{1}{L} \mathbb{E}[\mathrm{Tr}(G^n)] = \frac{1}{L} \sum_{i_1, \dots, i_n} \mathbb{E}[G_{i_1 i_2} \cdots G_{i_n i_1}]. \quad (20) $$

Here, we dropped the subscript $t$ of the measure $\mathbb{E}_t$ for simplicity. Whenever two indices $i_k$ in the above sum are equal, the expectation value factorizes according to Eq. (3). We therefore split this sum into sums where all indices are distinct—except for a given set of indices that are forced to be equal,

$$ \sum_{i_1, \dots, i_n} = \sum_{i_1, \dots, i_n \text{ distinct}} + \sum_{i_1 = i_2, \dots, i_n \text{ distinct}} + \dots + \sum_{i_1 = i_2 = \dots = i_n} \quad (21) $$

Such a splitting can be understood as a sum over partitions $\pi \in P(n)$ of the set $\{1, \dots, n\}$ into blocks $b \in \pi$ that group together all the indices $i_k$ that are forced to be equal. For example, the partitions corresponding to the three terms above are $\pi = \{\{1\}, \dots, \{n\}\}$, $\pi = \{\{1, 2\}, \{3\}, \dots, \{n\}\}$, and $\pi = \{\{1, 2, \dots, n\}\}$. The sum becomes

$$ \sum_{i_1, \dots, i_n} = \sum_{\pi \in P(n)} \sum_{\substack{i_1 \dots i_n \\ \text{distinct,} \\ \text{except } i_k = i_j \text{ whenever } k,l \\ \text{are in the same block of } \pi}}. \quad (22) $$

For $n=4$, two possible partitions are

$$ \pi_1 = \{\{1,3\}, \{2\}, \{4\}\} = \begin{tikzpicture}[baseline=(current bounding box.center), scale=0.7] 
    \draw (0,0) circle (1); 
    \node at (0,0) {$2$}; 
    \node at (0.5,0) {$1$}; 
    \node at (-0.5,0) {$4$}; 
    \draw (0,0) -- (-0.5,0); 
    \draw (0,0) -- (0.5,0); 
    \draw (0,0) -- (-1,0); 
    \draw (0,0) -- (1,0); 
    \node at (-0.5,0) {$3$}; 
    \node at (0.5,0) {$4$}; 
    \end{tikzpicture} \quad (23) $$

and

$$ \pi_2 = \{\{1,3\}, \{2,4\}\} = \begin{tikzpicture}[baseline=(current bounding box.center), scale=0.7] 
    \draw (0,0) circle (1); 
    \node at (0,0) {$2$}; 
    \node at (0.5,0) {$1$}; 
    \node at (-0.5,0) {$4$}; 
    \draw (0,0) -- (-0.5,0); 
    \draw (0,0) -- (0.5,0); 
    \draw (0,0) -- (-1,0); 
    \draw (0,0) -- (1,0); 
    \node at (-0.5,0) {$3$}; 
    \node at (0.5,0) {$4$}; 
    \end{tikzpicture} \quad (24) $$

Representing partitions by diagrams where nodes in the same block are connected by dashed lines very intuitively shows that $\pi_2$ is a crossing partition while $\pi_1$ is non-crossing. It turns out that in the sum in Eq. (20), terms corresponding to crossing partitions are of the order of $1/L$ and lower and therefore vanish for $L \to \infty$, while all noncrossing partitions are of order one. Here, we illustrate this fact for $n=4$ and the two examples above. The term corresponding to $\pi_1$ factorizes and becomes

$$ \frac{1}{L} \sum_{i_1=i_3, i_2=i_4} \mathbb{E}[G_{i_1 i_2} G_{i_2 i_1}] \mathbb{E}[G_{i_3 i_4} G_{i_4 i_3}] = \mathcal{O}(1). \quad (25) $$

Each of the two-loop expectation values scales as $1/L$, and each term in the sum is of order $1/L^3$. Since the sum carries over three indices running from 1 to $L$, this cancels, and the resulting scaling is of order one.

In contrast, the term corresponding to $\pi_2$ can factorize in two different ways, and it becomes

$$ \frac{1}{L} \sum_{i_1, i_2 \text{ distinct}} 2 \mathbb{E}[G_{i_1 i_2} G_{i_2 i_1}]^2 = \mathcal{O}(1/L). \quad (26) $$

The difference from $\pi_1$ is that now there are only two indices to sum over, and hence the scaling is of order $1/L$.

For the surviving noncrossing partitions, one can ask how the partition that determines which indices $i_k$ are equal is related to the product of loop expectation values of $G_{ij}$'s that appears after the factorization. If we modify Eq. (23) by connecting as many edges by solid lines as possible without crossing a dashed line, i.e.,

$$ \begin{tikzpicture}[baseline=(current bounding box.center), scale=0.7] 
    \draw (0,0) circle (1); 
    \node at (0,0) {$2$}; 
    \node at (0.5,0) {$1$}; 
    \node at (-0.5,0) {$4$}; 
    \draw (0,0) -- (-0.5,0); 
    \draw (0,0) -- (0.5,0); 
    \draw (0,0) -- (-1,0); 
    \draw (0,0) -- (1,0); 
    \node at (-0.5,0) {$3$}; 
    \node at (0.5,0) {$4$}; 
    \end{tikzpicture} \quad (27) $$

we see that, for the partition $\pi_1$, this product of loop expectation values, i.e., $\mathbb{E}[G_{i_1 i_2} G_{i_2 i_1}] \mathbb{E}[G_{i_3 i_4} G_{i_4 i_3}]$, is
determined by the solid lines. In fact, for every noncrossing partition $\pi$ of nodes, the solid lines define a unique noncrossing partition $\pi^*$ of edges, the dual partition (also called Kreweras complement). If we label an edge by the adjacent node with the lower number, we have, for example, $\pi_1^* = \{\{1, 2\}, \{3, 4\}\}$. Each block therein corresponds to a loop of $G_{ij}$'s. In this terminology, the expansion of $\phi(G^n)$ in Eq. (20) becomes

$$ \frac{1}{L} \sum_{\pi \in NC(n)} \sum_{\substack{i_1, \dots, i_n \\ \text{distinct,} \\ \text{except } i_j = i \\ \text{if } k, l \in b \text{ for } b \in \pi}} \prod_{b \in \pi^*} \mathbb{E}[G_{i_{b(1)}} i_{b(2)}} \dots G_{i_{b(|b|)}} i_{b(1)}], \quad (28) $$

where $NC(n)$ denotes all noncrossing partitions of the set $\{1, \dots, n\}$.

Instead of using this slightly awkward sum over indices $i_k$, we introduce a modified Kronecker delta

$$ \delta_{\pi} = \delta_{\pi}(i_1, \dots, i_n) = \prod_{b \in \pi} \delta_{i_{b(1)}, \dots, i_{b(|b|)}} \quad (29) $$

which sets all indices equal that belong to the same block in $\pi$. Next, we extend the sum to include the cases in which indices $i_k$ can become equal by replacing $\mathbb{E}[\dots]$ with $\mathbb{E}[\dots]^c$. This introduces only a subleading error of order $1/L$ because Eq. (4) ensures that connected expectation values do not become more dominant when some indices are equal.

Then, the expansion of $\phi(G^n)$ in Eq. (20) reads

$$ \frac{1}{L} \sum_{\pi \in NC(n)} \sum_{i_1, \dots, i_n} \delta_{\pi} \prod_{b \in \pi^*} \mathbb{E}[G_{i_{b(1)}} i_{b(2)}} \dots G_{i_{b(|b|)}} i_{b(1)}]^c. \quad (30) $$

In free probability, $\phi$ corresponds to the “expectation” value of the random matrix $G$, and we show that the moments of $G$ have a natural expansion in terms of a sum of noncrossing partitions. The terms in this sum are usually called free cumulants. However, comparing to the definition of free cumulants in Eq. (38) in the next section, one sees that the expansion above has a slightly different form due to the additional $\delta_{\pi}$. This issue is further discussed in Sec. IV F.

To conclude the argument, note that we could perform the whole derivation multiplying $\mathbb{E}[G_{i_1 i_2} \dots G_{i_n i_1}]$ in Eq. (20) with a smooth test function $h_{i_1, \dots, i_n}$ that would again appear in Eq. (30). By comparing the two equations, we would find that

$$ \mathbb{E}[G_{i_1 i_2} \dots G_{i_n i_1}] = \sum_{\pi \in NC(n)} \delta_{\pi^*}(i_1, \dots, i_n) \prod_{b \in \pi} \mathbb{E}[G_{i_{b(1)}} i_{b(2)}} \dots G_{i_{b(|b|)}} i_{b(1)}]^c, \quad (31) $$

an equation we will reuse in Sec. IV. Here, we interchange the role of $\pi$ and $\pi^*$, which is possible because, for

noncrossing partitions, they are in one-to-one correspondence. A full proof of this formula for any $n$ is given in Appendix E 1.

### C. Analogy with ETH

Recently, it has been observed in Ref. [36] that there is a very similar link to free probability in the context of the eigenstate thermalization hypothesis (ETH). Indeed, on the mathematical level, fluctuations of spatial coherences $G_{ij}$ in one-dimensional mesoscopic systems seem to behave in complete analogy to matrix elements $A_{ij} = \langle E_i | A | E_j \rangle$ of a local observable $A$ expressed in the energy eigenbasis $|E_i\rangle$ of a Hamiltonian $H$ in a closed system that obeys the ETH. In the context of the ETH, $A_{ij}$ is a random variable with respect to a fictitious ETH-random-matrix ensemble that captures its typical behavior.

Comparing to Ref. [39], which discusses higher-order fluctuations of $A_{ij}$ within the ETH, one realizes that the reasons for the emergence of the three properties (1)–(3) are analogous. In the ETH, they are the result of an average over small energy windows. In our context, they are the result of an average over small space windows, which we call ballistic cells. To complete the analogy, loop expectation values of $A_{ij}$ in the ETH scale with the density of states $e^{S(E_+)}$ at energy $E_+ = \frac{1}{2}(E_i + E_j)$. In our context, this corresponds to the “density of states” at sites $i$ and $j$, which is just a constant equal to the number of sites $L$ (because the physical length of our system is set to 1, for simplicity).

## III. INTRODUCTION TO FREE PROBABILITY

In classical probability, two variables are independent if (and only if) their moments factorize at all orders, $\mathbb{E}[X^n Y^m] = \mathbb{E}[X^n] \mathbb{E}[Y^m]$ for all $n, m \in \mathbb{N}$. One can therefore determine joint moments of any product of independent variables knowing only the moments of the individual independent variables. If instead $X$ and $Y$ are random matrices with independent entries, then it is less clear how we would achieve the factorization of, say, $\mathbb{E}[XYXY]$ into the moments of the independent variables $\mathbb{E}[X^2]$ and $\mathbb{E}[Y^2]$, on the level of matrices, since they do not commute, in general.

Free probability theory solves this issue by proposing an extension of the notion of independence for noncommutative random variables, called freeness. Free variables are not only required to be independent in the probabilistic sense but also to be algebraically independent in the sense that there are no algebraic relations between the variables. This is similar to generators in a free group, hence the name.

Given two noncommuting random variables $a$ and $b$ in some algebra $\mathcal{M}$ (e.g., algebra of large random matrices) and a linear functional $\phi: \mathcal{M} \to \mathbb{C}$ (that plays the role of the expectation value), then $a$ and $b$ are called free if for all
polynomials $P_1, \dots, P_l$ and $Q_1, \dots, Q_l$ with $\varphi(P_i(a)) = 0$ and $\varphi(Q_i(b)) = 0$, we have

$$ \varphi(P_1(a)Q_1(b)\cdots P_l(a)Q_l(b)) = 0. \quad (32) $$

The reason to evoke all possible polynomials in the definition is that any element in the subalgebras $\mathcal{A}$ and $\mathcal{B}$ generated by $a$ and $b$, respectively, can be written as $P_i(a)$ and $Q_i(b)$. Hence, freeness can also be understood as a statement about the subalgebras $\mathcal{A}$ and $\mathcal{B}$.

This definition implies, for example, that if $a$ and $b$ are free, then $\varphi(a^m b^n) = \varphi(a^m)\varphi(b^n)$. Additional structure occurs if one interchanges the order such that free variables are no longer grouped together. However, it always remains true that joint moments can be determined through the moments of the individual free variables therein. For example, $\varphi(abab) = \varphi(a^2)\varphi(b)^2 + \varphi(a)^2\varphi(b^2) - \varphi(a)^2\varphi(b)^2$.

The definition of freeness was proposed by Voiculescu in 1985, who founded the field of free probability theory while working on problems in operator algebras. More details can be found in his book [31], and a good introduction to the subject is provided in the lecture notes by Speicher [42] as well as in the book by Mingo and Speicher [35].

In the 1990s, Speicher proposed a complementary combinatorial approach to free probability by introducing what he called free cumulants. This is also the route we took in the last section to make a connection between coherent fluctuations and free probability theory. Before introducing free cumulants, let us review how cumulants in classical probability theory are defined and how they are related to partitions. This will then allow us to better appreciate why free cumulants are defined via noncrossing partitions.

### A. Classical cumulants and partitions

Let $\{X_1, \dots, X_N\}$ be a family of classical random variables with the moment-generating function

$$ Z[a, u] := \mathbb{E}[e^{u \sum_i a_i X_i}] \\ = \sum_{n \ge 0} \frac{u^n}{n!} \sum_{i_1 \cdots i_n} a_{i_1} \cdots a_{i_n} \mathbb{E}[X_{i_1} \cdots X_{i_n}], $$

where $u$ is an (optional) parameter that counts the order of the joint moment (or correlation function) $\mathbb{E}[X_{i_1} \cdots X_{i_n}]$. The joint cumulant (or connected correlation function) $\mathbb{E}[X_{i_1} \cdots X_{i_n}]^c$ is defined as the term proportional to $a_{i_1} \cdots a_{i_n}$ in the expansion of the cumulant-generating function $W[a, u] := \log Z[a, u]$,

$$ W[a, u] = \sum_{n \ge 0} \frac{u^n}{n!} \sum_{i_1 \cdots i_n} a_{i_1} \cdots a_{i_n} \mathbb{E}[X_{i_1}, \dots, X_{i_n}]^c. $$

In fact, cumulants and moments are related by a combinatorial formula. Expanding $Z[a, u]$ in terms of the cumulants and grouping together terms with the same power of $u$, one can derive that a moment $\mathbb{E}[X_{i_1} \cdots X_{i_n}]$ can be expressed as a sum of products of cumulants arranged according to partitions $\pi$ of the set $\{i_1, \dots, i_n\}$,

$$ \mathbb{E}[X_{i_1} \cdots X_{i_n}] = \sum_{\pi \in P(n)} \prod_{b \in \pi} \mathbb{E}[X_{i_{b(1)}} X_{i_{b(2)}} \cdots]^c, \quad (33) $$

where $b = \{b(1), b(2), \dots\}$ denotes the elements of a block of the partition $\pi$. The number of partitions of a set of $n$ elements is called the Bell number $B_n$, with recursion relations $B_{n+1} = \sum_{k=0}^{n} \binom{n}{k} B_k$ and $B_1 = 1$, $B_2 = 2$, $B_3 = 5$, $B_4 = 15$, $B_5 = 52$, etc.

Let us give an example for $n = 4$ and the set $\{1, 2, 3, 4\}$. The partition $\pi = \{\{1, 2\}, \{3, 4\}\}$ is represented by the diagram

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

This is in full analogy to Eqs. (23) and (24), except that here we use solid lines for better visibility since we do not discuss dual partitions in this section. The expansion of the moment $\mathbb{E}[X_1 X_2 X_3 X_4]$ in terms of products of cumulants $\prod_{b \in \pi} \mathbb{E}[X_{i_{b(1)}} X_{i_{b(2)}} \cdots]^c$, represented through the diagram of the corresponding partition $\pi$, becomes

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

where the subscript $\mathfrak{S}k$ suggests that, by cyclic permutation, there are in total $k$ such diagrams. Note that the last diagram (in a dotted box) corresponds to a crossing partition $\pi = \{\{1, 3\}, \{2, 4\}\}$. In free probability theory, these diagrams do not appear, as we will see below.

### 1. Classical cumulants of a single variable

The moment-cumulant relation (33) allows us to express the cumulants recursively through the moments. In the case of a single variable $X = X_1 = \cdots = X_N$, we illustrate how this can be done up to order four. Let us denote by $m_n = \mathbb{E}[X^n]$ and $c_n = \mathbb{E}[X^n]^c$ the moments and cumulants of this variable; then,
$$ 
\begin{aligned}
m_1 &= c_1, \\
m_2 &= c_2 + c_1^2, \\
m_3 &= c_3 + 3c_2c_1 + c_1^3, \\
m_4 &= c_4 + 4c_1c_3 + 3c_2^2 + 6c_1^2c_2 + c_1^4.
\end{aligned} 
\quad (36) $$

Note that the coefficients correspond exactly to the cyclic multiplicities of the diagrams. This can be solved recursively for $c_k$,

$$ 
\begin{aligned}
c_1 &= m_1, \\
c_2 &= m_2 - m_1^2, \\
c_3 &= m_3 - 3m_1m_2 + 2m_1^3, \\
c_4 &= m_4 - 4m_1m_3 + 12m_1^2m_2 - 3m_2^2 - 6m_1^4.
\end{aligned} 
\quad (37) $$

The generating function of cumulants is the logarithm of that of moments, as explained above. Below, we will see how this formula differs for the free cumulants.

## B. Free cumulants and noncrossing partitions

In the setting of noncommuting random variables $a_1, \dots, a_N$ in some algebra $\mathcal{M}$ with expectation value $\varphi: \mathcal{M} \to \mathbb{C}$, the free cumulants $\kappa_n$ are multilinear forms implicitly defined through moments by a sum over non-crossing partitions $\pi \in NC(n)$ of the set $\{i_1, \dots, i_n\}$,

$$ \varphi(a_{i_1} \cdots a_{i_n}) := \sum_{\pi \in NC(n)} \prod_{b \in \pi} \kappa_{|b|}(a_{i_{b(1)}}, a_{i_{b(2)}}, \dots), \quad (38) $$

where $|b|$ denotes the number of elements in the block $b = \{b(1), b(2), \dots\}$. For example, if we expand $\varphi(a_1 a_2 a_3 a_4)$ in terms of free cumulants, we could obtain all diagrams in Eq. (35), except the last one, which is the crossing. Therefore, the order of the arguments in $\kappa_n$ becomes important, even if the $a_i$ commute, hence the separation by the comma.

The free cumulants satisfy a number of properties that are analogous to properties of classical cumulants:

(i) $\kappa_n(a_1, \dots, a_n) = 0$ as soon as there appears a pair $a_i, a_k$ of free variables; see, e.g., Ref. [42].

(ii) As a result of multilinearity and the previous point, free cumulants of free variables $a$ and $b$ are additive, $\kappa_n(a + b, \dots, a + b) = \kappa_n(a, \dots, a) + \kappa_n(b, \dots, b)$; see Ref. [43].

(iii) Any variable $a$ whose free cumulants $\kappa_n(a, \dots, a)$ vanish for $n \ge 3$ is distributed according to Wigner's semicircle law of random matrix theory for the Gaussian unitary ensemble (GUE), so large GUE random matrices are analogous to Gaussian variables but from a free-probability point of view.

The number of noncrossing partitions of a set of $n$ elements is the Catalan number $C_n = (1/n+1)^{(2n)}$, with $C_1 = 1, C_2 = 2, C_3 = 5, C_4 = 14, C_5 = 42$, etc. The

formula (38) is triangular and can be inverted to express the free cumulants in terms of the moments.

## 1. Free cumulants of a single variable

In the case of a single variable $a = a_1 = \dots = a_N$, we denote by $\kappa_n := \kappa_n(a, \dots, a)$ and $m_n := \varphi(a^n)$ the free cumulant and the moment at order $n$, respectively. Then, we have

$$ 
\begin{aligned}
m_1 &= \kappa_1, \\
m_2 &= \kappa_2 + \kappa_1^2, \\
m_3 &= \kappa_3 + 3\kappa_2\kappa_1 + \kappa_1^3, \\
m_4 &= \kappa_4 + 4\kappa_1\kappa_3 + 2\kappa_2^2 + 6\kappa_1^2\kappa_2 + \kappa_1^4.
\end{aligned} 
\quad (39) $$

The equations can be solved for $\kappa_k$ recursively,

$$ 
\begin{aligned}
\kappa_1 &= m_1, \\
\kappa_2 &= m_2 - m_1^2, \\
\kappa_3 &= m_3 - 3m_1m_2 + 2m_1^3, \\
\kappa_4 &= m_4 - 4m_1m_3 + 10m_1^2m_2 - 2m_2^2 - 5m_1^4.
\end{aligned} 
\quad (40) $$

Note that the difference between standard and free cumulants only shows up at order four since here a crossing partition becomes possible for the first time. For a single variable, the relation between the moments and the free cumulants can be found by inverting the resolvent associated to the distribution; see Refs. [42,44].

## C. Free cumulants in random matrix theory

A relation between free probability theory and random matrices was first observed by Voiculescu in 1991 [33]. For example, he realized that matrices in the GUE with independent entries become free variables in the limit where the dimension of the matrix goes to infinity. Since then, many more connections between other random matrix ensembles and free probability have been found.

Here, we make one of these results more explicit, which applies to matrices that are rotated by random unitaries. Consider $N \times N$ random matrices of the form $X_A = U_N A_N U_N^\dagger$ and $Y_N = U_N B_N U_N^\dagger$, where $U_N$'s are chosen independently from the Haar distribution over the unitary group, and $A_N$ and $B_N$ are deterministic matrices with spectral densities $\mu_A$ and $\mu_B$. In other words, the moments $m_k := \lim_{N \to \infty} (1/N) \text{Tr}(A_N^k) = \int \lambda^k \mu_A(\lambda) d\lambda$ are all finite, and similarly for $B_N$. Thus, in the limit $N \to \infty$, the random matrices $X_N$ and $Y_N$ become free variables $a$ and $b$ (in some noncommutative probability space) with distribution $\mu_A$ and $\mu_B$ with respect to the measure $\varphi := (1/L) \mathbb{E} \text{Tr}$, where $\mathbb{E}$ is the Haar measure.

Furthermore, it is known that the classical cumulants of such matrices $X_N$ can be expressed as the free cumulants $\kappa_n := \kappa_n(a, \dots, a)$ of the spectral density $\mu_A$. In other words,
$$ \mathbb{E}[e^{N \mathrm{Tr}(P_N X_N)}] \asymp_{N \to \infty} e^{N \sum_{n=1}^{\frac{1}{N}} \kappa_n \mathrm{Tr}(P_N^n)}, \quad (41) $$

where $P_N$ is a sequence of matrices with fixed rank [such that $\mathrm{Tr}(P_N^k)$ does not scale with $N$], for instance, a rank-one projector.

This connection between Haar-distributed random matrices and free cumulants is adapted to the closed Q-SSEP in the steady state—which actually corresponds to an equilibrium situation due to the absence of current in the steady state. Indeed, since in the closed case the Q-SSEP dynamics is unitary and ergodic over the unitary group, its steady-state distribution is the one induced by the Haar measure on the orbit of the initial matrix of coherences $G_0$ [26]. In other words, the matrix of coherence $G$ is distributed as $UG_0U^\dagger$, with $U$ a Haar-distributed unitary $L \times L$ matrix. As a consequence of Eq. (41), the large deviation function for coherence fluctuations in the closed Q-SSEP is the generating function of the free cumulants of the spectral measure of the initial matrix of coherences $G_0$.

## D. Alternative derivation of free probability in mesoscopic systems

As an alternative to Sec. II B, one can derive the link between fluctuations in mesoscopic systems and free probability starting from the moment-cumulant formula (33) introduced earlier in this section. We also use this method in Appendix E 1, which contains the full proof of this link for any $n$.

We apply Eq. (33) to a loop of $G_{ij}$'s and obtain

$$ \mathbb{E}[G_{i_1 i_2} \dots G_{i_n i_1}] = \sum_{\pi \in P(n)} \prod_{b \in \pi} \mathbb{E}[G_{i_{b(1)} i_{b(1)+1}} \dots G_{i_{b(|b|)} i_{b(|b|)+1}}]^c, \quad (42) $$

where the subscript $t$ of $\mathbb{E}_t$ is dropped, for simplicity.

For $n=4$, two possible partitions, each with two blocks, are

$$ \pi_1 = \{\{1,2\}, \{3,4\}\} = \begin{tikzpicture}[baseline=-0.5ex, scale=0.5, every node/.style={scale=0.5}]
    \node[circle, fill, inner sep=1pt] (1) at (0,0) {};
    \node[circle, fill, inner sep=1pt] (2) at (1,0) {};
    \node[circle, fill, inner sep=1pt] (3) at (0,1) {};
    \node[circle, fill, inner sep=1pt] (4) at (1,1) {};
    \draw[->] (1) -- (2) node[left] {$2$};
    \draw[->] (1) -- (3) node[left] {$3$};
    \draw[->] (1) -- (4) node[left] {$4$};
    \node[below] at (2) {$1$};
    \node[below] at (3) {$2$};
    \node[below] at (4) {$3$};
    \node[above] at (2) {$4$};
    \end{tikzpicture} \quad (43) $$

and

$$ \pi_2 = \{\{1,3\}, \{2,4\}\} = \begin{tikzpicture}[baseline=-0.5ex, scale=0.5, every node/.style={scale=0.5}]
    \node[circle, fill, inner sep=1pt] (1) at (0,0) {};
    \node[circle, fill, inner sep=1pt] (2) at (1,0) {};
    \node[circle, fill, inner sep=1pt] (3) at (0,1) {};
    \node[circle, fill, inner sep=1pt] (4) at (1,1) {};
    \draw[->] (1) -- (2) node[left] {$2$};
    \draw[->] (1) -- (3) node[left] {$3$};
    \draw[->] (1) -- (4) node[left] {$4$};
    \node[below] at (2) {$1$};
    \node[below] at (3) {$2$};
    \node[below] at (4) {$3$};
    \node[above] at (2) {$4$};
    \end{tikzpicture} \quad (44) $$

We determine the contribution of the two partitions:

$$ \pi_1: \begin{aligned} & \mathbb{E}[G_{i_1 i_2} G_{i_2 i_3}]^c \mathbb{E}[G_{i_3 i_4} G_{i_4 i_1}]^c \\ & = \delta_{i_1 i_3} \mathbb{E}[G_{i_1 i_2} G_{i_2 i_1}]^c \mathbb{E}[G_{i_3 i_4} G_{i_4 i_3}]^c = \mathcal{O}(L^{-3}). \end{aligned} \quad (45) $$

In the second line, we use the $U(1)$ invariance of $\mathbb{E}$, which implies that the expression is nonzero only for $i_1 = i_3$. Finally, we evaluate the scaling of the expression with $L$. Here, we use $\mathbb{E}[G_{ij} G_{ji}]^c = \mathcal{O}(L^{-1})$ from Eq. (4), and the Kronecker delta becomes $\delta_{ij} \to \delta(x-y)/L$ for large $L$, where $x = i/L$ and $y = j/L$. Note that the resulting scaling with $\mathcal{O}(L^{-3})$ is the same as that for four loops, which also appears in the sum over partitions as the partition $\{\{1,2,3,4\}\}$ that contains only a single block. It turns out that all noncrossing partitions have the same scaling with $L$.

In contrast, crossing partitions are subleading. Evaluating $\pi_2$, one has

$$ \pi_2: \begin{aligned} & \mathbb{E}[G_{i_1 i_2} G_{i_3 i_4}]^c \mathbb{E}[G_{i_2 i_3} G_{i_4 i_1}]^c \\ & = \delta_{i_2 i_3} \delta_{i_1 i_4} \mathbb{E}[G_{i_1 i_2} G_{i_2 i_1}]^c \mathbb{E}[G_{i_2 i_2} G_{i_1 i_1}]^c = \mathcal{O}(L^{-4}). \end{aligned} \quad (46) $$

Here, we assume that no connected expectation value is more dominant than that of loops with the same number of points, i.e., $\mathbb{E}[G_{i_1 i_1} G_{i_2 i_2}]^c = \mathcal{O}(L^{-1})$, because the loop with the same number of points scales as $\mathcal{O}(L^{-1})$.

Note that the average over independent ballistic cells introduced in Sec. II A actually implies that $\mathbb{E}[G_{ii} G_{kk}]^c = 0$ if $i \neq k$. However, we anticipate that the toy model Q-SSEP predicts this term to scale with $L^{-2}$ if $i \neq k$. A possible conclusion is that even though the average over ballistic cells is a good approximation for leading-order terms, it is probably too crude to capture the complete dependence of correlation functions on the system size $L$.

In the case of noncrossing partitions, the delta functions that we introduce each time to “close” the blocks into loops can be elegantly described via the dual partition (Kreweras complement) $\pi^*$. This case is best understood diagrammatically. We associate the diagram

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

with $\pi_1$. The solid lines represent blocks in $\pi_1$ (here, understood as a partition of the edges), while the dashed lines define the blocks of a partition of the nodes, the dual partition $\pi_1^* = \{\{1,3\}, \{2\}, \{4\}\}$. If $\pi$ is noncrossing, then there is a unique way to connect a maximal number of nodes by dashed lines without crossing a solid line. We associate a Kronecker delta $\delta_{\pi^*(i_1, \dots, i_n)}$ with a dual partition $\pi^*$, equating all indices belonging to the same block of $\pi^*$, exactly as in Eq. (29) but with $\pi$ and $\pi^*$ interchanged. For example, $\delta_{\pi_1^*(i_1, \dots, i_4)} = \delta_{i_1 i_3}$.

This notation allows us to rewrite the moment-cumulant formula (42) as a sum over noncrossing partitions only.
$$g_t^s(x_1, \dots, x_n) = \sum_{\pi \in NC(n)} \delta_{\pi} \prod_{b \in \pi^*} g_t(x_{b(1)}, \dots, x_{b(|b|)}), \quad (56)$$

where

$$\delta_{\pi} \equiv \prod_{b \in \pi^*} \delta(x_{b(1)}, \dots, x_{b(|b|)}) \quad (57)$$

is defined in analogy to Eq. (29) but for continuous variables.

Finally, in Sec. IV F, we show how to construct a very simple steady-state solution for connected-loop expectation values. This solution exploits the insight below Eq. (30), where we noted that the presence of the $\delta_{\pi}$ prevents us from viewing connected-loop expectation values as the free cumulants of the measure $\mathbb{E}_t$. If we remove the $\delta_{\pi}$, then connected-loop expectation values are, by definition, the free cumulants of a new measure $\varphi_t$ that is necessarily different from $\mathbb{E}_t$. Surprisingly, $\varphi_t$ has a very simple steady-state solution from which the connected-loop expectation values $g_{\infty}$ can be determined recursively.

## B. Introduction to the model

The quantum symmetric simple exclusion process was first introduced in Refs. [25] (closed case) and [27] (open case) and was further elaborated in Refs. [26,28,30].

### 1. Definition

The model describes a one-dimensional chain with L sites on which spinless fermions hop to neighboring sites with random amplitudes. The boundary sites are coupled to particle reservoirs [45]. The bulk evolution of the system is stochastic and unitary. We describe it in terms of the system's density matrix $\rho_t$ as

$$\rho_t \rightarrow \rho_{t+dt} = e^{-idH_t} \rho_t e^{idH_t}, \quad (58)$$

where the Hamiltonian increment is defined as

$$dH_t = \sum_{j=1}^{L-1} c_{j+1}^{\dagger} c_j dW_t^j + c_j^{\dagger} c_{j+1} d\bar{W}_t^j,$$

with $c_j^{\dagger}(c_j)$ fermion creation (annihilation) operators at site j with the usual commutation relations $\{c_j^{\dagger}, c_k\} = \delta_{jk}$ and $\{c_j^{\dagger}, c_k^{\dagger}\} = \{c_j, c_k\} = 0$. The noisy Hamiltonian increment depends on a complex Brownian motion $dW_t^j$, which determines the random hopping amplitude along the edge $(j, j+1)$ at time t. There is one complex Brownian motion per edge, which means that $W_t^j$ at each instance t is a centered Gaussian random variable that has independent increments $dW_t^j = W_{t+dt}^j - W_t^j$ with variance $\mathbb{E}[dW_t^j d\bar{W}_t^k] = J \delta^{j,k} dt$ if $t = t'$, and zero otherwise. Note that J is a rate

parameter with dimension one over time, and we set $J = 1$ in the following.

To the stochastic and unitary bulk evolution in Eq. (58), we add a deterministic but dissipative Lindbladian evolution

$$\mathcal{L}_{\text{bdry}} = \alpha_1 \mathcal{L}_1^+ + \beta_1 \mathcal{L}_1^- + \alpha_L \mathcal{L}_L^+ + \beta_L \mathcal{L}_L^-,$$

which represents the interaction with the reservoirs. The operator $\mathcal{L}_j^+(\cdot) = c_j^{\dagger} \cdot c_j - \frac{1}{2} \{c_j c_j^{\dagger}, \cdot\}$ models particle injection and is multiplied by the injection rate $\alpha_j$, while $\mathcal{L}_j^-(\cdot) = c_j \cdot c_j^{\dagger} - \frac{1}{2} \{c_j^{\dagger} c_j, \cdot\}$ models particle extraction with rate $\beta_j$. For example, the density matrix of an isolated empty site $\tau_t = |0\rangle\langle 0|$ that evolves according to $\partial_t \tau_t = \alpha \mathcal{L}^+(\tau_t)$ will be occupied after a time interval $dt$ with probability $\alpha dt$. In other words,

$$|0\rangle\langle 0| \rightarrow \alpha dt |1\rangle\langle 1| + (1 - \alpha dt) |0\rangle\langle 0|.$$

The full evolution of the system's density matrix $\rho_t$ can be expressed as a stochastic differential equation (SDE) in the Itô convention (with Itô rules $dW_t^j d\bar{W}_t^k = \delta^{j,k} dt$) by expanding Eq. (58) up to order $dt$,

$$d\rho_t = -i[dH_t, \rho_t] - \frac{1}{2}[dH_t, [dH_t, \rho_t]] + \mathcal{L}_{\text{bdry}}(\rho_t)dt. \quad (59)$$

Occasionally, we refer to the “closed case,” by which we mean that there are no boundary reservoirs and the one-dimensional chain is closed periodically such that sites $1 \equiv L$ are identified.

### 2. Relation to the classical SSEP

The name of this model is inherited from the classical SSEP. The latter can be obtained from Q-SSEP in the special case where one is just interested in the mean density matrix $\bar{\rho}_t := \mathbb{E}[\rho_t]$, where the expectation $\mathbb{E}[\cdot \cdot \cdot]$ is taken with respect to the different realizations of the Brownian motions. This matrix evolves according to a Lindblad equation,

$$\partial_t \bar{\rho}_t = \mathcal{L}(\bar{\rho}_t) = \sum_{j=1}^{L-1} \mathcal{L}_j^{\text{edge}}(\bar{\rho}_t) + \mathcal{L}_{\text{bdry}}(\bar{\rho}_t),$$

with $\mathcal{L}_j^{\text{edge}}(\bar{\rho}) = l_j^- \bar{\rho} l_j^+ + l_j^+ \bar{\rho} l_j^- - \frac{1}{2} \{l_j^+ l_j^- + l_j^- l_j^+, \bar{\rho}\}$, where $l_j^+ = c_{j+1}^{\dagger} c_j$ and $l_j^- = c_j^{\dagger} c_{j+1}$. Writing $\bar{\rho}_t$ in the occupation number basis, this Lindbladian evolution preserves its diagonal structure, while off-diagonal elements vanish exponentially in time. The diagonal elements of $\bar{\rho}_t$ provide the probability that the system is in one of the $2^L$ states with well-defined occupation number $\hat{n}_i = c_i^{\dagger} c_i$ on each site. This corresponds to a configuration of the classical SSEP, where one specifies the number of particles $n_i = 0, 1$ on
each site. One can easily see that the Lindbladian evolution of the diagonal elements of $\bar{\rho}_t$ precisely corresponds to the Markov process of SSEP: During a time interval $dt$, a particle in the bulk can jump to the left or right neighboring site with probability $dt$ if the site is empty and particles are injected (extracted) at the boundaries with probability $\alpha_i dt$ ($\beta_i dt$). This correspondence can also be formulated via the moment-generating function,

$$ \langle e^{\sum_i a_i n_i} \rangle_{\text{ssep}} = \text{Tr}(\bar{\rho} e^{\sum_i a_i \hat{n}_i}), \quad (60) $$

illustrating the well-known correspondence between the Markov process and Lindbladian evolution. The relation between Q-SSEP and free cumulants we describe below implies a new relation between the classical SSEP and free probability [46].

### 3. Fluctuating coherences

If we go beyond the mean dynamics, then Q-SSEP has an additional structure, which describes purely quantum mechanical effects such as entanglement. This structure is inscribed into the coherences $G_{ij}(t) = \text{Tr}(\rho_t c_i^\dagger c_j)$. While

the coherences vanish in mean exponentially with time (we see that the mean corresponds to the classical SSEP, and there are no quantum correlations left), their fluctuations survive the long time limit, although they are subleading in the system size [27]. For example, at large system size $L \to \infty$ with $x = i/L \le y = j/L$ fixed, the connected quadratic fluctuation (or second cumulant) in the steady state is

$$ \mathbb{E}_{\infty}[G_{ij}G_{ji}]^c \sim \frac{1}{L}(\Delta n)^2 x(1-y), $$

where $\Delta n = n_a - n_b$ is the difference in the particle density between the boundary reservoirs; see Eq. (64). Here, we again adopt the point of view where the time dependence of $G_{ij}(t)$ is transferred to the measure $\mathbb{E}_t$ and the steady measure is denoted by $\mathbb{E}_\infty$.

Let us also note that since the Hamiltonian is quadratic in fermionic creation and annihilation operators, the coherences $G_{ij}$ completely characterize the state of the system according to Wick's theorem. From Eq. (59), one finds that their time evolution is given by the SDE

$$ \begin{aligned} dG_{ij} = & -i(G_{i,j-1}dW_t^{j-1} + G_{i,j+1}d\tilde{W}_t^j - G_{i-1,j}d\tilde{W}_t^{j-1} - G_{i+1,j}dW_t^i) + \delta_{ij}(G_{i+1,i+1} + G_{i-1,i-1})dt \\ & - 2G_{ij}dt - \sum_{p \in \{1,L\}} \left(\frac{1}{2}(\delta_{ip} + \delta_{jp})(\alpha_p + \beta_p - 1)G_{ij} - \delta_{pi}\delta_{pj}\alpha_p\right)dt. \end{aligned} \quad (61) $$

### 4. Three properties of $\mathbb{E}_t$

First, the measure $\mathbb{E}_t$ of Q-SSEP possesses a local $U(1)$ invariance, which can be seen as follows: The multiplication with local phases $G_{ij} \to \tilde{G}_{ij} = e^{-i\theta_i} G_{ij} e^{i\theta_j}$ leaves Eq. (61) invariant if the Brownian motions are multiplied by a phase $dW_t^j \to d\tilde{W}_t^j = e^{i(\theta_{j+1}-\theta_j)} dW_t^j$. Since $d\tilde{W}_t^j$ and $dW_t^j$ have the same distributions, $G$ and $\tilde{G}$ also have the same distribution if they agree at $t=0$, which means that initially $G$ is a diagonal matrix. In the scaling limit, this is always the case since off-diagonal terms vanish exponentially fast in system size (about $e^{-L^2t}$, where $t$ is the rescaled time).

Second, $n$-loop expectation values in Q-SSEP scale as $L^{-n+1}$ if indices are distinct because this scaling gives rise to the nontrivial equation (51).

Third, expectation values of products of loops factorize. This result is shown in Appendix G of Ref. [30], where the time-evolution equation of loop expectation values in Eq. (51) is derived (for the closed case). The time-evolution equation for the leading order of a product of loop expectation values allows for factorized solutions. Since the initial condition is always factorized (the noise has not yet correlated with the variables), the claim follows.

For Q-SSEP, properties (1)–(3) are satisfied, and as a result, in Q-SSEP, loop expectation values expand into noncrossing partitions [Eq. (31)] as claimed.

### C. Bulk vs boundary thermalization

When we casually talk about bulk and boundary thermalization, one should keep in mind that the open Q-SSEP does not actually “thermalize.” At large times, the system approaches a steady state [47], in which observables do not depend on time anymore. However, they are not described by a thermal density matrix since there is a steady current, so the system is out of equilibrium.

Here, we show that the density of fermions in the open Q-SSEP approaches its steady-state value much faster on the boundaries than in the bulk. In particular, the decay times scale with $O(L^2)$ in the bulk and $O(1)$ on the boundary. Because of the correspondence (60), this property also holds true for the classical SSEP [48].

Under the evolution outlined above, the mean fermion density $n_i(t) := \mathbb{E}_t[\text{Tr}(\rho_c c_i^\dagger c_i)]$ evolves according to

$$ \partial_t n_i(t) = \Delta n_i(t) + \sum_{p=1,L} \delta_{ip} (\alpha_p - (\alpha_p + \beta_p) n_p(t)), \quad (62) $$
where it is understood that the discrete Laplacian on the boundaries is truncated, i.e., $\Delta n_1 := n_2 - n_1$ and $\Delta n_L = n_{L-1} - n_L$. Here, we only discuss the special case, where the injection or extraction parameters are such that $\alpha_1 + \beta_1 = 1 = \alpha_L + \beta_L$. The general case is discussed in Appendix C, and it works analogously, with some help from numerical calculations. In the special case, Eq. (62) can be solved analytically,

$$n_j(t) = \sum_{k=1}^{L} \sin\left(\frac{jk\pi}{L+1}\right) \left(e^{(-2+2\cos(\frac{kg}{L+1}))t}(c_k - b_k) + b_k\right), \quad (63)$$

where the coefficients $c_k$ are determined by the initial condition and $b_k := (2/L + 1)[\alpha_1 \sin(k\pi/L + 1) + \alpha_L \sin(Lk\pi/L + 1)/2(1 - \cos(k\pi/L + 1))]$. This solution consists of two contributions: a time-dependent term that decays exponentially in time and a constant term that provides the steady-state value and can be simplified to $n_j(\infty) = [\alpha_1(L - j + 1) + \alpha_L j/L + 1]$. We study the time-scale with which the time-dependent term decays on the boundary compared to its decay in the bulk.

A bulk site is characterized by $j \sim aL$ with $a \sim \mathcal{O}(1)$. Because of the factor $e^{(-2+2\cos(k\pi/L+1))t} \approx e^{-\pi^2 k^2 t/L^2}$ (for large $L$), the biggest contribution to the sum comes from the term with $k=1$. Since we are in the bulk, its amplitude is finite, $\sin(j\pi/(L+1)) \sim \sin(a\pi)$. Therefore, we find a timescale of $t_{\text{decay}} \sim \mathcal{O}(L^2)$,

$$n_{\text{bulk}}(t) \sim \text{const} e^{-\pi^2 t/L^2} + n_{\text{bulk}}(\infty).$$

A site on the boundary is characterized by $j=1$ or $j=L$; therefore, the amplitude of the $k=1$ term, $\sin(j\pi/(L+1))$, will be zero for large $L$. Thus, a significant contribution only comes from terms where $k \sim bL$, with $b \sim \mathcal{O}(1)$, because then the amplitudes $\sin(jk\pi/(L+1)) \sim \sin(b\pi)$ or $\sim \sin(bL\pi)$ stay finite. The timescale with which these terms decay is of order one, $t_{\text{decay}} \sim \mathcal{O}(1)$,

$$n_{\text{bdry}}(t) \sim \text{const} e^{-b^2 \pi^2 t} + n_{\text{bdry}}(\infty).$$

Note that the value of the density on the boundaries after a time of $\mathcal{O}(1)$ is that of the steady state. For general injection or extraction parameters, these values are as in the classical SSEP,

$$n_1(\infty) = n_a := \frac{\alpha_1}{\alpha_1 + \beta_1}, \quad n_L(\infty) = n_b := \frac{\alpha_L}{\alpha_L + \beta_L}. \quad (64)$$

We conclude that timescales with which the boundary and bulk approach their steady-state values are separated by an order of $L^2$. In the scaling limit (49) where $t \to t/L^2$ denotes the rescaled time, this means that boundaries are

equal to the steady-state values at all times, while in the bulk it takes $t_{\text{decay}} \sim \mathcal{O}(1)$.

## D. Open boundary scaling limit

The basic question is whether the scaling limit (49) introduced for the closed Q-SSEP is also meaningful in the open case. In other words, does $g_t^s$ satisfy a nontrivial equation in this limit?

### 1. Density

We first answer this question on the level of the density $n_i(t) = \mathbb{E}_t[G_{ii}]$. It satisfies $L$ coupled differential equations in time given by Eq. (62). The aim is to represent these equations by a single partial differential equation in space and time for a continuous density $\rho(x, t) \approx n_{Lx}(L^2t)$. Since it is easily seen that the bulk satisfies a pure diffusion equation, the remaining question is what the correct boundary conditions are. The answer is provided by the observation made in the last section. In the scaling limit, the boundaries immediately take their steady-state values. Therefore,

$$\partial_t \rho(x, t) = \partial_x^2 \rho(x, t), \quad (65)$$

$$\rho(0, t) = n_a, \quad \rho(1, t) = n_b. \quad (66)$$

This result is confirmed numerically in Fig. (3), where a solution for the discrete density $n_i(t)$ is compared to the solution for the continuous density $\rho(x, t)$. In particular, the boundary values agree.

### 2. Higher-order fluctuations

The claim we make here is that the exact same equation [Eq. (51)] that holds in the closed case also holds in the

<figure><img src="image_8.png" /><figcaption>FIG. 3. Discrete fermion density $n_{Lx}(L^2t)$ for system sizes $L=24$ and $L=48$, together with the scaling limit $\rho(x, t)$ at $t=0.01$ as a function of space $x=i/L$. The extraction and injection rates are $\alpha_1 = \beta_1 = \alpha_L = \beta_L = 1$, and they do not fit the initial conditions. The agreement is very good.</figcaption></figure>
open case if we subject $g^s$ to boundary conditions that are equal to its steady-state value. These boundary conditions are known from Ref. [27] for the connected part: $g_\infty(x_1, ..., x_n) = 0$ for all $n \ge 2$ and $(x_1, ..., x_n) \in \partial[0, 1]^n$.

We test this claim on the level of the two-loop expectation values, for which Eq. (51) simplifies to

$$ (\partial_t - \Delta)g_t^s(x, y) = 2\partial_x\partial_y(\delta(x-y)g_t(x)g_t(y)), \quad (67) $$

and we identify one-loop values with the density, $g_t(x) = \rho(x, t)$. Thus, it follows that the time evolution of the connected correlation functions with appropriate boundary conditions is

$$ (\partial_t - \Delta)g_t(x, y) = 2\delta(x-y)\partial_x g_t(x)\partial_y g_t(y), \\ g_t(x, y) = 0 \quad \text{if } (x, y) \in \partial[0, 1]^2, \quad (68) $$

which is again a diffusion equation but with a new source term. To see this, note that the connected correlation function is defined by

$$ g_t(x, y) = \lim_{L \to \infty} L(\mathbb{E}_{L^2 t}[G_{ij}G_{ji}] - \delta_{ij}\mathbb{E}_{L^2 t}[G_{ii}]\mathbb{E}_{L^2 t}[G_{jj}]) \\ = g_t^s(x, y) - \delta(x-y)g_t(x)g_t(y). $$

The analytic solution of Eq. (68) (which is constructed in Appendix D) is compared to a numerical solution of the discrete evolution equations of $L(\mathbb{E}_{L^2 t}[G_{ij}G_{ji}] - \delta_{ij}\mathbb{E}_{L^2 t}[G_{ii}]^2)$ for given $L$, which can be derived from Eq. (61). Figures 4 and 5 show the result of this comparison for different system sizes $L$. The agreement is excellent.

Notice that the injection and extraction rates enter Eq. (68) only through the source term since the density $g_t(x)$ depends on these rates through its own boundary conditions. Apart from that, the equation for the connected two-loop expectation values $g_t(x, y)$ (and for all higher-order loops) never explicitly depend on the injection and

<figure><img src="image_9.png" /><figcaption>FIG. 4. Boundary conditions that fit the initial domain wall state ($n_a = 1, n_b = 0$).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 5. Boundary conditions that do not fit the initial domain wall state ($n_a = 1/2, n_b = 1/2$).</figcaption></figure>

extraction rates. As a consequence, without loss of generality, we use the convention

$$ n_a = 0, \quad n_b = 1 $$

throughout the paper. In particular, the density in the steady state $g_\infty(x) = n_a + x(n_b - n_a)$ reduces to $g_\infty(x) = x$ in this convention.

As a result, one can now derive the time-evolution equation of connected-loop expectation values in the scaling limit with open boundaries (see Appendix E 2). One finds

$$ (\partial_t - \Delta)g_t(x_1, \dots, x_n) = \sum_{i<j}^{n} 2\delta(x_i, x_j) \partial_i g_t(x_i, \dots, x_{j-1}) \partial_j g_t(x_j, \dots, x_{i-1}), \quad (69) $$

with boundary conditions

$$ g_t(x_1, \dots, x_n) = \begin{cases} n_a, & n_b \text{ for } n=1 \text{ and } x=0, 1 \\ 0, & n \ge 2 \text{ and some } x_i \in \{0, 1\}. \end{cases} \quad (70) $$

## E. Free cumulants in Q-SSEP

Comparing the moment-cumulant relation (31) [valid for any system satisfying Eqs. (1)–(3), in particular, for Q-SSEP] to the definition of free cumulants (38), one would be tempted to identify

$$ \mathbb{E}_t[G_{i_1 i_2} \cdots G_{i_n i_1}]^c \sim \kappa_n(G_{i_1 i_2}, \dots, G_{i_n i_1}), \quad (71) $$

with $\kappa_n$ the free cumulants of $\mathbb{E}_t$. However, this is not correct because of the presence of $\delta_\pi$. If we nonetheless insist on the identification (71), then the connected-loop expectation values $\mathbb{E}[G_{i_1 i_2} \cdots G_{i_n i_1}]^c$ are, by definition, the
free cumulants of a new measure $\varphi_t$ that is different from $E_t$. In other words,

$$ \varphi_t(G_{i_1i_2} \cdots G_{i_ni_1}) \\ := \sum_{\pi \in NC(n)} \prod_{b \in \pi} \mathbb{E}_t[G_{i_{b(1)}}i_{b(2)} \cdots G_{i_{b(|b|)}}i_{b(1)}]^c. \quad (72) $$

In terms of $g_t(x_1, \dots, x_n) \sim L^{n-1}\mathbb{E}_{L^{2t}}[G_{i_1i_2} \cdots G_{i_ni_1}]^c$, the new measure $\varphi_t$ defines a function (which we denote by the same name) of continuous positions and rescaled time

$$ \varphi_t(x_1, \dots, x_n) = \sum_{\pi \in NC(n)} \prod_{b \in \pi} g_t(x_{b(1)}, \dots, x_{b(|b|)}). \quad (73) $$

The time evolution of $\varphi_t$ is found to satisfy exactly the same equation as $g_t$ (see Appendix E 3),

$$ (\partial_t - \Delta)\varphi_t(x_1, \dots, x_n) \\ = \sum_{i<j}^{n} 2\delta(x_i, x_j) \partial_i \varphi_t(x_i, \dots, x_{j-1}) \partial_j \varphi_t(x_j, \dots, x_{i-1}). \quad (74) $$

However, the boundary conditions are different. If some $x_i \in \{0, 1\}$ lies on the boundary, then

$$ \varphi_t(x_1, \dots, x_i, \dots, x_n) = x_i \varphi_t(x_1, \dots, \hat{x}_i, \dots, x_n), \quad (75) $$

where the hat on $\hat{x}_i$ indicates that $x_i$ is missing from the set $\{x_1, \dots, x_n\}$.

## F. Steady-state solution—Inspired by free probability

The striking observation by Biane [32] was that $\varphi_t$, defined as a sum over noncrossing partitions of products of connected-loop expectation values, has a very simple solution in the steady state. In Appendix E 3, we show that for $t \to \infty$, our equation for $\varphi_t$ is indeed solved by

$$ \varphi_{\infty}(x_1, \dots, x_n) = \min(x_1, \dots, x_n). \quad (76) $$

As a consequence, we can use Eq. (73) to recursively reconstruct the connected correlation functions $g_\infty$ in the steady state. This works similarly to Eqs. (39) and (40). Denoting $\min(x_1, \dots, x_n) =: x_1 \wedge \cdots \wedge x_n$, we have

$$ g_\infty(x_1) = x_1, \\ g_\infty(x_1, x_2) = x_1 \wedge x_2 - x_1 x_2, \\ g_\infty(x_1, x_2, x_3) = x_1 \wedge x_2 \wedge x_3 - x_1(x_2 \wedge x_3)_{\oplus 3} + 2x_1x_2x_3, $$

while the four-point function $g_\infty(x_1, x_2, x_3, x_4)$ reads

$$ x_1 \wedge x_2 \wedge x_3 \wedge x_4 - x_1(x_2 \wedge x_3 \wedge x_4)_{\oplus 4} \\ - (x_1 \wedge x_2)(x_3 \wedge x_4)_{\oplus 2} + 2x_1x_2(x_3 \wedge x_4)_{\oplus 4} \\ + x_1x_3(x_2 \wedge x_4)_{\oplus 2} - 5x_1x_2x_3x_4, $$

where $\cdots_{\oplus q}$ denotes the sum of all terms obtained by $q$ successive cyclic permutations of the arguments of the term in question. Note that because of the absence of crossing partitions, $g_\infty(x_1, \dots, x_n)$ is no longer invariant under the permutation of its arguments for $n \ge 4$. In the example above, it is the term $(x_1 \wedge x_2)(x_3 \wedge x_4)_{\oplus 2}$ that prevents $g_\infty(x_1, \dots, x_4)$ from being invariant under the exchange of $x_2$ and $x_3$.

Indeed, ordering the variables as $0 \le x_1 \le x_2 \le x_3 \le x_4 \le 1$, we get $g_\infty(x_1, x_2) = x_1(1-x_2)$ and $g_\infty(x_1, x_2, x_3) = x_1(1-2x_2)(1-x_3)$ and

$$ g_\infty(x_1, x_2, x_3, x_4) = x_1(1 - 3x_2 - 2x_3 + 5x_2x_3)(1 - x_4), \\ g_\infty(x_1, x_3, x_2, x_4) = x_1(1 - 4x_2 - x_3 + 5x_2x_3)(1 - x_4), $$

whereas $g_\infty(x_1, x_3, x_4, x_2) = g_\infty(x_1, x_2, x_3, x_4)$, in agreement with Ref. [27].

The derivation of the steady-state solution presented here provides an alternative proof of Biane's formula (6.1) in Ref. [32] that relates the steady-state connected correlations of the open Q-SSEP to free cumulants of the measure $\varphi_\infty$. Our derivation extends this relation to finite times, though an explicit solution for $\varphi_t$ at finite times seems out of reach.

Note that the measure $\varphi_\infty$ can be realized as the Lebesgue measure on the interval $[0, 1]$ of the indicator function $\mathbb{I}_x := 1_{[0,x]}$, e.g.,

$$ \varphi_{\infty}(x, y) := \int_{0}^{1} \mathbb{I}_{x}(z) \mathbb{I}_{y}(z) dz = \min(x, y). \quad (77) $$

It is surprising that $\varphi_\infty$ has a realization in terms of commuting variables since free cumulants usually appear in a setting of noncommuting variables. At finite times, Eq. (74) suggests that $\varphi_t$ is not invariant under a permutation of its arguments.

## V. CONCLUSION

In this paper, we presented two main points: (1) a general argument as to why the fluctuations of spatial coherences in one-dimensional mesoscopic quantum systems could be well described by the framework of free probability theory and (2) a precise calculation that shows that the model Q-SSEP has a mathematical structure that fits into the framework of free probability (we used this structure to derive the time evolution of connected correlation functions). Specific to the open Q-SSEP is the observation that the density approaches its steady-state value much faster on the boundary than in the bulk, which we used to formulate
the correct boundary conditions for the time evolution of coherences in the scaling limit.

In both cases, the link to free probability can be reduced to three properties of the noise expectation value: (i) local $U(1)$ invariance, (ii) a large deviation scaling of correlation functions $\mathbb{E}[G_{i_1i_2} \cdots G_{i_ni_1}] \sim L^{-n+1}$, and (iii) the fact that expectation values of products of loops factorize. It is not surprising that the same three properties are responsible for the fact that a relation with free probability has been observed in the context of the ETH.

We should also stress that the link with free probability in the context of coherent fluctuations in mesoscopic systems is, in fact, not particular to the systems being out of equilibrium. Rather, this link emerges from a coarse-grained description under the assumption that, locally in space (i.e., within ballistic cells) and on timescales much shorter than the diffusion time, the system is ergodic. Such an assumption allows us to introduce the noise average as an average over all possible unitary transformations that the system could have undergone locally (i.e., within ballistic cells), and this noise average satisfies the three properties above.

However, we admit that the general picture developed in Sec. II is certainly oversimplified and calls for more details. First, the argument assumes a separation of timescales that gives rise to a crossover from ballistic to diffusive transport at some length scale $\ell$. One should provide criteria on the Hamiltonian for when this separation of timescales is satisfied. Second, one should explore the physical meaning of the length $\ell$ of ballistic cells. In the Introduction, we crudely argued that $\ell$ could be related to the mean free path of electrons in a disordered metal. Third, a better understanding of the perturbative expansion of an interacting Hamiltonian—the argument needed to derive the scaling of loop expectation values from classical MFT—will allow us to characterize the domain of validity of this theory, that is, for which class of systems it should be applicable.

To test the validity of this general picture, we plan to conduct numerical tests on more physical models. A first candidate would be a Floquet Heisenberg XXZ model with a staggered magnetic field, which breaks integrability but conserves the local $U(1)$ invariance. The time dependence of the Hamiltonian ensures that energy is not conserved, which would otherwise represent another conserved quantity.

Such numerical studies would also answer the question of whether or not Q-SSEP describes coherent fluctuations in a larger class of systems through the identification of sites in Q-SSEP with ballistic cells—an idea we developed in the Introduction. This interpretation is supported by a recent work [49], where the authors used a decomposition into equilibrated and statistically independent cells to characterize transport properties of quantum stochastic Hamiltonians. In particular, for the so-called dephasing model, i.e., free fermions with an independent Brownian

noise on each site, they found that the size of the cell scales as $1/\gamma$, where $\gamma$ is the strength of the noise. This result is interesting since the Q-SSEP can be obtained as a limit of the dephasing model for strong noise $\gamma \to \infty$ [25]. In this limit, the size of the cell becomes zero, which corresponds to the idea that for Q-SSEP, the mean free path $\ell$ has effectively shrunk to the lattice spacing $a_{uv}$, i.e., the ballistic cell contracts to a single site.

In spite of the progress reported here, we believe that it remains a challenge to construct a quantum mesoscopic fluctuation theory describing fluctuations of quantum coherences in generic, diffusive, many-body systems at coarse-grained mesoscopic scales.

## ACKNOWLEDGMENTS

This paper has been submitted simultaneously with “Eigenstate Thermalization Hypothesis and Free Probability” by Pappalardi, Foini, and Kurchan [36], where the relation between these two frameworks is discussed. The occurrence of free probability in both problems has a similar origin: the coarse graining at microscopic spatial or energy scales, and the unitary invariance at these microscopic scales. Thus, the use of free probability tools promises to be ubiquitous in chaotic or noisy many-body quantum systems. We thank Fabian Essler and Adam Nahum for numerous discussions on this topic. L. H. thanks Tony Jin for discussions about the general picture of fluctuating mesoscopic systems. This work was, in part, supported by CNRS, by the ENS, and by the ANR project “ESQuisses,” Contract No. ANR-20-CE47-0014-01.

## APPENDIX A: MEASURING COHERENCES

To gain a better understanding of the coherences $G_{ij} = \text{Tr}(\rho c_j^\dagger c_i)$, we outline an experiment to measure them, which was proposed in Ref. [18]. The setup is shown in Fig. 6. The idea is to probe the system at spatially separated places and let the two signals interfere before each output is measured separately. Let us outline the steps of the measurement protocols in detail.

(i) The total state of the system, left and right wires, is described by a state in the Hilbert space $\mathcal{H}_S \otimes \mathcal{H}_L \otimes \mathcal{H}_R$. Let us assume that the system is in a pure state and that, initially, the wires are empty and not yet coupled to the system,

$$|\psi^{(0)}\rangle = |\psi_S\rangle|0,0\rangle.$$

(ii) Now, we couple the two wires to the system. A very simple description of this coupling could be given by the unitary evolution with $U_{\text{int}} = e^{-i\lambda(c_L^\dagger c_i + c_R^\dagger c_j + \text{H.c.})}$, where $\lambda$ is the product of coupling strength and the time during which we allow the wires to couple to the system, and $c_L$ ($c_R$) are fermionic operators on the left (right) wire. If we tune the coupling strength
<figure><img src="image_11.png" /><figcaption>FIG. 6. Two wires attached to the system at sites i and j such that only one fermion can enter at a time. First, the fermions in the wire are allowed to interact via the beam splitter S. Then, their occupation numbers $n_L$ and $n_R$ are measured on each side. In the first measurement (a), one uses a symmetric beam splitter, which allows one to measure the imaginary part of $G_{ij}$. In the second measurement (b), one needs to use a beam splitter where the fermion that is transmitted from R to L accumulates a phase $\pi$, while it does not accumulate this phase when being transmitted in the other direction. In this way, one can measure the real part of $G_{ij}$.</figcaption></figure>

and duration such that $\lambda \ll 1$ is small, we can neglect $O(\lambda^2)$ terms and find the complete state to be

$$|\psi^{(1)}\rangle := U_{\text{int}}|\psi^{(0)}\rangle \\ = |\psi_S\rangle|0,0\rangle - i\lambda(c_i|\psi_S\rangle|1,0\rangle + c_j|\psi_S\rangle|0,1\rangle) \\ + O(\lambda^2).$$

(iii) Next, the fermions on the left and right wires interfere in a beam splitter. Written in the basis $\{|00\rangle, |01\rangle, |10\rangle, |11\rangle\}$, the beam splitter can, in general, be described by the scattering matrix

$$S = \begin{pmatrix} 1 & & & \\ & r' & t & \\ & & t' & r \\ & & & rr' - tt' \end{pmatrix},$$

where $r$ and $t$ ($r'$ and $t'$) are the reflection and transmission amplitudes for a fermion incident from the left (right) side. Though not important in our case—since the state $|1, 1\rangle$, where there is a fermion in each wire, is suppressed by $\lambda^2$—we give a short explanation for how to obtain the last entry $|1, 1\rangle \rightarrow (rr' - tt')|1, 1\rangle$. One has to take into account that the wave function is antisymmetric, $|1, 1\rangle = |\phi_L\rangle \otimes |\phi_R\rangle - |\phi_R\rangle \otimes |\phi_L\rangle$. Here, the position in the tensor product labels the fermion (say they are called 1 and 2), whereas $|\phi_L\rangle$ and $|\phi_R\rangle$ are single fermion states on the left and right wires. Then,

$$|1, 1\rangle \rightarrow (r|\phi_L\rangle + t|\phi_R\rangle) \otimes (r'|\phi_R\rangle + t'|\phi_L\rangle) \\ - (r'|\phi_R\rangle + t'|\phi_L\rangle) \otimes (r|\phi_L\rangle + t|\phi_R\rangle)$$

leads to the entry $(rr' - tt')$. Unitarity demands $|r|^2 + |t|^2 = 1$, $|r'|^2 + |t'|^2 = 1$ and $\bar{r}t' + \bar{t}r' = 0$

(the condition $|rr' - tt'|^2 = 1$ is then automatically fulfilled).

(a) Choosing a symmetric beam splitter, $r = r'$ and $t = t'$, allows one to measure the imaginary part of $G_{ij}$. Note that the unitary constraints can now be expressed as $r = \sin\theta e^{i\varphi}$ and $t = i \cos\theta e^{i\varphi}$, and we set the overall phase $\varphi = 0$ since this will not change the result. The state evolves to

$$|\psi^{(2,a)}\rangle = S^{(a)}|\psi^{(1)}\rangle \\ = |\psi_S\rangle|0,0\rangle \\ - i\lambda c_i|\psi_S\rangle(\sin\theta|1,0\rangle + i\cos\theta|0,1\rangle) \\ - i\lambda c_j|\psi_S\rangle(\sin\theta|0,1\rangle + i\cos\theta|1,0\rangle).$$

(b) If the beam splitter is symmetric—except for the fact that transmitted fermions incident from the right will accumulate an additional phase $\pi$, i.e., $r = r'$ and $t = -t'$—this allows one to measure the real part of $G_{ij}$. Note that the unitary constraints result in $r = \sin\theta e^{i\varphi}$ and $t = \cos\theta e^{i\varphi}$, and again we set $\varphi = 0$. The state evolves to

$$|\psi^{(2,b)}\rangle = S^{(b)}|\psi^{(1)}\rangle \\ = |\psi_S\rangle|0,0\rangle \\ - i\lambda c_i|\psi_S\rangle(\sin\theta|1,0\rangle + \cos\theta|0,1\rangle) \\ - i\lambda c_j|\psi_S\rangle(\sin\theta|0,1\rangle - \cos\theta|1,0\rangle).$$

(iv) Finally, we measure the particle number $n_L = c_L^\dagger c_L$ and $n_R = c_R^\dagger c_R$ on the left and right wires. Denoting averages with respect to the system $|\psi_S\rangle$ by $\langle...\rangle_S$ as in $G_{ij} = \langle c_j^\dagger c_i\rangle_S$, we find, for case (a),

$$\langle n_L \rangle^{(a)} = \lambda^2(\sin^2\theta \langle n_i \rangle_S + \cos^2\theta \langle n_j \rangle_S \\ - 2 \sin\theta \cos\theta \operatorname{Im}(G_{ij})),$$

$$\langle n_R \rangle^{(a)} = \lambda^2(\cos^2\theta \langle n_i \rangle_S + \sin\theta \langle n_j \rangle_S \\ + 2 \sin\theta \cos\theta \operatorname{Im}(G_{ij})).$$
Choosing an angle $\theta = \pi/4$ gives the imaginary part of $G_{ij}$,

$$2\lambda^2\mathrm{Im}(G_{ij}) = \langle n_R \rangle^{(a)} - \langle n_L \rangle^{(a)}.$$

For case (b), one obtains

$$\langle n_L \rangle^{(b)} = \lambda^2(\sin^2\theta\langle n_i \rangle_S + \cos^2\theta\langle n_j \rangle_S - 2\sin\theta\cos\theta\mathrm{Re}(G_{ij})),$$

$$\langle n_R \rangle^{(b)} = \lambda^2(\cos^2\theta\langle n_i \rangle_S + \sin\theta\langle n_j \rangle_S + 2\sin\theta\cos\theta\mathrm{Re}(G_{ij})).$$

Choosing the same angle $\theta = \pi/4$ gives the real part of $G_{ij}$,

$$2\lambda^2\mathrm{Re}(G_{ij}) = \langle n_R \rangle^{(b)} - \langle n_L \rangle^{(b)}.$$

Hence, both the real and imaginary parts of the coherence, as well as their statistical distribution, are experimentally measurable (at least in principle).

## APPENDIX B: SINGULAR BEHAVIOR OF NONCONNECTED CORRELATION FUNCTIONS

For the example of $n = 2$, we show here why Eq. (51) produces solutions $g_t^s$ that are singular at coinciding points, whereas solutions $g_t$ of Eq. (54) are continuous (but not differentiable). The corresponding equations for $n = 2$ are Eqs. (67) and (68). Namely,

$$(\partial_t - \Delta)g_t^s(x, y) = 2\partial_x\partial_y(\delta(x - y)g_t(x)g_t(y))$$

and

$$(\partial_t - \Delta)g_t(x, y) = 2\delta(x - y)\partial_x g_t(x)\partial_y g_t(y).$$

Integrating these equations across the diagonal line $\{x = y\}$ reveals that the van Neumann boundary conditions will be singular for $g_t^s$ and regular for $g_t$.

Let us perform this integration explicitly for Eq. (67). We rotate the variables $(x, y)$ by $\pi/4$ clockwise, $(v, u) = [(x - y/\sqrt{2}), (x + y/\sqrt{2})]$. Then, the derivative $\partial_v = (\partial_x - \partial_y/\sqrt{2})$ is orthonormal to the line $\{x = y\}$, and $\partial_u = (\partial_x + \partial_y/\sqrt{2})$ is parallel to the line. In these variables, the equation reads

$$(\partial_t - (\partial_u^2 + \partial_v^2))g_t^s \left( \frac{u+v}{\sqrt{2}}, \frac{u-v}{\sqrt{2}} \right) = (\partial_u^2 - \partial_v^2) \left( \delta(\sqrt{2}v)\rho_t \left( \frac{u+v}{\sqrt{2}} \right) \rho_t \left( \frac{u-v}{\sqrt{2}} \right) \right).$$

Integrating $\int_{-\epsilon}^{\epsilon} dv$ in this equation and keeping only terms of $O(1)$, we obtain

$$\partial_v g_t^s \left( \frac{u+v}{\sqrt{2}}, \frac{u-v}{\sqrt{2}} \right) \Big|_{-\epsilon}^{\epsilon} = \partial_v \left( \delta(\sqrt{2}v) \rho_t \left( \frac{u+v}{\sqrt{2}} \right) \rho_t \left( \frac{u-v}{\sqrt{2}} \right) \right) \Big|_{-\epsilon}^{\epsilon}.$$

The expression simplifies due to cyclic invariance, $g_t^s(x, y) = g_t^s(y, x)$,

$$2\partial_v g_t^s \left( \frac{u+v}{\sqrt{2}}, \frac{u-v}{\sqrt{2}} \right) \Big|_{v=\epsilon} = 2\sqrt{2}\delta'(\sqrt{2}v)\rho_t \left( \frac{u+v}{\sqrt{2}} \right) \rho_t \left( \frac{u-v}{\sqrt{2}} \right) \Big|_{v=\epsilon}.$$

Finally, we take $\epsilon \to 0^+$ and find the van Neumann boundary condition in both sectors $\{x > y\}$ and $\{x < y\}$,

$$\partial_v g_t^s|_{x=y^+} = -\partial_v g_t^s|_{x=y^-} = \sqrt{2}\delta'(0^+)\rho(x)^2.$$

The delta function causes this derivative to blow up; therefore, $g_t^s(x, y)$ is indeed singular at $x = y$. The same procedure for Eq. (68) yields

$$\partial_v g_t|_{x=y^+} = -\partial_v g_t|_{x=y^-} = -\frac{1}{\sqrt{2}}\rho_t'(x)^2, \quad (B1)$$

which is finite.

## APPENDIX C: SOLUTION FOR THE DISCRETE DENSITY

The time-evolution equation for the discrete density (62) can be rewritten as $\partial_t n = An + n$, with $n = (n_1, \dots, n_L)^T$, $b = (\alpha_1, \dots, \alpha_L)^T$, and

$$A = \begin{pmatrix} -1 - (\alpha_1 + \beta_1) & 1 & & & \\ & 1 & -2 & \dots & \\ & & \ddots & \ddots & \ddots & \\ & & & -2 & 1 & \\ & & & & 1 & -1 - (\alpha_L + \beta_L) \end{pmatrix}.$$

If $A$ is diagonalizable, $A = SDS^{-1}$, a general solution for an initial condition $n(0) = u$ is given by

$$\begin{aligned} n(t) &= \exp(At)u + (\exp(At) - 1)A^{-1}b \\ &= S\exp(Dt)(S^{-1}u + D^{-1}S^{-1}b) + A^{-1}b. \end{aligned} \quad (C1)$$

If $\alpha_1 + \beta_1 = 1 = \alpha_L + \beta_L$, $A$ is a so-called Tölpitz matrix and has eigenvalues $\lambda_k$ and normalized eigenvectors $v_k$,
$$ \lambda_k = -2 + 2 \cos(k), \\ (v_k)_j = \sqrt{2/(L+1)} \sin(jk), $$

with $k = (\pi/L + 1), (2\pi/L + 1), \dots, (L\pi/L + 1)$, from which Eq. (63) can be derived. For general injection or extraction rates, we obtain an ansatz

$$ \lambda_k = -2 + 2 \cos(k), \\ (v_k)_j = e^{ijk} + s_k e^{-ijk}, $$

where $s_k$ and $k$ are determined by the eigenvalue equations. Denoting $z := e^{ik}$, this leads to

$$ 0 = (z^{L+1} - z^{-(L+1)}) - (A+B)(z^L - z^{-L}) \\ + AB(z^{L-1} - z^{-(L-1)}), \\ s_k = -e^{2ik} \frac{A - e^{ik}}{A - e^{-ik}}, \qquad (C2) $$

where $A = 1 - (\alpha_1 + \beta_1)$ and $B = 1 - (\alpha_L + \beta_L)$. Note that even though this equation has $2L + 2$ solutions for $z$, only $L$ of them give rise to eigenvalues with linearly independent eigenvectors: If $k$ is a solution, then $-k$ is also a solution, but the corresponding eigenvectors are linearly dependent, $v_{-k} = s_k v_k$. Furthermore, $k = 0$ is a solution with a zero eigenvector. Numerical solutions of Eq. (C2) show that, for $A, B \in [1, -1]$, all solutions of $z$ lie on the unit circle; therefore (taking into account that $k \sim -k$ and $k \neq 0$), we have $k \in (0, \pi)$.

A site in the bulk, $j \sim aL$, approaches the steady-state value $(A^{-1}b)_j$ according to the time-dependent part of Eq. (C1). Again, since $\lambda_k < 0$, there is an exponential decay. Only the smallest solution for $k$, denoted by $k^*$, will contribute. In the limit of large system size $L$, one can check that $k^* \approx \pi/L$ is the smallest positive solution for $k$. Note that the amplitude $(v_{k^*})_j \approx e^{ia\pi} - e^{-ia\pi}$ of this term stays finite ($s_{k^*} \approx -1$). Then, the timescale with which the bulk approaches the steady-state value is, as before, $t_{\text{decay}} \sim \mathcal{O}(L^2)$.

At the boundary, the amplitude of the term corresponding to $k^*$ will be zero. In general, for $j = 1$, $(v_{k^*})_j = e^{ik^*} + s_{k^*}e^{-ik^*} = 0$. To have a nonzero amplitude, one has to consider terms where $k \sim bL$ with $b \sim \mathcal{O}(1)$, which leads to $t_{\text{decay}} \sim \mathcal{O}(1)$.

## APPENDIX D: ANALYTIC SOLUTION OF THE CONNECTED TWO-POINT FUNCTION

Here, we outline how to solve the connected two-point function analytically, which is needed for the comparison with the solution of the discrete equations in the scaling limit in Figs. 4 and 5.

### 1. Solution for the density

First, we construct an analytic solution of the density in Eq. (65) with the domain-wall initial condition $\rho(x, 0) = \Theta(1/2 - x)$ and boundary conditions $\rho(0, t) = n_a$, $\rho(1, t) = n_b$. We simplify the boundary conditions by subtracting the stationary solution, $\rho_\infty(x) = n_a + x(n_b - n_a)$. Then, we solve for $\tilde{\rho}(x, t) = \rho(x, t) - \rho_\infty(x)$, which has simpler boundary conditions:

$$ (\partial_t - \partial_x^2)\tilde{\rho} = 0, \\ \tilde{\rho}(x, 0) = \Theta(1/2 - x) - n_a - x(n_b - n_a), \\ \tilde{\rho}(0, t) = \tilde{\rho}(1, t) = 0. $$

We find the solution by an expansion in $\{\sin(n\pi x)\}_{n=1}^\infty$. These functions satisfy the correct boundary conditions and are orthogonal in the sense that $\int_0^1 \sin(n\pi x) \sin(m\pi x) dx = \frac{1}{2}\delta_{nm}$. Importantly, they form a complete basis of $L^2([0, 1])$, which justifies the expansion.

Taking into account the initial condition, we obtain

$$ \tilde{\rho}(x, t) = \sum_{n=1}^{\infty} c_n \sin(n\pi x) e^{-n^2\pi^2 t}, \\ c_n = \frac{2}{n\pi} (1 - (-1)^n (n_a - n_b) - 2n_a\delta_{n,\text{odd}} \\ - (-1)^{n/2}\delta_{n,\text{even}}). $$

In the special case where the boundary conditions match the initial conditions ($n_a = 1, n_b = 0$), one finds

$$ \rho(x, t) = 1 - x - \sum_{k=1}^{\infty} \frac{(-1)^k}{k\pi} \sin(2\pi kx) e^{-4\pi^2 k^2 t}. $$

### 2. Solution for the connected two-point function

In Appendix B, we showed that the connected two-point function $g_t(x, y)$ satisfies the van Neumann boundary conditions (B1) on the diagonal $\{x = y\}$. Here, we construct a solution of Eq. (68) on the lower triangle $T^+ = \{(x, y) \in [0, 1]^2 : x \ge y\}$: First, we identify a function that satisfies the boundary conditions,

$$ w(x, y, t) = y(1-x)\partial_x\rho(x, t)\partial_y\rho(y, t). $$

Note that for $t \to \infty$, $w$ becomes the correct stationary solution on $T^+$. Then, we solve for $f(x, y, t) := g_t(x, y) - w(x, y, t)$, which satisfies an inhomogeneous heat equation with homogeneous boundary conditions,
$$ (\partial_t - \Delta)f(x, y, t) = S(x, y, t) \\ := 2(1-x)\partial_x\rho(x, t)\partial_y^2\rho(y, t) \\ - 2y\partial_x^2\rho(x, t)\partial_y\rho(y, t), \\ f(x, y, 0) = -y(1-x)\delta(1/2-x)\delta(1/2-y), \\ f(x, 0, t) = 0 = g(1, 0, t) \quad (\text{Dirichlet condition}), \\ \partial_v f|_{x=y} = 0 \quad (\text{Neumann condition}). $$

with $\partial_v := (\partial_x - \partial_y)/\sqrt{2}$ as in Appendix B. This is solved, as before, by the method of eigenfunction expansion. Note that

$$ \psi_{nm}(x, y) := \sin(n\pi x) \sin(m\pi y) + \sin(m\pi x) \sin(n\pi y) $$

is a complete basis of $L^2(T^+)$ that satisfies the correct Dirichlet and Neumann boundary conditions. It satisfies

$$ \int_{T^+} \psi_{nm}(x, y) \psi_{kl}(x, y) dx dy = \frac{\delta_{nk} \delta_{ml} + \delta_{nl} \delta_{mk}}{4}. $$

We write $S(x, y, t) = \sum_{n \ge m \ge 1} \hat{S}_{nm}(t) \psi_{nm}(x, y)$ and $f(x, y, t) = \sum_{n \ge m \ge 1} \hat{f}_{nm}(t) \psi_{nm}(x, y)$, where the coefficients are given by

$$ \hat{S}_{nm}(t) = \begin{cases} 4 \int_{T^+} S \psi_{nm} & \text{if } n > m \\ 2 \int_{T^+} S \psi_{nn} & \text{if } n = m. \end{cases} $$

This leads to $\partial_t \hat{f}_{nm} + \pi^2(n^2 + m^2)\hat{f}_{nm} = \hat{S}_{nm}$, which is solved by

$$ \hat{f}_{nm}(t) = \underbrace{\hat{f}_{nm}(0) e^{-\pi^2(n^2+m^2)t}}_{= \hat{f}_{nm}^{\text{hom}}(t)} + \underbrace{\int_0^t e^{\pi^2(n^2+m^2)(\tau-t)} \hat{S}_{nm}(\tau)}_{= \hat{f}_{nm}^{\text{part}}(t)}. $$

We get a solution $f_{\text{hom}} = \sum_{n \ge m \ge 1} \hat{f}_{nm}^{\text{hom}}(t) \psi_{nm}$ of the homogeneous equation that satisfies the initial condition and a solution $f_{\text{part}} = \sum_{n \ge m \ge 1} \hat{f}_{nm}^{\text{part}}(t) \psi_{nm}$. One finds

$$ \begin{aligned} f_{\text{hom}}(x, y, t) &= \sum_{k,l \ge 0} (-1)^{k+l+1} e^{-\pi^2((2k+1)^2+(2l+1)^2)t} \\ & \quad \times \sin((2k+1)\pi x) \sin((2l+1)\pi y) \\ &= -\frac{1}{4} \vartheta_1(\pi x, e^{-4\pi^2 t}) \vartheta_1(\pi y, e^{-4\pi^2 t}). \end{aligned} $$

The precise expression for $\hat{S}_{nm}$ is rather complicated, and it is easier to solve for $f_{\text{part}}$ using the MATHEMATICA NDSolve function. The complete solution is then $g_t(x, y) = w(x, y, t) + f_{\text{hom}}(x, y, t) + f_{\text{part}}(x, y, t)$. By symmetry in x and y, i.e., $g_t(x, y) = g_t(y, x)$, this also determines a solution on $T^- := \{(x, y) \in [0, 1]^2 : x \ge y\}$.

# APPENDIX E: PROOFS

## 1. Expansion of loop expectation values into noncrossing partitions

The proof of Eq. (31) and its continuous version Eq. (56), which relate the correlation function $\mathbb{E}[G_{i_1 i_2} \cdots G_{i_n i_1}] := [n]$ to its connected part $\mathbb{E}[G_{i_1 i_2} \cdots G_{i_n i_1}]^c := [n]^c$ through a sum over noncrossing partitions (the notation $[n]$ and $[n]^c$ are only used in this subsection), is based on the three conditions (1)–(3) which can be stated more compactly as

(i) The measure $\mathbb{E}$ is $U(1)$ invariant (we will often denote $\mathbb{E}[\cdots] = [\cdots]$, where we dropped the subscript $t$ from $\mathbb{E}_t$).

(ii) The connected-loop expectation value of $n$ points scales as $[n]^c \sim (1/L^{n-1})$, and all other connected correlations with equal numbers are either of the same order or subleading.

Using the moment-cumulant formula (33), the essential point is to show that cumulants

$$ \mathbb{E}_{\pi}[G_{i_1 i_2} \cdots G_{i_n i_1}]^c := \prod_{b \in \pi} \mathbb{E}[G_{i_{b(1)}} i_{b(2)}} \cdots G_{i_{b(|b|)}} i_{b(1)}]^c \quad (\text{E1}) $$

corresponding to crossing partitions $\pi$ will be subleading compared to noncrossing partitions $\pi$ in the scaling limit.

### a. Scaling of noncrossing partitions

We start the proof by showing that if $\pi = \{b^{(1)}, \dots, b^{(m)}\}$ is a noncrossing partition of $\{12, 23, \dots, n1\}$, i.e., of the edges of a loop with $n$ nodes, then $[\pi]^c := \mathbb{E}_{\pi}[G_{i_1 i_2} \cdots G_{i_n i_1}]^c$ scales as $L^{-n+1}$ independently of the number of blocks $m$; therefore, it behaves in the same way as the single-loop expectation value $\mathbb{E}[G_{i_1 i_2} \cdots G_{i_n i_1}]^c$ with $n$ points.

First notice that if the blocks of $\pi$ are not nested into each other (e.g., $\pi = \{\{12, 23\}, \{34\}, \{45, 56, 61\}\}$), then, by $U(1)$ invariance, we have to connect starting and end points of each block by a Kronecker $\delta$. Once we arrive at the last block, this condition is already satisfied because of the other delta functions. We therefore need $m-1$ delta functions,

$$ \begin{aligned} [\pi]^c &= \underbrace{\delta \cdots \delta [b^{(1)}]^c \cdots [b^{(m)}]^c}_{m-1} \sim L^{-m+1} L^{-\sum_{i=1}^m |b^{(i)}|+m} \\ & \sim L^{-n+1}, \end{aligned} $$

where $\delta \sim L^{-1}$ in the scaling limit and the sum over the size $|b^{(i)}|$ of each block is equal to the total number of elements, $n$.

Next, assume that $\pi$ has some nested blocks (e.g., $\pi = \{\{12, 23\}, \{45\}, \{34, 56, 61\}\}$); then, treat each collection of nested blocks as a big block $B$ (e.g., $B = \{\{45\}, \{34, 56, 61\}\}$), such that the argument above applies to the non-nested blocks and the big blocks. Now, we can iterate the argument for each big block $B$ and for
possible collections of nested sub-blocks therein. In the end, $m-1$ delta functions are needed to close each block in $\pi$ to form a loop, which shows that $[\pi]^c \sim L^{-n+1}$ if $\pi$ is noncrossing.

## b. Scaling of crossing partitions

Now, assume that $\pi$ consists of a collection $B$ of noncrossing blocks with $|B|$ elements in total and a collection $C$ of crossing blocks that cannot be disentangled from one another with $|C|$ elements in total (e.g., $\pi = \{\{12, 34\}, \{23, 45\}, \{56\}, \{67\}\}, B = \{\{56\}, \{67\}\},$ and $C = \{\{12, 34\}, \{23, 45\}\}$). Then, we have $|B| + |C| = n$. If we treat $B$ as an independent partition of the cyclic set formed by all of its elements, then by the argument above, $[B]^c \sim L^{-|B|+1}$. Alternatively, we can view $B$ as the partition of the full loop from which we removed all edges that belong to the collection $C$ of crossing blocks, leaving us with a smaller loop of $|B|$ edges. If we assume that $[C]^c \sim L^{-|C|}$ (as will be shown below), then

$$[\pi]^c = \delta[B]^c [C]^c \sim L^{-1} L^{-|B|+1} L^{-|C|} \sim L^{-n},$$

where only a single delta function is needed to connect the collections $B$ and $C$ (because all the other delta functions are already included in $[B]^c$ and $[C]^c$). If $\pi$ has more than one collection $C$ of crossing blocks that cannot be disentangled, then it will scale with an even higher negative power of $L$. This shows that crossing partitions are subleading in the scaling limit compared to noncrossing partitions.

It remains to show that $[C]^c \sim L^{-|C|}$ (and probably even $[C]^c \sim L^{-|C|-1}$). The idea is to produce a collection of crossing blocks starting with noncrossing blocks and permuting its elements. Assume that $\{b^{(1)}, b^{(2)}\}$ is a collection of two noncrossing blocks

$$b^{(1)} = \{i_1 i_2, \dots, i_k i_{k+1}\}, \quad b^{(2)} = \{i_{k+1} i_{k+2}, \dots, i_n i_1\},$$

which would mean that there is only one delta function $\delta(i_l, i_{k+1})$ necessary for the product of these blocks to be nonzero. Then, construct a first crossing by inserting an element $i_l i_{l+1}$ from the “middle” of $b^{(1)}$ (and not from the boundary) into $b^{(2)}$,

$$b^{(1)} \rightarrow b^{(1)} = \{i_1 i_2, \dots, i_{l-1} i_l, i_{l+1} i_{l+2}, \dots, i_k i_{k+1}\},$$

$$b^{(2)} \rightarrow b^{(2)} = \{i_l i_{l+1}, i_{k+1} i_{k+2}, \dots, i_n i_1\},$$

which makes it necessary to insert a second delta function, $\delta(i_l, i_{l+1})$. In this way, one can continue to permute elements between the two blocks, create more crossings, and obtain $C = \{b^{(1)}, b^{(2)}\}$ of sizes $|b^{(1)}| = k'$ and $|b^{(2)}| = n - k'$. Each new crossing necessitates a new delta function. Note, however, that a new crossing is only created if one permutes an element (a) that is still in its

original block, (b) whose neighbors have not yet been permuted (we sort the elements $i_l i_{l+1}$ in ascending order with respect to the index $l$), and (c) that is not taken from the boundary of the block. Since no connected correlation function is more dominant than that of single loops, the two blocks, after an arbitrary permutation between them, will scale at most as $[b^{(1)}] \sim L^{-k'+1}$ and $[b^{(2)}] \sim L^{-(n-k')+1}$. We therefore have

$$[C]^c = \underbrace{\delta \cdots \delta}_{\#\text{crossings}+1} [b^{(1)}]^c [b^{(2)}]^c \sim L^{-n+1-\#\text{crossings}}.$$

As a consequence, in the case where $C$ consists of two crossing blocks, this shows that it scales at most with $[C]^c \sim L^{-|C|}$, where $|C| = n$ in our example.

If $C$ is a collection of $m$ crossing blocks, it is probably still true that the number of delta functions is equal to $\#\text{crossings} + 1$. However, here it is only necessary that $m$ crossing blocks cause at least $m$ delta functions. Then, the scaling is

$$[C]^c = \underbrace{\delta \cdots \delta}_{m} [b^{(1)}]^c \cdot [b^{(m)}]^c \sim L^{-|C|},$$

which is sufficient. Note that two crossing blocks cause at least two delta functions. However, each originally noncrossing block that we add to the collection of crossing blocks already comes with one delta function, even before we permute its elements with the other blocks. For example, if we have the crossing blocks $\{\{12, 34\}, \{23, 41\}\}$, which need two delta functions $\delta(2, 3)\delta(1, 4)$ and we add a block, then we have $\{\{12, 34\}, \{23, 45\}, \{56, 61\}\}$, which needs three delta functions, $\delta(2, 3)\delta(1, 4)\delta(4, 5)$. If we now cross the new block with the others, this cannot reduce the number of delta functions. Therefore, this argument shows that $m$ crossing blocks indeed cause $m$ delta functions and concludes the proof.

## 2. Time evolution of connected-loop expectation values

We show how Eq. (54) follows from Eq. (51) by induction over $n$. We start by writing the sum $2 \sum_{i<j} = \sum_{ij, j \neq i}$ in these equations as a sum over all ordered, noncrossing, and nonempty sets $r = \{r_1, r_2, \dots\}$ and $s = \{s_1, s_2, \dots\}$ with $i = r_1$ and $j = s_1$ such that $r \sqcup s = [n] = \{x_1, \dots, x_n\}$. Throughout this paper, we use the symbol $\sqcup$ to denote the union of ordered, noncrossing, and nonempty subsets. Furthermore, instead of $g_t(x_i, x_{i+1}, \dots, x_{j-1})$, we simply write $g(r)$ if $r = \{x_i, x_{i+1}, \dots, x_{j-1}\}$. Then, the two equations become
$$ \sum_{\pi \in NC([n]) \setminus [n]} \left( -\sum_{\substack{(x,y) \in d \\ d \in \pi^*}} \partial_x^g \partial_y^g D_\pi + \sum_{\substack{(x,y) \in b \\ b \in \pi}} \partial_x^g \partial_y^g D_{\pi \setminus b \cup b_1(x,y) \cup b_2(x,y)} \right) \\ = - \sum_{r \sqcup s = [n]} \partial_{r_1}^g \partial_{s_1}^g D_{\{r,s\}}. \qquad (E10) $$

Since $\pi \neq \{[n]\}$, one can always find a block $d \in \pi^*$ that consists of at least two nodes $(x, y)$. We can join the two blocks $b(x) \in \pi$ and $b(y) \in \pi$ to which the corresponding edges $x$ and $y$ belong. This forms a block $b$ of a new partition $\pi'$ that differs from $\pi$ only by this block, $\pi' \setminus b \cup b(x) \cup b(y) = \pi$. Note that if $\pi$ consists of only two blocks, then $\pi'$ will consist of a single block, so $\pi' = \{[n]\}$. Hence, we obtain a bijection

$$ \begin{gathered} \{(D_\pi, x, y)|\pi \in NC([n]) \setminus [n], (x, y) \in \text{dst.} d \in \pi^*\} \\ \leftrightarrow \\ \{\pi' \setminus b \cup b_1(x, y) \cup b_2(x, y)|\pi' \in NC([n]), (x, y) \in b \in \pi'\}. \end{gathered} $$

However, $\pi' = [n]$ is not available in the above sum; therefore, the terms

$$ - \sum_{\substack{(x,y) \in d \\ d \in \pi^*}} \partial_x^g \partial_y^g D_\pi $$

are uncanceled, where $\pi$ consists of only two blocks. Then, we can write $\pi = \{s, r\}$, where $r \sqcup s = [n]$, which establishes the above equality.

Equation (E8) then simplifies to

$$ \begin{aligned} (\partial_t - \Delta)g_t([n]) &= \sum_{r \sqcup s = [n]} \partial_{r_1}^g \partial_{s_1}^g D_{\{r,s\}} \\ &= \sum_{r \sqcup s = [n]} \delta(r_1, s_1) \partial_{r_1} g_t(r) \partial_{s_1} g_t(s), \end{aligned} \qquad (E11) $$

which is what we wanted to show.

*Proof of Eq. (E9).* Recall that our notation $(x, y) \in b$, with $b$ any block of a partition, assumes that $x \neq y$. Two identities for the $\delta$ functions of several variables that we will need are

$$ \sum_{x \in [n]} \partial_x \delta([n]) = 0. \qquad (E12) $$

$$ \Delta_{[n]} \delta([n]) = - \sum_{\substack{(x,y) \in [n] \\ x \neq y}} \partial_x \delta_y \delta([n]). \qquad (E13) $$

The subscript of the Laplacian denotes the set on which it acts. For example, $\Delta_b = \sum_{x \in b} \partial_x^2$ for some set $b \subset [n]$. For $\pi \in NC([n]) \setminus [n]$, one finds

$$ (\partial_t - \Delta) D_\pi = - \underbrace{\left( \Delta \prod_{d \in \pi^*} \delta(d) \right) \prod_{b \in \pi} g_t(b) - 2 \sum_{\substack{x \in [n] \\ x \neq y}} \left( \partial_x \prod_{d \in \pi^*} \delta(d) \right) \left( \partial_x \prod_{b \in \pi} g_t(b) \right)}_{(I)} + \underbrace{\prod_{d \in \pi^*} \delta(d) (\partial_t - \Delta) \prod_{b \in \pi} g_t(b)}_{(II)} $$

Using $\sum_{x \in [n]} = \sum_{e \in \pi^*} \sum_{x \in e}$,

$$ (I) = - \sum_{e \in \pi^*} \prod_{d \in \pi^* \setminus e} \delta(d) \left( (\Delta_e \delta(e)) \prod_{b \in \pi} g_t(b) + 2 \sum_{x \in e} \partial_x \delta(e) \partial_x \prod_{b \in \pi} g_t(b) \right). $$

With the help of the $\delta$-function identities (E12) and (E13), we obtain

$$ (I) = \sum_{e \in \pi^*} \prod_{d \in \pi^* \setminus e} \delta(d) \left( \sum_{(x,y) \in e} (\partial_x \partial_y \delta(e)) \prod_{b \in \pi} g_t(b) + 2 \sum_{(x,y) \in e} \partial_y \delta(e) \partial_x \prod_{b \in \pi} g_t(b) \right). $$

This expression can be written as a total derivative of $\partial_x \partial_y$ up to the missing term $\partial_x^g \partial_y^g$,

$$ (I) = \sum_{e \in \pi^*} \prod_{d \in \pi^* \setminus e} \delta(d) \sum_{(x,y) \in e} \left( \partial_x \partial_y \left( \delta(e) \prod_{b \in \pi} g_t(b) \right) - \delta(e) \partial_x \partial_y \prod_{b \in \pi} g_t(b) \right). $$

which, rewritten in terms of $D_\pi$, is the first term that appears in Eq. (E9),

$$ (I) = \sum_{e \in \pi^*} \sum_{(x,y) \in e} (\partial_x \partial_y - \partial_x^g \partial_y^g) D_\pi. $$

For the second term, we need to use Eq. (E3), which by assumption holds for all $k \le n - 1$. Since $\pi \in NC([n]) \setminus [n]$, the assumption applies,
$$ (II) = \prod_{d \in \pi^+} \delta(d) \sum_{c \in \pi} (\partial_t - \Delta_c) g_t(c) \prod_{b \in \pi \setminus c} g_t(b) \\ = \prod_{d \in \pi^+} \delta(d) \sum_{c \in \pi} \sum_{r \sqcup s = c} \delta(r_1, s_1) \partial_{r_1} g_t(r) \partial_{s_1} g_t(s) \prod_{b \in \pi \setminus c} g_t(b). $$

Instead of summing over $r \sqcup s = c$, we can also sum over $(x, y) \in c$ (where, as usual, $x \neq y$). In this case, $r$ and $s$ are the blocks $b_1(x, y)$ and $b_2(x, y)$ that result from cutting $c$ along the nodes $(x, y)$,

$$ (II) = \prod_{d \in \pi^+} \delta(d) \sum_{c \in \pi} \sum_{(x,y) \in c} \delta(x,y) \partial_x g_t(b_1(x,y)) \partial_y g_t(b_2(x,y)) \\ \times \prod_{b \in \pi \setminus c} g_t(b). $$

The expression can therefore be rewritten in terms of the partition $\pi \setminus c \cup b_1(x, y) \cup b_2(x, y)$, and it provides the second term in Eq. (E9),

$$ (II) = \sum_{c \in \pi} \sum_{(x,y) \in c} \partial_x^g \partial_y^g D_{\pi \setminus c \cup b_1(x,y) \cup b_2(x,y)}. $$

This concludes the derivation of the time-evolution equation of connected correlation functions.

### a. Time evolution of the new measure

Here, we present the proof of the evolution equation (74) of the new measure $\varphi_t$. We follow the explanation before Eq. (E3) to rewrite the evolution of the connected correlations $g_t$ as

$$ (\partial_t - \Delta)g_t([n]) = \sum_{r \sqcup s = [n]} \delta(r_1, s_1) \partial_{r_1} g_t(r) \partial_{s_1} g_t(s), \quad (E14) $$

where $r \sqcup s = [n]$ denotes the union of noncrossing subsets $r = \{r_1, r_2, \dots\}$ and $s = \{s_1, s_2, \dots\}$ of the set $[n] = \{x_1, \dots, x_n\}$. The proof of Eq. (74) is, by induction, over $n$. For $n=1$, $\varphi_t$ and $g_t$ are identical and therefore satisfy the same equation. Assume the formula holds for all $n \le k-1$ for some $k \in \mathbb{N}$.

We start by evaluating its left-hand side making use of Eq. (E14) (for better readability, we suppress the time argument),

$$ \begin{aligned} (\partial_t - \Delta)\varphi([n]) &= \sum_{\pi \in NC([n])} \sum_{c \in \pi} (\partial_t - \Delta_c)g(c)g_{\pi \setminus c} \\ &= \sum_{\pi \in NC([n])} \sum_{c \in \pi} \sum_{a \cup a' = c} \delta(a_1, a'_1) \partial_{a_1} g(a) \\ & \qquad \times \partial_{a'_1} g(a') g_{\pi \setminus c}. \end{aligned} \quad (E15) $$

Here, $\Delta_c = \sum_{x \in c} \Delta_x$, and $\pi \setminus c$ is the partition $\pi$ without the block $c$.

We continue with the right-hand side of Eq. (E14),

$$ \begin{aligned} & \sum_{r \sqcup s = [n]} \delta(r_1, s_1) \partial_{r_1} \varphi(r) \partial_{s_1} \varphi(s) \\ &= \sum_{r \sqcup s = [n]} \delta(r_1, s_1) \left( \partial_{r_1} \sum_{\rho \in NC(r)} g_{\rho} \right) \left( \partial_{s_1} \sum_{\sigma \in NC(s)} g_{\sigma} \right) \\ &= \sum_{\substack{r \sqcup s = [n] \\ \rho \in NC(r) \\ \sigma \in NC(s)}} \delta(r_1, s_1) \partial_{r_1} g(b(r_1)) \partial_{s_1} g(b(s_1)) g_{\rho \setminus b(r_1)} g_{\sigma \setminus b(s_1)}, \end{aligned} \quad (E16) $$

where the blocks $b(r_1) \in \rho$ and $b(s_1) \in \sigma$ are uniquely defined by the fact that $r_1 \in b(r_1)$ and $s_1 \in b(s_1)$.

Note that terms appearing in the sum in Eqs. (E15) and (E16) agree iff $a = b(r_1)$, $a' = b(s_1)$ and $\pi \setminus c = \rho \setminus b(r_1) \cup \sigma \setminus b(s_1)$. To show the equivalence between the two sums, we note that there is a bijection between the sets,

$$ \begin{aligned} & \{(a, a', \pi) | \pi \in NC([n]), a \cup a' =: c \in \pi\} \\ & \leftrightarrow \\ & \{(r_1, s_1, \rho, \sigma) | \rho \in NC(r), \sigma \in NC(s), r \sqcup s = [n]\}. \end{aligned} $$

For “→,” given $\pi \in NC([n])$, $a \cup a' =: c \in \pi$, this completely specifies terms in the sum in Eq. (E14). Define $r_1 = a_1$ and $s_1 = a'_1$. By the condition that $r \sqcup s = [n]$ is noncrossing, this uniquely defines the sets $r$ and $s$. Since $\pi$ is noncrossing, once we take away the “connecting block” $c$, $\pi \setminus c$ factorizes in a unique way into two noncrossing partitions $\tilde{\rho}$ of $r \setminus a$ and $\tilde{\sigma}$ of $s \setminus a'$. In other words, $\pi \setminus c = \tilde{\rho} \cup \tilde{\sigma}$. Therefore, we can define $\rho = \tilde{\rho} \cup \{a\}$ and $\sigma = \tilde{\sigma} \cup \{a'\}$ such that, indeed, $\pi \setminus c = \rho \setminus a \cup \sigma \setminus a'$. This process produces a corresponding term in the sum of Eq. (E16), which is completely specified by the data $(r_1, s_1, \rho, \sigma)$.

For “←,” given $\rho \in NC(r)$, $\sigma \in NC(s)$, $r \sqcup s = [n]$, we define $a = b(r_1)$, $a' = b(s_1)$, and $c = a \cup a'$. It remains to construct $\pi \in NC([n])$ such that $\pi \setminus c = \rho \setminus a \cup \sigma \setminus a'$. This is achieved by defining $\pi := \rho \cup \sigma \cup \{c\}$.

### 3. Steady-state solution

Here, we prove that $\varphi_{\infty}(x_1, \dots, x_n) = \min(x_1, \dots, x_n)$ is a steady-state solution of Eq. (74).

We represent the minimum by a sum of Heaviside functions $\Theta(\text{condition})$, which equal 1 if the condition is true and zero otherwise,

$$ \min(x_1, \dots, x_n) = \sum_{i=1}^{n} x_i \Theta(x_i < \{x_1, \dots, \hat{x}_i, \dots, x_n\}). $$

The hat on $\hat{x}_i$ suggests that $x_i$ is missing from the set $\{x_1, \dots, x_n\}$. We have $\partial_i \min(x_1, \dots, x_n) = \Theta(x_i < \{x_1, \dots, \hat{x}_i, \dots, x_n\})$. Furthermore, derivatives of the Heaviside function evaluate to
$$ \partial_i \Theta(x_i < \{x_1, \dots, \hat{x}_i, \dots, x_n\}) = -\sum_{j \neq i} \delta(x_i, x_j) \Theta(x_i < \{x_1, \dots, \hat{x}_i, \hat{x}_j, \dots, x_n\}), $$

$$ \partial_j \Theta(x_i < \{x_1, \dots, \hat{x}_i, \dots, x_n\}) = \delta(x_i, x_j) \Theta(x_j < \{x_1, \dots, \hat{x}_i, \hat{x}_j, \dots, x_n\}). $$

With these formulas, it is easy to check that

$$ - \Delta \min(x_1, \dots, x_n) = \sum_{i,j \atop i \neq j} \delta(x_i, x_j) \Theta(x_j < \{x_1, \dots, \hat{x}_i, \hat{x}_j, \dots, x_n\}). $$

Furthermore, we have

$$ \delta(x_i, x_j) \Theta(x_j < \{x_1, \dots, \hat{x}_i, \hat{x}_j, \dots, x_n\}) = \delta(x_i, x_j) \Theta(x_i < \{x_{i+1}, \dots, \hat{x}_j\}) \Theta(x_j < \{x_{j+1}, \dots, \hat{x}_{i-1}\}). $$

As a consequence,

$$ - \Delta \min(x_1, \dots, x_n) = \sum_{i,j \atop i \neq j} \delta(x_i, x_j) \partial_i \min(x_i, \dots, x_{j-1}) \partial_j \min(x_j, \dots, x_{i-1}). $$

This is all we needed to prove our claim.

[1] H. Touchette, *The Large Deviation Approach to Statistical Mechanics*, Phys. Rep. **478**, 1 (2009).

[2] J. R. Dorfman, T. R. Kirkpatrick, and J. V. Sengers, *Generic Long-Range Correlations in Molecular Fluids*, Annu. Rev. Phys. Chem. **45**, 213 (1994).

[3] C. Jarzynski, *Nonequilibrium Equality for Free Energy Differences*, Phys. Rev. Lett. **78**, 2690 (1997).

[4] G. E. Crooks, *Entropy Production Fluctuation Theorem and the Nonequilibrium Work Relation for Free Energy Differences*, Phys. Rev. E **60**, 2721 (1999).

[5] B. Derrida, *Non-equilibrium Steady States: Fluctuations and Large Deviations of the Density and of the Current*, J. Stat. Mech. **07023** (2007).

[6] 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).

[7] L. Bertini, A. D. Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, *Macroscopic Fluctuation Theory*, Rev. Mod. Phys. **87**, 593 (2015).

[8] B. Derrida, M. R. Evans, V. Hakim, and V. Pasquier, *Exact Solution of a 1D Asymmetric Exclusion Model Using a Matrix Formulation*, J. Phys. A **26**, 1493 (1993).

[9] B. Derrida, *An Exactly Soluble Non-equilibrium System: The Asymmetric Simple Exclusion Process*, Phys. Rep. **301**, 65 (1998).

[10] B. Derrida and M. R. Evans, *Bethe Ansatz Solution for a Defect Particle in the Asymmetric Exclusion Process*, J. Phys. A **32**, 4833 (1999).

[11] T. Bodineau and B. Derrida, *Current Fluctuations in Non-equilibrium Diffusive Systems: An Additivity Principle*, Phys. Rev. Lett. **92**, 180601 (2004).

[12] K. Mallick, *The Exclusion Process: A Paradigm for Non-equilibrium Behaviour*, Physica A (Amsterdam) **418**, 17 (2015).

[13] D. Bernard, *Can the Macroscopic Fluctuation Theory Be Quantized?*, J. Phys. A **54**, 433001 (2021).

[14] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, *Quantum Entanglement Growth under Random Unitary Dynamics*, Phys. Rev. X **7**, 031016 (2017).

[15] A. Nahum, S. Vijay, and J. Haah, *Operator Spreading in Random Unitary Circuits*, Phys. Rev. X **8**, 021014 (2018).

[16] M. P. A. Fisher, V. Khemani, A. Nahum, and S. Vijay, *Random Quantum Circuits* (to be published).

[17] A. C. Potter and R. Vasseur, *Entanglement Dynamics in Hybrid Quantum Circuits*, in *Entanglement in Spin Chains. Quantum Science and Technology*, edited by A. Bayat, S. Bose, and H. Johannesson (Springer, Cham, 2022), 10.1007/978-3-031-03998-0_9.

[18] M. J. Gullans and D. A. Huse, *Entanglement Structure of Current-Driven Diffusive Fermion Systems*, Phys. Rev. X **9**, 021007 (2019).

[19] A. H. Steinbach, J. M. Martinis, and M. H. Devoret, *Observation of Hot-Electron Shot Noise in a Metallic Resistor*, Phys. Rev. Lett. **76**, 3806 (1996).

[20] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, *Transport in Out-of-Equilibrium XXZ Chains: Exact Profiles of Charges and Currents*, Phys. Rev. Lett. **117**, 207201 (2016).

[21] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, *Emergent Hydrodynamics in Integrable Quantum Systems Out of Equilibrium*, Phys. Rev. X **6**, 041065 (2016).

[22] B. Doyon, *Lecture Notes on Generalised Hydrodynamics*, SciPost Phys. Lect. Notes 18 (2020), https://scipost.org/series/collection/2018-08_integrability/.

[23] V. Alba, B. Bertini, M. Fagotti, L. Piroli, and P. Ruggiero, *Generalized-Hydrodynamic Approach to Inhomogeneous Quenches: Correlations, Entanglement and Quantum Effects*, J. Stat. Mech. (2021) 114004.

[24] P. Ruggiero, P. Calabrese, B. Doyon, and J. Dubail, *Quantum Generalized Hydrodynamics*, Phys. Rev. Lett. **124**, 140603 (2020).

[25] M. Bauer, D. Bernard, and T. Jin, *Stochastic Dissipative Quantum Spin Chains (I): Quantum Fluctuating Discrete Hydrodynamics*, SciPost Phys. **3**, 033 (2017).

[26] M. Bauer, D. Bernard, and T. Jin, *Equilibrium Fluctuations in Maximally Noisy Extended Quantum Systems*, SciPost Phys. **6**, 45 (2019).

[27] D. Bernard and T. Jin, *Open Quantum Symmetric Simple Exclusion Process*, Phys. Rev. Lett. **123**, 080601 (2019).

[28] D. Bernard and T. Jin, *Solution to the Quantum Symmetric Simple Exclusion Process: The Continuous Case*, Commun. Math. Phys. **384**, 1141 (2021).

[29] D. Bernard and L. Piroli, *Entanglement Distribution in the Quantum Symmetric Simple Exclusion Process*, Phys. Rev. E **104**, 014146 (2021).

[30] D. Bernard, F. H. L. Essler, L. Hruza, and M. Medenjak, *Dynamics of Fluctuations in Quantum Simple Exclusion Processes*, SciPost Phys. **12**, 42 (2022).
[31] D. V. Voiculescu, *Free Probability Theory*, Vol. 12 (American Mathematical Society, 1997), 10.1090/fic/012.

[32] P. Biane, *Combinatorics of the Quantum Symmetric Simple Exclusion Process, Associahedra and Free Cumulants*, arXiv:2111.12403.

[33] D. Voiculescu, *Limit Laws for Random Matrices and Free Products*, Inventiones Math. **104**, 201 (1991).

[34] A. Guionnet, *Uses of Free Probability in Random Matrix Theory*, Proceedings of the International Conference on Random Matrices in Physics and Math. **9** (2009), 10.1142/9789814304634_0005.

[35] J. A. Mingo and R. Speicher, *Free Probability and Random Matrices*, Vol. 35 (Springer, New York, 2017).

[36] S. Pappalardi, L. Foini, and J. Kurchan, *Eigenstate Thermalization Hypothesis and Free Probability*, Phys. Rev. Lett. **129**, 170603 (2022).

[37] M. Bauer, D. Bernard, and T. Jin, *Universal Fluctuations around Typicality for Quantum Ergodic Systems*, Phys. Rev. E **101**, 012115 (2020).

[38] T. Jin, *Equilibration of Quantum Cat States*, SciPost Phys. **9**, 004 (2020).

[39] L. Foini and J. Kurchan, *Eigenstate Thermalization Hypothesis and Out of Time Order Correlators*, Phys. Rev. E **99**, 042139 (2019).

[40] It is a bit unsatisfying that *U* depends on the indices of the specific average we consider (here, it would be $E_i[G_{ij}]$). Intuitively, a general *u* should be a unitary band matrix with bandwidth of the order of *ℓ*. However, since band matrices do not form a group (the bandwidth grows under multiplication), one cannot take Haar averages on

them, and the construction of the noise would not be well defined.

[41] There is a small subtlety: The first term consists of $\ell$ terms of the form $G_{ii'}^2$ when $i' = j'$. These also scale as $O(1)$, but they are suppressed by the $1/\ell^2$ factor. Their total contribution is $1/\ell$, which goes to zero.

[42] R. Speicher, *Lecture Notes on "Free Probability Theory"* (2019).

[43] J. Novak and P. Sniady, *What Is a Free Cumulant?*, Not. Am. Math. Soc. **58**, 300 (2009), https://scholar.google.com/citations?view_op view_citation&hl=id&user=HHE8yGwAAAAAJ&citation_for_view=HHE8yGwAAAAAJ:2osOgNQ5qMEC.

[44] P. Biane, *Free Probability for Probabilists*, arXiv:math/9809193.

[45] The model can be defined on any graph, but we are dealing with Q-SSEP on a line.

[46] M. Bauer, D. Bernard, P. Biane, and L. Hruza, *Bernoulli Variables, Classical Exclusion Processes and Free Probability*, arXiv:2211.01710.

[47] To be precise, if one is looking at the mean behavior, then the system approaches a mean steady state, but if one is looking at fluctuations, the system is then approaching a steady distribution of states.

[48] Although we expect this result to be known, we did not find any paper on the classical SSEP where this is shown explicitly.

[49] T. Jin, J. a. S. Ferreira, M. Filippone, and T. Giamarchi, *Exact Description of Quantum Stochastic Models as Quantum Resistors*, Phys. Rev. Res. **4**, 013109 (2022).