Entanglement renormalization circuits for $2d$ Gaussian Fermion States

The $1d$ GMPS [4] and GMERA procedure begins by examining the correlation matrix restricted to a sub-block of $B$ nearby sites, $C|_{B}$. Diagonalizing $C|_{B}$ yields a set of sub-block eigenvector, eigenvalue pairs: $\phi_{\alpha},0\leq\eta_{\alpha}\leq 1$. The degree to which block eigenorbital $\alpha$ is localized within the block is characterized by the Shannon entropy of its eigenvalue:

For a correlation matrix representing a locally-entangled state such as a ground-state of a local Hamiltonian, for a sufficiently large radius $r$, some eigenvalues will be very close to 0 or 1, corresponding to $S(\eta)\approx 0$, due to the limited entanglement structure. These eigenmodes, collectively referred to as frozen modes, as their occupation number has almost no quantum fluctuations, and they are essentially uncorrelated with the rest of the system outside the block in question.

Through a sequence of nearest-neighbor rotation gates, one can transform the basis to isolate these modes at the left boundary of the block. These rotations provide a local-circuit approximation to the the unitary that rotates an eigen-orbital of the full correlation matrix $C$ onto a single site, with the accuracy of the approximation controlled by the inter-block entanglement $S(\eta)$. The remaining modes with intermediate eigenvalues that are not close to $0,1$ are termed courier modes as they carry the quantum correlations between different blocks cannot be accurately disentangled by local operations acting within the block.

A particularly attractive feature of states prepared by sequential or MERA-type circuits is that their observables can be sampled using using a small number of qubits by leveraging qubit recycling techniques to opportunistically measuring qubits that have completed their role in the circuit, and resetting them to reuse in other parts of the circuit. For the $1d$ GMPS procedure, qubit reuse enables implementation with number of qubits scaling as the block size $B$ independent of system size [16]. Moreover, $1d$ MERA circuits (including GMERA) can be recast as seqential circuits, and implemented with number of physical qubits scaling as $\sim\log(L)$, i.e. logarithmically in the length of the system, $L$. These techniques can dramatically increase the size and complexity of systems that can be modeled with a given amount of quantum memory.

While the GMPS and GMERA methods work well in quasi-$1d$ geometries, $1d$ ground-states can generally be efficiently simulated by classical MPS methods, and quantum advantages on ground-state simulation problems likely reside in at least two spatial dimensions (as opposed to quantum dynamics simulation, which are classically difficult in any dimension). In the remainder of this work, we focus on generalizing the $1d$ GMERA methods to $2d$ GMERA methods. To start, we illustrate why a direct generalization of the $1d$ techniques generally leads to an uncontrolled build up of entanglement causing the compression scheme to fail. To address this issue, we next introduce an extra Wannierization procedure, to attempt to localize the entanglement produced within each local block during the GMERA procedure.

The applicability of Gaussian Matrix Product State (GMPS) methods can be best understood by first examining quasi-$1d$ systems, where the width of the system is much smaller than its length. In this geometry, when we process each block sequentially, the spatial spread of courier modes remains confined by the narrow width of the system. Crucially, each subsequent block completely encompasses the courier modes generated by the previous block as shown in Fig. 3(a), allowing for consistent and effective mode distillation throughout the procedure. This complete coverage ensures that we can maintain control over the entanglement structure as we progress along the quasi-$1d$ system.

A straightforward generalization of this approach fails for $2d$ systems. As depicted in Fig. 3, distilling orbitals in one block create courier modes with extended entanglement, which cannot be removed by subsequent blocks. Iterating the block-distillation procedure then inevitably leads to uncontrolled build-up of entanglement and a failure of the GMPS procedure. Below, we introduce a modified version of this scheme that introduces a Wannierization step to maintain the locality of courier modes to avoid this entanglement build-up. The modification is inspired by a Zipper Entanglement Renormalization procedure (ZER) [5], a recently-introduced holographic method that uses continuous-time Hamiltonian evolution to produce a quasi-local MERA-like unitary to disentangle Gaussian fermion states. Our $2d$ GMERA can be loosly thought of as a strictly local version of ZER that enables its implementation with gate-based quantum circuits.

For concreteness, we illustrate the technique first for a honeycomb, as illustrated in Fig. 4 for block radius $r=4$, and then comment on generalizations to a general $2d$ lattice. The honeycomb lattice 2d GMERA procedure begins with an enlarged unit cell of size $\frac{r}{2}a\times\frac{r}{2}a$ where $a$ is the lattice vector and $r$ is the block radius (Fig. 4(a)(i)). Each GMERA step consists of three sequential disentangling layers (Fig. 4(a)(ii)), with each layer processing a set of non-overlapping hexagonal blocks that are shifted relative to each other.

The procedure for each disentangling layer works as follows:

The procedure described above can be generalized to arbitrary lattices by adapting the following components:

1. 1.

   Block Structure: Define appropriate blocks with radius $r$ that cover the entire lattice. The number and arrangement of these blocks will depend on the lattice geometry.

2. 2.

   Layer Organization: Determine the minimum number of disentangling layers needed to ensure entanglement within the same length scale are considered. For the honeycomb lattice, three layers are sufficient, but other lattices may require different numbers.

3. 3.

   Overlapping Regions: Identify the junction points where blocks from different disentangling layers meet. These regions will serve as the centers for Wannierization of courier modes.

4. 4.

   Coarse-Graining Structure: Ensure that the coarse-grained lattice formed by the courier modes preserves the geometric structure of the original lattice, allowing for iterative application of the GMERA procedure.

The key insight of the 2d GMERA approach is the use of Wannierization to maintain the locality of courier modes in the overlapping region among different disentanling steps, preventing the uncontrolled spread of entanglement that would otherwise occur in direct applications of 1d techniques to 2d systems. By processing blocks in layers and carefully localizing the courier modes to overlapping regions, we can systematically disentangle the state while preserving its local entanglement structure.

We demonstrate the $2d$ GMERA procedure for a honeycomb Haldane model [6] for a single spin species of electrons, whose phase diagram encompasses multiple paradigmatic phases including both trivial and topological (Chern) gapped insulators trivial insulator, separated by quantum critical Dirac semimetal:

where $c_{i}^{\dagger}$ creates an electron on site $i$, and the remaining parameters are detailed below.

The nearest neighbor hopping, $t$, alone, gives a gapless Dirac semimetal (DSM) with two Dirac cones. The DSM is a quantum critical point with long-range entanglement, and infinite correlation length corresponding to scale-invariant correlations that decay algebraically with distance. These critical correlations and entanglement are expected to be captured by a (G)MERA wavefunction that is asymptotically self-similar.

Turning on a staggered potential, $V_{A}$, with opposite sign on the $A$ ($B$) sublattice, gives each Dirac cone opposite sign masses. This results in a topologically trivial insulator, with only short-range entanglement, that should be accurately described by a $2d$ PEPS wavefunction. Alternatively, we anticipate that this short-range entangled state can be accurately approximated by terminating the GMERA recursion after a finite number of RG steps.

When $V_{A}=0$, the (imaginary) second neighbor hopping, $t_{H}$, whose sign $s_{ij}=\pm 1$ is assigned based on whether $ij$ goes along or against the direction of the arrows shown in Fig. 5(a), gives the Dirac cones the same sign mass resulting in a Chern insulator with Chern number $C={\rm sign}(t_{H})$. Despite being short-range correlated, the Chern insulator has long-range entanglement and cannot be prepared from the trivial insulator by a short depth circuit. Moreover, there is an obstacle to preparing a Chern insulator with a scale-invariant MERA [8], which produces either strictly vanishing or asymptotically infinite correlation length, in contrast to the Chern insulator’s finite but non-vanishing correlation length. Our $2d$ GMERA procedure bypasses this no-go constraint by breaking the self-similarity assumption. Though it is geometrically self-similar, preserving the structure of the honeycomb, the unitary operations are scale dependent (different for each RG step), being generated from the physical correlation matrix at each scale.

For general parameter values, Fig. 5(b) shows a schematic phase diagram. To isolate the physics of each representative phase, in the following numerical simulations, we focus our attention on the red and blue shaded line segments where at most one of $t_{H}$ or $V_{A}$ is non-zero.

To benchmark the performance of our GMERA algorithm, we analyze the error in approximating the correlation matrix as a function of block radius $r$ for different phases and system sizes. We define the (root-mean-square, RMS) error-per-site as:

where $C_{ij}$ and $C_{ij}^{\text{approx}}$ are the exact and GMERA-approximated correlation matrices respectively, and $N$ is the total number of sites. As shown in Fig. 5(ci,ii), the error decreases exponentially with the block radius $r$, as $\epsilon\sim e^{-\alpha r}$ where the decay rate $\alpha$ depends on the parameter values $t_{H},v_{A}$, but is independent of system size (as expected for short-range correlated states).

The effectiveness of GMERA in compressing states with different correlation lengths is evident in the systematic downward shift of these lines. For the trivial insulator phase (Fig. 5(ci)), increasing the on-site potential $V_{A}$ from $0$ to $0.4$ results in shorter correlation lengths and correspondingly lower error at fixed block radius. Similar behavior is observed in the Chern insulator (Fig. 5(cii)), though we remark that the correlation length exhibits a non-mononotonic dependence on $t_{H}$ for $t_{H}>0$ This behavior aligns with intuitive expectation that states with shorter correlation lengths require less quantum information to be carried by the courier modes between blocks, enabling more efficient compression at fixed block size.

To further understand the relationship between compression efficiency and correlation length, we analyze how the fitted slope $\alpha$ from the error scaling $\text{error-per-site}\approx\exp(-\alpha r)$ varies with the inverse correlation length $1/\xi$. As shown in Fig. 5(d), we plot $\alpha$ against $1/\xi$ for both trivial and Chern insulator phases, with the correlation length $\xi$ extracted from the exponential decay of the real-space correlation functions. For the trivial insulator phase (blue points, $V_{A}$ varying from $0.05$ to $0.4$), the magnitude of the slope $|\alpha|$ increases approximately linearly with $1/\xi$, indicating that states with shorter correlation lengths can be compressed more efficiently. A similar linear trend is observed in the Chern insulator phase though with a different slope (red points, $t_{H}$ varying from $0.06$ to $0.2$). Crucially, we note that the parameter $\alpha$ remains finite at the gapless/critical DSM point ($V_{A},t_{H}\rightarrow 0$), showing that the GMERA can efficiently capture scale-invariant patterns of long-range entanglement.

Finally, in all cases the decay rate, $\alpha$, shows little dependence on system size beyond $L\sim 24$. Though it is impossible to rule out very slow (e.g. logarithmic) drifts with $L$, with any finite size numerics, our data strongly suggests that $\alpha$ remains constant throughout the phase diagram in the thermodynamic limit ($L\rightarrow\infty$).

For short-range correlated states, it may not be necessary to continue the GMERA RG recursion indefinitely. Truncating the RG recursion of the GMERA procedure after a fixed number of RG steps (and distilling all remaining orbitals in the final step) generates a Gaussian PEPS (GPEPS). Examining the accuracy of this resulting GPEPS versus the number of RG steps before termination (Fig. 6) reveals fundamental distinctions between the entanglement structure in the trivial and topological phases, despite both phases having short-range correlations.

A key step in implementing the $2d$ GMERA is efficiently compiling the block disentanglers acting on $B$ modes into sequences of two-body fermionic rotations. The block disentangler is a unitary matrix $\tilde{v}^{(i)}=(v^{(i)}_{f}|\tilde{v}^{(i)}_{c}|v^{(i)}_{e})$ that separates the modes into filled, courier, and empty sectors. This $B\times B$ single-particle unitary can be decomposed into $O(B^{2})$ $2\times 2$ single-particle rotations by rotating $\tilde{v}^{(i)}$ into identity column by column - each column requires at most $B$ $2\times 2$ rotations, leading to $O(B^{2})$ rotations in total. Each $2\times 2$ rotation can then be converted into a two-body fermionic rotation of the form $e^{-i\theta_{ij}(c^{\dagger}_{i}c_{j}+\text{h.c.})}$ acting on a many-body fermionic Hilbert space.

In a quantum processor built from fermionic modes, these fermionic rotations could be directly implemented. Implementing the Gaussian fermion rotations of the form $e^{-i\theta_{ij}(c_{i}^{\dagger}c_{j}^{\vphantom{\dagger}}+h.c.)}$ on a conventional, qubit-based quantum processor requires encoding the fermionic bilinears into qubits using a fermion-to-qubit encoding (F2QE) scheme. F2QE adds additional circuit complexity, which depends both on the encoding scheme, dimensionality, and connectivity of the fermion interactions. Below, we briefly reviewing the standard $1d$ F2QE based on the Jordan Wigner (JW) transformation [20] and show that its overhead scales only with the block radius, $r$, and not the system size for for $1d$ GMPS and GMERA methods.

The $2d$ GMERA procedure can be implemented in a qubit-efficient manner by recasting it as a sequential circuit and leveraging mid-circuit measurement and reset (MCMR) techniques [16], following similar strategies developed for holographic simulation of $1d$ MERA circuit [9]. Rather than implementing the full $2d$ GMERA circuit in parallel, which would require $O(L^{2})$ qubits for an $L\times L$ system, we can implement it sequentially by preparing one row of blocks with radius $r$ at a time. This sequential approach dramatically reduces the required number of qubits while preserving the ability to measure all desired observables.

Using MCMR techniques, where physical qubits are measured and reset after each row operation, the total number of required qubits scales as $O(r^{2}L\log L)$. This includes $O(r^{2}L)$ physical qubits to encode one row of fermions with radius $r$ and O(log L) bond qubits to mediate correlations between rows. Details of the implementation of qubit reuse MCMR can be found in Appendix B.3.

Our numerical simulations place empirical upper bounds on the quantum resources required to prepare certain types of area-law states. Namely, for gapped phases of correlated electrons that are adiabatically connected to a non-interacting electron state the GMERA procedure can be augmented by a constant depth adiabatic evolution [3] to approximate the interacting ground-state. This class of short-range correlated states that can be reliably prepared in this manner include (correlated) trivial and Chern insulators, but not, for example, fractional Chern insulators. For $2d$ systems with short-range interactions, this adiabatic dressing can be implemented with a fermion-to-qubit encoding with constant overhead, though simulating long-range (e.g. Coulomb) interactions may add additional polynomial (in system size) overhead.

The results in Fig. 5 suggest that the error in approximating Gaussian area-law entangled states by $2d$ as a GMERA decays exponentially in the block radius, $r\sim\sqrt{B}$. Equivalently, to achieve a desired target error-per-site, $\epsilon$, one needs to consider $r\gtrsim\log 1/\epsilon$. The number of two-qubit gates, $G(B)$, required to implement a general Gaussian unitary within a block scales as: $G_{\rm block}\sim B^{2}r\sim r^{5}\sim\log^{5}1/\epsilon$ where the $B^{2}$ factor comes from the number of fermion rotations to perform the distillation step, and the factor of $r$ comes from F2QE. Since the GMERA circuit has logarithmic depth, the gate-count for the entire GMERA circuit in system size $L\times L$ is:

where the factor of $(L/r)^{2}$ comes from the number of blocks in each RG step. This scaling applies for the long-range entangled Chern insulator and Dirac semimetal phases. For a short-range entangled trivial insulator with correlation length $\xi$, one can truncate the MERA RG procedure after $O(\log\xi/B)$ RG steps, removing the $\log L$ factor from the required gate-count.

Finally, we note that with qubit reuse techniques [16], the number of qubits, $Q$, required to implement the GMERA circuit can be reduced from $L^{2}$ (without qubit reuse) to either $O(L\log L\log^{2}(1/\epsilon))$ (for Chern and DSM) or $O(L\log^{2}(1/\epsilon))$ for a trivial insulator (truncating GMERA to GPEPS reduce the a $\log L$ factor).

We have introduced a $2d$ GMERA procedure to compress a $2d$ Gaussian state, represented by its correlation matrix, into a logarithmic-depth quantum circuit. This approach extends previous work by Fishman and White for $1d$ GMPS and GMERA compression schemes, by introducing a crucial Wannierization procedure to limit the build-up of higher-dimensional correlations that cause a direct extension of $1d$ techniques to fail. The approach can also be considered as a local, circuitization of the recently-introduced continuous time zipper entanglement renormalization (ZER) method [5].

By numerically implementing this procedure for a Haldane model on the honeycomb, empirically reveal that the approximation errors in the $2d$ GMERA procedure decrease exponentially in the size of the local blocks used in the procedure.

Compared to other recent approaches like isometric Gaussian fermionic TNS (isoGfTNS)[26], the $2d$ GMERA offers distinct advantages in circuit depth for large systems. While uni-isoGfTNS requires $O(2L)$ circuit depth and alt-isoGfTNS demands $O(L^{2})$ depth, the $2d$ GMERA approach scales as $O((\log L)\log^{5}(1/\epsilon))$, providing an exponential advantage for large system sizes at fixed target error $\epsilon$. This favorable scaling makes $2d$ GMERA particularly suitable for quantum simulation of large-scale systems.

A number of targets remain for future inquiry. Examples include generalizing the approach to systems without translational invariance (for example to simulate disordered systems or edge states of topological insulators), or to metallic states with Fermi surfaces that violate the area-law entanglement (requiring either block size to scale at least logarithmically in system size, or to consider a more complicated branching MERA type network [27, 28]).

For near-term quantum devices, an important direction is optimizing the GMERA implementation for smaller-scale demonstrations that can be realized with current NISQ hardware. This includes developing architecture-specific compilation strategies - for example, exploiting the all-to-all connectivity of trapped-ion systems and neutral atom arrays to minimize circuit depth and gate count. Other key considerations include exploring error mitigation techniques suited to the GMERA structure, and identifying the minimal system sizes needed to demonstrate key qualitative features of different phases. For example, one could focus on distinguishing between trivial and topological phases using small honeycomb patches that are accessible with current quantum processors, while carefully accounting for hardware connectivity constraints and gate error rates.

Since Gaussian states can be efficiently simulated classically, the main practical utility of an efficient Gaussian fermion state preparation procedure is to serve as an input to other variational or non-variational algorithms for including correlation effects. Highly entangled states in $2d$ or $3d$ is a ripe target for achieving practical quantum advantage. For gapped systems, correlations can, in-principle be built in by adiabatic dressing as described in the previous section. For gapless systems, whose gap vanishes as a polynomial in system size, other approaches may be required. Previous investigations in $1d$ systems indicate that using a Gaussian mean-field/Hartree-Fock as a starting point for variational optimization can dramatically reduce the number of variational parameters required to accurately capture correlated ground-states [2]. It would be interesting to assess whether the same type of variational improvements can be achieved in $2d$ (a challenging task for classical simulations that may require a quantum processor).