Entanglement entropy from nonequilibrium work

TL;DR

This paper introduces a novel fluctuation theorem-based method to estimate Rènyi entanglement entropy in quantum many-body systems, enabling efficient analysis of large-scale simulations and universal properties.

Contribution

It develops a new framework linking nonequilibrium work fluctuation theorems to entanglement entropy estimation, improving computational efficiency and scalability.

Findings

Efficient estimation of entanglement entropy using nonequilibrium work methods.

Application to large-scale quantum Monte Carlo simulations of SU(N) spin models.

Access to unprecedented system sizes up to 192 x 96 spins.

Abstract

The R\'enyi entanglement entropy in quantum many-body systems can be viewed as the difference in free energy between partition functions with different trace topologies. We introduce an external field $\lambda$ that controls the partition function topology, allowing us to define a notion of nonequilibrium work as $\lambda$ is varied smoothly. Nonequilibrium fluctuation theorems of the work provide us with statistically exact estimates of the R\'enyi entanglement entropy. This framework also naturally leads to the idea of using quench functions with spatially smooth profiles, providing us a way to average over lattice scale features of the entanglement entropy while preserving long distance universal information. We use these ideas to extract universal information from quantum Monte Carlo simulations of SU(N) spin models in one and two dimensions. The vast gain in efficiency of this method allows us to access unprecedented system sizes up to 192 x 96 spins for the square lattice Heisenberg antiferromagnet.

Equations

$$S_A^{(n)}=\frac{1}{1-n}\ln (\mathrm{Tr}{\rho }_A^n),$$

$$S_A^{(n)}=\frac{1}{1-n}\ln \left(\frac{Z_A^{(n)}}{Z_{\text{o}}^{(n)}}\right),$$

$${\mathcal{Z}}_A^{(n)}(\lambda )=\sum _{B\subseteq A}{\lambda }^{N_B}(1-\lambda {)}^{N_A-N_B}{Z}_B^{(n)},$$

$$S_A^{(n)}=\frac{1}{1-n}{\int }_0^{1}d\lambda \frac{\partial \ln {\mathcal{Z}}_A^{(n)}(\lambda )}{\partial \lambda }.$$

$$W_A^{(n)}=-\frac{1}{\beta }{\int }_{t_i}^{t_f}dt\frac{d\lambda }{dt}\frac{\partial \ln g_A(\lambda (t),N_B(t))}{\partial \lambda }$$

$$S^{(n)}_{A}=\frac{1}{1-n}\ln\left(\Big{\langle}e^{-\beta W^{(n)}_{A}}\Big{\rangle}\right),$$

$$\lambda (x,l(t),\delta )=\frac{1}{1+e^{(x-l(t)+\frac{1}{2})/\delta }}.$$

$$H=\frac{J}{N}\sum _i\sum _{\alpha ,\beta =1}^N|{\alpha }_i{\beta }_{i+1}⟩⟨{\beta }_i{\alpha }_{i+1}|.$$

$$S^{(n)}(l)=\frac{c}{6}\left(1+\frac{1}{n}\right)\log \left(\frac{L}{\pi }\sin \left(\frac{\pi l}{L}\right)\right)+b.$$

$$H=-\frac{J_1}{N}\sum _{⟨ij⟩}\sum _{\alpha ,\beta =1}^N|{\alpha }_i{\alpha }_j⟩⟨{\beta }_i{\beta }_j|-\frac{J_2}{N}\sum _{⟨⟨ij⟩⟩}\sum _{\alpha ,\beta =1}^N|{\alpha }_i{\beta }_j⟩⟨{\beta }_i{\alpha }_j|.$$

$$S^{(n)}(L/2)=aL+\frac{N_g}{2}\log (L)+b.$$

$$P_{\text{join}}=\min \left\{\frac{\lambda }{1-\lambda },1\right\}{P}_{\text{split}}=\min \left\{\frac{1-\lambda }{\lambda },1\right\}.$$

$$\begin{array}{rl}\Delta \ln g_A(\lambda ,N_B)=& (N_A-N_B(t_m))\ln \left(\frac{1-\lambda (t_{m+1})}{1-\lambda (t_m)}\right)\\ & +N_B(t_m)\ln \left(\frac{\lambda (t_{m+1})}{\lambda (t_m)}\right).\end{array}$$

$$\frac{Z_{A}}{Z_{\o}}=\bigg{\langle}\prod_{m}\frac{g_{A}(\lambda(t_{m+1}),N_{B}(t_{m}))}{g_{A}(\lambda(t_{m}),N_{B}(t_{m}))}\bigg{\rangle},$$

$${\mathcal{Z}}^{(n)}(\lambda )=\sum _{B\subseteq \mathcal{L}}\left(\prod _{x\in B}\lambda (x)\prod _{x\in \bar{B}}[1-\lambda (x)]\right)Z_B^{(n)}.$$

$$g(\lambda (x),B)=\prod _{x\in B}\lambda (x)\prod _{x\in \bar{B}}[1-\lambda (x)],$$

$$W^{(n)}=-\frac{1}{\beta }{\int }_{t_i}^{t_f}dt\frac{d\lambda }{dt}\frac{\partial \ln g(\lambda (x,t),B(t))}{\partial \lambda }.$$

$$\Delta \ln g=\ln \left(\frac{g(\lambda (x,t_{m+1}),B(t_m))}{g(\lambda (x,t_m),B(t_m))}\right),$$

$$\Delta \ln g=\ln \left(\prod _{x\in B(t_m)}\frac{\lambda (x,t_{m+1})}{\lambda (x,t_m)}\prod _{x\in \bar{B}(t_m)}\frac{1-\lambda (x,t_{m+1})}{1-\lambda (x,t_m)}\right).$$

$$|\psi ⟩=\underset{m\to \infty }{\lim }{H}^m|{\psi }_{\text{trial}}⟩.$$

$$S_A^{(2)}=-\ln \left(\frac{⟨\psi |{\text{Swap}}_A|\psi ⟩⟨\psi |{\text{Swap}}_A|\psi ⟩}{⟨\psi |\psi ⟩⟨\psi |\psi ⟩}\right),$$

Full text

Entanglement entropy is a quantity of basic importance in the characterization of quantum many-body wavefunctions. It signals the departure from a simple tensor product of subsystem wavefunctions for a spatial bipartition. Though in general the entanglement entropy depends only on the boundary size of the bipartition–the “area law”–and captures local correlations, violations to the area law can contain long distance universal information [Amico et al., 2008; Eisert et al., 2010; Laflorencie, 2016]. For instance, one dimensional critical systems show a logarithmic growth of the entanglement entropy with subsystem size where the prefactor depends on the central charge of the underlying conformal field theory [Holzhey et al., 1994; Vidal et al., 2003; Calabrese and Cardy, 2004]. Examples of universal features in higher dimensions include logarithmic contributions in symmetry broken phases with gapless Goldstone modes [Metlitski and Grover, 2011; Castro-Alvaredo and Doyon, 2012; Song et al., 2011; Ju et al., 2012; Luitz et al., 2015; Laflorencie et al., 2015; Rademaker, 2015; Kallin et al., 2011; Kulchytskyy et al., 2015; Luitz and Laflorencie, 2017] and a universal negative constant in topological phases [Kitaev and Preskill, 2006; Levin and Wen, 2006; Grover et al., 2013].

Low entanglement methods such as the density matrix renomalization group (DMRG) most naturally allow for the computation of entanglement entropy, which has been enormously successful in characterizing critical spin chains and identifying 2D gapped topological phases of matter [Laflorencie, 2016]. However, 2D gapless systems, long-range interactions and 3D systems still pose significant challenges for the DMRG method. This motivates the development of entanglement entropy methods in quantum Monte Carlo (QMC) that, barring a sign problem, do not face the same difficulties.

In fact, QMC simulations have made remarkable progress in this direction. The ability to interpret the Rényi entanglement entropies in terms of replica partition functions has enabled the introduction of various QMC estimators [Melko et al., 2010; Hastings et al., 2010; Humeniuk and Roscilde, 2012; Inglis and Melko, 2013]. While this has allowed for the investigation of many physical systems, the calculations remain costly due to the need to perform independent simulations either as a function of temperature or to incrementally compute the entanglement entropy for the region of interest. These inefficiencies are compounded in two or more dimensions where the area law dominates and universal features are subtle. Though improved estimators have been introduced [Luitz et al., 2014,Kulchytskyy et al., 2015], these represent either incremental advancements or are restricted to specific models. Simulations in this realm have yet to reach the large scale typically enjoyed by other more common QMC measurements.

At the same time, the ability to view the Rényi entanglement entropy as a difference of free energies begs the question of whether one can apply nonequilibrium work relations [Jarzynski, 1997, Crooks, 1999] that have been widely used in the molecular dynamics community [Ritort, 2006]. Indeed, this idea was recently explored in the context of classical path integral Monte Carlo [Alba, 2017], paving the way for the extension to QMC that we present here.

We will first present the basic idea for a fixed subsystem size using replica partition functions. We then show how this leads to the idea of smooth, space dependent quench functions that generalize the Rényi entanglement entropy while preserving its universal features and allows for computations as a function of subsystem size. These are used to extract universal features of SU($N$) magnets in one and two dimensions. Finally we extend the method to $T=0$ projector QMC simulations in the valence bond basis [Sandvik, 2005, Sandvik and Evertz, 2010], enabling ground state entanglement entropy measurements of the 2D Heisenberg model on truly large lattices.

General method: The Rényi entanglement entropy is given by

[TABLE]

where $n$ is the Rényi index ($n=2$ throughout this work) and $\rho_{A}=\mathrm{Tr}_{\bar{A}}\rho$ is the reduced density matrix of a subsystem $A$. We express the density matrix as $\rho=e^{-\beta H}/Z$ with the partition function $Z=\mathrm{Tr}e^{-\beta H}$. In what follows $\beta$ will be made sufficiently large so that only the ground state contributes.

$S^{(n)}_{A}$ can be conveniently re-expressed as [Calabrese and Cardy, 2004]

[TABLE]

where the replica partition functions are defined as $Z^{(n)}_{A}=\mathrm{Tr}\left(\left(\mathrm{Tr}_{\bar{A}}e^{-\beta H}\right)^{n}\right)$ and $Z^{(n)}_{\o}=\left(\mathrm{Tr}e^{-\beta H}\right)^{n}=Z^{n}$. Notice that $Z^{(n)}_{A}$ and $Z^{(n)}_{\o}$ are defined in exactly the same way, where $\o$ means that there are no sites in the $A$ subsystem. We now wish to define a function $\mathcal{Z}^{(n)}_{A}(\lambda)$ such that $\mathcal{Z}^{(n)}_{A}(0)=Z^{(n)}_{\o}$ and $\mathcal{Z}^{(n)}_{A}(1)=Z^{(n)}_{A}$. This can be written explicitly as

[TABLE]

where here $B$ denotes all possible subsets of $A$ including the empty set $\o$ and $A$ itself. $N_{A}$ and $N_{B}$ denote the number of sites contained in the sets $A$ and $B$, respectively. Here $\lambda$ controls the probability for an individual site in the $A$ subsystem to be traced only once. To ease the notation in what follows we define $g_{A}(\lambda,N_{B})=\lambda^{N_{B}}(1-\lambda)^{N_{A}-N_{B}}$.

The entanglement entropy can then be computed as:

[TABLE]

If one wanted to numerically compute this quantity using quantum Monte Carlo techniques, independent equilibrium simulations on a fine grid of points between $\lambda=0$ and $\lambda=1$ would need to be performed. Each simulation would measure the equilibrium average $\partial\ln\mathcal{Z}^{(n)}_{A}(\lambda)/\partial\lambda=\langle\frac{N_{B}}{\lambda}-\frac{N_{A}-N_{B}}{1-\lambda}\rangle_{\lambda}$ and the resulting curve from all the simulations would need to be numerically integrated.

This approach has obvious practical limitations. It is much more desirable to estimate the entanglement entropy directly from a single simulation. Fortunately, this can be done by defining a nonequilibrium process in which $\lambda$ is varied smoothly from 0 to 1:

[TABLE]

Here $\lambda(t_{i})=0$, $\lambda(t_{f})=1$ and $W^{(n)}_{A}$ is the total work done throughout the process. We emphasize here that the work is a random variable that will follow a distribution. One instance of the work is given by the sum of the increments $\partial\ln g_{A}(\lambda(t),N_{B}(t))$ along a nonequilibrium path in the configuration space of $\mathcal{Z}^{(n)}_{A}(\lambda)$ as $\lambda$ is varied from 0 to 1. The average $\beta\langle W^{(n)}_{A}\rangle/(n-1)$ then approaches $S^{(n)}_{A}$ as $t_{f}-t_{i}\to\infty$. If the quench time is finite then the average will overestimate $S^{(n)}_{A}$ due to the nonequilibrium entropy production associated with irreversibility.

Again the situation is undesirable, since we would like to accurately estimate the entanglement entropy even when the quench time is finite. Here we can greatly prosper from a well known nonequilibrium fluctuation theorem of work, Jarzynski’s equality [Jarzynski, 1997], which in this context reads:

[TABLE]

Remarkably this equality holds true regardless of the quench time.

Quench protocols: In the previous section we took the external field $\lambda$ to be constant in space and varying in time between 0 and 1. This is depicted in panel $(a)$ of Fig. 1. This provides us with a very efficient way of computing the entanglement entropy of large regions, which we will demonstrate in the final section using $T=0$ projector QMC simulations. At this point, however, we would like to introduce another quench protocol that can compute the entanglement entropy as a function of subsystem size in a single nonequilibrium simulation.

This can be achieved by making $\lambda$ a function of both time and space in such a way that the subsystem partition “slides” across the lattice, in which case Eq. (3) can be easily adapted (see supplemental material). This naturally leads us to the idea that $\lambda$ should be a smooth function of space so that individual sites are quenched gradually while the partition moves. A convenient choice is given by

[TABLE]

Here $x$ is the spatial coordinate (from 0 to $L-1$) along the chain and $\lambda$ depends on time through the parameter $l(t)$ which represents the center of one boundary between the $A$ subsystem and the rest of the chain (see panel $(b)$ and $(c)$ of Fig. 1). A shift of $\tfrac{1}{2}$ is introduced so that $l(t)$ is centered between sites of the chain. We have also introduced a parameter $\delta$ that controls the width of the boundary. As we will discuss momentarily, $\delta$ allows us to suppress lattice scale features of the entanglement entropy while preserving the universal features.

The quench protocol as defined by Eq (7) involves varying $l(t)$ linearly in time from $l(t_{i})=-p\delta$ to $l(t_{f})=L-1+p\delta$, where the constant $p$ is used to ensure that $\lambda(x,l(t_{i}),\delta)$ and $\lambda(x,l(t_{f}),\delta)$ are approximately 0 and 1 everywhere, respectively. As $l(t)$ is varied along the length of the chain, the total work is computed as a function of time. Since $l(t)$ represents the current location of the subsystem boundary, this allows us to compute the entanglement entropy as a function of subsystem size in a single nonequilibrium simulation. Furthermore, since $l(t)$ is real valued, the entanglement entropy is computed as a continuous curve that can be sampled up to the number of time steps used in the nonequilibrium process.

We now wish to discuss the role of the parameter $\delta$, which is a central idea that emerges from our new method. When $\delta$ is small enough so that the width of the subsystem boundary is much less than the lattice spacing, our protocol allows for the exact computation of $S^{(2)}_{A}$ as a function of subsystem size. The entanglement entropy thus defined contains strong lattice-scale features that can often obscure universal properties that one wishes to extract. By increasing $\delta$ so that the subsystem partition is spread over a few lattice sites we effectively average over features below a certain scale. As we intuitively expect, the presence of $\delta$ will only affect non-universal features of the entanglement entropy, namely the area law contribution.

SU(N) chain: As a first test case for our new method, we start with the SU($N$) Heisenberg antiferromagnet on a periodic one dimensional chain [Sutherland, 1975, Affleck, 1988]. The Hamiltonian can be written simply as a nearest neighbor permutation:

[TABLE]

The ground state of this model is critical, and belongs to the family of conformally invariant Wess-Zumino-Witten nonlinear sigma models with central charge $c=N-1$ [Affleck, 1988]. This system gives rise to the celebrated log violation of the area law, which is given by the following form:

[TABLE]

The Rényi entanglement entropies of this system also show oscillations as a function of subsystem size [D’Emidio et al., 2015]. In Fig. 2 we use our newly developed nonequilibrium method in the framework of the stochastic series expansion QMC algorithm [Sandvik, 2010], combined the with quench protocol in Eqn. (7) to show that lattice scale oscillations of the entanglement entropy can be suppressed by smoothing out the subsystem boundary using the parameter $\delta$. Crucially the central charge is insensitive to the presence of $\delta$, which only affects the nonuniversal constant (area law) contribution (see supplemental materials). As expected, we find perfect agreement with the central charge $c=N-1$.

2D staggered SU(N) Model: Recently, the entanglement entropy has been used to detect Goldstone modes in magnetically ordered systems [Metlitski and Grover, 2011, Kulchytskyy et al., 2015]. To this end, we select a family of two dimensional SU($N$) symmetric models that support gapless spin wave excitations. The Hamiltonian is given by [Kaul et al., 2013]

[TABLE]

Here we consider periodic square lattices, using the fundamental representation of SU($N$) on one sublattice and the conjugate to the fundamental on the other sublattice. $J_{1}>0$ and $J_{2}>0$ are the nearest and next-nearest neighbor couplings, respectively.

This model is an SU($N$) generalization of the spin-$\tfrac{1}{2}$ Heisenberg antiferromagnet with a ferromagnetic next-nearest neighbor interaction. When $J_{2}=0$ the ground state is magnetically ordered with gapless Goldstone modes for $2\leq N\leq 4$ and forms an valence bond solid for $N\geq 5$ [Kaul and Sandvik, 2012]. When $J_{2}>0$ magnetic order is enhanced which greatly facilitates the ability to accurately extract the number of Goldstone modes by fitting to our entanglement entropy data.

The presence of Goldstone modes produces a log contribution to the entanglement entropy [Metlitski and Grover, 2011]. In the simplest case, we can take an $L$ x $L$ cylinder and cut it in half. The resulting entanglement entropy as a function of the linear system size $L$ is given by

[TABLE]

Here $N_{g}$ is the number of Goldstone modes and b contains a geometrical constant.

We again apply our nonequilibrium method in combination with the quench protocol in Eq. (7), where now the quench function is taken to be spatially constant in the $y$-direction and the partition is swept along the $x$-direction. Fig. (3) shows the resulting curves for system sizes $L=6,8,..,20$ for SU(2). We take the midpoint of our data and perform a fit to Eq. (11). This procedure is performed for SU(2), SU(3) and SU(4) with $J_{2}/J_{1}=1,2,3.5$, respectively. The inset of Fig. (3) shows the entanglement entropy at the center cut with the area law piece subtracted. On semi-log axes, the log contribution manifests itself as a straight line. Within one percent accuracy we find the number of Goldstone modes to be $N_{g}=2,4,6$ for $N=2,3,4$ respectively.

Projector QMC simulations: We finally wish to extend our method to $T=0$ projector QMC simulations in the valence bond basis [Sandvik, 2005, Sandvik and Evertz, 2010]. The motivation here is that one can much more easily reach the groundstate as compared with the finite temperature stochastic series expansion algorithm. Furthermore, we directly compute the half-system second Rényi entanglement entropy by performing the quench analogous to panel (a) of Fig. 1. These efficiency gains allow us to access unprecedented system sizes for the square lattice spin-$\tfrac{1}{2}$ Heisenberg model, shown in Fig. 4. Finally, we note that the convergence of $N_{g}$ by fitting to the form in Eq. (11) is very slow when a next-nearest-neighbor coupling is absent.

Conclusion: We have introduced a highly efficient nonequilibrium method for computing the Rényi entanglement entropy in the context of statistically exact QMC simulations. Remarkably, in all of our simulations we observe that the QMC statistical error remains constant when the number of measurement sweeps is proportional to the number of sites in the entangling region (see supplemental material). Room for improvement is also still possible since a linear quench protocol for a fixed region is most certainly not optimal.

Acknowledgements: We thank Frédéric Mila for his encouragement and support for this project. We also thank Ribhu K. Kaul for suggesting the 2D staggered SU(N) model as a test case for this method. SSE simulations were carried out on the Fidis cluster at EPFL and projector simulations on Comet at SDSC with the Extreme Science and Engineering Discovery Environment (XSEDE) Allocation No. TG- DMR130040.

Note added: Recently, another important work on com- puting entanglement entropy using thermodynamic concepts via a field theoretic mapping has also been introduced [Mendes-Santos et al., 2020].

I Supplemental material

I.1 QMC sampling and measurement

Our task is to perform nonequilibrium simulations in the space of configurations contained in $\mathcal{Z}^{(n)}_{A}(\lambda)$ and add up all of the increments $\partial\ln g_{A}(\lambda(t),N_{B}(t))$ along a path between $\lambda=0$ and $\lambda=1$. We must first describe how to stochastically sample the $\mathcal{Z}^{(n)}_{A}(\lambda)$ configurations in equilibrium. We employ the stochastic series expansion QMC algorithm [Sandvik, 2010] using the replica trick [Hastings et al., 2010], where one samples $n$ independent copies of the partition function in a single simulation. These are the configurations of $Z^{(n)}_{\o}$, which are one subset of $\mathcal{Z}^{(n)}_{A}(\lambda)$. In order to transition to other subsets, spins in different (independent) traces need to be joined into a single trace. This can be done if the spin states in the different traces match each other. The weight of the configuration will also change upon joining spins, since this changes the prefactor $g_{A}(\lambda,N_{B})=\lambda^{N_{B}}(1-\lambda)^{N_{A}-N_{B}}$. We take the acceptance probability of joining or splitting spins (with the requirement that the spin states must match) to be the ratio of the weights of the new and old configurations. This gives the the acceptance probabilities

[TABLE]

Each site in the $A$ subsystem is given the opportunity to join or split (depending on its current trace topology). Once this rewiring is done, the connectivity of the configuration is fixed and one may perform a standard QMC update.

We now need to know how to measure the work increments as we move along a nonequilibrium trajectory in the space of configurations contained in $\mathcal{Z}^{(n)}_{A}(\lambda)$. The increment between time $t_{m}$ and $t_{m+1}$ is given by

[TABLE]

In other words, the increment is calculated by fixing the configuration and computing the change in $\ln g_{A}$ as $\lambda$ is incremented. This is in direct analogy to the way in which nonequilibrium work increments are computed in classical systems, except there one computes the change in energy of a configuration.

At this point two comments are in order. Firstly, one can avoid numerical rounding errors incurred by summing the $\log$ increments by instead taking a product of the arguments. Secondly, the factor of $-1/\beta$ in the definition of the work (Eq. (5)) completely drops out of the formulation when one uses Jarzynski’s equality (Eq. (6)). It is only included to give $W^{(n)}_{A}$ the proper units of work. With these ideas in mind, we can succinctly write the Jarzynski estimator for the partition function ratio as follows:

[TABLE]

where here we have suppressed the Rényi index.

Now that we have the ingredients necessary to update and perform work measurements on our configurations, we will briefly outline the the main steps of the algorithm. First we equilibrate the system in the $Z^{(n)}_{\o}$ ensemble of configurations (i.e. $\lambda=0$). The equilibrated configuration is then saved to a file, and the nonequilibrium process begins. First we measure the work increment by Eq. (13) (or better to multiply the $\lambda$ factors from the first time step in Eq. (14)), then the current value of $\lambda$ is incremented and each spin in the $A$ subsystem is given the opportunity to change its trace topology according the the probabilities in Eq. (12). The trace topology is then held fixed and a regular QMC update is performed. Then the cycle begins again by computing the next increment and so on. Throughout the entire nonequilibrium process all of the work increments are accumulated and saved, in practice sampling them at regular intervals.

In order to initiate a new nonequilibrium process, we first read in the saved configuration that was equilibrated at the beginning, then we re-equilibrate for some smaller number of QMC steps (always staying in the $Z^{(n)}_{\o}$ ensemble) and again save the configuration to a file. This new starting configuration is used for the next nonequilibrium process.

We note that if the quenches are fast and the entangling region is large, it can happen that at the end of the quench when $\lambda=1$ spins have not entirely joined and the ratio formula in Eq. (14) gives zero. This should be avoided by increasing the quench time, but (so as not to lose the work measurement) one can multiply by one (instead of zero) for these spins. Although we have tried to avoided it, this does not seem to have any detectable effect on the accuracy of the entanglement entropy measurements, even for reasonably short quench times.

I.2 Space dependent quench functions

In order to clearly see how to treat space dependent quench functions, we write the formula for $\mathcal{Z}^{(n)}(\lambda)$ as:

[TABLE]

Here $B$ is now summed over all proper subsets of the entire lattice. One can now interpolate between $Z^{(n)}_{\o}$ and $Z^{(n)}_{A}$ by using the space dependent quench function $\lambda(x,t)=\lambda(t)\chi_{A}(x)$, where the indicator function $\chi_{A}(x)=1$ if $x\in A$ (zero otherwise), and $\lambda(t)=(t-t_{i})/(t_{f}-t_{i})$. This is the same function appearing in panel (a) of Fig. 1, but formulating things in this way allows us to treat arbitrary quench functions.

Computing the work increments in this case is no more complicated. First we define the $g$ function as

[TABLE]

and the dynamical work is

[TABLE]

The work increment accumulated between two adjacent time steps $t_{m}$ and $t_{m+1}$ is given by

[TABLE]

or

[TABLE]

The total work and Jarzynski estimator for the Rényi entanglement entropy follow naturally, and the space dependent $\lambda$ is used in the joining and splitting probabilities.

I.3 QMC versus ED

When using any numerical method, it is always important to compare with exact results. We have gone to great lengths to check the accuracy of our method for small systems that can be diagonalized exactly. Fig. 5 shows a comparison of our QMC method with exact diagonalization of an $\mathrm{SU}(2)$ $L=16$ chain. We have used the quench function in Eq. (7) for different values of $\delta$. The QMC data is the colored curves with shading as the error bar (this can be seen in the zoomed inset), and the black curves are exact results obtained by diagonalizing the reduced density matrices for all possible bipartitions and weighting them with the appropriate factors of $\lambda(x,l,\delta)$. We find perfect agreement between the QMC and ED. We also see that when $\delta$ is made small enough the exact Rényi entanglement entropy of a block subsystem (the black dots) is produced.

We also provide the same type of comparison in Fig. 6 for our two dimensional model on an SU(2) $L=4$ square lattice with $J_{2}/J_{1}=2$ using the same quench function as in the the main text. Again we find perfect agreement within the error bars.

I.4 Smooth quench fuctions and universal information

We now wish to show that our smooth quench function preserves the universal features of the entanglement entropy. In Fig. 7 we use our newly developed nonequilibrium method combined with the quench protocol in Eq. (7) to show that lattice scale oscillations of the entanglement entropy can be suppressed by smoothing out the subsystem boundary using the parameter $\delta$. Crucially the central charge is insensitive to the presence of $\delta$, which only affects the nonuniversal constant (area law) contribution.

I.5 Efficiency of Method

It is clear from the results presented in the main text that this method is very efficient. But here we would like to give a direct comparison of our nonequilibrium quenches with the extended ensemble method [Humeniuk and Roscilde, 2012] combined with the increment trick [Hastings et al., 2010]. In Fig. 8 we show the statistical error (QMC error bar) as a function of subsystem size $L_{A}=L/2$ for $\mathrm{SU}(N)$ chains. Here we compute the second Rényi entanglement entropy of half of the chain in two ways: (1) by directly quenching the entire half chain using our nonequilibrium method and (2) by using the extended ensemble method with increments of one site at a time. For the extended ensemble method, we compute each increment in a separate (equilibrium) simulation with 28 x 10 binned measurements each consisting of 10,000 measurement sweeps. For our nonequilibrium method we quench the entire half chain at once with 28 x 10 work realizations each consisting of $L_{A}$ $\mathrm{x}$ 10,000 nonequilibrium time steps. As such, the total number of measurement sweeps used for each method is identical. Very remarkably, in the nonequilibrium case the error bar stays flat as long as the number of nonequilibrium time steps is increased in proportion to the number of sites quenched. Significant computational resources are also saved on equilibration, since the independent equilibrium increments need to be separately equilibrated.

I.6 $T=0$ converged data for 2D model

Here we would like to provide a more detailed view of the data presented in the main paper for the two dimensional SU($N$) model. In order to obtain $T=0$ converged data, we have had to set $\beta=L^{2}$ with $J_{1}=1$. These extremely low temperatures prohibit us from simulating systems much larger than $L=20$. In Fig. 9 we again show the raw data for $\mathrm{SU}(2)$ appearing in the main text, this time with subfigures that show (from left to right) the $\log$ contribution on linear axes, the QMC data with the fit subtracted, and the fitted value of $N_{g}$ as a function of the smallest system size used. We show the same type of plot for $\mathrm{SU}(3)$ and $\mathrm{SU}(4)$ in Fig. 10 and Fig. 11, respectively.

I.7 Finite-$T$ data for 2D model

We have also collected data at finite temperatures, this time only scaling $\beta$ proportional to $L$ and not $L^{2}$. Remarkably, the finite temperature effects do not influence the extraction of $N_{g}$ as we show in Fig. 12 for $\mathrm{SU}(3)$ and Fig. 13 for $\mathrm{SU}(4)$. These results at finite temperature arguably produce even more high quality fits than the data converged at $T=0$.

Finally we include data for SU(4) at a lower value of $J_{2}/J_{1}$ to illustrate the finite size effects in this case. Fig. 14 shows finite temperature data collected with $\beta=3L$ and $J_{2}/J_{1}=2$. We can see a systematic dependence of $N_{g}$ on $L_{\mathrm{min}}$ in our numerical fits, with a drift toward the true value $N_{g}=6$.

I.8 Projector QMC simulations

We now briefly outline how the nonequilibrium Rényi entanglement entropy measurements can be extended to projector QMC simulations in the valence bond basis [Sandvik, 2005, Sandvik and Evertz, 2010]. Here instead of using Monte Carlo techniques to sample the quantum partition function at finite temperature, one instead chooses a trial wave function that is projected into the ground state by acting with powers of the Hamiltonian. The main advantage (aside from some additional technical simplifications) is that very good trial wavefunctions can be used, so that the number of operators needed to reach the ground state is far fewer compared to the stochastic series expansion at low temperatures.

In the projector method we use the fact that the unnormalized ground state $|\psi\rangle$ can be written as

[TABLE]

One can then write the density matrix (dividing by the norm) and use it to construct the Rényi entanglement entropy for a region $A$. Just as with the replica partitions functions, the Rényi entanglement entropy will take the form of a log of a ratio. In this case, however, it will be a ratio of wave function overlaps. For the second Rényi entropy of a region $A$, we can write it as follows:

[TABLE]

where $\text{Swap}_{A}$ means that the spin degrees of freedom have been exchanged in the region $A$ between the two overlaps in the numerator.

We now find ourselves in exactly the same situation as with the finite temperature simulations, we just need to vary an external field from 0 to 1 in such a way that it takes us from configurations contained in the denominator of Eq. (21) to configurations contained in the numerator. To achieve this we can couple the external field $\lambda$ to the single site swap operator on each site of the $A$ subsystem. We use the same prefactor $g_{A}$ as before, and the measurement process is identical. To be clear, during the simulation single site swap operators are inserted or removed throughout the $A$ subsystem with the probabilities in Eq. (12) whenever the spin degrees of freedom match in the center of both wave function overlaps.

For the 2D Heisenberg antiferromagnet we choose to work in the combined basis [Sandvik and Evertz, 2010] of $S^{z}$ spin values and valence bond (singlet) coverings of the square lattice. Working with valence bonds allows for the construction of good trial states for the Heisenberg ground state and using the $S^{z}$ spin values enables efficient loop updates of the configurations. We take our trial states to be valence bond coverings with a $1/r^{3}$ potential [Sandvik and Evertz, 2010].

Firstly in order to verify our $T=0$ projector code, we make a comparison with the second Rényi entanglement entropy obtained by exact diagonalization on a 4 x 4 (periodic) lattice, which is shown in Fig. 15.

In order to ensure that our measurements are consistent (especially on very large system sizes), we have performed separate nonequilibrium work calculations in both the forward and reverse directions. We find perfect agreement in all cases. The forward and reverse measurements can be combined so that the final result is obtained at each point along the trajectory ($\lambda$ between 0 and 1). The average of all of these points can be used as the final average and a conservative estimate of the error bar can be taken as the average error bar of all the points plus the standard deviation of the points from the mean. This is illustrated in Fig. 16 for a 24 x 24 system that has been cut in half, which we compare with the same measurement reported in [Kallin et al., 2011].

Finally, in Fig. 17 we show the convergence of our large scale entanglement data as a function of the projection power $m$. Here $m$ needs to be scaled with the number of sites, and the proportionality constant increases with the system size. We find that the entanglement monotonically approaches the value in the ground state.