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

Quantum coherences characterise 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 neighbouring 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 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.

From the point of view of statistical physics, a system is completely characterised 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 is 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 and proportional to the appropriate thermodynamic potential (e.g. the entropy or the free energy), see e.g. [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 arXiv:2204.11680v4 [cond-mat.stat-mech] 31 Mar 2023 variables in out-of-equilibrium systems, which show rich and new features. For example, density correlations in non-equilibrium steady states extend to macroscopic distances, as has been experimentally observed by Dorfmann [2]. Furthermore, fluctuations relations that go beyond the linear response regime, have been shown to be generically applicable to out-of-equilibrium systems. [3,4]. And the large deviation principle introduced above in the equilibrium context has been realised to hold also to the non-equilibrium setting -with the need to formulate appropriate out-of-equilibrium thermodynamic potentials, see e.g. [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-ofequilibrium by coupling the boundaries to reservoirs at different chemical potentials. MFT allows 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, that 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 models are exactly solvable [8][9][10][11] (for a review, see [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 = Tr(ρc † i c j ), the off-diagonal elements of the density matrix ρ when expressed in the particle number basis on each site i, j = 1, · · · , L. A big difference to 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 from partial traces of powers of ρ. Fluctuations of the entanglement entropy have already been studied in the context of unitary random circuits [14,15] (see [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. I) 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 FIG. 1. A schematic representation of a mesoscopic system coupled to two reservoirs of different chemical potentials µa and µ b . Here denotes the ballistic length, above which the transport is diffusive, and L φ the coherence length, below which interference effects can be observed.
of the order of the coherence length L φ of the dynamical degrees of freedom (e.g. for electrons in a disordered metal of size L φ ≈ 1µm [19]). Another important length scale is the mean free path , 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 , 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 L (we identified the physical length with the number of lattice sites by setting the lattice constant to one, a uv = 1). At length above we expect diffusive transport, and at length below L φ we observe quantum mechanical interference. This is the mesoscopic regime we are interested in.
The decomposition of a many-body quantum systems into thermodynamically large cells finds many examples in the literature. Notably, it appears as the main feature of Generalized Hydrodynamics (GHD) which has allowed to study out-of-equilibrium dynamics in integrable systems [20,21] (see [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 generalised Gibbs ensemble, it looses 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 [24]. Furthermore, GHD is mostly concerned with ballistic transport due to the integrable nature of the systems to which it applies. It therefore does not seem to be a good candidate for a quantum version of MFT.
A toy model that allows to study fluctuations in the mesoscopic regime is the quantum symmetric simple exclusion process (Q-SSEP) [25][26][27][28][29][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 non-equilibrium steady state [27], although they are sub-leading in the system size. It is in this sense that Q-SSEP describes mesoscopic trans-port. Furthermore, fluctuations of coherences satisfy a large deviation principle in analogy to macroscopic variables in MFT. This can be seen as a 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 non-commutative random variables pioneered by Dan Voiculescu in the 90's [31]. More precisely, connected fluctuations of coherences in Q-SSEP correspond to free cumulants which generalise 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 sums over non-crossing partitions. This is a structure of which we make repeated use in this paper: Firstly, 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 [30]). And secondly, it also allows 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 mathematician Philippe Biane [32]. Here we generalise his result to all times 0 < t < ∞ 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 characterise 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 realisation of its concepts through large random matrices [33][34][35]. For a generic mesoscopic system, the coarse-grained view of Fig. I in terms of ballistic cells, which undergo a very rapid (unitary) evolution -much faster than the system's evolution as a whole, can be modelled by large random matrices whose size scales with and whose distribution is invariant under the rapid evolution of the ballistic cell, i.e. by the some unitary group (for the moment we leave open the question which degrees of freedom inside the ballistic cells are modelled precisely). In so-far it is not surprising that fluctuations on mesoscopic scales, i.e. between ballistic cells, show signs of free probability. We believe that its 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 section III we give a small introduction to the subject of free probability and in particular explain the relevance of crossing and non-crossing partitions. Section IV deals with the toy model Q-SSEP. In 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, sections IV E explains how to identify free cumulants within Q-SSEP. And IV F explain how the relation to free cumulants can be used to find a simple steady state solution. Details and proofs are given in a few appendices. In particular, 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 was observed that the mathematical structure of the toy model Q-SSEP has a relation to free probability. Then we realised that one can reduce this relation to three simple conditions which apply to generic mesoscopic systems.
Before stating them some explanation is necessary. In the spirit of Fig. I we need to assume that the system exhibits a separation of time scales responsible for fast ballistic transport on scales of the mean free path (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. Mesoscopic observables become therefore random variables with a (yet to be defined) expectation value E t . If such observables do not depend on long range coherences, such as the local particle numbern i := c † i c i (c † i is a fermionic creation operator on site i = 1, · · · , L) or the associated particle current, then we assume that they are well described by classical MFT, with the re- In contrast to this, we are interested in the statistical properties of coherences G ij (t) = Tr(ρ t c † i c j ), a purely quantum mechanical property without classical analogue -and therefore outside the scope of classical MFT. In appendix A we outline a procedure how one could in principle experimentally measure G ij , following an idea in [18].
The three conditions sufficient for the link with free probability concern the statistical properties of coher-ences w.r.t. to the noise expectation value E t in the limit 1 L. Here we adopt the view that G ij is a random variable with time dependent probability distribution, hence the subscript t.
1. Local U (1)-invariance: The expectation value is invariant under G ij → e −iθi G ij e iθj which is a multiplication with local phases. In other words, unless {i 1 , . . . , i n } is a permutation of {j n , . . . , i n } we have 2. Expectation values of "loops" with distinct indices, 3. Expectation values of products of "loops" factorise at leading order. This holds in particular if i 1 = j 1 .
We speak about "loops" because we can associate diagrams to the expectation values of the various products of G ij by connecting nodes i and j by a directed edge. The diagram associated to (2) is indeed a loop, In fact, the local U (1) invariance implies that a diagram built in this way is non-zero 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 factorisation (3) also implies that connected expectation values scale as even if indices become equal. If for example i 1 = i k and no other indices coincide, 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 . What remains, must then scale at least one order of magnitude lower in L. This explains (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 (2) scales as most with 1/L for any choice of i and j. This illustrates nicely that a loops 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 E t and shows how to derive the three conditions from this definition. The consecutive subsection 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 ρ. 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 ∂ t n + ∂ x j = 0. For the following discussion it is not important if the system is open or closed.
The separation of time scales between fast and slow degrees of freedom ensures that ballistic cells don't exchange particles with each other during times smaller than the typical time scale t of ballistic cells. In other words, for t < t , each ballistic cell evolves as a "closed" system and its dynamics is described by a unitary and particle preserving transformations, ρ → U (i) ρ U (i) † . Here U (i) acts only on the Hilbert space H (i) ≡ C (2 ) of the cell around i. Particle preservation means that with N (i) = i ∼in i is the total number operator in the cell around site i. The sum carries over all sites i within distance /2 from i. As a consequence, U (i) decomposes into subspaces Λ p (p-th exterior product of C ) of dimension p where the number of fermions is fixed to p. That is, U (i) consists of blocks U (i) p on the diagonal. To construct the average over ballistic cells, we make the ergodic hypothesis that within a time interval t the ballistic cell has undergone all possible unitary and particle preserving transformation. Then time averages over intervals t can be replaced by a uniform Haar average [. . . ] 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-manybody system can be characterised 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 [38] and [39] (the latter in the context of ETH, though both talk about very similar things). 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.
The cells around sites i and j in the mesoscopic system above correspond 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 conserves the number of particles inside each cell.
If the separation of two sites i and j is larger than , 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 where denotes the Haar average. Note that by the cyclic property of the trace, we could have also chosen to average ρ t instead of c † i c j which makes clear that the average carries over all possible evolutions of the cell.
The time average over t promotes G ij (t) to a random variable that is only sensitive to the long (with respect to t ) time behaviour of the system. Therefore we also expect the right hand side not to depend on the time at which ρ t is evaluated, as long as t ∈ [t, t + t ]. Here we chose to evaluate at time t.
For quadratic fluctuations this generalises to and similarly for higher order fluctuations. a. U(1)-invariance. The averages appearing in (6) and (8) 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( i θ ini ), which are generated by the particle number operatorn i = c † i c i . The creation and annihilation operators c † i and c i carry the charge +1 and −1 with respect to this U (1) transformation. That is Then it is easy to see that any combination of operators inside an average [· · · ] U ⊗··· must be composed of an equal number of c † i '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 (1). For example, (8) is non-zero only if i = l and j = k.
Evaluating averages such as (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 sites in the cells. This is shown explicitly in the remainder of this subsection. The discussion can be skipped without impact on the comprehension of the following.
Unitaries in the one-particle sector are of the form where c Here u := e M * = 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 (comment [40]). For example, evaluating (6) using (11), 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 non-random as a result of the law of large numbers. Evaluating (8), using appropriate Haar averages one finds for any choice of i and j, up to corrections in 1/ , one finds b. Scaling with system-size. The scaling of loop expectation values (2) with system-size L can be derived from the assumption that the system satisfies the classical macroscopic fluctuation theory (MFT). In particular, this means, that the particle density n i at site i satisfies a large deviation principle, Prob(n 1 , · · · , n L ) ∼ e −LF (n1,··· ,n L ) , where F is the so-called quasi potential [7]. For a system in equilibrium with a heat bath, F would be the difference in free energy between the density profile {n 1 , · · · , n L } in question and the equilibrium density profile. One of the great insights of MFT is that the large deviation principle also extends out-of-equilibrium systems if they are diffusive.
Denoting the noise average at time t within classical MFT by · · · t the large deviation principle implies that connected density correlations scale as where all indices are different.
In the quantum description, we replace n i →n i = c † i c i and · · · t → E[Tr(ρ t · · · )] such that where E is the average over ballistic cells introduced above. If we assume the Hamiltonian to be of the form where H 0 is quadratic, and furthermore we initialise the system in a Gaussian state of the form , then we can apply perturbation theory to decompose (15) into a sum of Feynman diagrams.
More precisely, we employ the interaction picture where one divides the full time evolution U t = W t U (0) t into a quadratic evolution U (0) t = e −iH0t , followed by an evolution with the time ordered exponential The density matrix at time t is preserves its Gaussian form. Now (15) becomes which can be expanded using Wick's theorem and the resulting contractions can be denoted by Feynman diagrams.
For n = 2 we would get where the first diagram is the 1-particle irreducible scattering of two fermions created and annihilated at sites i and j, the second diagram is the multiplication of two propagators (Green's functions) between sites i and j, i.e. G ij G ji = Tr(ρ t c † i c j )Tr(ρ t c † j c i ), and the last diagram is the amplitude to stay on the sites i and j, i.e. the multiplication of the densities G ii G jj .
When taking the expectation value E, these three diagrams should be equal to the corresponding value within MFT, n i n j t = n i n j c t + n i t n j t .
Since the second term n i t n j t = E t [G ii ]E t [G jj ] (to leading order) is equal to the third diagram, the first term n i n j c t ∼ L −1 must equal the first plus the second diagram in (17). If we assume that the expectation value E of these two diagrams has the same scaling with L, we can conclude that the expectation value of the second diagram satisfies the scaling we claimed in the beginning of the section, that is E[G ij G ji ] ∼ L −1 . This argument can be repeated for any n in a similar fashion if one assumes that the expectation value of 1-particle irreducible diagrams with more legs has the same scaling as the expectation values of the multiplication of several crossing propagators. We admit that this assumption is a weak point in our argument and one should start here if one is interested to understand for which class of Hamiltonians this general picture applies. In particular, our argument would not apply if perturbation theory breaks down, signalling a possible phase transition, or in presence of elementary excitation of a non-perturbative nature.
c. Factorisation of products of loops. The factorisation in (3) for the case i 1 = j 1 is trivial, because the average E t treats each ballistic cell independently by construction. It remains to treat the case i 1 = j 1 . We illustrate this case by an example of a 2-loop.
Assuming the scaling from the last section, a look at (13) reveals that for i = j the leading contribution comes from the second terms which is O(1) while the first terms is O(1/L) (comment [41]). That is, to leading order, and restricted to the one-particle sector, we have, In the last equality we used (12). This calculation shows that for the given example the expectation value of products of loops (here G 2 ii ) 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 (3).

B. Link with free probability
The link with free probability occurs when expressing loop expectation values E[G i1i2 ...G ini1 ] in terms of their connected parts. We will find that this expansion involves a sum over non-crossing partitions only, see (31). Such non-crossing partitions are at the heart of a mathematical theory of non-commuting random variables that is called free probability theory and introduced in the next section. By means of non-crossing partitions one can define socalled 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, (20) Here we dropped the subscript t of the measure E t for simplicity. Whenever two indices i k in the above sum are equal, the expectation value factorises according to (3). We therefore split this sum into sums where all indices are distinct -except for a given set of indices that are imposed to be equal, Such a splitting can be understood as a sum over partitions π ∈ P (n) of the set {1, · · · , n} into blocks b ∈ π that group together all the indices i k that are imposed to be equal. For example, the partitions corresponding to the three terms above are π = {{1}, · · · , {n}}, π = {{1, 2}, {3}, · · · , {n}} and π = {{1, 2, · · · , n}}. The sum becomes i1,...,in = π∈P (n) i1,...,indistinct, except i k =i l whenever k,l are in the same block of π .
For n = 4 two possible partitions are for example and Representing partitions by diagrams where nodes in the same block are connected by dashed lines shows very intuitively that π 2 is a crossing partition while π 1 is noncrossing. It turns out that in the sum (20) terms corresponding to crossing partitions are of the order of 1/L and lower and therefore vanish for L → ∞, while all non-crossing partitions are of order one. Here we illustrate this fact for n = 4 and the two examples above.
The term corresponding to π 1 factorises and becomes Each of the 2-loop expectation values scale 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 to this the term corresponding to π 2 can factorise in two different ways and becomes The difference to π 1 is that now there are only two indices to sum over and hence the scaling is of order 1/L. For the surviving non-crossing 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 factorisation. If we modify (23) by connecting as many edges by solid lines as possible without crossing a dashed line, i.e.
we see that for the partition π 1 this product of loop ex- , is determined by the solid lines. In fact, for every non-crossing partition π of nodes, the solid lines define a unique noncrossing partition π * 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 π * 1 = {{1, 2}, {3, 4}}. Each block therein corresponds to a loop of G ij 's. In this terminology the expansion of (28) where N C(n) denotes all non-crossing partitions of the set {1, · · · , n}.
Instead of using this slightly akward sum over indices i k we introduce a modified Kronecker-delta that sets all indices to be equal that belong to the same block in π. Next, we extend the sum to include the cases in which indices i k can become equal by replacing E[· · · ] with E[· · · ] c . This introduces only a sub-leading error of order 1/L, because (4) ensures that connected expectation values do not get more dominant when some indices become equal. Then the expansion of φ(G n ) in (20) reads In free probability φ corresponds to the "expectation" value of the random matrix G and we have shown that the moments of G have a natural expansion in terms of a sum of non-crossing partitions. The terms in this sum are usually called free cumulants. However, comparing to the definition of free cumulants (38) in the next section, one sees that the expansion above has a slightly different form due to the additional δ π . This issue is further discussed in section IV F To conclude the argument, note that we could have done the whole derivation multiplying E[G i1i2 · · · G ini1 ] in (20) with a smooth test function h i1,··· ,in that would again appear in (30). By comparing the two equations we would find that an equation we will reuse in section IV. Here we interchanged the role of π and π * which is possible because, for non-crossing 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
Simultaneously to our work it has been observed in [36] that there is a very similar link to free probability in the context of the eigenstate thermalisation 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 = E i |A|E j of a local observable A expressed in the energy eigenbasis |E i of a hamiltonian H in a closed system that obeys ETH. In the context of ETH, A ij is a random variable with respect to an fictitious ETH-random-matrix-ensemble that captures its typical behaviour.
Comparing to [39], which discusses higher order fluctuation of A ij within ETH, one realises that the emergence of the three properties (1)-(3) have analogous reasons. In 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 called ballistic cells. To complete the analogy, loop expectation values of A ij in ETH scale with the density of states e S(E+) at energy . In our context, this corresponds to the "density of states" at sites i and j, which is just a constant and equal to the number of sites L (because the physical length of our system is set to one for simplicity).

III. INTRODUCTION TO FREE PROBABILITY
In classical probability, two variables are independent if (and only if) their moments factorise at all orders, for all n, m ∈ 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 factorisation of, say, E[XY XY ] into the moments of the independent variables E[X 2 ] and E[Y 2 ], on the level of matrices, since they don't 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 similarly to generators in a free group, hence the name.
Given two non-commuting random variables a and b in some algebra M (e.g. algebra of large random matrices) and a linear functional ϕ : M → C (that plays the role of the expectation value), then a and b are called free if for all polynomials P 1 , · · · , P l and Q 1 , · · · , Q l with ϕ(P i (a)) = 0 and ϕ(Q i (b)) = 0 we have The reason to evoke all possible polynomials in the definition is that any element in the subalgebras A and 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 A and B. This definition implies for example that if a and b are free, then ϕ(a m b n ) = ϕ(a m )ϕ(b n ). Additional structure occurs if one interchanges the order such that free variables are no longer grouped together. But it remains always true that joint moments can be determined through the moments of the individual free variables therein. For The definition of freeness was proposed by Dan 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 provide the lecture notes by Roland Speicher [42] as well as the book by Mingo and Speicher [35].
In the 1990's, 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 allow us then to better appreciate why free cumulants are defined via noncrossing partitions.

A. Classical cumulants and partitions
Let {X 1 , · · · , X N } be a family of classical random variables with moment-generating-function where u is an (optional) parameter that counts the order of the joint moment (or correlation function) 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 E[X i1 · · · X in ] can be expressed as a sum of products of cumulants arranged according to partitions π of the set {i 1 , · · · , i n }, This is in full analogy to (23) and (24), except that here we used solid lines for better visibility, since we don't talk about dual partitions in this section. The expansion of the moment E[X 1 X 2 X 3 X 4 ] in terms of products of represented through the diagram of the corresponding partition π, becomes where the subscript k suggests that by cyclic permutation there are in total k such diagram. Note the fact, that the last diagram (in a dotted box) corresponds to a crossing partition π = {{1, 3}, {2, 4}}. In free probability theory these diagrams do not appear as we will see below.
a. 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 = · · · = X N we illustrate how this can be done up to order four. Let us denote by m n = E[X n ] and c n = E[X n ] c the moments and cumulants of this variable, then Note that the coefficients correspond exactly to the cyclic multiplicities of the diagrams. This can be solved recursively for c k , 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 non-crossing partitions
In the setting of non-commuting random variables a 1 , · · · , a N in some algebra M with "expectation value" ϕ : M → C, the free cumulants κ n are multilinear forms implicitly defined through moments by a sum over noncrossing partitions π ∈ N C(n) of the set {i 1 , · · · , i n }, ϕ(a i1 · · · a in ) =: (38) where |b| denotes the number of elements in the block b = {b(1), b(2), · · · }. For example, if we would expand ϕ(a 1 a 2 a 3 a 4 ) in terms of free cumulants we could get all diagrams in (35), except the last one which is crossing. Therefore the order of the arguments in κ n becomes important, even if the a i were to commute, hence the separation by the comma.
The free cumulants satisfy a number of properties that are analogous to properties of classical cumulants: • κ n (a 1 , · · · , a n ) = 0 as soon as there appears a pair a i , a k of free variables, see e.g. [42] • As a result of multiliniarity and the last bullet point, free cumulants of free variables a and b are additive, κ n (a + b, · · · , a + b) = κ n (a, · · · , a) + κ n (b, · · · , b), see [43] • Any variable a whose free cumulants κ n (a, · · · , a) vanish for n ≥ 3 is distributed according Wigner's semi-circle law of random matrix theory for the Gaussian unitary ensemble (GUE), so that large GUE random matrices are analogous of Gaussian variables but from a free probability point of view.
The number of non-crossing partitions of a set of n elements is the Catalan number C n = 1 n+1 2n n , with C 1 = 1, C 2 = 2, C 3 = 5, C 4 = 14 and C 5 = 42, etc. The formula (38) is triangular and can be inverted to express the free cumulants in terms of the moments.
a. Free cumulants of a single variable. In the case of a single variable a = a 1 = ... = a N , we denote by κ n := κ n (a, · · · , a) and m n := ϕ(a n ) the free cumulant and the moment at order n, respectively. Then we have The equations can be solved for κ k recursively, Note, that the difference between standard and free cumulants only shows up at order 4 since here a crossingpartition become 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 [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 realised that matrices in the Gaussian unitary ensemble (GUE) with independent entries become free variables in the limit where the dimension of the matrix goes to infinity. Since then many more connection between other random matrix ensembles and free probability have been found.
Here we will make one of these results more explicit, which applies to matrices that are rotated by random unitaries. Consider N × N random matrices of the form where U N 's are choses independently from the Haar distribution over the unitary group and A N and B N are deterministic matrices with spectral densities µ A and µ B . That is, the moments m k := lim N →∞ 1 N Tr(A k N ) = λ k µ A (λ)dλ are all finite, and similarly for B N . The statement is that in the limit N → ∞ the random matrices X N and Y N become free variables a and b (in some non-commutative probability space) with distribution µ A and µ B with respect to the measure ϕ := 1 L E Tr, where E is the Haar measure. It is furthermore known that the classical cumulants of such matrices X N can be expressed as the free cumulants κ n ≡ κ n (a, · · · , a) of the spectral density µ A . That is, where P N is a sequence of matrices with fixed rank (such that Tr(P k N ) 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]. That is, the matrix of coherence G is distributed as U G 0 U † with U a Haar distributed unitary L × L matrix. As a consequence of (41) above, 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 section II B, one can derivation the link between fluctuations in mesoscopic systems and free probability starting from the moment-cumulant formula (33) introduced earlier in this section. This is also the way we chose in appendix E 1 which contains the full proof of this link for any n.
We apply (33) to a loop of G ij 's and obtain where the subscript t of E t is dropped for simplicity. For n = 4 two possible partitions, each with two blocks, are and We determine the contribution of the two partitions: In the second line we used the U (1) invariance of E that implies that the expression is non-zero only for i 1 = i 3 . Finally we evaluated the scaling of the expression with L. Here we used that E[G ij G ji ] c = O(L −1 ) by (4) and that the Kronecker delta becomes δ ij → δ(x − y)/L for large L, where x = i/L and y = j/L. Note that the resulting scaling with O(L −3 ) is the same as that of a 4-loop, 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 non-crossing partitions have the same scaling with L.
In contrast to this, crossing partitions turn out to be sub-leading. Evaluating π 2 one has Here we assumed that no connected expectation value is more dominant than that of loops with the same number of points, i.e.
, because the loop with the same number of points scales as O(L −1 ). Note that the average over independent ballistic cells introduced in section II A actually implies that E[G ii G kk ] c = 0 if i = k. However, we anticipate that the toy model Q-SSEP predicts this term to scale with L −2 if i = 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 non-crossing 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) π * . This is best understood diagrammatically. To π 1 we associate the diagramm The solid lines represent blocks in π 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 If π is non-crossing then there is a unique way to connect a maximal number of nodes by dashed lines without crossing a solid line. To a dual partition π * we associate a Kronecker delta δ π * (i 1 , ..., i n ), that equates all indices belonging to the same block of π * , exactly as in (29) but with π and π * interchanged. For example δ π * 1 (i 1 , ..., i 4 ) = δ i1i3 . This notation allows us to rewrite the momentcumulant-formula (42) as a sum over non-crossing partitions only. Denoting the set of non-crossing partitions of {1, . . . , n} by N C(n), this is (48) This is exactly the same formula as (31).

A. Summary of results
The quantum symmetric simple exclusion process (Q-SSEP) describes a one-dimensional chain with L sites on which spinless fermions can hop to neighbouring sites with random amplitudes. In the closed Q-SSEP this chain is closed periodically, while in the open Q-SSEP the chain is coupled to two reservoirs that can inject and extract fermions on the boundary with different rates. It is therefore possible to have a steady current of fermions through the chain keeping the system out-of-equilibrium. The model allows to study the time evolution of spatial coherences G ij (t) = Tr(ρ t c † i c j ) which fluctuate due to the random hopping. The measure E t of this randomness converges to a unique steady measure E ∞ after long times which has been characterised in [27].
The main results on the open Q-SSEP in this paper are equations for the time evolution of the n-loop expectation value E t [G i1i2 G i2i3 · · · G ini1 ] and its connected part, in a non-trivial scaling limit where L becomes large. In deriving them, we make use of non-crossing partitions which appear in the same way as in section II B because Q-SSEP satisfies the three conditions (1)-(3) (see section IV B).
It has been shown in [30], that the closed Q-SSEP has a meaningful scaling limit for L → ∞ if positions i k and time t are scaled diffusively according to In this limit one defines the expectation value of on nloops which satisfies Pictorially, this equation reduces the evolution of an n-loop expectation value to diffusion sourced by the product of two loop expectation values with less nodes that emerge by pinching the original loop along the nodes x i and x j , We therefore say that the evolution for g s t has a triangular structure.
In section IV D we will show that this equation is also valid for the open Q-SSEP, if we require g s t (x 1 , · · · , x n ) on the boundary x i ∈ ∂[0, 1] to be equal to its steady state value at all times. Due to the triangular structure of the equation this behaviour can be traced back to the fermion density, i.e. the 1-loop g s t (x), which approaches its steady state value on the boundary immediately, in contrast to its evolution in the bulk which happens on a much slower time scale. This is explained in section IV C.
As discussed below (3), loop expectation values jump by an order of L if two indices become equal. This is reflected by the fact that solutions g s t (x 1 , · · · , x n ) of (51) become singular whenever two of its arguments x i = x j are equal, hence the choice for the superscript "s". We show this in appendix B for the example of n = 2.
In appendix E 2 we derive the time evolution of connected loop expectation values in the scaling limit. Connected n-loop expectation values are defined as (53) We find that they satisfy with boundary conditions g t (x 1 , · · · , x n ) = n a , n b for n = 1 and x = 0, 1 0 for n ≥ 2 and some x i ∈ {0, 1}, (55) where only the one-point function (fermion density) depends on particle density of the two reservoirs n a and n b . This equation produces indeed solutions that are continuous if two arguments become equal, see appendix B. The equation differs from (51) by the relative position of the derivatives and the delta function.
The derivation of equation (54) crucially depends on the fact that loop expectation values (moments) and their connected parts (cumulants) are related by a sum over non-crossing partitions such as in (31). In the scaling limit this relation becomes is defined in analogy to (29), but for continuous variables. Finally, in section IV F we show how to construct a very simple steady state solution for connected loop expectation values. This exploits the insight below (30) where we noted that the presence of the δ π prevents us from viewing connected loop expectation values as the free cumulants of the measure E t . If we remove the δ π , then connected loop expectation values are, by definition, the free cumulants of a new measure ϕ t that is necessarily different from E t . Surprisingly ϕ t has a very simple steady state solution from which the connected loop expectation values g ∞ can be determined recursively.

B. Introduction to the model
The quantum symmetric simple exclusion process was first introduced in [25] (closed case) and [27] (open case) and was further elaborated in [26,28,30].
a. Definition. The model describes a onedimensional chain with L sites on which spinless fermions hop to neighbouring sites with random amplitudes. The boundary sites are coupled to particle reservoirs (comment [45]). The bulk evolution of the system is stochastic and unitary. We describe it in terms of the systems density matrix ρ t as where the Hamiltonian increment is defined as with c † j (c j ) fermion creation (annihilation) operators at site j with the usual commutation relations {c † j , c k } = δ jk and {c † j , c † k } = {c j , c k } = 0. The noisy Hamiltonian increment depends on a complex Brownian motion dW j t which determines the random hopping amplitude along the edge (j, j + 1) at time t. There is one complex Brownian motion per edge. This means that W j t at each instance t is a centred Gaussian random variable that has independent increments dW j t = W j t+dt − W j t with variance E[dW j t dW k t ] = Jδ j,k dt if t = t , and zero otherwise. J is a rate parameter with units one over time and we set J = 1 in the following.
To the stochastic and unitary bulk evolution (58), we add a deterministic but dissipative Lindbladian evolution that represents the interaction with the reservoirs. The •} models particle injection and is multiplied by the injection rate α j while •} models particle extraction with rate β j . For example, the density matrix of an isolated empty site τ t = |0 0| that evolves according to ∂ t τ t = αL + (τ t ) will be occupied after a time interval dt with probability αdt. That is, The full evolution of the systems density matrix ρ t can be expressed as a stochastic differential equation (SDE) in Itô convention (with Itô rules dW j t dW k t = δ j,k dt) by expanding (58) up to order dt, Occasionally we will 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 ≡ L are identified.
b. Relation to the classical SSEP. The name of this model is inherited from the classical symmetric simple exclusion process (SSEP). The latter can be obtained from Q-SSEP in the special case where one is just interested in the mean density matrixρ t := E[ρ t ], where the expectation E[· · · ] is taken with respect to the different realisations of the Brownian motions. This matrix evolves according to a Lindblad equation, Writingρ 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ρ t provide the probability that the system is in one of the 2 L states with well defined occupation numbern i = c † i 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 see easily that the Lindbladian evolution of the diagonal elements ofρ 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 neighbouring site with probability dt if the site is empty and particles get injected (extracted) at the boundaries with probability α i dt (β i dt). This correspondence can also be formulated via the moment generating function, and it illustrates the well known correspondence between Markov processes and Lindbladian evolutions. The relation between Q-SSEP and free cumulants we shall describe below implies a new relation between the classical SSEP and free probability [46]. c. 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) = Tr(ρ t c † i c j ). While the coherences vanish in mean exponentially with time (we saw 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 sub-leading in the system size [27]. For example, at large system size L → ∞ with x = i/L ≤ y = j/L fixed, the connected quadratic fluctuation (or 2nd cumulant) in the steady state is where ∆n = n a − n b is the difference in the particle density between the boundary reservoirs, see (64). Here we adopted again the point of view where the time dependence of G ij (t) is transfered to the measure E t and the steady measure is denoted by E ∞ . Let us also note that since the Hamiltonian is quadratic in fermionic creation and annihilation operators, the coherences G ij completely characterise the state of the system due to Wick's theorem. From (59) one finds that their time evolution is given by the SDE d. Three properties of E t . Firstly, the measure E t of Q-SSEP possesses a local U (1) invariance which can bee seen as follows: The multiplication with local phases, G ij →G ij = e −iθi G ij e iθj leaves (61) invariant if also the Brownian motions are multiplied by a phase dW j t → dW j t = e i(θj+1−θj ) dW j t . Since dW j t and dW j t have the same distributions, also G andG 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 (∼ e −L 2 t where t is the rescaled time).
Secondly, n-loop expectation values in Q-SSEP scale as L −n+1 if indices are distinct because this scaling gives rise to the non-trivial equation (51).
Thirdly, expectation values of products of loops factorise. This is shown in appendix G of [30] where the time evolution equation of loop expectation values (51) is derived (for the closed case). The reason is that the time evolution equation for the leading order of a product of loop expectation values allows for factorised solutions. Since the initial condition is always factorised (the noise has not yet correlated the variables) the claim follows.
This shows that for Q-SSEP properties (1)-(3) are satisfied and as a result Q-SSEP has the expansion of loop expectation values into non-crossing partitions (31) as claimed.

C. Bulk vs. boundary thermalisation
When we casually talk about bulk and boundary thermalisation, one should keep in mind that the open Q-SSEP does not actually "thermalise". At large times, the system approaches a steady state (comment [47]), in which observables do not depend on time any more longer. But they are not described by a thermal density matrix since there is a steady current, so that 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. Due to the correspondence (60) this property also hold true for the classical SSEP (comment [48]).
Under the evolution outlined above the mean fermion density n i (t) := E t [Tr(ρc † i c i )] evolves according to where it is understood that the discrete Laplacian on the boundaries is truncated, i.e. ∆n 1 := n 2 − n 1 and ∆n L = n L−1 − n L . Here we only discuss the special case, where the injection/extraction parameters are such that α 1 + β 1 = 1 = α L + β L . The general case is discussed in Appendix C and works analogously with a little bit of help of numerical calculations. In the special case, (62) can be solved analytically where the coefficients c k are determined by the initial condition and b k : . 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 (∞) = α1(L−j+1)+α 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 characterised by j ∼ aL with a ∼ O(1). Due to the factor e (−2+2 cos( kπ L+1 ))t ≈ e −π 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π/(L + 1)) ∼ sin(aπ). Therefore we find a time scale of t decay ∼ O(L 2 ), n bulk (t) ∼ const.e −π 2 t/L 2 + n bulk (∞).
A site on the boundary is characterized by j = 1 or j = L and therefore the amplitude of the k = 1 term, sin(jπ/(L + 1)), will be zero for large L. A significant contribution therefore only comes from terms where k ∼ bL with b ∼ O(1), because then the amplitudes sin(jkπ/(L + 1)) ∼ sin(bπ) or ∼ sin(bLπ) stay finite. The time scale with which these terms decay is of order one, t decay ∼ O(1), n bdry (t) ∼ const. e −b 2 π 2 t + n bdry (∞).
Note that the value of the density on the boundaries after a time of O(1) is that of the steady state. For general injection/extraction parameters these values are as in the classical SSEP, We conclude that time scales 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 → 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 decay ∼ O(1) of time.

D. Open boundary scaling limit
The basic question is whether the scaling limit (49) introduced for the closed Q-SSEP is meaningful also in the open case. In other words, does g s t satisfy a nontrivial equation in this limit? a.
Density We first answer this question on the level of the density n i (t) = E t [G ii ]. It satisfies L coupled differential equations in time given by (62). The aim is to represent these by a single partial differential equation in space and time for a continuous density ρ(x, t) ≈ n Lx (L 2 t). Since it is easily seen that the bulk satisfies a pure diffusion equation the remaining question is what are the correct boundary conditions. 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, This is confirmed numerically in Figure (3), where a solution for the discrete density n i (t) is compared to the solution for the continuous density ρ(x, t). In particular the boundary values agree. b. Higher order fluctuations. The claim we make here is that the exact same equation (51) that holds in the closed case also holds in the open case if we subject g s to boundary conditions that are equal to its steady state value. These are known from [27] for its connected part: g ∞ (x 1 , · · · , x n ) = 0 for all n ≥ 2 and (x 1 , · · · , x n ) ∈ ∂[0, 1] n .
We test this claim on the level of the 2-loop expectation values, for which (51) simplifies to and we identified 1-loops with the density, g t (x) = ρ(x, t).
From this it follows that the time evolution of the connected correlation functions with appropriate boundary conditions is which is a again a diffusion equation, but with a new source term. To see this, note that the connected correlation function is defined by The analytic solution of (68) (which is constructed in appendix D) is compared to a numerical solution of the discrete evolution equations of L( for given L -which can be derived from (61). Figures 5 and 4 show the result of this comparison for different system sizes L. The agreement is excellent. Notice that the injection and extraction rates enter (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 2-loop expectation values g t (x, y) (and for all higher order loops) never explicitly depend on the injection and extraction rates. As a consequence, without loss of generality, we use the convention n a = 0, n b = 1, in the following, through out the paper. In particular, the density in the steady state g ∞ (x) = n a + x(n b − n a ) reduces to g ∞ (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, with boundary conditions g t (x 1 , · · · , x n ) = n a , n b for n = 1 and x = 0, 1 0 for n ≥ 2 and some x i ∈ {0, 1}. (70)

E. Free cumulants in Q-SSEP
Comparing the moment-cumulant relation (31) (valid for any system satisfying (1)-(3), in particular for Q-SSEP) to the definition of free cumulants (38) one would be tempted to identify with κ n the free cumulants of E t . However, this is not correct due to the presence of δ π . If we nontheless insist on the identification (71) then the connected loop expectation values E[G i1i2 · · · G ini1 ] c are, by definition, the free cumulants of a new measure ϕ t that is different from E t . That is, In terms of g t (x 1 , · · · , x n ) ∼ L n−1 E L 2 t [G i1i2 · · · G ini1 ] c the new measure ϕ t defines a function (which we denote by the same name) of continuous positions and rescaled time The time evolution of ϕ t is found to satisfy exactly the same equation as g t (see appendix E 3) However the boundary conditions are different. If some x i ∈ {0, 1} lies on the boundary, then where the hat onx i indicates that x i is missing from the set {x 1 , · · · , x n }.

F. Steady state solution -inspired by free probability
The striking observation by Biane [32] was that ϕ t , defined as a sum over non-crossing partitions of products of connected loop expectation values, has a very simple solution in the steady state. In appendix E 4 we show that for t → ∞, our equation for ϕ t is indeed solved by ϕ ∞ (x 1 , · · · , x n ) = min(x 1 , · · · , x n ) As a consequence, we can use (73) to recursively reconstruct the connected correlations functions g ∞ in the steady state. This works similarly to (39) and (40). Denoting min(x 1 , · · · , x n ) =: x 1 ∧ · · · ∧ x n , we have while the four point function g ∞ (x 1 , x 2 , x 3 , x 4 ) reads where · · · q denotes the sum of all terms obtained by q successive cyclic permutations of the arguments of the term in question. Note that due to the absence of crossing partitions, g ∞ (x 1 , · · · , x n ) is no longer invariant under the permutation of its arguments for n ≥ 4. In the example above it is the term (x 1 ∧ x 2 )(x 3 ∧ x 4 ) 2 which prevents g ∞ (x 1 , · · · , x 4 ) from being invariant under the exchange of x 2 and x 3 .
The derivation of the steady state solution presented here provides an alternative proof of Biane's formula (6.1) in [32] that relates the steady state connected correlations of the open Q-SSEP to free cumulants of the measure ϕ ∞ . Our derivation extends this relation to finite times, though an explicit solution for ϕ t at finite times seems out of reach.
Note that the measure ϕ ∞ can be realized as the Lebesgue measure on the interval [0, 1] of the indicator function I x := 1 [0,x] , e.g.
It is surprising that ϕ ∞ has a realization in terms of commuting variables, since free cumulants usually appear in a setting of non-commuting variables. At finite times (74) suggests that ϕ t is not invariant under a permutation of its arguments.

V. CONCLUSION
In this paper we presented two main points. (1) A general argument why the fluctuations of spatial coherences in one dimensional mesoscopic quantum systems could be well described by the framework of free probability theory. (2) A precise calculation that shows that the model Q-SSEP has a mathematical structure that fits into the framework of free probability and 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 E[G i1i2 · · · G ini1 ] ∼ L −n+1 , and (iii) the fact that expectation values of products of loops factorise. 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 eigenstate thermalisation hypothesis (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 outof-equilibrium. Rather, this link emerges from a coarsegrained description under the assumption that, locally in space (i.e. within ballistic cells) and on time scales much shorter than the diffusion time, the system is ergodic. Such an assumption allows to introduce the noise average as an average over all possible unitary transformation 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 section II is certainly oversimplified and calls for more details. Firstly, the argument assumes a separation of time scales which gives rise to a crossover from ballistic to diffusive transport at some length scale . One should provide criteria on the Hamiltonian for when this separation of time scales is satisfied. Secondly, one should explore the physical meaning of the length of ballistic cells. In the introduction we crudely argued that could be related to the mean free path of electrons in a disordered metal. And thirdly, 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 characterise 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 timedependence of the Hamiltonian ensures that energy is not conserved which would otherwise represent another conserved quantity.
Such numerical studies would also answer the question if 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 use a decomposition into equilibrated and statistically independent cells to characterise 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 find that the size of the cell scales as 1/γ where γ 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 γ → ∞ [25]. In this limit, the size of the cell becomes zero which corresponds to the idea that for Q-SSEP the mean free path has effectively been 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. and J. Kurchan [36], where the relation between these two frameworks is discussed. The occurrence of free probability in both problems has a similar origin: the coarsegraining at microscopic either 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. LH 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 number ANR-20-CE47-0014-01.

Appendix A: Measuring coherences
To gain a better understanding of the coherences G ij = Tr(ρc † j c i ) we outline an experiment to measure them that was proposed in [18]. The setup is shown in Figure (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. 6. Two wires are 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 number nL and nR is measured on each side. In the first measurement (a) one uses a symmetric beam splitter, which allows to measure the imaginary part of Gij. 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 π, while it does not accumulate this phase when being transmitted in the other direction. In this way one can measure the real part of Gij.
• The total state of system, left and right wire is described by a state in the Hilbert space H S ⊗ H L ⊗ 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, • 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 int = e −iλ(c † L ci+c † R cj +h.c.) , where λ 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 and duration such that λ 1 is small, we can neglect O(λ 2 ) terms and find the comp/lete state to be • Next, the fermions in the left and right wire interfere in a beam splitter. Written in the basis {|00 , |01 , |10 , |11 } the beam splitter can in general be described by the scattering matrix 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 where there is a fermion in each wire is suppressed by λ 2 , lets give a short explanation for how to to obtain last entry |1, 1 → (rr − tt )|1, 1 . One has take into account that the wave function is antisymmetric, |1, 1 = |φ L ⊗ |φ R − |φ R ⊗ |φ L . Here the position in the tensor product labels the fermion (say they are called 1 and 2), whereas |φ L and |φ R are single fermion states in the left and right wire. Then leads to the entry (rr − tt ). Unitarity demands |r| 2 + |t| 2 = 1, |r | 2 + |t | 2 = 1 andrt +tr = 0 (the condition |rr − tt | 2 = 1 is then automatically fulfilled).
-(a): Choosing a symmetric beam splitter, r = r and t = t , allows to measure the imaginary part of G ij . Note that the unitary constraints can now be expressed as r = sin θe iϕ and t = i cos θe iϕ and we set the overall phase ϕ = 0 since this will not change the result. The state evolves to If the beam splitter is symmetric, except that transmitted fermions incident from the right will accumulate an additional phase π, i.e. r = r and t = −t , this allows to measure the real part of G ij . Note that the unitary constraints result it r = sin θe iϕ and t = cos θe iϕ and again we set ϕ = 0. The state evolves to • Finally, we measure the particle number n L = c † L c L and n R = c † R c R in the left and right wire. Denoting averages w.r.t the system |ψ S by ... S as in G ij = c † j c i S we find for case (a) n L (a) = λ 2 sin 2 θ n i S + cos 2 θ n j S − 2 sin θ cos θ (G ij ) n R (a) = λ 2 cos 2 θ n i S + sin θ n j S + 2 sin θ cos θ (G ij ) .
Choosing the same angle θ = π/4 gives the real part of G ij , Hence, both the real and imaginary parts of the coherence, as well as their statistical distribution, are experimentally measurable (at least in principle).
Integrating these equations across the diagonal line {x = y} reveals that the van Neumann boundary conditions will be singular for g s t and regular for g t . Let us do this explicitly for (67). We rotate the variables (x, y) by π/4 clockwise, (v, u) = ( x−y √ 2 , x+y √ 2 ). Then the derivative ∂ v = ∂x−∂y √ 2 is orthonormal to the line {x = y} and ∂ u = ∂x+∂y √ 2 is parallel to the line. In these variables the equation reads Integrating − dv this equation and keeping only terms of O(1), The expression simplifies due to cyclic invariances, g s t (x, y) = g s t (y, x), Finally we take → 0 + and find the van Neumann boundary condition in both sectors {x > y} and {x < y}, The delta function will make this derivative blow up and therefore g s t (x, y) is indeed singular at x = y. The same procedure for (68) yields which is finite.
Appendix C: Solution for the discrete density The time evolution equation for the discrete density (62) can be rewritten as ∂ t n = An + n, with n = (n 1 , ..., n L ) T , b = (α 1 , ..., α L ) T and If A is diagonalisable, A = SDS −1 , a general solution for an initial condition n(0) = u is given by If α 1 + β 1 = 1 = α L + β L , A has is a so-called Tölpitz matrix and has eigenvalues λ k and normalised eigenvectors v k , with k = π L+1 , 2π L+1 , ..., Lπ L+1 , from which (63) can be derived. For general injection/extraction rates we make an ansatz where s k and k are determined by the eigenvalue equations. Denoting z := e ik this leads to where A = 1 − (α 1 + β 1 ) and B = 1 − (α L + β L ). Note that even though this equation has 2L + 2 solutions for z, only L on them give rise to eigenvalues with linearly independent eigenvectors: If k is a solution then also −k is a solution, but the corresponding eigenvectors are linearly dependent, v −k = s k v k . Furthermore k = 0 is a solution with zero eigenvector. Numerical solutions of C2 show that for A, B ∈ [1, −1] all solutions of z lie on the unit circle and therefore (taking into account that k ∼ −k and k = 0) we have k ∈ (0, π).
A site in the bulk, j ∼ aL, approaches the steady value (A −1 b) j according to the time dependent part of (C1). Again, since λ 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 * ≈ π/L is the smallest positive solution for k. Note that the amplitude (v k * ) j ≈ e iaπ − e −iaπ of this term stays finite (s k * ≈ −1). Then, the time scale with which the bulk approaches the steady value is, as before, t decay ∼ O(L 2 ).
At the boundary, the amplitude of the term corresponding to k * will be zero. E.g. for j = 1, (v k * ) 1 = e ik * + s k * e −ik * = 0. To have a non-zero amplitude, one has to consider terms where k ∼ bL with b ∼ O(1), which leads to t decay ∼ O(1).

Appendix D: Analytic solution of the connected 2-point function
Here we outline how to solve the connected 2-point function analytically, which was needed for the comparison with the solution of the discrete equations in the scaling limit in Fig. (4) and (5).
We find the solution by an expansion in {sin(nπx)} ∞ n=1 . These functions satisfy the correct boundary conditions and are orthogonal in the sense that 1 0 sin(nπx) sin(mπx)dx = 1 2 δ nm . Importantly, they form a complete basis of L 2 ([0, 1]) which justifies the expansion.
Taking into account the initial condition, this leads tõ In the special case where the boundary conditions match the initial conditions (n a = 1, n b = 0) one finds b. Solution for the connected two point function. In appendix B we saw 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 (68) on the lower triangle T + = {(x, y) ∈ [0, 1] 2 : x ≥ y}: First we identify a function that satisfies the boundary conditions, w(x, y, t) = y(1 − x)∂ x ρ(x, t)∂ y ρ(y, t).
Note that for t → ∞, 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 a inhomogeneous heat equation with homogeneous boundary conditions . This is solved, as before, by the method of eigenfunction expansion. Note that ψ nm (x, y) := sin(nπx) sin(mπy) + sin(mπx) sin(nπy) is a complete basis of L 2 (T + ) that satisfies the correct Dirichlet and Neumann boundary conditions. It satisfies We write S(x, y, t) = n≥m≥1Ŝ nm (t)ψ nm (x, y) and f (x, y, t) = n≥m≥1f nm (t)ψ nm (x, y), where the coefficients are given byŜ This leads to the ∂ tfnm + π 2 (n 2 + m 2 )f nm =Ŝ nm which is solved bŷ .
The precise expression forŜ nm is rather complicated, and it is easier to solve for f part using the Mathematica NDSolve function. The complete solution is then g t (x, y) = w(x, y, t) + f hom (x, y, t) + f 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) ∈ [0, 1] 2 : x ≥ y}.
a. Scaling of non-crossing partitions. We start the proof by showing that if π = {b (1) , · · · , b (m) } is a non-crossing partition of {12, 23 · · · , n1}, i.e. of the edges of a loop with n nodes, then [π] c := E π [G i1i2 · · · G ini1 ] c scales as L −n+1 independently of the number of blocks m -and therefore behaves in the same way as the single loop expectation value E[G i1i2 · · · G ini1 ] c with n points.
First notice, that if the blocks of π are not nested into each other (e.g π = {{12, 23}, {34}, {45, 56, 61}}) then, by U (1) invariance, we have to connect starting and endpoints of each block by a Kronecker δ. Once we arrive at the last block, this conditions is already satisfied due to the other delta functions. We therefore need m − 1 delta functions, where δ ∼ 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 π has some nested blocks (e.g. π = {{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 possible collections of nested sub-blocks therein. In the end we are left with m − 1 delta functions needed to close each block in π to form a loop. This shows that [π] c ∼ L −n+1 if π is non-crossing.
b. Scaling of crossing partitions. Now assume that π consists of a collection B of non-crossing blocks with |B| elements in total and a collection C of crossing blocks that cannot be disentangled from another with |C| elements in total (e.g. 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 ∼ L −|C| (as will be shown below) then where only a single delta function is needed to connect the collection B and C (because all the other delta functions are already included in [B] c and [C] c . If π has more than one collection C of crossing blocks that cannot be disentangled, then it will scale with even higher negative power of L. This shows that crossing partitions are sub-leading in the scaling limit compared to non-crossing partitions. It remains to show that [C] c ∼ L −|C| (probably it is even true that [C] c ∼ L −|C|−1 ). The idea is to produce a collection of crossing blocks starting with non-crossing blocks and permuting its elements. Assume that {b (1) , b (2) } is a collection of two non-crossing blocks which would mean that there is only one delta function δ(i 1 , i k+1 ) necessary for the product of these blocks to be non-zero. 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) , which makes it necessary to insert a second delta function, δ(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 that (a) is still in its original block, (b) whose neighbours 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) 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) ] ∼ L −k +1 and [b (2) ] ∼ L −(n−k )+1 . We therefore have As a consequence, in the case where C consists of two crossing blocks, this shows that it scales at most with [C] c ∼ 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 #crossings + 1. However, all we need here is that m crossing blocks will cause at least m delta functions. Then the scaling is is sufficient. We saw that 2 crossing blocks cause at least 2 delta functions. But each originally non-crossing 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 2 delta functions δ(2, 3)δ(1, 4) and we add a block, we have {{12, 34}, {23, 45}, {56, 61}, which needs 3 delta functions, δ(2, 3)δ(1, 4)δ(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.

Time evolution of connected loop expectation values
We will show how (54) follows from (51) by induction over n. We start by writing the sum 2 i<j = ij,j =i in these equations as a sum over all ordered, non-crossing and non-empty sets r = {r 1 , r 2 , ...} and s = {s 1 , s 2 , ...} with i = r 1 and j = s 1 such that r s = [n] ≡ {x 1 , ..., x n }. Throughout this paper we use the symbol to denote the union of ordered, non-crossing and non-empty subsets. Furthermore, instead of g t (x i , x i+1 , ..., x j−1 ) we simply write g(r) if r = {x i , x i+1 , ..., x j−1 }. Then the two equations become, and δ(r 1 , s 1 )∂ r1 g t (r)∂ s1 g t (s).
For n = 1 the two equations are identical. Let's assume that (E3) holds for all k ≤ n − 1 for some n ∈ N. We will use (E2) to show that it holds also for k = n.
"←": One starts with (D π , x, y). Since the two marked nodes (x, y) ∈ d ∈ π * belong to a block of the dual partition, one can cut the diagram at these two nodes without breaking any block b ∈ π. The cut therefore defines independent non-crossing partitions ρ and σ on the two non-crossing subsets r s = [n] where r 1 = x and s 1 = y.
With the help of (E6) and (E5) we can rewrite (E2) and isolate the connected correlation function g([n]), which is the term we are aiming for, The only information which is missing is the action of (∂ t − ∆) on D π . We claim that (∂ t − ∆)D π = (x,y)∈d d∈π * (∂ x ∂ y − ∂ g x ∂ g y )D π + (x,y)∈b b∈π ∂ g x ∂ g y D π\b∪b1(x,y)∪b2(x,y) , the derivation of which comes at the end of this section. Let us explain the arising terms: • The first term is a sum over all blocks d of the dual partition π * , from which we choose all possible tuples (x, y) with x = y. The symbol ∂ g x means that the derivative only acts on the g's and not the delta functions that appear in D π , while ∂ x acts on both.
• The second term is a sum over all blocks b ∈ π, from which we choose all possible tuples of edges (x, y) with x = y. By the convention (E4), we denote the neighbouring nodes by the same name. This allows us to cut the block b of edges along the nodes (x, y). The two resulting blocks are denoted by b 1 (x, y) and b 2 (x, y). The partition π\b ∪ b 1 (x, y) ∪ b 2 (x, y) is the one where b was removed and replaced by the two blocks b 1 (x, y) and b 2 (x, y).
Note that the term with ∂ x ∂ y is cancelled once we plug (E9) into (E8 Let us explain this: Since π = {[n]} one can always find a block d ∈ π * that consists of at least two nodes (x, y). We can join the two blocks b(x) ∈ π and b(y) ∈ π to which the corresponding edges x and y belong. This forms a block b of a new partition π which differs from π only by this block, π \b ∪ b(x) ∪ b(y) = π. Note, that if π consists of only two blocks, then π will consist of a single block, so π = {[n]}. There is hence a bijection {(D π , x, y) | π ∈ N C([n])\[n], (x, y) ∈ d st. d ∈ π * } ↔ {π \b ∪ b 1 (x, y) ∪ b 2 (x, y) | π ∈ N C[n], (x, y) ∈ b ∈ π }. But π = [n] is not available in the above sum, therefore leaving all those terms δ(r 1 , s 1 )∂ r1 g t (r)∂ s1 g t (s), which is what we wanted to show. Proof of (E9). Recall that our notation (x, y) ∈ b, for b any block of a partition, assumes that x = y. Two identities about δ-functions of several variables that we will need are This expression can be written as a total derivative of ∂ x ∂ y up to the missing term ∂ g x ∂ g y , which, rewritten in terms of D π , is the first term that appears in (E9), For the second term we need to use (E3), which by assumption holds for all k ≤ n − 1. Since π ∈ N C([n])\[n] the assumption applies, c∈π r s=c δ(r 1 , s 1 )∂ r1 g t (r)∂ s1 g t (s) b∈π\c g t (b).
Instead of summing over r s = c, we can also sum over (x, y) ∈ c (where as usual x = 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), c∈π (x,y)∈c δ(x, y)∂ x g t (b 1 (x, y))∂ y g t (b 2 (x, y)) b∈π\c g t (b).
This concludes the derivation of the time evolution equation of connected correlation functions.

Time evolution of the new measure
Here, we present the proof of the evolution equation (74) of the new measure ϕ t . We follow the explanation before (E3) to rewrite the evolution of the connected correlations g t as (∂ t − ∆)g t ([n]) = r s= [n] δ(r 1 , s 1 )∂ r1 g t (r)∂ s1 g t (s), where r s = [n] denotes the union of non-crossing subsets r = {r 1 , r 2 , · · · } and s = {s 1 , s 2 , · · · } of the set [n] = {x 1 , · · · , x n }. The proof of (74) is by induction over n. For n = 1, ϕ t and g t are identical and therefore satisfy the same equation. Assume the formula holds for all n ≤ k − 1 for some k ∈ N. We start by evaluating its left hand side making use of (E14) (for better readability we suppress write the time argument) = π∈N C([n]) c∈π a∪a =c δ(a 1 , a 1 )∂ a1 g t (a)∂ a 1 g t (a )g π\c .
Here ∆ c = x∈c ∆ x and π\c is the partition π without the block c.