Measuring Rényi entanglement entropy with high efficiency and precision in quantum Monte Carlo simulations

The study investigates how to efficiently and precisely compute Rényi entanglement entropy (EE) in interacting quantum many-body systems, where EE encodes universal information about phases and phase transitions. In (2+1)d, EE generally follows an area law with universal subleading terms: logarithmic corrections whose coefficients depend on geometry and intrinsic properties (e.g., CFT data), and topological entanglement entropy (TEE) that equals the logarithm of the total quantum dimension in gapped topological phases. Obtaining these universal coefficients requires finite-size scaling and has been challenging in Quantum Monte Carlo (QMC) due to unstable estimators, enlarged configuration spaces with replicas, and high computational cost—especially in frustrated or multi-spin models. The authors aim to overcome these obstacles by combining nonequilibrium measurements based on Jarzynski’s equality with an increment (divide-and-conquer) protocol to create a nonequilibrium increment method that improves speed, precision, and scalability for extracting universal EE contributions across diverse (2+1)d systems.

Prior work established universal EE terms in critical and gapped phases and in systems with spontaneous symmetry breaking. QMC has enabled Rényi EE measurements, but system sizes and precision were limited by estimator instability and cost, particularly in frustrated and fermionic systems. Approaches include equilibrium replica-based estimators, Wang–Landau and linked-cluster methods, and nonequilibrium Jarzynski-based estimators, as well as the increment trick for equilibrium measurements. For TEE, planar constructions (Kitaev–Preskill/Levin–Wen) cancel leading area terms but suffer from severe finite-size effects on lattices, often underestimating the expected 2γ on torus geometries. DMRG studies on cylinders frequently identify minimum entropy states (MES) and observe γ = ln 2 in Z₂ topological phases. However, controlled extraction of 2γ on torus geometries via QMC for frustrated models (e.g., Kagome spin liquids) has remained elusive due to limited accessible system sizes and temperatures.

The method computes the 2nd Rényi entanglement entropy S^(2) via a nonequilibrium increment protocol within finite-temperature stochastic series expansion (SSE) QMC. Using the replica trick, one defines Z^(2) as the partition function of two replicas with the entanglement region A glued across replicas, while the complement Ā remains independent. Introduce a λ-parametrized family Z^(2)(λ) that interpolates between unglued (λ=0) and fully glued (λ=1) ensembles by weighting configurations with g_λ(A,N_B)=λ^{N_B}(1−λ)^{N_A−N_B}, where N_B counts sites of A connected across replicas. The Rényi entropy equals S_R = (1/β) ∫_0^1 dλ ∂_λ ln Z^(2)(λ) = ln ⟨e^{−β W}⟩, with W the nonequilibrium work accumulated along a protocol λ:0→1, using Jarzynski’s equality ΔF=−β^{-1} ln⟨e^{−β W}⟩. To enhance accuracy and reduce bias from finite quench time, the λ-interval is split into N increments of width h=1/N, inserting identities Z^(2)(kh)/Z^(2)(kh), and performing independent nonequilibrium paths on each sub-interval λ∈[kh,(k+1)h]. The total entropy is S_R = (1/β) Σ_{k=0}^{N−1} ln ⟨e^{−β W_k}⟩, where W_k is the work for the k-th segment. Algorithmic steps: (1) Thermalize a configuration for Z^(2)(0) and duplicate to create two replicas. (2) Launch N parallel processes; process k starts at λ_k=kh. (3) For sites in A, topology updates (joining/leaving the glued geometry) are proposed with probabilities P_join = min{1, λ/(λ−λ_k)} and P_leave = min{1, (1−λ)/(1−λ_k)}, followed by standard MC updates on each replica. (4) Thermalize at fixed λ=λ_k. (5) Begin nonequilibrium ramp: increase λ by Δλ each step, record g_{λ(t_m)}(A,N_B(t_m)) (or equivalently accumulate Δ ln g), perform MC sweeps, and repeat until λ reaches (k+1)h. (6) Compute ⟨e^{−β W_k}⟩ over paths for each k, and sum contributions over k to obtain S_R. This divide-and-conquer protocol effectively lengthens the total quench time by factor N without increasing wall-clock time, improving convergence and reducing systematic bias. Implementation details: SSE QMC in S^z basis at inverse temperature β chosen proportional to system size to approximate ground-state properties (e.g., 1/T=L or fixed below anyon gaps). Entanglement cuts avoid corners when targeting area law + TEE, or include controlled corners to extract corner logarithms. The method applies to complex Hamiltonians (dimerized Heisenberg, Z₂ gauge models, Kagome multi-spin interactions) and large lattices. Convergence is benchmarked by varying total quench time and number of increments N, showing faster convergence and smaller error bars for larger N at fixed computational cost.

- Heisenberg antiferromagnet with spontaneous SU(2) symmetry breaking (Goldstone modes): For a torus L×(L/2) with cylindrical bipartition (no corners), S^(2)(L) fits S_2(L)=0.092(1) L + 1.00(9) ln(L) − 1.63(3) over L ∈ [40,160], yielding logarithmic coefficient s=1.00(9), matching the prediction s=N_G(d−1)/2=1 (N_G=2 in d=2). This identifies the expected Goldstone contribution with high precision using standard SSE (1/T=L) without specialized bases. - O(3) quantum critical point (columnar-dimer J1–J2 Heisenberg model on square lattice): With a square entangling region (four corners, boundary length 2L) and careful placement to avoid crossing strong J2 bonds, the fit gives S^(2)(L)=0.167(2) L − 0.081(4) ln(L) − 0.124(7), extracting the universal corner log coefficient s=0.081(4), consistent with O(3) CFT expectations and improving on previous system sizes and precision. - Z₂ topological order, checkerboard toric-code-like toy model: DMRG on cylinders (L_x=72, L_y=6–12) yields γ=0.69(1)≈ln 2; QMC on tori with varying L shows Rényi EE extrapolates to 2γ=1.35(8), consistent with 2 ln 2. DMRG diagnostics of the deconfined–confined transition locate h_c≈0.333 via plaquette energy slope and entropy peaks, consistent with literature. - Minimum entropy state (MES) selection: For bipartition of a torus into two cylinders (two noncontractible boundaries), the expected TEE is 2γ=2 ln 2 only for MES ground states. Evidence indicates the QMC sampling, like DMRG on cylinders, preferentially accesses MES, enabling recovery of the full 2γ in torus geometries. - Kagome Z₂ quantum spin liquid (BFG model): For J_⊥/J_z=0.0625 (inside QSL) at T^{-1}=480 (below anyon gaps), with torus geometries of aspect ratios L_x/L_y=2 and 1 and cylindrical bipartitions without corners, S^(2)(L)=2α L_y − 2γ shows linear scaling with converged intercept 2γ=1.4(2)≈2 ln 2, across L_y=12–16 (type-I) and L_y=6–10 (type-II). Access to sufficiently large sizes (L_y≳12) is crucial to reach the asymptotic TEE, explaining prior failures using Levin–Wen constructions. - Algorithmic performance: On a 32×16 Heisenberg system with A=16×16, the original single-path (N=1) nonequilibrium method converges to S^(2)=4.84(1) only at long quench times; the increment method (N=24) achieves faster convergence and smaller error bars at comparable computational cost by effectively extending total quench time via parallel increments.

The results demonstrate that the nonequilibrium increment method enables precise extraction of universal EE corrections—logarithmic terms linked to Goldstone modes and corner contributions at CFT criticality, and TEE in Z₂ topological phases—from large-scale QMC simulations across challenging (2+1)d models. By stabilizing Jarzynski-based estimators through incremental, parallelized quenches, the method overcomes limitations of equilibrium estimators and previous nonequilibrium protocols, allowing accurate finite-size scaling beyond leading area-law terms. The recovery of 2γ=2 ln 2 for torus bipartitions indicates effective sampling of minimum entropy states and provides unambiguous evidence of Z₂ topological order, including in the frustrated Kagome BFG model where prior QMC methods struggled. These advances broaden the scope of entanglement-based diagnostics in numerical many-body studies, offering a practical pathway to access universal data (CFT coefficients, quantum dimensions) and to inform experimental efforts in identifying topological and highly entangled phases. The approach is adaptable to other complex systems, including interacting fermionic models and gauge-matter theories, and complements newer probes (e.g., disorder operators, higher-form symmetry diagnostics) that capture nonlocal order and emergent phenomena.

The study introduces a nonequilibrium increment method for measuring Rényi entanglement entropy in QMC with markedly improved efficiency and precision. Applied to several prototypical (2+1)d systems, it (i) verifies the universal logarithmic coefficient from Goldstone modes in the Heisenberg antiferromagnet, (ii) determines the universal corner log coefficient at an O(3) quantum critical point, and (iii) unambiguously extracts the full topological entanglement entropy 2γ=2 ln 2 in Z₂ topological phases, including the Kagome Z₂ quantum spin liquid. The method’s divide-and-conquer nonequilibrium protocol yields faster convergence and reduced bias at fixed computational resources, enabling large system sizes and controlled finite-size scaling. Future directions include extending the approach to interacting fermionic systems (e.g., Gross-Neveu and QED₃ universality classes), non-Fermi liquids and quantum critical metals, and integrating with other nonlocal probes (disorder operators, higher-form symmetry diagnostics) to deepen understanding of strongly entangled quantum matter.

- Nonequilibrium protocol sensitivity: Accurate Jarzynski estimates require sufficiently slow quenches; too-fast ramps bias results due to incomplete topology changes and integral discretization. The increment method alleviates but does not eliminate this, necessitating convergence checks versus quench time and number of increments N. - Ground-state sector selection (MES): Extracted TEE on torus depends on the ground state; while evidence suggests QMC sampling finds MES yielding 2γ, this may be model- and protocol-dependent. - Finite-size and temperature constraints: Reliable extraction of universal terms requires large boundary lengths (e.g., L_y≳12 in Kagome QSL) and temperatures below excitation gaps (anyon gaps in topological phases). Smaller systems or higher temperatures can obscure asymptotic scaling. - Geometry-dependent corrections: Corners and strong-bond inhomogeneities can introduce additional finite-size effects; careful choice of entanglement cuts is needed to isolate targeted universal terms. - Applicability to sign-problem systems: While demonstrated for sign-free bosonic/spin models, extension to generic fermionic/sign-problem systems may require additional algorithmic advances and higher computational cost.