Infinite Projected Entangled-Pair States (iPEPS)

• iPEPS is a variational tensor-network ansatz that represents infinite 2D quantum many-body states using a periodically repeated unit cell of local tensors.
• Its contraction techniques, such as CTMRG, approximate the infinite network’s environment to accurately extract physical observables like energy and correlation lengths.
• Optimization methods including imaginary-time evolution and variational minimization systematically improve accuracy by controlling the virtual bond dimension D.

Infinite projected entangled-pair states (iPEPS) are a variational tensor-network ansatz for two-dimensional quantum many-body wavefunctions defined directly in the thermodynamic limit. The construction replaces an infinite lattice state by a periodic tiling of a finite unit cell of local tensors whose virtual indices are contracted according to the lattice geometry; the expressive power is controlled systematically by the virtual bond dimension $D$, with $D=1$ yielding a product state and larger $D$ admitting progressively more entanglement (Corboz, 2015).

1. Tensor-network ansatz and state representation

On the square lattice, the basic local object is a rank-5 tensor $A$ with one physical index and four virtual indices connecting nearest neighbors. In the notation used across the literature, a site tensor may be written as $A^{s}_{lurd}$, where $s$ has local dimension $d$ and $(l,u,r,d)$ are virtual indices of dimension $D$. For a translationally invariant phase, a single tensor can be repeated on every site; for states with spontaneous symmetry breaking, a larger periodic supercell is used, such as two tensors for Néel order or larger cells for stripe states (Corboz, 2015).

The wavefunction is obtained by contracting all neighboring virtual indices across the infinite lattice,

$$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$

or, in periodic-unit-cell notation, by repeating a finite set of tensors $D=1$0 over the plane (Corboz, 2015). The number of variational parameters per tensor scales as $D=1$1. This finite-parameter description is the central reason iPEPS can work directly at infinite system size while remaining systematically improvable.

The same formalism extends beyond the square lattice by coarse graining or embedding. Honeycomb, kagome, triangular, square-kagome, ruby, and triangle-honeycomb lattices can all be represented within a square-lattice CTMRG framework by blocking several microscopic degrees of freedom into an enlarged on-site physical space, or by using iPESS-like decompositions and then mapping to square-lattice contractions (Naumann et al., 2023). In three dimensions, the analogous infinite PEPS construction on the cubic lattice uses rank-7 tensors with one physical and six virtual indices, again repeated periodically in the thermodynamic limit (Vlaar et al., 2021).

A particularly important generalization is the translationally invariant representation of incommensurate order. In the spiral construction, one keeps a single-site translationally invariant iPEPS $D=1$2 and applies a position-dependent physical rotation,

$D=1$3

so that the spiral wavelength is encoded in a continuous wavevector $D=1$4 rather than in an enlarged tensor unit cell (Hasik et al., 2023). This decouples the spatial period of the order from the tensor-network unit cell.

2. Contraction, environments, and transfer matrices

Expectation values in iPEPS require contraction of an infinite two-dimensional tensor network, which is not performed exactly in practice. The dominant approach in the cited works is corner transfer matrix renormalization group (CTMRG), or closely related CTM variants, where the infinite environment of a local region is approximated by four corner tensors and four edge tensors with an environment bond dimension $D=1$5 controlling the contraction accuracy (Corboz, 2015).

The contraction is typically formulated in the double layer, obtained by contracting bra and ket tensors over physical indices. For a local tensor $D=1$6, the associated double-layer tensor has effective virtual dimension $D=1$7 on each leg. In standard single-site translationally invariant CTMRG, the environment consists of four corners $D=1$8 and four edges $D=1$9; for larger unit cells, separate environments are maintained for each inequivalent tensor position (Hasik et al., 2023). Directional CTMRG updates iteratively absorb rows or columns of double-layer tensors and truncate the enlarged boundary spaces back to $D$0 by projectors constructed from singular-value decompositions of local environment patches (Naumann et al., 2023).

The transfer-matrix viewpoint is central. Row-to-row or column-to-column transfer matrices built from the double-layer network determine long-distance properties, and the leading correlation length is extracted from the dominant eigenvalues,

$D$1

or equivalently with the largest and second-largest eigenvalues under the indexing convention used in a given paper (Hasik et al., 2023). This quantity is both a physical diagnostic and a control parameter for finite-correlation-length scaling.

For standard ground-state calculations, $D$2 is increased until contraction errors are negligible compared to finite-$D$3 errors; in the 2D Hubbard study, values up to several hundreds were sufficient for this purpose (Corboz, 2015). In the anisotropic triangular-lattice spiral simulations, optimization used $D$4 up to $D$5, while observables were measured at larger $D$6 until stable (Hasik et al., 2023). The recent GPU-oriented implementation Ace-TN keeps the same CTMRG structure but accelerates projector construction and contractions, using randomized SVD and simple multi-GPU parallelization for directional moves (Richards et al., 18 Mar 2025).

A recurrent practical subtlety is gauge freedom. Unlike matrix product states, PEPS do not admit an exact global canonical form. Nevertheless, local gauge fixing, bond balancing, and phase fixing in the SVDs used to construct CTM projectors substantially improve numerical conditioning and fixed-point convergence (Phien et al., 2015). This is not merely a technical detail: ill-conditioned norm tensors can destabilize optimization and may even generate spurious minima when $D$7 is too small (Hasik et al., 2022).

3. Optimization schemes and tensor updates

The oldest widely used optimization route is imaginary-time evolution. For a local Hamiltonian $D$8 with nearest-neighbor terms, one Trotterizes $D$9 into small local gates,

$A$0

with a second-order decomposition implemented by reversing the gate order on alternating steps (Corboz, 2015). Acting with a two-site gate enlarges the corresponding virtual bond, after which one truncates back to the target bond dimension $A$1.

The main truncation strategies differ by how much of the global environment they use. In the simple update, local diagonal bond-weight matrices approximate the environment, leading to low cost but limited accuracy in strongly entangled or critical regimes (Phien et al., 2015). In the full update, the truncation is determined by minimizing a cost function with the full CTM environment. In the Hubbard formulation, the local truncation minimizes

$A$2

or equivalently the associated quadratic form $A$3 built from overlaps of the post-gate and truncated states (Corboz, 2015). This environment-aware truncation substantially improves fidelity to the evolved state.

Several refinements address the computational bottlenecks of the full update. Fast full update evolves the tensors and the environment simultaneously, performing only one CTM iteration after each gate instead of fully reconverging the environment; this yields major speedups and improved stability, especially when combined with local gauge fixing of the norm tensor (Phien et al., 2015). Loop update and full loop update replace bondwise truncation by cyclic optimal truncation around a plaquette, explicitly targeting redundant internal correlations generated on closed loops; benchmark results show improved critical behavior relative to conventional simple and full updates (Zheng et al., 2019).

A second major optimization family is direct variational minimization of the Rayleigh quotient

$A$4

In the CTM-based variational scheme, local updates solve generalized eigenvalue problems of the form

$A$5

for effective Hamiltonian and norm matrices obtained by summing Hamiltonian contributions over the infinite lattice with CTM environments (Corboz, 2016). This global energy-minimization strategy was shown to give substantially more accurate energies than the then-standard imaginary-time methods at similar computational cost.

Automatic differentiation has generalized this approach further. In AD-based iPEPS, the entire energy evaluation—including CTMRG—is treated as a differentiable program, and gradients are propagated through the converged CTM fixed point rather than through an unrolled sequence of CTM iterations. The review literature emphasizes backward-mode AD, fixed-point differentiation, and optimizers such as nonlinear conjugate gradient and L-BFGS as the current state of the art for variational ground-state simulations (Naumann et al., 2023).

4. Extrapolation, finite-correlation-length scaling, and critical behavior

Finite $A$6 and finite $A$7 introduce distinct cutoffs, and a central methodological question is how to extrapolate finite-iPEPS data toward the exact thermodynamic answer. A widely used older heuristic was to fit observables as functions of $A$8, but several papers show that this can be irregular and unreliable, especially in competitive or critical regimes (Corboz, 2015).

A more controlled alternative is the truncation-error extrapolation introduced for full-update imaginary-time evolution. Defining

$A$9

where $A^{s}_{lurd}$0 is the full-update truncation cost and $A^{s}_{lurd}$1 the Trotter step, one finds empirically that the variational energy often depends much more smoothly on $A^{s}_{lurd}$2 than on $A^{s}_{lurd}$3 (Corboz, 2015). In the 2D Hubbard benchmarks, third-order polynomial fits in $A^{s}_{lurd}$4 produced near-exact energies at $A^{s}_{lurd}$5 and excellent agreement with AFQMC at half filling, while also resolving the competition between vertical and diagonal stripe states at $A^{s}_{lurd}$6, $A^{s}_{lurd}$7.

A different, now standard, route is finite-correlation-length scaling (FCLS). For Lorentz-invariant critical points, the optimal finite-$A^{s}_{lurd}$8 iPEPS approximation exhibits a finite correlation length even at criticality; this provides the two-dimensional analogue of one-dimensional finite-entanglement scaling (Corboz et al., 2018). In this framework, the finite system size $A^{s}_{lurd}$9 of finite-size scaling is replaced by the intrinsic iPEPS correlation length $s$0. For an order parameter $s$1 near criticality,

$s$2

and at criticality $s$3 (Corboz et al., 2018). The same logic yields practical crossing and collapse procedures that avoid higher moments of the order parameter. Applied to interacting spinless fermions on the honeycomb lattice, this gave $s$4, $s$5, and $s$6, in agreement with QMC (Corboz et al., 2018).

In ordered Goldstone phases, FCLS takes a particularly simple form. For the anisotropic triangular Heisenberg model, the energy and magnetization were extrapolated as

$s$7

with $s$8 in the Néel phase and both $s$9 and $d$0 examined in the spiral phase (Hasik et al., 2023). For coupled spin-$d$1 ladders, separate transfer-matrix correlation lengths $d$2 and $d$3 were used in an anisotropic FCLS analysis, which remained accurate down to a narrow interval around the critical point (Hasik et al., 2022).

Imaginary-time evolution itself can also be used to extract universal critical dynamics before convergence to the ground state. In the 2D transverse-field Ising model, a fully polarized initial state gives

$d$4

while a weakly magnetized initial product state yields the critical initial-slip law

$d$5

Using full-update iPEPS, the reported extrapolations gave $d$6 and $d$7, with the measured slopes close to the 3D Ising value $d$8 (Lin et al., 14 Nov 2025). This establishes that universal short-time imaginary-time scaling can be resolved at moderate bond dimension.

These extrapolation methods also delimit what finite-$d$9 iPEPS cannot yet do cleanly. Near continuous transitions, optimization becomes harder, finite-$(l,u,r,d)$0 effects can generate shallow spurious minima, and the accessible $(l,u,r,d)$1 may remain below the physical correlation length in a narrow critical window (Hasik et al., 2022). For chiral phases, finite-$(l,u,r,d)$2 PEPS are additionally constrained by no-go theorems, which reappear numerically as small long-distance “gossamer” tails in otherwise exponentially decaying correlators (Chen et al., 23 Dec 2025).

5. Physical applications across ordered, frustrated, and fermionic systems

iPEPS was developed as a ground-state ansatz, and its core successes remain benchmark calculations for strongly correlated two-dimensional Hamiltonians. In the single-band 2D Hubbard model, iPEPS on the infinite square lattice gave lower variational energies than fixed-node Monte Carlo in the doped, strongly correlated regime at $(l,u,r,d)$3, $(l,u,r,d)$4, and $(l,u,r,d)$5 (Corboz, 2015). At half filling, the method reproduced AFQMC benchmarks with relative energy errors of order $(l,u,r,d)$6 at $(l,u,r,d)$7 and about $(l,u,r,d)$8 at $(l,u,r,d)$9 already before extrapolation, and the $D$0-based analysis strongly favored a vertical superconducting stripe over a diagonal stripe in the doped case (Corboz, 2015).

For frustrated magnets with incommensurate order, the translationally invariant spiral construction removes the commensurability constraints of finite clusters. On the anisotropic triangular lattice, the optimized wavevector $D$1 evolves continuously from Néel order at $D$2 to $D$3 order at $D$4, while the phase diagram inferred from FCLS places a nonmagnetic regime between the ordered phases, with a quantum spin liquid region between $D$5 and $D$6 (Hasik et al., 2023). The same work showed that a single-site spiral iPEPS at the isotropic triangular point matches a conventional 3-site unit-cell iPEPS within $D$7 in energy and magnetization at $D$8.

For coupled two-leg ladders, U(1)-symmetric iPEPS captures the continuous transition between a Néel phase and a gapped paramagnet, with energy extrapolations consistent with QMC down to a narrow interval near the critical point $D$9 (Hasik et al., 2022). The earlier study of the same model identified a separate, more algorithmic issue: generic bond dimensions can be incompatible with emergent multiplet structures on virtual bonds, leading to spurious finite magnetization in the paramagnetic phase unless selected values of $$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$0 are used (Hasik et al., 2019). This does not invalidate iPEPS, but it shows that symmetry implementation and bond-sector design can be decisive in practice.

Beyond two spatial dimensions, the iPEPS formalism extends to infinite 3D lattices. For the cubic-lattice Heisenberg model, a boundary-iPEPS-plus-CTMRG contraction scheme and finite-cluster contractions yielded competitive results, including an extrapolated energy per site $$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$1, consistent with QMC $$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$2 (Vlaar et al., 2021). For the 3D Bose-Hubbard model, the same framework produced Mott-lobe boundaries close to QMC and competitive with B-DMFT (Vlaar et al., 2021).

The scope of iPEPS now includes fermionic topological phases. U(1)-symmetric fermionic iPEPS with swap-gate statistics were optimized for a $$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$3 fractional Chern insulator on the honeycomb Haldane model, with bond dimensions up to $$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$4 (Chen et al., 23 Dec 2025). That work identified a critical bond dimension $$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$5 above which the variational energy, charge homogeneity, and edge entanglement spectrum all become consistent with the target FCI phase. The equal-time Green’s function showed rapid exponential decay with $$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$6 plus a small long-distance tail, and the pair-correlation function matched the continuum $$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$7 Laughlin form with an effective magnetic length $$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$8 (Chen et al., 23 Dec 2025).

6. Dynamical, topological, and computational extensions

Although iPEPS originated as a ground-state method, several works have extended it to real-time, thermal, and open-system dynamics. A first-principles variational compression algorithm showed that real- and imaginary-time evolution can be simulated by applying small Trotter steps, producing an enlarged-bond iPEPS, and variationally compressing it back to the target $$|\Psi(A)\rangle = \sum_{\{s_i\}} \mathrm{Contraction}\!\left(\prod_i A^{s_i}_{l_i u_i r_i d_i}\right)|\{s_i\}\rangle,$$9 using CTMRG overlaps (Czarnik et al., 2018). A subsequent efficient formulation replaced the global fidelity objective by a local fidelity estimator and introduced unitary disentanglers on ancillas for thermal-state purifications, leading to accurate estimates of the critical temperature of the 2D transverse-field Ising model: about $D=1$00 accuracy at $D=1$01 and $D=1$02 at $D=1$03 (Czarnik et al., 2018).

Real-time iPEPS can also address disordered and spectroscopic problems. Disorder averaging for binary fields can be restored to translational invariance by coupling each physical spin to an ancilla spin prepared in $D=1$04, so that the disorder average becomes a clean expectation value in an enlarged Hilbert space (Hubig et al., 2018). Using this trick, short-time dynamics of a disordered square-lattice Heisenberg model showed a clear slowdown at strong disorder, with a percolation-based interpretation in the $D=1$05 limit (Hubig et al., 2018). Spectral functions are accessed by applying a local operator in a large unit cell, evolving the perturbed state in real time with fast-full update, measuring non-equal-time correlators by CTMRG at each step, and Fourier transforming them. In the 2D transverse-field Ising model this reproduced the main features of the dynamical structure factor, including the magnon line and higher-energy two-magnon structures, already at moderate $D=1$06 and unit-cell size (Espinoza et al., 2024).

Topological order can be extracted directly from iPEPS transfer matrices. For Abelian phases, infinite MPO projectors acting on virtual indices resolve torus states with definite anyon flux and permit direct extraction of modular $D=1$07 and $D=1$08 matrices from overlap data on rotated tori (Francuz et al., 2019). The method was then generalized to non-Abelian topological order by identifying MPO symmetries and their fusion rules, constructing central idempotents that project onto definite anyon sectors, and computing Renyi-2 topological entanglement entropies and modular matrices directly in the infinite-size limit (Francuz et al., 2020). For exact Fibonacci and Ising string-net tensors, the reported $D=1$09 and $D=1$10 matrices were recovered with machine-precision accuracy (Francuz et al., 2020).

At the same time, there are limits to what bulk iPEPS data alone can certify. In the anisotropic triangular study, the absence of magnetic and dimer order, short-ranged correlations, and competitive energies supported a quantum spin liquid, but the authors emphasized that iPEPS without explicit topological diagnostics cannot determine the topological order or definitively distinguish a gapped from a gapless QSL (Hasik et al., 2023). This is an important corrective to the common misconception that low energy alone fixes the phase.

The methodological ecosystem has also broadened. Automatic-differentiation-based variational iPEPS is now supported by libraries such as peps-torch, quimb, PEPSKit, variPEPS, and YASTN, with block-sparse symmetry support and fixed-point differentiation through CTMRG (Naumann et al., 2023). Ace-TN provides a GPU-first implementation of CTMRG and full or fast-full update, with randomized SVD and multi-GPU execution; the reported speedups reach about $D=1$11 for full update on a single GPU versus CPU at larger $D=1$12, and about $D=1$13–$D=1$14 on $D=1$15A100 versus a 20-core CPU in contraction-dominated regimes (Richards et al., 18 Mar 2025).

Taken together, these developments define iPEPS less as a single algorithm than as a family of thermodynamic-limit tensor-network methods. The unifying core remains the same: a finite set of local tensors, an approximate but systematically improvable contraction of the infinite network, and a bond dimension $D=1$16 that organizes the approximation. What has changed is the range of questions that can now be addressed with that core—from variational energies and order parameters to critical exponents, relaxation dynamics, spectral functions, anyon data, and fermionic chiral topological phases.