Variational Quantum Simulation of Valence-Bond Solids

We introduce a hybrid quantum-classical variational algorithm to simulate ground-state phase diagrams of frustrated quantum spin models in the thermodynamic limit. The method is based on a cluster-Gutzwiller ansatz where the wave function of the cluster is provided by a parameterized quantum circuit whose key ingredient is a two-qubit real XY gate allowing to efficiently generate valence-bonds on nearest-neighbor qubits. Additional tunable single-qubit Z- and two-qubit ZZ-rotation gates allow the description of magnetically ordered and paramagnetic phases while restricting the variational optimization to the U(1) subspace. We benchmark the method against the J1-J2 Heisenberg model on the square lattice and uncover its phase diagram, which hosts long-range ordered Neel and columnar anti-ferromagnetic phases, as well as an intermediate valence-bond solid phase characterized by a periodic pattern of 2x2 strongly-correlated plaquettes. Our results show that the convergence of the algorithm is guided by the onset of long-range order, opening a promising route to synthetically realize frustrated quantum magnets and their quantum phase transition to paramagnetic valence-bond solids with currently developed superconducting circuit devices.

We introduce a hybrid quantum-classical variational algorithm to simulate groundstate phase diagrams of frustrated quantum spin models in the thermodynamic limit. The method is based on a cluster-Gutzwiller ansatz where the wave function of the cluster is provided by a parameterized quantum circuit whose key ingredient is a two-qubit real XY gate allowing to efficiently generate valence-bonds on nearest-neighbor qubits. Additional tunable single-qubit Z-and two-qubit ZZrotation gates allow the description of magnetically ordered and paramagnetic phases while restricting the variational optimization to the U(1) subspace. We benchmark the method against the J 1 -J 2 Heisenberg model on the square lattice and uncover its phase diagram, which hosts longrange ordered Néel and columnar antiferromagnetic phases, as well as an intermediate valence-bond solid phase characterized by a periodic pattern of 2×2 strongly-correlated plaquettes.
Our results show that the convergence of the algorithm is guided by the onset of long-range order, opening a promising route to synthetically realize frustrated quantum magnets and their quantum phase transition to paramagnetic valence-bond solids with currently developed superconducting circuit devices.
Hybrid quantum-classical variational algorithms, so-called variational quantum algorithms (VQA), are at the center of current research for their potentialities in providing useful applications of currently developed noisy intermediate scale quantum (NISQ) devices [1]. They consist in a generic feedback loop where the NISQ device provides a quantum state via a parameter-Daniel Huerga: daniel.huerga@ubc.ca ized quantum circuit (PQC) that is tuned by a classical computer so as to optimize a certain objective function encoding the problem of interest [2][3][4]. The variational quantum eigensolver [5] was one of the first VQA proposed to approximate the ground-state and energy of finite stronglycorrelated fermionic Hamiltonians as an alternative to the phase estimation algorithm [6], which provides the exact ground-state solution but requires coherence times on the quantum devices unreachable with current technology.
The initially expected unleashed potentialities of VQA towards providing quantum advantage on problems including machine learning, optimization, and the simulation of strongly-correlated electron systems -a foundational motivation driving research in quantum computation [7][8][9][10][11][12][13]-have been narrowed down due to the identification of various limitations. Specifically, the optimization landscape has been shown to be plagued with the so-called barren plateaux [14], large flat regions hindering optimization that may appear irrespective of the optimization routine used [15] and induced by noise [16]. Additionally, suboptimal minima have been argued to render the classical optimization problem NP-hard [17]. In spite of these limitations, local objective functions, e.g. the energy of a local Hamiltonian, may still be efficiently optimized for shallow inexpressive PQCs in certain regimes [16,18,19].
In particular, two-dimensional (2D) frustrated quantum magnets, characterized by the impossibility of finding a spin arrangement satisfying all local energy constraints simultaneously in any locally-rotated basis [20], provide a natural testbed arena for VQAs and NISQ devices. They pose a challenge to state-of-the-art classical numerical methods [21][22][23] and at the same time host a plethora of phases and phenomena of both fundamental and applied interest. In particular, quantum paramagnetic phases either breaking or preserving translational symmetry -so-called valence-bond solids (VBS) and quantum spin-liquids (QSL), respectively-, have important implications for layered materials [24][25][26] and quantum computation [27][28][29][30][31]. Recent developed hybrid approaches have mainly focused on the simulation of 1D lattices [32][33][34][35][36], while approaches to 2D have been more scarce and limited to finite systems [36,37].
Here, we introduce a cluster-Gutzwiller VQA to simulate ground-state phase diagrams of frustrated 2D quantum spin models in the thermodynamic limit. We build upon the grounds of hierarchical mean-field theory (HMFT) [38], an algebraic framework based on the use of clusters for which a Gutzwiller ansatz represents the lowest order approximation. Furnished with a scaling analysis, it allows to uncover ground-state phase diagrams in the thermodynamic limit characterized by co-existence and competition of different long-range orders (LROs) and quantum paramagnetic phases [39][40][41][42][43][44]. Aiming at overcoming the scaling limitations of HMFT, here we present its quantum-assisted approach, dubbed Q-HMFT, where the cluster wave function is generated via a PQC whose central element is a parameterized real two-qubit XY gate that efficiently generates valence-bonds on nearest-neighbor (NN) qubits. Respecting the planar square connectivities of currently developed chips [45,46], we provide a systematic PQC construction resulting in shallower depths than other commonly used PQC ansatze, which either admix symmetry sectors [47,48] or require higher chip connectivities [11]. Information about the thermodynamic limit is implemented at the objective function level (i.e. the energy) through the self-consistent mean-field embedding concomitant to the cluster-Gutzwiller ansatz. Such mean-field embedding push the convergence of Q-HMFT within LRO phases, while paramagnetic phases are accessed by smoothly tuning the system through a quantum phase transition. We benchmark numerically Q-HMFT on the paradigmatic antiferromagnetic J 1 -J 2 model on the square lattice [49,50], which hosts Néel and columnar anti-ferromagnet (CAF) LROs, as well as a long-dabeted intermediate paramagnetic phase [41,[51][52][53][54][55][56][57][58][59][60][61][62][63][64]. Consistently with previous classical results from HMFT [41] and other stateof-the-art algorithms [53,54,58,61], we find that the Néel melts onto a VBS characterized by a 2×2 plaquette ordering. Convergence results from up to a 4×4-qubit cluster suggest the potential scalability of Q-HMFT, making it a promising route to quantum simulate frustrated magnets and valence-bond solid phases with currently developed superconducting qubit devices beyond classical capabilities.
1 Quantum-assisted hierarchical meanfield Let us consider a generic translational invariant 2-body local Hamiltonian in the infinite lattice, refers to the S=1/2 SU(2) spin operators at site i, and J i,j to the two-spin interaction strength. We tile the lattice with N -site equivalent clusters preserving as much as possible the original symmetries of the lattice. We classify the Hamiltonian terms on those acting within clusters, h R = (i,j)∈R h i,j , and those acting on two different clusters, h I R,R = i∈R,j∈R h i,j , where R labels the position of the cluster in the tiled lattice, or superlattice. We define a uniform cluster-Gutzwiller ansatz as the product state of N -site clusters in the infinite lattice, where we have restricted all clusters R to be in a same normalized state, |ψ(θ) , parameterized by θ={θ k }. Considering the wave function (1), the energy per spin reduces to the contribution of a single cluster and its mean-field embedding [41][42][43], where · θ refers to the expectation value with the parameterized cluster wave function, |ψ(θ) . The first summand in (2) accounts for the intra-cluster contributions, while the second runs over all twospin inter-cluster interactions with the neighboring clusters, (the 1/2 preventing double counting) and provides information about the thermodynamic limit by allowing the explicit breakdown of symmetries and concomitant onset of LROs. Generalization to n-body interactions appearing in e.g. ring-exchange models is straightforward [39,43,65]. The optimal energy E=E(θ ) with The wave function of the cluster is provided by a PQC whose parameters {θ} are tuned so as to optimize the energy density in the thermodynamic limit (2), which is an upper bound to the exact, E 0 . parameters θ is obtained upon minimization of the variational energy (2). Generically, the exact limit is recovered at N →∞. Thus, by increasing the cluster size, we assess the stability of the phase diagram. However, exact diagonalization of the clusters is typically limited to N 30 with classical computing means, which conforms the major bottleneck of the method [40-42, 44, 65].
In order to overcome such limitation, here we consider a PQC generating the cluster wave function, where |ψ 0 is an initial easy-to-prepare state, here onwards fixed to |ψ 0 = |0101 . . . 01 for reasons that will be clear below, with |0 . = |↑ , |1 . = |↓ , and U(θ) a digitalized unitary transformation, an hermitian operator acting on either one or two qubits (see Fig.1).
For SU(2) Hamiltonians on a bipartite lattice, we may restrict the variational search to the U(1) subspace of null total magnetization, j∈ S z j =0. Within this subspace, the selfconsistent mean-field embedding is restricted to the z-magnetization, i.e. S x j = S y j =0. This restriction allows to reduce the variational search and the PQC depth while still allowing for the description of symmetry breaking (long-range ordered magnetic) solutions, as well as SU (2) preserving (paramagnetic) ones. Within the twoqubit S z =0 subspace, spanned by |0 . =|10 and |1 . =|01 , any transformation can be generated by the XY-Heisenberg, In particular, a rotation generated by the DM interaction, V ij = S y ij , leads to a purely real XY-gate, which efficiently transforms an uncorrelated Néel into a singlet, i.e. a valence-bond, (4) can be considered as a generalization of a Givens rotations [66] to Hilbert spaces, and its use has been proposed for the digital preparation of fermionic Slater determinants [12,13]. Recently, a controlled version of it has been argued to be universal [67]. From the experimental standpoint, recent implementations have shown fidelities ∼ 97% [68].
Considering the planar square connectivities of currently developed superconducting quantum chips [45,46], we use L×L qubit clusters with even L. We heuristically construct the PQC (3) as it follows. First, we apply a set of REAL-XY gates (4) to all pairs of NN qubits in an order that (i) packs as much as possible twoqubit gates per layer and (ii) favors the C 4 symmetry of the square lattice. Second, in order to tune the sign-structure of the resulting wave function and increase the amount of correlations, we add two-qubit parameterized ZZ- in the same order described before. Finally, we add a layer of single-qubit parameterized Z-rotation gates, U Z j (θ)=e −iθZ j . Both ZZ-and Z-rotations are realizable experimentally [69]. Such a XY-ZZ-Z macro-layer, which minimizes the circuit depth within the considered gateset, is recursively applied a few times m to improve accuracy. Specifically, as the product of g-type two-qubit gates applied to NN dimers in a columnar pattern along α={x, y}, where r k =(x k , y k ) labels the position of the qubit considering the bottom-left corner as the origin, andê α is the unit vector. That is, 0 ≤ x k ≤ L − 2 (for even x k ) and 0 ≤ y k ≤ L − 1 for U g x , and vice versa for U g y . The operator T α U g α refers to the translation of U g α byê α , and U Z = j U Z r j (θ m j ) is the last single-qubit layer. For L=2, no translation operators are needed to uncover all NNs, thus d=5m (see Fig.1). For L≥4, the PQC (5) has n=(5L 2 − 4L)m variational parameters and a depth d=9m (see Supplementary Material). As shown in the following, a small m=2 is sufficient to exactly reproduce classical HMFT with 2×2 clusters, and provides with a good quantitative approximation to HMFT-4×4 results.
We perform ideal noiseless classical simulations of Q-HMFT on the J 1 -J 2 model (6) with 2×2 and 4×4 clusters, and m=2 and 4 macro-layers. We compute energy, magnetization and dimer observables, as well as provide details on its convergence, and compare with the classical HMFT results with the same clusters, which provide with a satisfactory extrapolation to the N →∞ limit, E ∞ (J 2 = 0) −0.64 [41], when compared to that obtained with state-of-the-art computational techniques, E ∞ (J 2 = 0) −0.67 [56,71]. As a classical optimizer, we use the gradientbased L-BFGS-B algorithm [72,73], where the gradient of the energy is computed via two-point first order approximation, ∂ θ k E = (E θ k +δ − E θ k )/|δ|, for δ = 10 −10 along k direction in variational space. We uncover the phase diagram by first obtaining an optimal solution out of various random PQC initializations at each extreme (i.e. J 2 =0 and J 2 =1), and then smoothly tuning the Hamiltonian parameters on both directions. The optimal parameters of one value of J 2 are used as initial parameters in the optimization of the immediate following point in the phase diagram. We compute the energy (2) and its derivatives to detect quantum phase transitions, as well as magnetic and dimer orders to characterize the phases.
In Fig. 2, we show the optimal energy per spin and first derivative with respect to J 2 as computed with Q-HMFT with m=2 and m=4 macrolayers, together with HMFT results. For 2×2 clusters, Q-HMFT with m=2 and fixing equal parameters within the XY and ZZ layers, respectively (i.e. with 12 variational parameters), reproduces exactly classical HMFT results and describes a Néel phase and its melting through second order quantum phase transition onto a plaquette-VBS at J c 1 To characterize the phases, first we inspect the magnetization, where r i = (r x i , r y i ) refers to the position of site i in the infinite lattice and k to a vector in the first Brillouin zone. Néel and CAF LROs are signaled by finite M z (π, π) and M z (0, π) (equivalently, M z (π, 0)), respectively. In Fig. 3 (top) we show Q-HMFT results on the magnetization together with HMFT-4×4 for comparison. For 2×2 clusters, Q-HMFT (m=2) reproduces exactly HMFT results (not shown for clarity purposes) showing the continuous vanishing of the Néel order parameter, an intermediate paramagnetic region, and a following discontinuous onset of CAF order, consistent with the description of second and first order phase transitions, respectively. For 4×4 clusters, Q-HMFT tends towards HMFT results upon increasing m, essentially providing the same phase diagram as with (Q-)HMFT-2×2. In particular, the vanishing of the Néel order parameter is discontinuous with m=2 (consistent with the first order transition seen in the energy), but such discontinuity decreases upon increasing m=4, suggesting a continuous vanishing for m>4 as found with HMFT. (π, π), 4×4 (0, π), 4×4 (π, π), Q-2×2 (2) (0, π), Q-2×2 (2) (π, π), Q-4×4 (2) (0, π), Q-4×4 (2) (π, π), Q-4×4 (4) (0, π), Q-4×4 (4) In order to characterize the possible onset of VBS order in the intermediate paramagnetic phase, we define the dimer observable along α = x, y allowing to distinguish between dimer or plaquette order [53], where i, j α refers to nearest-neighbor bonds along α. Specifically, D x =D y signals the onset of dimer-VBS along one direction, while D x =D y provides a signature of a C 4 preserving plaquette-VBS. In Fig. 3 (bottom) we show dimerization for Q-HMFT and its classical counterpart. For 2×2, Q-HMFT reproduces exactly HMFT results (not shown for clarity purposes) showing a nonzero D x =D y within the Néel phase that increases in absolute value and reaches a plateaux indicating the plaquette-VBS at J c 1 2 . Upon further tuning J 2 , both dimer observables show a discontinuity at J c 2 2 consistent with the first order transition found upon inspecting the energy. For 4×4 clusters, HMFT essentially provides the same picture, but with a reduced value of the dimerization throughout the phase diagram. Q-HMFT approximates quantitatively well HMFT results within the LRO phases, but shows a slight breakdown of the C 4 symmetry, i.e. D y <D x , within the intermediate plaquette-VBS. Nevertheless, such difference is reduced upon increasing m, tending towards the HMFT-4×4 results.
Last, we analyze the scalability of Q-HMFT by inspecting the variance with respect to the first variational parameter, Var[∂ θ 1 1 E θ ], estimated over 100 random initializations of the PQC, and the convergence of Q-HMFT-4×4 with m=2 macrolayers. As shown in Fig.4, the variance is strongly suppressed upon increasing the cluster size from 2×2 to 4×4 well within all three phases (i.e. at J 2 =0, 0.5, and 1), . Nevertheless, within the LRO phases, i.e. at J 2 =0 (Néel) and J 2 =1 (CAF), the convergence is pushed by the nonzero mean-fields concomitant to the onset of LRO. The system may explore suboptimal minima, but in most of the cases some extra tens of iterations permit the system to escape from them and approximate the classical HMFT result. On the contrary, in the quantum paramagnetic phase (J 2 =0.5), which is characterized by a null mean-field embedding both in (Q-)HMFT for the J 1 -J 2 model, the system explores a manifold of states characterized by different dimer orderings and magnetizations without arriving to converge (see Supplementary Material).

Conclusion and outlook
We have presented an hybrid algorithm to quantum simulate 2D quantum frustrated magnets in infinite lattices. Based on the cluster-Gutzwiller ansatz, we present a quantum-assisted approach to HMFT [41][42][43][44], dubbed Q-HMFT, where the wave function of a finite cluster is provided by a U(1) symmetry preserving PQC that respects the native NN square connectivities of currently developed superconducting quantum circuits, while the thermodynamic limit is accounted for by a self-consistent mean-field embedding introduced at the objective function level, i.e. the energy. The main ingredient of the PQC is a parameterized REAL-XY gate performing generalized Givens rotations in the two-qubit odd-parity subspace and efficiently generating valence-bonds. Addition of parameterized ZZ-and Z-rotations furnish a fixed-depth macro-layer for general even L×L-qubit clusters that is recursively applied a few m times to increase accuracy.
We provide noiseless benchmark numerical results on the paradigmatic frustrated J 1 -J 2 Heisenberg model with L=2 and 4 clusters, and m=2 and 4 macro-layers, and show that m=2 is sufficient to approach quantitatively well previous classical HMFT results [41]. In particular, m=2 exactly reproduces HMFT-2×2, describing a phase diagram hosting Néel and columnar antiferromagnetic (CAF) phases, as well as an intermediate quantum paramagnetic plaquette-VBS phase, in accordance with results from state-ofthe-art classical numerical techniques [53,54,61]. With 4×4 clusters, Q-HMFT provides an excellent approximation to HMFT within LRO phases, while the plaquette-VBS phase shows a slight dimerization -and thus breakdown of C 4that is decreased when increasing the number of macro-layers to m=4. We conjecture that HMFT results with generic L≥4 should be recovered for m≥6, for which numerical simulations require the use of more sophisticated optimization techniques [74][75][76].
Interestingly, although the variance of the gradient is strongly suppressed with system size, what is commonly considered an indicator of the emergence of barren-plateaux [14], the onset of magnetization (through nonzero mean-field embedding) pushes the convergence of the algorithm within Néel and CAF phases when randomly initialized, suggesting an eventual avoidance of barren-plateaux. On the contrary, the quantum paramagnetic VBS phase (characterized by null mean-fields) is accessed by smoothly tuning the system across the Néel-VBS quantum phase transition, as randomly initializing the PQC leads the algorithm to get lost exploring a manifold of states with different orders. These results highlight the key role of the mean-field embedding on the convergence of the algorithm, and suggest the potential scalability of Q-HMFT to large cluster sizes inaccessible by classical HMFT, i.e. even L>4, for which numerical simulation would require highly involved classical algorithms [77][78][79][80]. Further investigation on the efect of noise on the emergence of barren-plateaux [16] would be required to ultimately determine the scalability of the method.
Q-HMFT may be implemented with other technologies, in particular those hosting native XY gates [33] capable to realize different 2D lattices [81]. Moreover, the simulation of 2D Hamiltonians hosting exact VBS ground-states [65,82] via Q-HMFT offers a means for benchmarking quantum devices [83]. In addition, the symmetryguided construction of the PQC makes it suitable for developing error mitigation strategies [84], or the description of low-energy excitations over the ground-state [42,85]. From the experimental standpoint, we expect these results to motivate further development and refinement of the family of parameterized XY gates [86][87][88] as key elements for variational quantum algorithms.
A 2×2-and 4×4-cluster-Gutzwiller Given a tiling of the lattice and the cluster-Gutzwiller ansatz, the energy per spin of a translational invariant Hamiltoinan reduces to the contribution of the terms acting within the cluster and those connecting different clusters. The variational parameters determining the wave function of the cluster are determined through the Rayleigh-Schrödinger variational principle. When considering a uniform ansatz (i.e. translational invariant in the coarsed lattice), such energy minimization reduces to minimizing the energy of a single cluster in a self-consistent mean-field embedding (see e.g. [43] of main text). For the specific case of the 2×2 cluster used in this work (Fig. 5), the energy per spin of the J 1 -J 2 Hamiltonian is where all expectation values are taken with respect to the wave function of the cluster |ψ(θ) . Equivalently, the energy per spin obtained with the 4×4 cluster-Guzwiller ansatz, If the parameterization of the wave function is lineal, i.e. |ψ(θ) = θ {S z j } |{S z j } with θ ∈ C uncovers the whole Hilbert space of the cluster. Equating the derivative of the energy to zero, ∂ θ k E(θ) = 0, leads to a non-lineal set of equations that can be cast in matrix form, and solved by iteratively performing exact diagonalization of the L×L cluster with open boundary conditions and a set of self-consistent mean-fields, ψ(θ)| S j |ψ(θ) , acting on its boundaries. LRO is generically identified by a non-zero mean-field embedding, i.e. S j =0, while quantum paramagnetic phases are generically identified by optimal solutions where the mean-field embeddings is null, i.e. S j =0.

C Convergence
Although the variance of the gradient of the energy is strongly suppressed well within the three main phases identified in the J 1 -J 2 Heisenberg AFM, i.e. at J 2 =0, 0.5, 1, the noiseless numerical simulations of Q-HMFT with m=2 show that the convergence of the Q-HMFT algorithm is pushed by the onset of long-range order (see Fig. 4 in main text), eventually escaping from noise-free barren-plateaux. On the contrary, within quantum paramagnetic phases, the system does not arrive to converge and explores states closely in energy but with different orders. Within the HMFT framework, quantum paramagnetic phases are generically identified by converged solutions with null mean-field embeddings (see e.g. Refs. [41,43,44] of main text). In general, in the classical cluster-Gutzwiller algorithm, the system is able to find optimal quantum paramagnetic solutions, even if the initial seed breaks the symmetries (in this case SU (2)). However, in Q-HMFT the system explores a manifold of states close in energy but with different orders, not being able to stabilize a quantum paramagnetic state with well-defined dimerization (see Fig. 7).