# Entanglement distribution in the quantum symmetric simple exclusion process

Denis Bernard, Lorenzo Piroli

► To cite this version:

Denis Bernard, Lorenzo Piroli. Entanglement distribution in the quantum symmetric simple exclusion process. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics [2001-2015], 2021, 104 (1), pp.014146. 10.1103/PhysRevE.104.014146. hal-03157129

HAL Id: hal-03157129

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

Submitted on 26 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Entanglement distribution in the quantum symmetric simple exclusion process

Denis Bernard¹ and Lorenzo Piroli²,³

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

²Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Straße 1, 85748 Garching, Germany

³Munich Center for Quantum Science and Technology, Schellingstraße 4, 80799 München, Germany

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

(Received 17 February 2021; accepted 7 July 2021; published 29 July 2021)

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 we focus on the late-time distribution of Rényi-*q* entropies for a subsystem of size *l*. 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 *l*/*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.

DOI: 10.1103/PhysRevE.104.014146

## I. INTRODUCTION

Many physical phenomena admit a description in terms of random variables, whose dynamics is dictated by stochastic processes. While they have been traditionally introduced for open systems, where randomness is acquired through the interaction with the environment [1], stochastic processes have recently received renewed attention in connection with investigations of typical features of isolated many-body systems. This trend was driven by the study of random unitary circuits [2], which proved to be ideal toy models to investigate aspects of the dynamics that are notoriously hard to tackle, including entanglement growth [2–8], operator spreading [9–16], dynamical correlations [16–18], and scrambling of quantum information [19–21]. Similar ideas were also explored in the context of continuous-time Hamiltonian dynamics [22–27] and stochastic conformal field theories [28].

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–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 generator, 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–43], quantum coherent effects have been 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 large-deviation 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 a 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 entanglement-related quantities in the Q-SSEP and its generalizations.

The rest of this article is organized as follows. In Sec. II we introduce the Q-SSEP and review previous results on

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Open access publication funded by the Max Planck Society.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Pictorial representation of the Q-SSEP. M fermionic particles, initially in a pure product state, hop with random amplitudes between neighboring sites. We focus on the probability distribution of the entanglement of a subsystem $\mathcal{A}_\ell$, containing $\ell$ neighboring sites.</figcaption></figure>

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

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

where $c_j$ and $c_j^\dagger$ are canonical fermionic operators, with $\{c_j, c_k^\dagger\} = \delta_{j,k}$, and $W_t^j$ and $\bar{W}_t^j$ are pairs of complex-conjugated Brownian motions. They satisfy $dW_t^j d\bar{W}_t^k = \delta_{j,k}dt$ and $d\bar{W}_t^j dW_t^k = dW_t^k d\bar{W}_t^j = 0$ for $t \neq 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 $|\Psi(0)\rangle = c_1^\dagger \cdots c_M^\dagger |0\rangle$, where $|0\rangle$ is the vacuum. We consider the entanglement of a subsystem $\mathcal{A}_\ell = \{1, \dots, \ell\}$, as measured by the Rényi-$q$ entropies

$$S_q(t) = (1-q)^{-1} \ln \text{tr}[\rho_{\ell}^q(t)], \quad (2)$$

where $\rho_\ell(t) = \text{tr}_{L\setminus\ell}|\Psi(t)\rangle\langle\Psi(t)|$. Clearly, $s_q(t) = S_q(t)/\ell$ is a stochastic variable distributed according to some probability density $p_{q,t}(s)$, with $0 \le s \le \ln 2$. Our goal is to compute the full distribution of $s_q(t)$ for large times, namely, $p_q(s) = \lim_{t\to\infty} p_{q,t}(s)$, in the limit of large $L$, $\ell$, and $M$, where we fix the ratios $\xi = \ell/L$ and $m = M/L$.

As a preliminary observation, note that the initial state $|\Psi(0)\rangle$ satisfies Wick's theorem, and its density matrix is completely specified by its covariance matrix $(G_0)_{i,j} := \langle\Psi(0)|c_i^\dagger c_j|\Psi(0)\rangle$. Since the Hamiltonian is quadratic, this remains true for the evolved state $|\Psi(t)\rangle$, 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^{(\ell)}$ the matrix obtained by selecting the first $\ell$ rows and columns of a matrix $A$, we have

$$S_q(t) = (1-q)^{-1} \sum_{j=1}^{\ell} \ln[\lambda_j^q + (1-\lambda_j)^q], \quad (3)$$

where $\{\lambda_j\}_{j=1}^\ell$ are the eigenvalues of $G_t^{(\ell)}$, satisfying $0 \le \lambda_j \le 1$.

In order to make progress, we use that the density $p_q(s)$ satisfies a large-deviation principle; in particular, we prove that, for $\ell < L$, $\ln p_q(s) \sim -\ell^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:

$$I_q(s) = -\inf_w \{ws - f_q(w)\}, \quad (4)$$

where we introduced $f_q(w) = -\ell^{-2} \ln \mathcal{F}_q(w)$, with $\mathcal{F}_q(w) = \lim_{t\to\infty} \mathbb{E}_t[e^{-w\ell S_q(t)}]$. Writing $S_q(t) = S_{q,\ell}[G_t]$, where we defined the function

$$S_{q,\ell}[G_t] = (1-q)^{-1} \text{tr} \ln [(G_t^{(\ell)})^q + (1-G_t^{(\ell)})^q], \quad (5)$$

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

$$\mathcal{F}_q(w) = \int_{U(L)} d\eta(V) \exp(-w\ell S_{q,\ell}[G_V]), \quad (6)$$

with $G_V = V^\dagger G_0 V$ and where $d\eta(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–55]. While we 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[\{\lambda_j\}]$ on the set of eigenvalues of $G_V^{(\ell)}$, which allows us to express Eq. (6) in the form

$$\mathcal{F}_q(w) = \int \left( \prod_{j=1}^{\ell} d\lambda_j \right) P[\{\lambda_j\}] e^{-w\ell S_q[\{\lambda_j\}]}, \quad (7)$$

where $S_q[\{\lambda_j\}] = (1-q)^{-1} \sum_j \ln[\lambda_j^q + (1-\lambda_j)^q]$. For the initial state chosen, simple manipulations give $G_V^{(\ell)} = V_{M,\ell}^\dagger V_{M,\ell}$, where $V_{M,\ell}$ is the $M \times \ell$ submatrix containing the first $M$ rows and $\ell$ columns of $V$. Thus, in order to evaluate Eq. (7), we need the probability distribution induced on the eigenvalues of $V_{M,\ell}^\dagger V_{M,\ell}$, when $V$ is sampled from the Haar invariant measure. It turns out that the latter is known in RMT [56] and takes the form

$$P[\{\lambda_i\}] = \frac{1}{N} \prod_{j<k} |\lambda_j - \lambda_k|^2 \prod_{i=1}^{\ell} \lambda_i^{M-\ell} (1-\lambda_i)^{L-\ell-M}, \quad (8)$$

where $N$ is a normalization constant. This distribution defines the so-called $\beta$-Jacobi ensemble (with $\beta = 2$) and has been recently exploited for the computation of averaged subsystem entanglement in the context of random noninteracting
fermionic ensembles [57–60] (see also Ref. [61]). Note that the distribution depends on both $\ell$ and $M$. In the following, we may restrict to $\ell \le L/2$, since the entanglement for pure states is symmetric under $\ell \mapsto L - \ell$. Furthermore, we may also choose $\ell \le M$ [62].

Following Refs. [51–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 of 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

$$ \mathcal{F}_q(w) = \frac{1}{N} \int_0^1 \left( \prod_{j=1}^{\ell} d\lambda_j \right) e^{-\ell^2 E_w[\{\lambda_j\}]}, \quad (9) $$

with

$$ E_w[\{\lambda_i\}] = -\frac{2}{\ell^2} \sum_{i<j} \ln|\lambda_i - \lambda_j| - \frac{(M-\ell)}{\ell^2} \sum_i \ln \lambda_i - \frac{(L-M-\ell)}{\ell^2} \sum_i \ln(1-\lambda_i) + wS_q[\{\lambda_i\}]/\ell. $$

Within the CG approach, the function $E_w[\{\lambda_i\}]$ is interpreted as the energy of a gas of charged particles with coordinates $\lambda_j \in [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 large-$\ell$ limit, the configuration of the Coulomb charges may be described in terms of the normalized density $\rho(\lambda, \ell) = \ell^{-1} \sum_j \delta(\lambda - \lambda_i)$, and the multiple integral in Eq. (9) can be replaced by an integral over all possible density functions $\rho(\lambda)$, i.e.,

$$ \mathcal{F}_q(w) = \frac{\int \mathcal{D}\rho e^{-\ell^2 E_w[\rho]}}{\int \mathcal{D}\rho e^{-\ell^2 E_0[\rho]}}, \quad (10) $$

where the denominator corresponds to the normalization constant $N$. To the leading order in $\ell$ [66], $E_w[\rho]$ reads

$$ E_w[\rho] = -\int_0^1 d\lambda \int_0^1 d\mu \rho(\lambda)\rho(\mu) \ln|\lambda - \mu| + \int_0^1 d\lambda \rho(\lambda) V_w(\lambda) + u \left\{ \int_0^1 d\lambda \rho(\lambda) - 1 \right\}, \quad (11) $$

where we introduce the Lagrange multiplier $u$ enforcing normalization and the effective potential

$$ V_w(\lambda) = -\left(\frac{m}{\xi} - 1\right) \ln \lambda - \left(\frac{1-m}{\xi} - 1\right) \ln(1-\lambda) + \frac{w}{1-q} \ln[\lambda^q + (1-\lambda)^q], \quad (12) $$

where $m$ and $\xi$ are the density of fermions and the rescaled interval length introduced before. The functional integrals in Eq. (10) may be evaluated by the saddle-point method. This yields $\int \mathcal{D}\rho e^{-\ell^2 E_w[\rho]} \sim e^{-\ell^2 E_w[\rho^*]}$, where $\rho^*_w(\lambda)$ is the “optimal” charge density, minimizing $E_w[\rho^*_w(\lambda)]$

<figure><img src="image_4.png" /><figcaption>FIG. 2. (a) Shifted effective potential $V_w(\lambda)$, defined in Eq. (12), for different values of $w$ and $m = 0.5$, $\xi = 0.25$, and $q = 2$. (b) Optimal charge distribution for $q = 2$ and $m = \xi = 0.5$. The plot shows $\rho^*_w(\lambda)$ for three values of $w$, each corresponding to a different regime.</figcaption></figure>

satisfying

$$ \int d\lambda \rho_w^*(\lambda) \ln|\mu - \lambda| = (1/2)V(\mu) + u/2. \quad (13) $$

Differentiating the last equation with respect to $\mu$, we arrive at

$$ \int d\lambda \frac{\rho_w^*(\lambda)}{\mu - \lambda} = \frac{1}{2} V_w'(\mu), \quad (14) $$

where $f$ denotes the principal-value integral. Equation (14) can be formally solved using the so-called Tricomi's formula [68,69] or a resolvent method [70,71], which both yield integral representations of the solution which can, in general, be evaluated numerically. Plugging

$$ f_q(w) = E_w[\rho_w^*] - E_0[\rho_0^*], \quad (15) $$

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 \le \xi \le m \le 1/2$, the effective potential (12) is always bounded from below. Furthermore, for $w$ negative and with large absolute value, $V_w(\lambda)$ has a single local minimum close to $\lambda = 1/2$ (cf. Fig. 2). Recalling that $\rho_w^*(\lambda)$ describes the distribution of charges subject to the external potential $V_w(\lambda)$, we expect $\rho_w^*(\lambda)$ to develop an increasingly sharp peak around this point. This is consistent with our intuition based on the quantum problem: for $w \to -\infty$, maximal entanglement entropies are favored in the average corresponding to $\mathcal{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 $\lambda = 1/2$. For $w$ very large, instead, $V_w(\lambda)$ develops a local maximum close to $\lambda = 1/2$, and the Coulomb charges are pushed at the boundaries of [0,1], eventually depleting its central region. Accordingly, we expect $\rho_w^*(\lambda)$ to become peaked around $\lambda = 0$ and $\lambda = 1$ and vanish in the neighborhood of $\lambda = 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 \to \pm\infty$ 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 $\xi$ and $m$, they are very cumbersome for general $q$ and $\xi$, $m < 1/2$. For this reason, here we report only the case $q = 2$ and $\xi = m = 1/2$. Furthermore, we only present the final result of our analysis, while all the details of our derivations will be reported elsewhere [72].

In general, we find that $\rho_w^*(\lambda)$ 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 [69]. In our case, for $m = \xi = 1/2$ and $q = 2$, there are three phases, separated by the points $w_1^* = -2 - \sqrt{2}$ and $w_2^* = 1 + \sqrt{2}$. The first one is characterized by states with large entanglement and corresponds to $-\infty < w < w_1^*$. In this case, $\rho_w^*(\lambda)$ has nonzero support over the interval $J^1 = [\nu_-, \nu_+] \subset [0, 1]$, with $\nu_\pm = (1 \pm \nu)/2$ and $\nu = -\sqrt{-2w - 1/(w+1)}$. It reads

$$ \rho_w^*(\lambda) = -(2/\pi)(w+1) \frac{\sqrt{\lambda - \nu_-} \sqrt{\nu_+ - \lambda}}{\lambda^2 + (1-\lambda)^2}. \quad (16) $$

As expected, $\rho_w^*(\lambda)$ becomes a $\delta$ function peaked around $\lambda = 1/2$ for $w \to -\infty$. Next, for $w_1^* < w < w_2^*$, we enter a transition regime: $\rho_w^*(\lambda)$ has support over the whole interval $J^{II} = (0, 1)$ and develops two integral singularities at its boundaries. It reads

$$ \rho_w^*(\lambda) = \frac{1}{\sqrt{\pi^2 \lambda (1-\lambda)}} g(\lambda), \quad (17) $$

with

$$ g(\lambda) = 1 + w\{1 - 2^{-1/2}[\lambda^2 + (\lambda - 1)^2]^{-1}\}. \quad (18) $$

Note that for $w = 0$ we recover the spectral density of the Jacobi ensemble (see, e.g., Refs. [56,73,74]). As $w$ varies from $w_1^*$ to $w_2^*$, the charge density decreases at the center of the interval, eventually vanishing in $\lambda = 1/2$ at $w = w_2^*$. Here, we enter the third phase, spanning $w_2^* < w < \infty$, which is that of low-entangled states. In this regime, $\rho_w^*(\lambda)$ has nonvanishing support over $J^{III} = (0, \nu_-) \cup (\nu_+, 1)$, with $\nu_\pm = (1 \pm \nu)/2$ and

$$ \nu = \frac{\sqrt{(w-1)^2 - 2}}{w+1}. \quad (19) $$

It has the form

$$ \rho_w^*(\lambda) = \frac{|\lambda - 1/2| [2 + wh(\lambda)]}{\pi \sqrt{(1-\lambda)\lambda} \sqrt{(1-2\lambda)^2 - \nu(w)^2}}, \quad (20) $$

with $h(\lambda) = 2 - \sqrt{2 + 2\nu(w)^2}/[1 - 2(1 - \lambda)\lambda]$. Importantly, we see that as $w \to \infty$ the support of $\rho_w^*(\lambda)$ localizes around 0 and 1, yielding vanishing entanglement. We plot the optimal density $\rho_w^*(\lambda)$ in Fig. 2, for three values of $w$ corresponding to the phases discussed above.

Let us also mention how this picture is modified when $\xi < 1/2$ (and $m = 1/2$). In this case, the potential $V_w(\lambda)$ is divergent at $\lambda = 0$ and 1, and the support of the optimal charge $\rho_w^*(\lambda)$ is strictly contained in $[0,1]$. Accordingly, we find that phases I and II merge, so that $\rho_w^*(\lambda)$ only displays two phases,

<figure><img src="image_5.png" /><figcaption>FIG. 3. (a) Rate function for the Rényi-2 entropy, for $m = 1/2$ and different values of $\xi$. For $\xi < m$ and $\xi = m$, respectively, two and three phases appear, which correspond to different colors. (b) Analytic predictions for the derivative $I_2'(s)$ (solid lines), against numerical data from Monte Carlo simulations for $\ell = 32$, $L = \ell/\xi$, and $M = L/2$ (dots). The numerical error is not visible at the scales of the plot.</figcaption></figure>

separated by the point

$$ w^*(\xi) = \frac{1 + 2\sqrt{2}\sqrt{(1-\xi)\xi}}{2\xi}. \quad (21) $$

The qualitative features of the optimal distributions remain the same, although they do not display singularities at the boundaries of their support for $\xi \neq 1/2$.

From the knowledge of $\rho_w^*(\lambda)$, we can compute the rate function $I_q(s)$. First, it is convenient to rewrite the Legendre transform (4) as

$$ I_q(s) = -w_s s + f_q(w_s), \quad (22) $$

where $w_s$ satisfies $df_q(w_s)/dw = s$. Using Eq. (15) and the fact that $\rho_w^*(\lambda)$ is the saddle point of $E_w[\rho]$, this condition is equivalent to $S_q[\rho_w^*] = s$, where

$$ S_q[\rho] = \int d\lambda \, \rho(\lambda) (1-q)^{-1} \ln[\lambda^q + (1-\lambda)^q]. \quad (23) $$

From Eq. (22) we see that $I_q(s)$ can be computed by evaluating numerically simple integrals [75]. We followed this procedure to generate plots of the function $I_q(s)$ for different values of $\xi$, as reported in Fig. 3. As a general feature, we see that the rate function develops singularities at $s = 0$ and $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[\rho_w^*] = s$ and express $w_s$ as a function of $s$. While this is difficult for general values of $q$, $\xi$, and $m$, due to the complicated form of $\rho_w^*(\lambda)$, it may be done in some cases. In particular, fully analytic results can be obtained for $q = 2$ and $\xi = m = 1/2$. In this case, $I_2(s)$ can be written explicitly in phase II, displaying the simple form

$$ I_2(s) = \frac{(s - \bar{s})^2}{2\gamma}, \quad (24) $$

where $\bar{s} = 3 \ln 2 - 2 \ln(1 + \sqrt{2})$ is the average Rényi-2 entropy, while $\gamma \simeq 0.06$ is a numerical constant. Hence, for $w \in (w_1^*, w_2^*)$ the probability density for the Rényi-2 entropy is simply Gaussian. In phases I and III, instead, a large-$w$ expansion reveals that $I_2(s)$ develops logarithmic singularities for $s \to 0$ and $s \to \ln 2$: we find

$$ I_2(s) = -\frac{1}{2} \ln |s - \bar{s}| + O(|s - \bar{s}|), \quad (25) $$

with $\bar{s} = 0$ and $\bar{s} = \ln 2$, respectively.
We have tested our predictions against Monte Carlo sim-
ulations [76], numerically constructing a histogram of the
probability $p_q(s)$ based on a sampling of Eq. (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 imple-
ment 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 Ref. [53], this
method gives us access to the derivative of the rate function
$I'_2(s)$ for finite systems [72]. The numerical data obtained
using this method are reported in Fig. 3 for the case $q =
m^{-1} = 2$ and different values of $\xi$. The plot shows excellent
agreement with our predictions, revealing that finite-size ef-
fects 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 differ-
ent 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 bound-
aries [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 Ref. [34]. These
questions are left for future work.

ACKNOWLEDGMENTS

We are very grateful to Bruno Bertini and Satya N. Majumdar for reading the manuscript and for valuable comments. D.B. acknowledges useful discussions with M. Bauer and J.-B. Zuber. D.B. acknowledges support from the CNRS, while L.P. acknowledges support from the Alexander von Humboldt Foundation.

[1] H.-P. Breuer and F. Petruccione, *The Theory of Open Quantum Systems* (Clarendon Press, Oxford, 2002).

[2] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Phys. Rev. X 7, 031016 (2017).

[3] T. Rakovszky, F. Pollmann, and C. W. von Keyserlingk, Phys. Rev. Lett. 122, 250602 (2019).

[4] T. Zhou and A. Nahum, Phys. Rev. B 99, 174205 (2019).

[5] M. J. Gullans and D. A. Huse, Phys. Rev. X 9, 021007 (2019).

[6] M. Znidaric, Commun. Phys. 3, 100 (2020).

[7] Y. Huang, IOP SciNotes 1, 035205 (2020).

[8] T. Zhou and A. Nahum, Phys. Rev. X 10, 031066 (2020).

[9] A. Nahum, S. Vijay, and J. Haah, Phys. Rev. X 8, 021014 (2018).

[10] C. W. von Keyserlingk, T. Rakovszky, F. Pollmann, and S. L. Sondhi, Phys. Rev. X 8, 021013 (2018).

[11] A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev. X 8, 041019 (2018).

[12] C. Sünderhauf, D. Pérez-García, D. A. Huse, N. Schuch, and J. I. Cirac, Phys. Rev. B 98, 134204 (2018).

[13] T. Rakovszky, F. Pollmann, and C. W. von Keyserlingk, Phys. Rev. X 8, 031058 (2018).

[14] V. Khemani, A. Vishwanath, and D. A. Huse, Phys. Rev. X 8, 031057 (2018).

[15] N. Hunter-Jones, arXiv:1812.08219.

[16] A. J. Friedman, A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev. Lett. 123, 210603 (2019).

[17] P. Kos, B. Bertini, and T. Prosen, Phys. Rev. X 11, 011022 (2021).

[18] P. Kos, B. Bertini, and T. Prosen, Phys. Rev. Lett. 126, 190601 (2021).

[19] P. Hosur, X.-L. Qi, D. A. Roberts, and B. Yoshida, J. High Energy Phys. 02 (2016) 004.

[20] B. Bertini and L. Piroli, Phys. Rev. B 102, 064305 (2020).

[21] L. Piroli, C. Sünderhauf, and X.-L. Qi, J. High Energy Phys. 04 (2020) 063.

[22] M. Bauer, D. Bernard, and T. Jin, SciPost Phys. 3, 033 (2017).

[23] E. Onorati, O. Buerschaper, M. Kliesch, W. Brown, A. Werner, and J. Eisert, Commun. Math. Phys. 355, 905 (2017).

[24] M. Knap, Phys. Rev. B 98, 184416 (2018).

[25] D. A. Rowlands and A. Lamacraft, Phys. Rev. B 98, 195125 (2018).

[26] T. Zhou and X. Chen, Phys. Rev. E 99, 052212 (2019).

[27] C. Sünderhauf, L. Piroli, X.-L. Qi, N. Schuch, and J. I. Cirac, J. High Energy Phys. 11 (2019) 038.

[28] D. Bernard and P. Le Doussal, Europhys. Lett. 131, 10007 (2020).

[29] J. Cotler, N. Hunter-Jones, and D. Ranard, arXiv:2010.11922.

[30] F. Carollo, R. L. Jack, and J. P. Garrahan, Phys. Rev. Lett. 122, 130605 (2019).

[31] F. Carollo and C. Pérez-Espigares, Phys. Rev. E 102, 030104(R) (2020).

[32] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Phys. Rev. Lett. 94, 030601 (2005).

[33] L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Rev. Mod. Phys. 87, 593 (2015).

[34] M. Bauer, D. Bernard, and T. Jin, SciPost Phys. 6, 45 (2019).

[35] C. Kipnis, S. Olla, and S. Varadhan, Commun. Pure Appl. Math. 42, 115 (1989).

[36] G. Eyink, J. L. Lebowitz, and H. Spohn, Commun. Math. Phys. 140, 119 (1991).

[37] B. Derrida, J. Stat. Mech. (2007) P07023.

[38] B. Derrida, J. Stat. Mech. (2011) P01030.

[39] K. Mallick, Phys. A (Amsterdam, Neth.) 418, 17 (2015).

[40] D. Bernard and T. Jin, Phys. Rev. Lett. 123, 080601 (2019).

[41] F. H. L. Essler and L. Piroli, Phys. Rev. E 102, 062210 (2020).

[42] R. Frassek, C. Giardina, and J. Kurchan, SciPost Phys. 9, 54 (2020).

[43] R. Frassek, C. Giardina, and J. Kurchan, SciPost Phys. 10, 135 (2021).

[44] D. Bernard and T. Jin, Commun. Math. Phys. 384, 1141 (2021).

[45] T. Jin, A. Krajenbrink, and D. Bernard, Phys. Rev. Lett. 125, 040603 (2020).
[46] B. Oksendal, *Stochastic Differential Equations: An Introduction with Applications* (Springer-Verlag, Berlin, 2003).

[47] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett. **90**, 227902 (2003).

[48] H. Touchette, Phys. Rep. **478**, 1 (2009).

[49] O. Giraud, J. Phys. A: Math. Theor. **40**, 2793 (2007).

[50] P. Facchi, U. Marzolino, G. Parisi, S. Pascazio, and A. Scardicchio, Phys. Rev. Lett. **101**, 050502 (2008).

[51] C. Nadal, S. N. Majumdar, and M. Vergassola, Phys. Rev. Lett. **104**, 110501 (2010).

[52] A. De Pasquale, P. Facchi, G. Parisi, S. Pascazio, and A. Scardicchio, Phys. Rev. A **81**, 052324 (2010).

[53] C. Nadal, S. N. Majumdar, and M. Vergassola, J. Stat. Phys. **142**, 403 (2011).

[54] P. Facchi, G. Florio, G. Parisi, S. Pascazio, and K. Yuasa, Phys. Rev. A **87**, 052324 (2013).

[55] P. Facchi, G. Parisi, S. Pascazio, A. Scardicchio, and K. Yuasa, J. Phys. A: Math. Theor. **52**, 414002 (2019).

[56] P. J. Forrester, *Log-Gases and Random Matrices* (Princeton University, Princeton, NJ, 2010).

[57] C. Liu, X. Chen, and L. Balents, Phys. Rev. B **97**, 245126 (2018).

[58] P. Zhang, C. Liu, and X. Chen, SciPost Phys. **8**, 94 (2020).

[59] P. Lydzba, M. Rigol, and L. Vidmar, Phys. Rev. Lett. **125**, 180604 (2020).

[60] P. Lydzba, M. Rigol, and L. Vidmar, Phys. Rev. B **103**, 104206 (2021).

[61] E. Bianchi, L. Hackl, and M. Kieburg, Phys. Rev. B **103**, L241118 (2021).

[62] Indeed, if $\ell > M$, we can use that $V_{M,\ell}^\dagger V_{M,\ell}$ has $\ell - M$ zero eigenvalues and that its nonvanishing eigenvalues coincide with

those of $V_{M,\ell} V_{M,\ell}^\dagger$. In this way, we exchanged the role of $\ell$ and $M$: redefining formally $\ell' = M$, $M' = \ell$, and $V' = V^\dagger$, we reduced to the case $\ell' < M'$ treated in the main text.

[63] P. Vivo, S. N. Majumdar, and O. Bohigas, Phys. Rev. Lett. **101**, 216809 (2008).

[64] P. Vivo, S. N. Majumdar, and O. Bohigas, Phys. Rev. B **81**, 104202 (2010).

[65] K. Damle, S. N. Majumdar, V. Tripathi, and P. Vivo, Phys. Rev. Lett. **107**, 177206 (2011).

[66] In replacing the multiple integral with the functional integral, one needs to take into account the Jacobian $J[\rho]$ associated with the change of coordinates. It can be argued that $J[\rho] \sim e^{O(\ell)}$, and it may thus be neglected at the leading order in $\ell$ [67].

[67] D. S. Dean and S. N. Majumdar, Phys. Rev. E **77**, 041108 (2008).

[68] F. G. Tricomi, *Integral Equations* (Dover, New York, 1985).

[69] S. N. Majumdar and G. Schehr, J. Stat. Mech. (2014) P01012.

[70] A. C. Pipkin, *A Course on Integral Equations*, Texts in Applied Mathematics Vol. 9 (Springer-Verlag, New York, 1991).

[71] F. D. Gakhov, *Boundary Value Problems* (Elsevier, Amsterdam, 2014).

[72] Denis Bernard and Lorenzo Piroli, (in preparation).

[73] P. J. Forrester, J. Phys. A: Math. Theor. **45**, 145201 (2012).

[74] H. M. Ramli, E. Katzav, and I. P. Castillo, J. Phys. A: Math. Theor. **45**, 465005 (2012).

[75] In principle, numerical computation of $I_q(s)$ involves a double integral. However, one can get rid of the latter by using the saddle-point equation $\int d\lambda \rho_w^*(\lambda) \ln|\mu - \lambda| = (1/2)V(\mu) + u/2$, as done in Ref. [53]. See Ref. [72] for further details.

[76] W. Krauth, *Statistical Mechanics: Algorithms and Computations*, Oxford Master Series in Physics Vol. 13 (Oxford University, Oxford, 2006).