Controlled and scalable variance computation for iPEPS ground states and steady states: cell-size and boundary-dimension requirements, long-distance contraction drift, and direct Liouvillian-norm minimization

The large-cell corner transfer matrix renormalization group (LC-CTMRG) method recently made it possible to compute the energy variance $\mathrm{Var}(E)=\langle H^{2}\rangle-\langle H\rangle^{2}$ of infinite projected entangled-pair states (iPEPS) accurately enough to extrapolate ground-state energies of two-dimensional quantum systems to the exact zero-variance limit, and extended the same machinery to the quantity $\epsilon=\langle\langle\mathcal{L}^{\dagger}\mathcal{L}\rangle\rangle$ that measures how close an approximate steady state of an open quantum system is to the true one. Four prerequisites for turning this framework into a routinely reliable and scalable tool remain open: (1) the resource requirements — the minimal cell size $L(D)$, set by the bond-dimension-dependent effective correlation length, and the contraction boundary dimension $\chi(D,L)$ — have not been systematically characterized beyond the models benchmarked so far; (2) small $\chi$ induces an artificial additive drift of summed correlators at large cell sizes, and no systematic strategy exists to detect and quantify this long-distance contraction noise; (3) direct variational optimization of the steady state by minimizing $\epsilon$ via automatic differentiation is considered computationally daunting because of the many terms in the cell, and no scalable implementation exists; (4) the Liouvillian correlator $\langle\langle\mathcal{L}^{\dagger}(r)\mathcal{L}_{0}\rangle\rangle$ need not vanish at large distance for approximate steady states, so the finite-cell truncation justification available for ground-state variance has no counterpart for open systems, and no extrapolation protocol in $\epsilon$ has been established. Resolving these questions is what stands between the variance and $\epsilon$ as heuristically used benchmarks and as controlled accuracy certificates for two-dimensional tensor-network simulation.

Background

Infinite projected entangled-pair states (iPEPS) are variational tensor-network wave functions that describe two-dimensional quantum lattice systems directly in the thermodynamic limit: a unit cell of rank-5 tensors with auxiliary bond dimension $D$ is repeated periodically over the lattice, and $D$ controls the variational accuracy. A central practical difficulty is extrapolating observables to the exact infinite-$D$ limit: the energy typically depends non-monotonically on $D$, so naive $1/D$ extrapolations are unreliable, and correlation-length-based extrapolations apply only to gapless Lorentz-invariant states. A route used for decades in variational and quantum Monte Carlo communities is to extrapolate in the energy variance $\mathrm{Var}(E)=\langle H^{2}\rangle-\langle H\rangle^{2}$, which is zero if and only if the state is an exact eigenstate; the energy is then fitted as a function of the variance and extrapolated to the zero-variance limit. For iPEPS, however, $\langle H^{2}\rangle$ requires correlators between every pair of Hamiltonian terms in the infinite lattice, which is why the earliest frameworks — the $H$-environment scheme of Corboz and the gradient-based scheme of Vanderstraeten et al. — were complex and converged slowly with respect to the contraction accuracy.

Contracting a two-dimensional tensor network requires an approximate scheme; the standard choice is the corner transfer matrix renormalization group (CTMRG), an iterative algorithm that encodes the environment of the infinite network in corner and edge tensors truncated to a boundary dimension $\chi$, with accuracy improving as $\chi$ grows. The large-cell (LC-)CTMRG method of Cortés Estay, Kamar, and Corboz computes the variance by enlarging the unit cell to $L\times L$, applying one Hamiltonian term $H_{0}$ at the center, and contracting the resulting overlap network with CTMRG sweeps moving outward from the center, so that each correlator $\langle H_{b}H_{0}\rangle$ is evaluated term by term (an idea introduced for two-point functions by Arias Espinoza and Corboz in the context of spectral functions). Because the CTMRG projectors 'see' the perturbation $H_{0}$, the variance converges at much smaller $\chi$ than before (around $\chi=30$ instead of $\chi>200$ at $D=4$), making variance extrapolation competitive: for the square-lattice Heisenberg antiferromagnet it reproduces the quantum Monte Carlo energy, and it underpinned a high-accuracy determination of the plaquette-to-spin-liquid transition in the Shastry-Sutherland model. Two finite parameters must nevertheless be chosen by hand: the cell size $L$, which must exceed the effective correlation length $\xi(D)$ that grows with $D$, and the boundary dimension $\chi$; moreover, small $\chi$ introduces a subtle artifact — truncation errors spoil the cancellation of the connected correlator (the correlator with its disconnected part subtracted) at long distances, producing an additive drift of the summed variance as $L$ grows, which the authors could only recommend to monitor and check.

The same machinery extends to open quantum systems. A Markovian environment is described by a Lindblad master equation $\dot{\rho}=\mathcal{L}[\rho]$ with a Liouvillian superoperator $\mathcal{L}$ (a linear map on density matrices built from local Hamiltonian and jump-operator terms); the steady state satisfies $\mathcal{L}|\rho_{s}\rangle\rangle=0$. In two dimensions the density matrix can be vectorized and represented as an infinite projected entangled-pair operator (iPEPO), and steady states are usually reached by long-time evolution using simple- or full-update schemes (Kshetrimayum, Weimer, and Orús; Czarnik, Dziarmaga, and Corboz; Mc Keever and Szymańska). The quantity $\epsilon=\langle\langle\rho|\mathcal{L}^{\dagger}\mathcal{L}|\rho\rangle\rangle$ per site is the Liouvillian analogue of the variance: it vanishes exactly on the steady state, can be computed with the LC-CTMRG machinery, and provides both a quality certificate for approximate steady states and a diagnostic to locate first-order transitions via hysteresis. For approximate steady states, however, the Liouvillian correlator $\langle\langle\mathcal{L}^{\dagger}(r)\mathcal{L}_{0}\rangle\rangle$ need not decay at large distance $r$, because $\langle\langle\mathcal{L}\rangle\rangle\neq 0$ — unlike the Hamiltonian case, where connected correlators of an exact eigenstate cluster — so the finite-$L$ truncation argument that justifies cutting off the variance sum at some cell size fails, and no extrapolation protocol in $\epsilon$ exists. In one dimension, variational steady states have been obtained since 2015 by minimizing a norm of the master equation (Weimer) or, equivalently, by computing the matrix-product-operator ground state of $\mathcal{L}^{\dagger}\mathcal{L}$ (Cui, Cirac, and Bañuls); the analogous direct minimization in two dimensions is what the LC-CTMRG paper leaves open, since gradients of $\epsilon$ — taken with the automatic-differentiation machinery used for ground-state energies (Liao et al.) — must be evaluated through the large-cell contraction for the many terms in the cell.

Problem Statement

The overarching objective is to turn the large-cell CTMRG variance framework into a controlled and scalable accuracy certificate for two-dimensional tensor-network simulation, valid both for ground states and for steady states of open systems. The source paper (Cortés Estay, Kamar, and Corboz, Phys. Rev. B 113, 245114 (2026)) explicitly raises direct steady-state optimization by minimizing $\epsilon$ (component (3) below) as a future direction; the following related open questions remain:

(1) Resource requirements. The cell size $L$ must exceed the effective correlation length $\xi(D)$, which grows with the bond dimension $D$; the source paper observes that the larger $D$ is, the larger the required cell size to reach convergence, and the longer the correlation length, the larger the required cell size to capture all relevant terms contributing to the variance. Determine, for at least two distinct model classes (one gapless and one gapped or frustrated), the minimal cell size $L(D)$ needed to capture the variance-contributing correlators to within a stated target accuracy of the variance (fixed here as a relative error of at most $10^{-3}$ in the variance density), and the contraction boundary dimension $\chi(D,L)$ required for the same target accuracy; both finite parameters currently must be chosen by hand, with their convergence checked case by case.

(2) Long-distance contraction drift. Small $\chi$ induces an artificial additive drift of summed correlators at large cell sizes: contraction errors lead to imperfect cancellation of the connected correlator at long distance, resulting in an additive contribution with increasing $L$, so that currently the convergence as a function of $L$ can only be monitored and checked for drifts by hand, especially when relatively small values of $\chi$ are used. Develop systematic strategies to detect and quantify this long-distance contraction noise at the stated variance tolerance, going beyond case-by-case manual monitoring.

(3) Direct steady-state optimization by minimizing $\epsilon=\langle\langle\mathcal{L}^{\dagger}\mathcal{L}\rangle\rangle$. This is the future direction explicitly raised by the source paper: it would be interesting to perform the variational optimization of the steady state directly by minimizing $\epsilon$, which is technically possible via automatic differentiation of the large-cell contraction, by analogy with ground-state energy minimization, but whose computational cost would be considerably higher because of the many terms involved in the cell; the paper notes that the gradient for each of these terms could be parallelized, making the approach in principle feasible. Making this direct minimization computationally practical and scalable in two dimensions remains open.

(4) Open-system correlator structure. As noted in a footnote of the source paper, the Liouvillian correlator $\langle\langle\mathcal{L}^{\dagger}(r)\mathcal{L}_{0}\rangle\rangle$ does not vanish at large distance for approximate steady states, because typically $\langle\langle\mathcal{L}\rangle\rangle\neq 0$ — unlike the Hamiltonian case, where connected correlators of an exact eigenstate cluster — so finite-cell estimates of $\epsilon$ lack the truncation justification available for the ground-state variance, and no extrapolation protocol in $\epsilon$ exists. How this affects the reliability of $\epsilon$ as an accuracy certificate for two-dimensional open systems remains open.

The four components are linked by the common objective stated above: components (1) and (2) concern the ground-state half of the framework, components (3) and (4) the open-system half, and together they determine whether variance-based accuracy assessment can be used routinely. All four components are required; each can be reported independently, but the task is resolved only when all four are addressed, and a component counts as addressed when its claim is supported by a concrete worked model or protocol.

References

1. Emilio Cortés Estay, Naushad A. Kamar, and Philippe Corboz, 'Accurate computation of the energy variance for ground states and ⟨⟨L†L⟩⟩ for steady states using iPEPS', Phys. Rev. B 113, 245114 (2026), DOI: 10.1103/fkyw-c5tn, arXiv:2511.22669, https://arxiv.org/abs/2511.22669
2. Philippe Corboz, Yining Zhang, Boris Ponsioen, and Frédéric Mila, 'Quantum spin liquid phase in the Shastry-Sutherland model revealed by high-precision infinite projected entangled-pair states', Phys. Rev. Lett. 136, 186701 (2026), DOI: 10.1103/67xm-5d8v, arXiv:2502.14091, https://arxiv.org/abs/2502.14091
3. Philipp Schmoll, Jan Naumann, Erik Lennart Weerda, Jens Eisert, and Yasir Iqbal, 'An extended non-magnetic phase in the spin-1/2 Heisenberg antiferromagnet from the ruby to the maple-leaf lattice', arXiv:2407.07145 (2024, v3 June 2026), https://arxiv.org/abs/2407.07145
4. Philippe Corboz, 'Variational optimization with infinite projected entangled-pair states', Phys. Rev. B 94, 035133 (2016), DOI: 10.1103/PhysRevB.94.035133, arXiv:1605.03006, https://arxiv.org/abs/1605.03006
5. Laurens Vanderstraeten, Jutho Haegeman, Philippe Corboz, and Frank Verstraete, 'Gradient methods for variational optimization of projected entangled-pair states', Phys. Rev. B 94, 155123 (2016), DOI: 10.1103/PhysRevB.94.155123, arXiv:1606.09170, https://arxiv.org/abs/1606.09170
6. Juan Diego Arias Espinoza and Philippe Corboz, 'Spectral functions with infinite projected entangled-pair states', Phys. Rev. B 110, 094314 (2024), DOI: 10.1103/PhysRevB.110.094314, arXiv:2405.10628, https://arxiv.org/abs/2405.10628
7. Hai-Jun Liao, Jin-Guo Liu, Lei Wang, and Tao Xiang, 'Differentiable programming tensor networks', Phys. Rev. X 9, 031041 (2019), DOI: 10.1103/PhysRevX.9.031041, arXiv:1903.09650, https://arxiv.org/abs/1903.09650
8. Hendrik Weimer, 'Variational principle for steady states of dissipative quantum many-body systems', Phys. Rev. Lett. 114, 040402 (2015), DOI: 10.1103/PhysRevLett.114.040402, arXiv:1409.8307, https://arxiv.org/abs/1409.8307
9. Jian Cui, J. Ignacio Cirac, and Mari Carmen Bañuls, 'Variational matrix product operators for the steady state of dissipative quantum systems', Phys. Rev. Lett. 114, 220601 (2015), DOI: 10.1103/PhysRevLett.114.220601, arXiv:1501.06786, https://arxiv.org/abs/1501.06786
10. Augustine Kshetrimayum, Hendrik Weimer, and Román Orús, 'A simple tensor network algorithm for two-dimensional steady states', Nat. Commun. 8, 1291 (2017), DOI: 10.1038/s41467-017-01511-6, arXiv:1612.00656, https://arxiv.org/abs/1612.00656
11. C. Mc Keever and M. H. Szymańska, 'Stable iPEPO tensor-network algorithm for dynamics of two-dimensional open quantum lattice models', Phys. Rev. X 11, 021035 (2021), DOI: 10.1103/PhysRevX.11.021035, arXiv:2012.12233, https://arxiv.org/abs/2012.12233
12. Dainius Kilda, Alberto Biella, Marco Schiró, Rosario Fazio, and Jonathan Keeling, 'On the stability of the infinite projected entangled pair operator ansatz for driven-dissipative 2D lattices', SciPost Phys. Core 4, 005 (2021), DOI: 10.21468/SciPostPhysCore.4.1.005, arXiv:2012.03095, https://arxiv.org/abs/2012.03095
13. Piotr Czarnik, Jacek Dziarmaga, and Philippe Corboz, 'Time evolution of an infinite projected entangled pair state: An efficient algorithm', Phys. Rev. B 99, 035115 (2019), DOI: 10.1103/PhysRevB.99.035115, arXiv:1811.05497, https://arxiv.org/abs/1811.05497
14. Yining Zhang, Qi Yang, and Philippe Corboz, 'Accelerating two-dimensional tensor network contractions using QR decompositions', Phys. Rev. B 113, L201106 (2026), DOI: 10.1103/ptls-kr9z, arXiv:2505.00494, https://arxiv.org/abs/2505.00494
15. Qi Yang and Philippe Corboz, 'Efficient iPEPS simulation on the honeycomb lattice via QR-based corner transfer matrix renormalization group', Phys. Rev. B 113, 085109 (2026), DOI: 10.1103/9gmp-byx8, arXiv:2509.05090, https://arxiv.org/abs/2509.05090
16. Yining Zhang and Philippe Corboz, 'Fast two-dimensional tensor-network contraction via subspace iteration', arXiv:2607.15158 (2026), https://arxiv.org/abs/2607.15158