Entanglement distribution in the Quantum Symmetric Simple Exclusion Process

We study the probability distribution of entanglement in the Quantum Symmetric Simple Exclusion Process, a model of fermions hopping with random Brownian amplitudes between neighboring sites. We consider a protocol where the system is initialized in a pure product state of $M$ particles, and focus on the late-time distribution of R\'enyi-$q$ entropies for a subsystem of size $\ell$. By means of a Coulomb gas approach from Random Matrix Theory, we compute analytically the large-deviation function of the entropy in the thermodynamic limit. For $q>1$, we show that, depending on the value of the ratio $\ell/M$, the entropy distribution displays either two or three distinct regimes, ranging from low- to high-entanglement. These are connected by points where the probability density features singularities in its third derivative, which can be understood in terms of a transition in the corresponding charge density of the Coulomb gas. Our analytic results are supported by numerical Monte Carlo simulations.

The relevance of stochastic models for generic systems relies on the assumption that the properties of individual random realizations are close to the averaged ones. While this is often a natural expectation, it is typically difficult to obtain quantitative results on the full probability distribution of coherent phenomena such as quantum entanglement [4,8,[29][30][31]. At the same time, understanding the nature of fluctuations is clearly an important task, and a necessary step towards the generalization of powerful methods developed for classical stochastic systems, such as the well-established Macroscopic Fluctuation Theory [32,33].
Here, we initiate a series of investigations aimed at understanding entanglement fluctuations in a prototypical model for quantum many-body stochastic dynamics: the Quantum Simple Symmetric Exclusion Process (Q-SSEP), cf. Fig. 1. This model, recently introduced in Refs. [22,34], describes fermions hopping with random amplitudes between neighboring sites, and is particularly useful from the theoretical point of view. On the one hand, given the quadratic form of the Hamiltonian gener- ator, it allows us to employ analytic techniques which are not available in other models. On the other hand, while its mean dynamics reduces to the classical SSEP [35][36][37][38][39][40][41][42][43], quantum coherent effects were shown to display a rich phenomenology in this system and its generalizations [40,44,45], making it an ideal toy model to build a quantitative understanding of quantum fluctuations.
We focus on the simplest setting where the system is initialized in a pure product state, and compute the largedeviation function for the Rényi-q entropy of subsystems at late times. Using the Coulomb gas (CG) approach from Random Matrix Theory (RMT), we find that it displays distinct phases, with two of them corresponding to states approaching either a pure state or the maximally mixed one (defining regimes of low and high entanglement, respectively). These regimes are separated by critical points where the probability density features singularities in its third derivative and which can be understood in terms of a transition in the corresponding charge density of the CG. Our results are supported by numerical Monte Carlo simulations, and open the way towards further studies of fluctuations of entanglementrelated quantities in the Q-SSEP, and its generalizations.
The rest of this article is organized as follows. In Sec. II we introduced the Q-SSEP and review previous results on the characterization of the stationary state approached at late times. In Sec. III we lay out the Coulomb-gas approach to the computation of the large-deviation function. We derive a set of equations whose exact solution is presented in Sec. IV. Finally, our conclusions are reported in Sec. V.

II. THE MODEL
We consider a chain of L sites with periodic boundary conditions. The Q-SSEP is formally defined by the Hamiltonian generator where c j , c † j are canonical fermionic operators, with {c j , c † k } = δ j,k , and W j t ,W j t are pairs of complex conjugated Brownian motions. They satisfy dW j t dW k t = δ j,k dt, and dW j t dW k t = dW k t dW j t = 0 for t = t , where we used the standard notation in Itô calculus [46]. The system is initialized in a pure product state of M particles. Late-time properties turn out to be independent of the specific initial state chosen, but for concreteness we may take |Ψ(0) = c † 1 · · · c † M |0 , where |0 is the vacuum. We consider the entanglement of a subsystem A = {1, . . . , }, as measured by the Rényi-q entropies where ρ (t) = tr L\ |Ψ(t) Ψ(t)|. Clearly, s q (t) = S q (t)/ is a stochastic variable distributed according to some probability density p q,t (s), with 0 ≤ s ≤ ln 2. Our goal is to compute the full distribution of s q (t) for large times, namely p q (s) = lim t→∞ p q,t (s), in the limit of large L, , M , where we fix the ratios ξ = /L, m = M/L. As a preliminary observation, note that the initial state |Ψ(0) satisfies Wick's theorem, and its density matrix is completely specified by its covariance matrix (G 0 ) i,j := Ψ(0)|c † i c j |Ψ(0) . Since the Hamiltonian is quadratic, this remains true for the evolved state |Ψ(t) , and the system is effectively described by the evolved covariance matrix G t . The latter also fully determines the value of the Rényi entropies S q (t) [47]: denoting by A ( ) the matrix obtained by selecting the first rows and columns of a matrix A, we have where {λ j } j=1 are the eigenvalues of G In order to make progress, we use that the density p q (s) satisfies a large-deviation principle; in particular, we will prove that, for < L, ln p q (s) ∼ − 2 I q (s), for some rate function I q (s). In this situation, the Gärtner-Ellis theorem applies [48], stating that I q (s) can be computed from the knowledge of the cumulant generating function by a Legendre transform: where we can make use of a result derived in Ref. [34], relating large-time expectation values to averages over the unitary group U (L) equipped with the Haar invariant measure. Explicitly, we obtain with G V = V † G 0 V , and where dη(V ) denotes the Haar measure over U (L). It follows from Eqs. (4) and (6) that the problem is reduced to computing the distribution of the subsystem entanglement for a random pure fermionic Gaussian state. It is important to stress that this is different from the analogous problem for Haar random states sampled over the whole many-body space (having dimension 2 L ). In that case several exact results were obtained for the full probability distribution of entanglement [49][50][51][52][53][54][55]. While we will employ similar techniques, qualitative and quantitative differences arise in our case.

III. THE COULOMB GAS APPROACH
The Haar measure over U (L), induces a probability distribution P [{λ j }] on the set of eigenvalues of G ( ) V , which allows us to express Eq. (6) in the form where is the M × sub-matrix containing the first M rows and columns of V . Thus, in order to evaluate (7), we need the probability distribution induced on the eigenvalues of V † M, V M, , when V is sampled from the Haar invariant measure. It turns out that the latter is known in RMT [56], and takes the form has been recently exploited for the computation of averaged subsystem entanglement in the context of random non-interacting fermionic ensembles [57][58][59][60], see also [61]. Note that the distribution depends on both and M . In the following, we may restrict to ≤ L/2, since the entanglement for pure states is symmetric under → L − . Furthermore, we may also choose ≤ M [62].
Following Refs. [51][52][53][54][55], we use Eq. (8) as the starting point of our computations, which are based on the CG approach. This is a method routinely applied in RMT, consisting in a mapping between random matrix eigenvalues and repulsive point charges [56]. The CG analysis of the Jacobi ensemble has been already employed in physical problems, e.g. to study the conductance and the shot noise power for a mesoscopic cavity with two leads [63,64], or to compute the so-called Andreev conductance of a metal-superconductor interface [65]. In order to see how it works, we rewrite with Within the CG approach, the function E w [{λ i }] is interpreted as the energy of a gas of charged particles with coordinates λ j ∈ [0, 1], which are subject to an external potential. The integral (9) is then understood as a thermal partition function for the CG. In the largelimit, the configuration of the Coulomb charges may be described in terms of the normalized density ρ(λ, ) = where the denominator corresponds to the normalization constant N . To the leading order in [66], where we introduced the Lagrange multiplier u enforcing normalization, and the effective potential where m, ξ are the density of fermions and the rescaled interval length introduced before. The functional integrals in (10) may be evaluated by the saddle-point method.

This yields
Differentiating the last equation with respect to µ we arrive at where − denotes the principal-value integral. Eq. (14) can be formally solved using the so-called Tricomi's formula [67,68], or a resolvent method [69,70], which both yield integral representations of the solution which can in general be evaluated numerically. Plugging which is derived from the saddle-point method, into Eq. (4), this finally allows us to obtain a numerical value for I q (s). In fact, we find that Eq. (14) can be solved analytically for all integers q > 1. Before discussing the mathematical details, however, it is interesting to observe that its qualitative features can be understood based on the analysis of the CG picture, as we now briefly discuss. First of all, we note that for 0 ≤ ξ ≤ m ≤ 1/2, the effective potential (12) is always bounded from below. Furthermore, for w negative and with large absolute value, V w (λ) has a single local minimum close to λ = 1/2, cf. Fig. 2. Recalling that ρ * w (λ) describes the distribution of charges subject to the external potential V w (λ), we expect ρ * w (λ) to develop an increasingly sharp peak around this point. This is consistent with our intuition based on the quantum problem: for w → −∞, maximal entanglement entropies are favored in the average corresponding to F q (w), and the most significant states in the average approach the maximally mixed one, i.e. all the eigenvalues of the covariance matrix should be close to λ = 1/2. For w very large, instead, V w (λ) develops a local maximum close to λ = 1/2, and the Coulomb charges are pushed at the boundaries of [0, 1], eventually depleting its central region. Accordingly, we expect ρ * w (λ) to become peaked around λ = 0 and λ = 1, and vanishing in a neighborhood of λ = 1/2. In terms of the quantum problem, this means that all eigenvalues of the covariance matrix are close to 0 or 1, i.e. the entanglement vanishes, and we approach a pure state. We will see that the two limits w → ±∞ correspond to different phases of the rate function.

IV. THE EXACT SOLUTION
We now present our analytic solution to Eq. (14). While we were able to obtain explicit expressions for all integers q > 1, and arbitrary ξ and m, they are very cumbersome for general q and ξ, m < 1/2. For this reason, here we report only the case q = 2 and ξ = m = 1/2. Furthermore, we will only present the final result of our analysis, while all the details of our derivations will be reported elsewhere [71].
In general, we find that ρ * w (λ) displays either two or three distinct phases as a function of w, separated by points where I q (s) develops a discontinuity in its third derivative. Similar kinds of "third order phase transitions" are ubiquitous in RMT, appearing in a wide variety of contexts [68]. In our case, for m = ξ = 1/2, and q = 2, there are three phases, separated by the points w * 1 = −2 − √ 2 and w * 2 = 1 + √ 2. The first one is characterized by states with large entanglement, and corresponds to −∞ < w < w * 1 . In this case, ρ * w (λ) has nonzero support over the interval J I = [ν − , ν + ] ⊂ [0, 1], with ν ± = (1 ± ν)/2 and ν = − √ −2w − 1/(w + 1). It reads As expected, ρ * w (λ) becomes a delta-function peaked around λ = 1/2 for w → −∞. Next, for w * 1 < w < w * 2 , we enter a transition regime: ρ * w (λ) has support over the whole interval J II = (0, 1), and develops two integral singularities at its boundaries. It reads with Note that for w = 0 we recover the spectral density of the Jacobi ensemble, see e.g. [56,72,73]. As w varies from w * 1 to w * 2 , the charge density decreases at the center of the interval, eventually vanishing in λ = 1/2 at w = w * 2 . Here, we enter the third phase, spanning w * 2 < w < ∞, which is that of low-entangled states. In this regime, ρ * w (λ) has non-vanishing support over J III = (0, ν − ) ∪ (ν + , 1), with ν ± = (1 ± ν)/2 and It has the form with Importantly, we see that as w → ∞ the support of ρ * w (λ) localizes around 0 and 1, yielding vanishing entanglement. We plot the optimal density ρ * w (λ) in Fig. 2, for three values of w corresponding to the phases discussed above.
Let us also mention how this picture is modified when ξ < 1/2 (and m = 1/2). In this case, the potential V w (λ) is divergent at λ = 0, 1, and the support of the optimal charge ρ * w (λ) is strictly contained in [0, 1]. Accordingly, we find that phases I and II merge, so that ρ * w (λ) only displays two phases, separated by the point The qualitative features of the optimal distributions remain the same, although they do not display singularities at the boundaries of their support for ξ = 1/2. From the knowledge of ρ * w (λ), we can compute the rate function I q (s). First, it is convenient to rewrite the Legendre transform (4) as where w s satisfies df q (w s )/dw = s. Using Eq. (15), and the fact that ρ * w (λ) is the saddle point of E w [ρ], this condition is equivalent to S q [ρ * ws ] = s, where From Eq. (22) we see that I q (s) can be computed by evaluating numerically simple integrals [74]. We followed this procedure to generate plots of the function I q (s) for different values of ξ, as reported in Fig. 3. As a general feature, we see that the rate function develops singularities at s = 0, s = ln 2. We also note that we may read off the average value for the entropy, corresponding to the minimum of I q (s).
To obtain an analytic form for I q (s), one should invert the relation S q [ρ * ws ] = s, and express w s as a function of s. While this is difficult for general values of q, ξ and m, due to the complicated form of ρ * w (λ), it may be done in some cases. In particular, fully analytic results can be obtained for q = 2, ξ = m = 1/2. In this case, I 2 (s) can be written explicitly in phase II, displaying the simple form wheres = 3 ln 2 − 2 ln(1 + √ 2) is the average Rényi-2 entropy, while γ 0.06 is a numerical constant. Hence, for w ∈ (w * 1 , w * 2 ) the probability density for the Rényi-2 entropy is simply Gaussian. In phase I and III, instead, a large-w expansion reveals that I 2 (s) develops logarithmic singularities for s → 0 and s → ln 2: we find withs = 0,s = ln 2, respectively. We have tested our predictions against Monte Carlo simulations [75], numerically constructing a histogram of the probability p q (s) based on a sampling of (8). Since the distribution of the Rényi entropies is highly peaked around its average, a standard Metropolis approach is not adequate to efficiently explore a wide range of its values, and we implemented the numerical scheme introduced in Ref. [53], where one forces the Metropolis algorithm to explore regions of large values of the Rényi entropy. As explained in [53], this method gives us access to the derivative of the rate function I 2 (s) for finite systems [71]. The numerical data obtained using this method are reported in Fig. 3 for the case q = m −1 = 2, and different values of ξ. The plot shows excellent agreement with our predictions, revealing that finite-size effects are very small for the set of parameters considered.

V. CONCLUSIONS
We have computed the large-deviation function for the entanglement of subsystems in the steady state of the Q-SSEP. We have shown that its distribution is characterized by different phases connected by points where the probability density features singularities in its third derivative. Our work raises several questions. First, it would be interesting to understand how our predictions are modified for suitable generalizations of the model, such as the Q-SSEP with dissipative boundaries [40,44], or its "asymmetric" version [45]. Furthermore, a natural direction to explore pertains to the dynamics of entanglement, which should be in principle accessible from the stochastic equations of motion studied in [34]. These questions are left for future work.