Out-of-equilibrium protocol for Rényi entropies via the Jarzynski equality
In recent years entanglement measures, such as the von Neumann and the Rényi entropies, provided a unique opportunity to access elusive feature of quantum many-body systems. However, extracting entanglement properties analytically, experimentally, or in numerical simulations can be a formidable task. Here, by combining the replica trick and the Jarzynski equality we devise a new effective out-of-equilibrium protocol for measuring the equilibrium Rényi entropies. The key idea is to perform a quench in the geometry of the replicas. The Rényi entropies are obtained as the exponential average of the work performed during the quench. We illustrate an application of the method in classical Monte Carlo simulations, although it could be useful in different contexts, such as in quantum Monte Carlo, or experimentally in cold-atom systems. The method is most effective in the quasi-static regime, i.e., for a slow quench. As a benchmark, we compute the Rényi entropies in the Ising universality class in dimensions. We find perfect agreement with the well-known Conformal Field Theory (CFT) predictions.
In recent years entanglement measures have arisen as new diagnostic tools to unveil universal behaviors in quantum many-body systems. Arguably, the most popular and useful ones are the Rényi entropies and the von Neumann entropy amico-2008 (); calabrese-2009 (); eisert-2010 (); laflorencie-2016 () (entanglement entropies). Given a system in a pure state and a bipartition into an interval and its complement (see Figure 1), the Rényi entropies for part are defined as
with the reduced density matrix of , and its -th moment. The limit defines the von Nemann entropy . Due to being non-local, extracting analytically, experimentally, or even in numerical simulations can be a challenging task, execpt for free-fermion and free-boson models, for which the entropies can be obtained exactly in arbitrary dimensions eisler-2009 ().
A large class of effective measurement protocols for the Rényi entropies are based on the replica trick. The key observation is that for a generic model can be obtained as calabrese-2004 ()
where , and is the partition function on the so-called -sheets Riemann surface (see Figure 2 (b)), which is defined by “gluing” together independent replicas of the model through part . Importantly, the replica trick lies at the heart of all the known methods for measuring entanglement in cold-atom experiments cardy-2011 (); abanin-2012 (); daley-2012 (); islam-2015 (); kaufman-2016 (); pichler-2016 (). For instance, very recently (see Ref. islam-2015, ) has been succesfully measured in ultra-cold bosonic systems from the interference between two replicas of a many-body state. Moreover, the ratio of partition functions in (2) can be sampled using classical buividovich-2008 (); caraglio-2008 (); alba-2010 (); gliozzi-2010 (); alba-2011 (); alba-2013 () and quantum Monte Carlo techniques hastings-2010 (); melko-2010 (); singh-2011 (); isakov-2011 (); humeniuk-2012 (); kaul-2013 (); inglis-2013 (); iaconis-2013 (), providing an efficient method to calculate Rényi entropies. Extensions of these techniques for systems in the continuum herdman-2014 (), and interacting fermions zhang-2011 (); tubman-2012 (); mcminis-2013 (); grover-2013 (); assaad-2014 (); broecker-2014 (); wang-2014 (); shao-2015 (); drut-2015 (); drut-2016 (); porter-2016 () are also available. Monte Carlo methods work effectively in any dimension and at any temperature, for sign-problem-free models. Oppositely, the Density Matrix Renormalization Group uli-2005 (); uli-2011 () (DMRG) provides the most effective way to access the full spectrum of for one-dimensional systems, whereas is less effective in higher dimensions. A severe issue of all the replica-trick-based protocols is that as the size of increases, is dominated by rare configurations. The commonly used strategy to mitigate this issue is based on the increment trick hastings-2010 (); alba-2011 (); humeniuk-2012 (). This consists in splitting the ratio in (2) as a product of intermediate terms, which have to be measured separately. Their number typically grows as , which is the main drawback of the protocol.
Here we propose a new out-of-equilibrium framework for measuring the Rényi entropies. Similar protocols to access entanglement in cold-atom experiments have been explored in Ref. cardy-2011, and Ref. abanin-2012, . Our approach combines the replica trick (2) and the Jarzynski equality jarzynski-1997 (). Crucially, the latter allows to relate the ratio of partition functions corresponding to two equilibrium thermodynamic states to the exponential average of the work performed during an arbitrary far-from-equilibrium process connecting them. The idea of the method is to modify (i.e., quenching) the geometry of the replicas, gradually driving the system from the geometry with independent replicas to that of the -sheets Riemann surface. Using the Jarzynski equality, the Rényi entropies are then extracted from the statistics of the work performed during the quench, similar to Ref. cardy-2011, . Some applications of the Jarzyinski equality to detect entanglement have been also presented in Ref. hide-2010, . The efficiency of the method depends dramatically on the rate at which the geometry is modified. Precisely, the number of independent protocol realizations needed increases upon increasing . However, in the quasi-static regime, i.e., large , can be extracted from a single realization. In this regime the Rényi entropies depend only on the average work and the standard deviation of the work fluctuations, reflecting that the work distribution function becomes gaussian.
Here we illustrate the approach in the framework of classical Monte Carlo simulations. Specifically, we focus on the Ising universality class in dimensions, although the method works in any dimension, and it could be extended to quantum Monte Carlo. Remarkably, even for moderately large system sizes the Rényi entropies can be extracted in a single simulation, in contrast with methods using the increment trick, which require typically independent simulations.
Ii Method: Quenching the replica geometry
The replica geometries used in the method are illustrated in Figure 2. Panel 2 (a) shows a single replica (sheet), where and are interpreted as the spatial and imaginary time directions, respectively. The single-sheet partition function is obtained as the path integral , where is a field and is the euclidean action of the model. We impose the periodic boundary conditions and , with and the number of sites along the and direction, respectively. For simplicity, we assume interaction only between nearest-neighbor sites , i.e., , with the interaction strength. In the replica trick (2) one has to consider replicas of the model ( sheets). The partition function on independent sheets is , where now denote fields living on the replica . We now consider the situation with coupled sheets. First, on each sheet a branch cut lying along the spatial direction is introduced (as in Figure 2 (a)). The replicas are coupled through (for see Figure 2 (b)). The partition function on the coupled sheets is , where
Here denotes links crossing the cut, and is the coupling between fields next to cut and living on replicas and . The replicas are coupled in a cyclic fashion, meaning that the replica indices are defined . The partition function on the -sheets Riemann surface (cf. (2)) corresponds to .
Any ratio can be calculated using the Jarzynski equality jarzynski-1997 (). Specifically, let us consider a system at equilibrium at an initial time . Let be its partition function. Now let us imagine that the system is driven to a new equilibrium state at , which is described by , using an arbitrary out-of-equilibrium protocol. For an Hamiltonian system this could be a quench, in which some parameters of the Hamiltonian are varied with time. The Jarzynski equality states that
Here is the infinitesimal work performed between time and , is the inverse temperature, and denotes the average over different realizations of the quench protocol. The Jarzynski equality has been verified in several systems, and it is routinely used to extract free energy differences in out-of-equilibrium experiments hummer-2001 (); liphardt-2002 (); douarche-2005 (); collin-2005 (); bustamante-2005 (); blickle-2006 (); harris-2007 (); saira-2012 (); jarzynski-2011 (); seifert-2012 (); shuoming-2015 (), and in Monte Carlo simulations (see for instance Ref. caselle-2016 ()).
The idea of our method is to use (4) performing a quench of the couplings (cf. (3)) from the starting configuration with ( independent replicas) and final one with (coupled replicas). Any quench protocol is expected to give the same result. We choose the ramp protocol
Here is the Monte Carlo time, and is the quenching rate. Importantly, has to be large enough to ensure initial thermal equilibrium.
where is the integrated work performed during the Monte Carlo history, i.e., , and is the change in energy between two consecutive Monte Carlo steps and , calculated using the field configurations at , i.e.,
where is the work variance. Eq. (8) is derived assuming that the work distribution function is gaussian, and it is also known as Callen-Welton fluctuation dissipation relation callen-1951 (). The second term in (8) corresponds to the work dissipated during the quench. Interestingly, Eq. (6) implies that the number of independent simulations needed to obtain a reasonable estimate of increases exponentially with (see Ref. jarzynski-2006, for a rigorous result). On the other hand from (8), one has that in the quasi-static limit , implying that can be extracted from a single protocol realization.
Iii Numerical checks in the Ising universality class
We now provide numerical evidence supporting the correctness and efficiency of our Monte Carlo method. We consider the two-dimensional classical critical Ising model (cf. (10)). The model has a critical point at . Here we consider very elongated lattices with (anisotropic limit). In this limit universal properties are the same as in the one-dimensional quantum critical Ising chain at zero temperature, which is defined by the Hamiltonian , with the Pauli matrices acting on site of the chain. The critical behavior of both models is described by a Conformal Field Theory yellow-book () (CFT) with central charge .
In our Monte Carlo simulations we employed the Swendsen-Wang method swendsen-1987 (), although any other type of update can be used. In the Swendsen-Wang update, at any Monte Carlo time one assigns to each link of the lattice (here the connected replicas, see Figure 2 (b)) an auxiliary activation variable. Links connecting pairs of aligned spins and not crossing the branch cut are then activated with probability , while links connecting aligned spins around the branch cut are activated with probability , with given in (5). In our Monte Carlo simulations we fixed and . Then, all the different clusters of spins are identified using the rule that pairs of spins connected by an activated link are in the same cluster. Finally, all the spins belonging to the same cluster are flipped with probability .
We should mention that for models that can be mapped to the random cluster model, à la Fortuin Kasteleyn, it is usually convenient to write Monte Carlo observables in terms of cluster-related variables. Typically, at least for equilibrium simulations, this leads to a considerable speedup. This can be done for the Rényi entropy estimator (6), as we discuss in A. Surprisingly, our results suggest that the cluster estimator performs worse than the standard local estimator (6). This is in contrast with cluster estimators for the Rényi entropies in equilibrium Monte Carlo simulations, which outperform local estimators caraglio-2008 ().
The results for are illustrated in Figure 3. Panel (a) shows as obtained using (6) and plotted versus . The symbols are Monte Carlo data for system sizes with , averaged over independent realizations of the driving protocol. In the scaling limit , the behavior of is obtained from CFT as calabrese-2004 ()
where is the central charge and a non universal constant. The dashed-dotted lines in the Figure are obtained from (9) by fitting , after fixing . The good agreement with the data confirms the validity of the method. This is also clear from the perfect linear behavior in panel (b), where we plot versus . The inset in Figure 3 (a) plots as obtained from a single realization of the driving protocol, i.e., a single Monte Carlo simulation. The dashed-dotted line is the same as in the main Figure. The good agreement with the data suggests that the protocol (5) with is already close to the quasi-static regime, at least for .
It is also interesting to investigate the behavior of as a function of the value of during the Monte Carlo dynamics. This is discussed in Figure 4. The data are for a single realization of the driving protocol. At , one has , reflecting that at the early stage of the Monte Carlo the two replicas are disconnected. Interestingly, exhibits a maximum around . The expected symmetry , which reflects that the zero-temperature one-dimensional system is in a pure state, is observed only at (see arrows in the Figure).
We now discuss the fluctuations of between different realizations of the driving protocol. Figure 5 plots as a function of for independent realizations of the driving protocol. The inset (a) shows a zoom around the region . Fluctuations between different realizations of the dynamics are of the order of a few percents. The behavior as a function of of the statistical error associated with the fluctuations between different realizations is illustrated in inset (b). Here is calculated as the standard deviation of the fluctuations, considering a sample of independent simulations. Clearly, increases mildly as a function of . The dashed-dotted line in the inset is a fit to . Finally, it is interesting to investigate the dependence of on . Figure 6 shows as a function of for several values of . Most notably, for any one has that the symmetry is not present. We should mention that this resembles the behavior of when the total system is not pure, for instance at finite temperature. Moreover, the data for suggest a volume-law behavior (see inset in Figure 6), again as in finite temperature systems.
We presented a novel out-of-equilibrium framework for measuring the Rényi entropies by combining the Jarzynski equality and the replica-trick. As an application, we presented a new classical Monte Carlo method to measure the Rényi entropies.
This work opens numerous research directions. First, it would be useful to implement the approach in the framework of quantum Monte Carlo. Furthermore, it would interesting to understand analytically the behavior of as a function of (see Figure 6). An intriguing possibility is that could be interpreted as a finite-temperature for the one-dimensional system. This could provide an alternative way to obtain the finite-temperature Rényi entropies. Interestingly, our approach could be used for the moments of the partially transposed reduced density matrix alba-2013 (); chung-2014 (). These are the main ingredients to construct the logarithmic negativity vidal-2002 (); plenio-2005 (); calabrese-2012 (), which is a good entanglement measure for mixed states. Another interesting direction is to consider a protocol in which the length of is also quenched. This would allow to obtain the entropies for different subsystem sizes in a single simulation. Finally, it is important to understand whether our framework could provide a viable alternative to measure Rényi entropies in cold-atom experiments or in NMR quantum simulators fan-2016 (); li-2016 ().
I am greatly indebted with Claudio Bonati for drawing to my attention Ref. caselle-2016, and for useful discussions. I thank Pasquale Calabrese and Paola Ruggiero for reading the manuscript and for useful suggestions. I acknowledge support from the ERC under the Starting Grant 279391 EDEQS. This project has received funding fromt the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowoska-Curie grant agreement No 702612.
- (1) L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Rev. Mod. Phys. 80, 517 (2008).
- (2) P. Calabrese, J. Cardy, and B. Doyon Eds, J. Phys. A 42, 500301 (2009).
- (3) J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. 82, 277 (2010).
- (4) N. Laflorencie, arXiv:1512.03388.
- (5) V. Eisler and I. Peschel, J. Phys. A: Math. Theor. 42 504003 (2009).
- (6) P. Calabrese and J. Cardy, J. Stat. Mech.: Theor. Exp. P06002 (2004).
- (7) J. Cardy, Phys. Rev. Lett. 106, 150404 (2011).
- (8) D. A. Abanin and E. Demler, Phys. Rev. Lett. 109, 020504 (2012).
- (9) A. J. Daley, H. Pichler, J. Schachenmayer, and P. Zoller, Phys. Rev. Lett. 109, 020505 (2012).
- (10) R. Islam, R. Ma, P. M. Preiss, M. E. Tai, A. Lukin, M. Rispoli, and M. Greiner, Nature 528, 77 (2015).
- (11) A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, arXiv:1603.04409
- (12) H. Pichler, G. Zhu, A. Seif, P. Zoller, and M. Hafezi, arXiv:1605.08624.
- (13) P. V. Buividovich and M. I. Polikarpov, Nucl. Phys. B 802, 458 (2008).
- (14) M. Caraglio and F. Gliozzi, JHEP 11, 076 (2008).
- (15) V. Alba, L. Tagliacozzo, and P. Calabrese, Phys. Rev. B 81, 060411 (2010).
- (16) F. Gliozzi and L. Tagliacozzo, J. Stat. Mech.: Theor. Exp., P01002 (2010).
- (17) V. Alba, L. Tagliacozzo, and P. Calabrese, J. Stat. Mech.: Theor. Exp. ,P06012 (2011).
- (18) V. Alba, J. Stat. Mech.: Theor. Exp. P05013 (2013).
- (19) M. B. Hastings, I. Gonzalez, A. B. Kallin, and R. G. Melko, Phys. Rev. Lett. 104, 157201 (2010).
- (20) R. G. Melko, A. B. Kallin, and M. B. Hastings, Phys. Rev. B 82, 100409 (2010).
- (21) R. R. P. Singh, M. B. Hastings, A. B. Kallin, and R. G. Melko, Phys. Rev. Lett. 106, 135701 (2011).
- (22) S. V. Isakov, M. B. Hastings, and R. G. Melko, Nature Physics 7, 772 (2011).
- (23) S. Humeniuk and T. Roscilde, Phys. Rev. B 86, 235116 (2012).
- (24) R. K. Kaul, R. G. Melko, and A. W. Sandvik, Ann. Rev. Cond. Matt. Phys. 4, 179 (2013).
- (25) S. Inglis and R. G. Melko, Phys. Rev. E 87, 013306 (2013).
- (26) J. Iaconis, S. Inglis, A. B. Kallin, and R. G. Melko, Phys. Rev. B 87, 195134 (2013).
- (27) C. M. Herdman, P. N. Roy, R. G. Melko, and A. Del Maestro, Phys. Rev. B 89, 140501 (2014).
- (28) Y. Zhang, T. Grover, and A. Vishwanath, Phys. Rev. Lett. 107, 067202 (2011).
- (29) N. M. Tubman and J. McMinis, arXiv:1204.4731.
- (30) J. McMinis and N. M. Tubman, Phys. Rev. B 87, 081108(R) (2013).
- (31) T. Grover, Phys. Rev. Lett. 111, 130402 (2013).
- (32) F. F. Assaad, T. C. Lang, and F. P. Toldin, Phys. Rev. B 89, 125121 (2014).
- (33) P. Broecker and S. Trebst, J. Stat. Mech. P08015 (2014).
- (34) L. Wang and M. Troyer, Phys. Rev. Lett. 113, 110401 (2014).
- (35) J. Shao, E. Kim, F. D. M. Haldane, E. H. Rezayi, Phys. Rev. Lett. 114, 206402 (2015).
- (36) J. E. Drut and W. J. Porter, Phys. Rev. B 92, 125126 (2015).
- (37) J. E. Drut and W. J. Porter, Phys. Rev. E 93, 043301 (2016).
- (38) W. J. Porter and J. E. Drut, arXiv:1605.07085.
- (39) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
- (40) U. Schollwöck, Annals of Physics 326, 96 (2011).
- (41) C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
- (42) J. Hide and V. Vedral, Phys. Rev. A 81, 062303 (2010).
- (43) G. Hummer and A. Szabo, Proc. Natl. Acad. Sci. USA 98, 3658 (2001).
- (44) J. Liphardt, S. Dumont, S. B. Smith, I. Tinoco Jr., and C. Bustamante, Science 296, 1832 (2002).
- (45) F. Douarche, S. Ciliberto, A. Petrosyan, and I. Rabbiosi, Europhys. Lett. 70, 593 (2005).
- (46) D. Collin, F. Ritort, C. Jarzynski, S. B. Smith, I. Tinoco Jr, and C. Bustamante, Nature 437, 231 (2005).
- (47) C. Bustamante, J. Liphardt, and F. Ritort, Phys. Today 58, 43 (2005).
- (48) V. Blickle, T. Speck, L. Helden, U. Seifert, and C. Bechinger, Phys. Rev. Lett. 96, 070603 (2006).
- (49) N. C. Harris, Y. Song, and C.-H. Kiang, Phys. Rev. Lett. 99, 068101 (2007).
- (50) O.-P Saira, Y. Yoon, T. Tanttu, M. Möttönen, D. V. Averin, and J. P. Pekola, Phys. Rev. Lett. 109, 180601 (2012).
- (51) C. Jarzynski, Annu. Rev. Cond. Matt. Phys. 2, 329 (2011).
- (52) U. Seifert, Rep. Prog. Phys. 75, 126001 (2012).
- (53) A. Shuoming, J.-N. Zhang, M. Um, D. Lv, Y. Lu, J. Zhang, Z.-Q. Yin, H. T. Quan, and K. Kim, Nat. Phys. 11, 193 (2015).
- (54) M. Caselle, G. Costagliola, A. Nada, M. Panero, and A. Toniato, Phys. Rev. D 94, 034503 (2016).
- (55) H. B. Callen and T. A. Welton, Phys. Rev. 83, 34 (1951).
- (56) C. Jarzynski, Phys. Rev. E 73, 046105 (2006).
- (57) P. Di Francesco, P. Mathieu, and D. Senechal, Conformal Field Theory (Springer-Verlag, New York, 1997).
- (58) R. H. Swendsen and N.-S. Wang, Phys. Rev. Lett. 58, 86 (1987).
- (59) C.-M. Chung, V. Alba, L. Bonnes, P. Chen, A. M. Läuchli, Phys. Rev. B 90, 064401 (2014).
- (60) G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 (2002).
- (61) M. B. Plenio, Phys. Rev. Lett. 95, 090503 (2005); J. Eisert, quant-ph/0610253.
- (62) P. Calabrese, J. Cardy, and E. Tonni, Phys. Rev. Lett. 109, 130502 (2012).
- (63) R. Fan, P. Zhang, H. Shen, H. Zhai, arXiv:1608.01914.
- (64) J. Li, R. Fan, H. Wang, B. Ye, B. Zeng, H. Zhai, X. Peng, J. Du, arXiv:1609.01246
- (65) P. W. Kastaleyn and C. M. Fortuin, J. Phys. Soc. Jpn. (Suppl.) 26, 11 (1969).
Appendix A Cluster estimator for the Rényi entropies
For statistical models that are mappable to the random cluster model, i.e., whose partition function admits a representation á la Fortuin Kasteleyn fk-1969 (), it is possible to rewrite (7) in terms of cluster variables (see Ref. caraglio-2008, and alba-2013, for similar results in equilibrium Monte Carlo simulations). Typically, this leads to improved estimators with much smaller variance, as compared with the standard one expressed in terms of spin variables. The physical reason is that each fixed cluster configuration corresponds to many spin configurations. To be concrete, we focus on the two-dimensional Ising model on a square lattice, which is defined by the Hamiltonian
where and is the coupling strength. The partition function of the Ising model at inverse temperature can be expressed using the Fortuin Kasteleyn representation as
Here the sum is over the set of subgraphs on the lattice, i.e., arbitrary sets of lattice sites, denotes the (activated) links connecting points in (it is for the other links), while is the number of connected components of (clusters). In (11), we defined .
To calculate , one first considers a stack of disconnected replicas of the model and of (11). Each configuration of activated links (cluster configuration) on the disconnected sheets is also a valid one on the -sheets Riemann surface (see Figure 2), in which the copies are connected. The only difference is that now activated links across the branch cut connect sites on different replicas. Importantly, the same configuration of activated links corresponds to different numbers of clusters on the independent sheets and the Riemann surface. The ratio of partition functions in (2) is written as
where, given a fixed set of activated links, is the total number of clusters on the -sheets Riemann surface with cut , while is the total number of clusters for the model on the independent sheets. In (12) the statistical average is performed using the model defined on the independent sheets. Interestingly, from formula (12) one has that the entropies depend only on the number of clusters.
In our out-of-equilibrium setup the Monte Carlo dynamics happens on the coupled sheets (see Figure 2). Clearly, the Fortuin-Kasteleyn representation (11) still holds. However, in contrast to (11), the probability of activating the links across the branch cuts are dynamical variables (cf. (5)). On the other hand, since the change in energy (cf. (7)) is calculated at fixed fields configuration, the number of clusters does not change between two consecutive Monte Carlo updates. Thus, one has
In the definition of the product is over all the links across the branch cut, which connect sites both on the same sheet and on different ones, , with as defined in (5), and . For small , which is always the case in the simulations, from (14) one obtains
The strategy to use (14) and (13) in Monte Carlo simulations is straightforward using the Swendsen-Wang algorithm. Between two successive update steps of the Monte Carlo one has to measure both and , which are then substituted in (14). Note that in (14), the sum is over the links crossing the branch cuts, and only if the link is active. Thus, the infinitesimal work obtained from (14) is integrated on the Monte Carlo history, to obtain the exponent in (13).
The validity and efficiency of (14) is investigated in Figure 7. Panels (a) and (b) show the Rényi entropy calculated using (6) and (14), respectively. The data are for the critical Ising model on the lattice with and , and for a subystem with sites. In the Figure, is plotted versus the inter-replica coupling . In both panels the data are obtained from simulations with Monte Carlo steps. The different symbols correspond to different realization of the out-of-equilibrium protocol. Surprisingly, the data obtained using the cluter estimator (14) (panel(b)) show much larger fluctuations between different realizations as compared with the local estimator (6) (panel (a)). This is in costrast with equilibrium Monte Carlo simulations for calculating the Rényi entropies, where cluster estimators performs better than local ones caraglio-2008 ()