Two-Dimensional Z2 Lattice Gauge Theory on a Near-Term Quantum Simulator Variational Quantum Optimization, Confinement, and Topological Order

We propose an implementation of a two-dimensional Z 2 lattice gauge theory model on a shallow quantum circuit, involving a number of single-and two-qubit gates comparable to what can be achieved with present-day and near-future technologies. The ground-state preparation is numerically analyzed on a small lattice with a variational quantum algorithm, which requires a small number of parameters to reach high ﬁdelities and can be eﬃciently scaled up on larger systems. Despite the reduced size of the lattice we consider, a transition between conﬁned and deconﬁned regimes can be detected by measuring expectation values of Wilson loop operators or the topological entropy. Moreover, if periodic boundary conditions are implemented, the same optimal solution is transferable among all four diﬀerent topological sectors, without any need for further optimization on the variational parameters. Our work shows that variational quantum algorithms provide a useful technique to be added in the growing toolbox for digital simulations of lattice gauge theories.


I. INTRODUCTION
Platforms for quantum computation are undergoing steady development and, nowadays, several architectures for the implementation of quantum circuits up to a few tens of qubits are available.In most of these platforms, qubits are considerably noisy, with coherence times that enable the application of only a few layers of unitary gates before losing quantum correlations.Therefore, for a successful use of these quantum computers the task of efficiently initializing the system in a target state becomes of the uttermost importance.
Variational quantum algorithms (VQAs) [1,2] are a promising technique to exploit noisy quantum hardware with shallow circuits.Among such methods, the quantum approximate optimization algorithm [3][4][5] (QAOA) is a popular strategy for preparing the ground state of a target many-body Hamiltonian H targ .Starting from a product state, usually the ground state of a simple mixing Hamiltonian H mix , it applies a series of unitary evolution operators generated alternately by H targ and H mix .The corresponding evolution times are treated as variational parameters to be optimized via a classical minimization of the energy.QAOA has been studied extensively for solving classical optimization problems [3,5,6], ground state approximation of quantum spin systems [4,[7][8][9][10], and has been shown to be computationally universal [11,12].Moreover, it has been successfully implemented on a trapped-ions quantum simulator [13] and on superconducting quantum circuits [14,15].
When considering current quantum computers based on superconducting qubits, the typical platforms are constituted by qubits arranged in two-dimensional ar-rays [16,17].These systems allow us to manipulate the qubits via single and two-qubit gates and, in most cases, two-qubit gates are local, i.e. they can be applied on neighbouring qubits only.These platforms are usually initialized starting from simple product states, while the efficient creation of complex entangled states is a nontrivial challenge.Despite the limitations of the noisy hardware currently available, these architectures open the path to various applications.One of the most important is the implementation of error-correcting codes and, to this purpose, considerable attention has been devoted to the realization of surface codes [18], thus to the preparation via local gates of states displaying topological order.In this context, the first experimental realization of the ground state of a surface code has been successfully achieved [19] and it allowed for the study of the main topological properties of its anyonic excitations.Another strategic application is the simulation of interacting quantum many-body problems that are intrinsically hard to simulate with classical computers [20].
In this framework, lattice gauge theories (LGT) emerged as a paradigmatic research subject.They constitute the backbone of particle physics, and many of their important features display a non-perturbative nature, requiring therefore advanced numerical techniques to be studied.Furthermore, some of the simplest twodimensional (2D) LGTs share the same topological properties of surface codes which, indeed, can be seen as the extreme deconfined limit of systems with Z 2 gauge symmetry.In the last decade, the application of quantum technologies to LGT became a lively field of research [21][22][23][24][25][26], progressing both on the development of several technologies and algorithms to tackle the complexity of LGTs and on the study of LGTs themselves.
One of the tasks which can be addressed through quantum simulation of LGTs is the study of their static properties.To this purpose, a key step is an efficient initialization of their ground states, allowing the investigation of their phase diagrams at low temperatures.
In this work, we explore the possibility of studying the ground state of a 2D pure lattice gauge theory within the framework of quantum circuits and digital quantum simulations.Indeed, the recent developments of quantum simulations provide complementary approaches to other quantum many-body approximation techniques such as Tensor Networks [23,[27][28][29][30][31], which are challenging to implement in 2D with current technologies, in particular for what regards quantum dynamics.We shall focus on the 2D Z 2 LGT, which is known to display a confinement-deconfinement phase transition between a trivial (confined) phase and a topologically ordered (deconfined) phase, matching the topological features of surface codes.As we show, most of the interesting ground state properties linked to topological order, which are usually described in the thermodynamic limit, can be characterized even with small lattices.
Our approach to realize the ground states of a 2D LGT is based on single and two-qubit gates only, and it is amenable to experimental implementations with limited resources, in particular with a number of qubits and fanout already accessible to present-day hardware or nearfuture devices.We apply QAOA to prepare the ground state at arbitrary values of the coupling and we show that the algorithm reaches high fidelities within a very small number of variational parameters, corresponding to a shallow quantum circuit.To reliably find optimal or quasi-optimal minima, we employ a two-step local optimization procedure [32], which provides regular schedules that can be efficiently transferred to larger systems.
Targeting the ground state in the confined phase, where there is no long-range entanglement, can always be performed efficiently and our numerical simulations suggest that QAOA can be scaled up to larger sizes without increasing the circuit depth.Concerning the preparation of states in the deconfined phase, instead, it is known that topologically ordered states cannot be obtained exactly with circuits of fixed depth for growing system size.In particular, for the ground states of the toric code, the required circuit depth scales linearly with the system width [19,[33][34][35][36].This is also a general property of QAOA, where long-range correlation and perfect control on the system is attained only with an extensive number of layers [4,37,38].
As a consequence, when targeting states in different phases, we compare two strategies: either we apply directly the QAOA evolution on a trivial product state or we first build exactly the toric code state and then apply the variational circuit from that starting point.The two approaches offer optimal results for targeting states in the confined or deconfined phase, respectively.They also display similar performances for the small system sizes we consider, except for the required overhead of the second approach.
The ground states prepared with QAOA are then used to characterize the crossover from the confinedtopologically trivial phase to the deconfined-topologically ordered one.In particular, we focus on the behavior of Wilson loop operators and of the topological entropy.We also discuss the possibility of exploring the ground state degeneracy when the lattice has periodic boundary conditions.Remarkably, all these indicators of a topological phase transition display very small deviations from their expected behavior in the thermodynamic limit despite the reduced size of our lattices.The successful implementation of 2D Z 2 LGT and the correct description of this non-trivial crossover -with very limited resources in terms of qubit numbers and circuit depth -provides a proof of principle of the feasibility of quantum simulations of deconfined and topological phases of lattice gauge theories in general.
The rest of this article is organized as follows.In Sec.II we describe the Z 2 lattice gauge theory model, introducing the main properties that characterize its topological order.In Sec.III we describe the implementation of QAOA for this model and how to represent the relevant unitary operators in terms of single and twoqubit gates.The main numerical analysis is presented in Sec.IV, where we focus on the performance of QAOA, its scalability to larger systems, and the characterization of the topologically ordered phase.Finally, we summarize our main results and some future perspectives in Sec.V.

II. Z2 LATTICE GAUGE THEORY
In this work, we consider a pure Z 2 gauge theory model on a regular square lattice.The discretized gauge fields are represented by qubits on the links of the lattice.Using a lattice of size L × L, there are 2L 2 qubits if periodic boundary conditions are imposed.The Hamiltonian we use is the sum of two competing terms which represent "electric" and "magnetic" noncommuting contributions.Their structure comes from an analogy with the QED Hamiltonian, where both space and the gauge group U (1) are discretized: the real space becomes a lattice, and U (1) is discretized to Z n .Here we focus on the smallest discrete group n = 2, which is naturally encoded in terms of qubits.The electric contribution to the Hamiltonian is where the index l runs over all the links in the lattice and the Pauli matrices are denoted σα l , with α = x, y, z.This specific choice of H E is motivated by the QED analogy, since the electric field enters the Hamiltonian via E 2 and our term is positive definite.A spin in the eigenstate

σx
l |+ l = |+ l brings no contribution to the electric energy and corresponds to a vanishing electric field.The state |− l indicates instead the presence of a Z 2 electric excitation on the link l, with energy cost assigned by H E .The magnetic term reads where p labels the plaquettes of the lattice and the plaquette operator B p involves the product of the four spin variables σz around the four links p 1 , • • • p 4 of the p-th plaquette (see Fig. 1).In particular, B p = −1 represents a magnetic flux through the p-th plaquette, and the interaction term h H B assigns an energy 2h to each of these excitations.The electric Hamiltonian H E effectively provides a kinetic energy to the magnetic fluxes.The local gauge constraint is the analog of Gauss's law and it selects the physically relevant sector of the Hilbert space.
For each vertex v of the lattice, physical states must be left invariant by gauge transformations, thus satisfying where the star operator A v is the product of the spin operators σx on the four links connected to the vertex v, as represented in Fig. 1.
The Hamiltonian in Eq. ( 1) has two well known limits, for h → 0 and h → ∞.When only H E is present (h → 0), the electric ground state is a trivial product state with all spins aligned along the x direction |Ω E = l |+ l , which satisfies the local gauge constraints in Eq. ( 4) and corresponds to the absence of any electric field excitation.In the opposite limit h → ∞ only the magnetic term remains, the system behaves like a surface code [18,39,40] and displays topological order.In this case, the number of ground states depends on the boundary conditions.The ground states of H B are the simultaneous eigenstates of all plaquette and star operators with eigenvalue 1 and correspond to the absence of magnetic fluxes.
In the case of open and smooth boundaries [18], there is a single magnetic ground state which can be expressed as an equal amplitude linear superposition of all possible contractible electric flux loops Γ: Here N is a normalization factor and W Γ the Wilson loop operator associated to a closed path Γ, defined as the product of σz matrices on the links belonging to Γ: Since σz |± = |∓ , Wilson loops applied to the electric ground state |Ω E create closed lines of electric field excitations.The second form of |Ω B in Eq. ( 5) expresses it as the normalized product of the projectors on the eigenstate of each plaquette operator with eigenvalue 1.
The magnetic coupling h drives the system across a topological phase transition, occurring at h c , between the electric and the magnetic phases, which are distinguished by different behaviors of the expectation values of the Wilson loop operators.From the definitions of the limiting ground states |Ω E , |Ω B , it follows that in the two limits h = 0 and h → ∞ we have for all paths Γ At a finite value of h, the expectation value of Wilson loops on the ground state decreases exponentially in the size of Γ, with a leading contribution given by [41,42] where A Γ and P Γ are the area enclosed by the loop Γ and its perimeter, respectively, while χ and δ are two positive functions.For h < h c , the system is in a phase dominated by the electric term H E , χ(h) > 0 and W Γ decays with an "area law".This means that large loops of electric excitations are strongly suppressed, which is a signature of confinement [43,44].In the opposite deconfined phase, where h > h c and the dominant term is H B , χ(h) → 0, and the behavior of large Wilson loops follows a "perimeter law".When periodic boundary conditions (PBC) are considered in both directions, the Hamiltonian acquires an extra Z 2 × Z 2 symmetry related to non-contractible 't Hooft loops.Consider a closed path C in the dual lattice.C crosses orthogonally a sequence of links of the direct lattice that we denote schematically by ∩ C. The 't Hooft loop operator commutes with the Hamiltonian (1) for any closed loop C.However, if C is contractible, τ C can always be expressed as a product of star operators A v , so that [τ C , Ĥ] = 0 does not provide any additional information that is not already contained in the gauge-invariance of the Hamiltonian.Considering PBC, the lattice becomes a torus and there are indeed two inequivalent non-contractible loops, whose corresponding 't Hooft operators τ h and τ v provide new symmetries.Fig. 2 shows examples of 't Hooft loop operators.
In the trivial limit h → 0, the expectation values of τ h and τ v on the electric ground state |Ω E are +1.Switching on strings of electric excitations, one sees that τ h and τ v measure the winding number of the electric field lines in the vertical and horizontal directions, respectively.With PBC, we can construct two non-contractible Wilson loop operators W h and W v , which commute with H B and A v , but vary the electric field winding numbers since they respectively anticommute with τ v and τ h : In the topological limit h → ∞, when [W h , Ĥ] = [W v , Ĥ] = 0, we get four degenerate ground states characterized by different eigenvalues of the 't Hooft loops, corresponding in the basis |τ h , τ v to: At finite values of h > h c , the perfect degeneracy between these four states gets lifted by an energy splitting vanishing exponentially with L.
The phase transition occurs at h c = 3.04438 (2) [45] and it can be understood by considering the duality between the Hamiltonian in Eq. ( 1) and the 2-dimensional quantum Ising model with transverse field (2D-TFIM), whenever the lattice has PBC, the gauge-symmetry A v |ψ phys = |ψ phys is imposed at each vertex, and the Hilbert space is restricted to the τ h = τ v = +1 sector [46].Indeed, we can define new Pauli spin variables X p and Z p on the dual lattice, where p denotes the plaquette centers, by identifying where p and p are neighboring plaquettes and l(p, p ) is the link shared by p, p .With the mapping in Eq. ( 10), the Hamiltonian becomes a transverse-field Ising model on the dual lattice: One can check that the algebra generated by the new operators is the same as the original one, confirming the unitary equivalence of the two models.Notice that the number of degrees of freedom is now halved: L 2 qubits (one for each plaquette), instead of 2L 2 (one for each link).This is an effect of the gauge symmetries, which are now automatically incorporated in the model.Finally, it is important to mention the fact that this duality fixes the global Z 2 symmetry of the Ising model: in the original representation, the product of all plaquette operators is the identity p B p = 1.In the 2D-TFIM, this is reflected in the condition p σx p |ψ = |ψ , which means that the physical states must be invariant under a global spin flip.
In this work, we shall employ both the original formulation of Eq. ( 1) and its dual model in Eq. (11).The dual Ising model will be exploited to speed up our numerical analysis of the Z 2 LGT and, in particular, to verify the scalability of QAOA between different system sizes.

III. QAOA AND CIRCUIT IMPLEMENTATION
A. Ground state preparation with QAOA To prepare the ground state of the LGT Hamiltonian in Eq. (1), we use the Quantum Approximate Optimization Algorithm (QAOA) [3].Although QAOA was initially proposed as a tool to look for approximate solutions to classical combinatorial optimization problems, it can be easily generalized to construct the ground state of manybody quantum Hamiltonians.Considering the two terms H B and H E in the LGT Hamiltonian, QAOA consists in constructing the following variational Ansatz where β = β 1 , . . ., β P and γ = γ 1 , . . ., γ P are 2P free real parameters, and the unitary operators U (γ m , β m ), for m = 1 For a given choice of the coupling h, which identifies a target Hamiltonian H targ (h) = H E + h H B , an approximation of the associated ground state is found using a classical minimization of the variational energy in this 2P-dimensional energy landscape.
The optimal energy at the global minimum E P (γ * , β * ) is a monotonically decreasing function of P.However, determining exactly the global minimum is, in general, not a trivial task [5], since local optimization routines tend to get trapped into one of the many local minima of the 2P-dimensional search space.We will discuss below an effective strategy to search for optimal (or quasi-optimal) solutions by a two-step QAOA procedure that starts from a linear schedule for the parameters, in the spirit of a digitized Quantum Annealing [48].

B. Circuit implementation of the QAOA Ansatz
The QAOA variational wavefunction in Eq. ( 12) is obtained by applying P layers of local unitary operators, by alternating the time evolutions generated by plaquette and electric field interactions.In what follows, we are going to describe how to implement the operations involved in each layer of the variational circuit by using only single and two-qubit gates.Since we shall focus on a single application of the unitary operations e −iβm H E and e −iγm H B , the index m will be dropped from the parameters.
The electric term of Eq. ( 2) is a sum of single-qubit operators, therefore the evolution it generates can be realized as a product of single-qubit rotations around the x axis by the angle β, up to an irrelevant global phase.The computational basis we adopt hereafter is the σz eigenbasis.Therefore, we employ Hadamard gates to diagonalize H E , and we reproduce the electric evolution during a single QAOA step by simultaneously applying operators U p (β) = e iβ σz to all qubits, i.e. a global rotation of angle β around the z-axis [49].A schematic representation of the single-qubit gates required is sketched in Fig. 3(b).
The implementation of the time evolution associated with plaquette operators is less trivial [50][51][52], but it can be realized in a local way as a combination of single-and two-qubit gates.Fig. 3(a) shows that a single-plaquette unitary operator e iγBp is obtained by a suitable combination of CNOT gates and a single-qubit rotation U p (γ) applied to the fourth qubit of the plaquette.The fourth qubit is the target of all CNOTs and it is restored to its initial logical state by the last three CNOTs, such that U p (γ) successfully applies the phase γB p only.An alternative technique based on ancillary qubits is presented in Refs.[53][54][55].
For a lattice composed of several plaquettes, the circuits in Fig. 3(a) cannot be simultaneously run on all of them, since two neighboring plaquettes share a qubit: as shown in Fig. 3(c), the qubit 4 not only acts as the target for the left plaquette but also as one of the controls for the right plaquette.The time evolution of the plaquette operators, however, can still be efficiently parallelized.For the sake of simplicity we shall consider first systems with an even number of columns (for an even number of rows the situation is formally equivalent).In this case, we can decompose the whole lattice into sets of two neighboring horizontal plaquettes, each set with the same structure as depicted in Fig. 3(c).
We focus on a single set, which corresponds to our basic unit.We show in Fig. 4 the corresponding quantum circuit, that will be run in parallel for all such sets.Neighboring pairs of plaquettes share qubits at their boundary: this can be understood by ideally replicating the pair in Fig. 3(c), to build a lattice.For instance, the qubits 3 and 6 of our set also correspond to the qubits 2 and 5 of the set above the one in exam.Similarly, qubits 2 and 5 are homologous of qubits 3 and 6 for the set below, whereas qubit 7 matches qubit 1 of the plaquette pair lying on the right of the one depicted, and so on.
The algorithm defined in Fig. 4 performs the rotation of both plaquettes in 12 steps and it can be run in parallel for all plaquette pairs.This procedure is based on applying the phase gates U p (γ) on qubits 4 and 7 (thus on all the qubits of the lattice lying on the vertical links) to implement the plaquette rotation.The gates partially depicted in green are related to the simultaneous realization of the same algorithm for the neighboring plaquettes.All the qubits lying on the vertical links are required to be connected via CNOTs to the four neighboring qubits, whereas the horizontal qubits are just connected to two neighbors each.As a result of the previous scheme, each Figure 4. Algorithm to implement the plaquette operator e iγB on the two plaquettes depicted in Fig. 3(c).The labeling of the qubit lines emphasizes that all the boundary qubits are shared with the neighboring plaquette pairs.The partially depicted green gates (dashed lines) are related to the simultaneous implementation of the same algorithm on a couple of neighboring plaquettes.They connect the displayed qubits with qubits belonging to the surrounding lattice sites, based on suitable translations of the two-plaquette unit.
of the P steps of the QAOA can be realized with a circuit of depth 13 on systems with open boundaries, or systems with closed boundaries and an even number of rows or columns.In this work, we focus most of the numerical investigation on a lattice with 3 × 3 plaquettes and periodic boundary conditions.Its implementation on actual quantum hardware requires some additional care due to the boundary conditions and, in that case, each QAOA step can be realized with a circuit of depth 18 (see Appendix A).
To evaluate the scaling of the preparation of a target ground state of the Z 2 LGT, it is also important to consider the preparation of the initial states |ψ 0 , which will subsequently be modified by the QAOA layers.The electric ground state |Ω E is a trivial product state and it can be prepared by applying Hadamard gates on all qubits, to rotate them into the eigenstate |+ of the σx operator.
The toric code (magnetic) ground state |Ω B displays instead topological order and long-range entanglement [34].To initialize this state, we follow the technique adopted in Ref. [36]: assuming that each qubit is initially in the eigenstate | ↑ of σz , in each plaquette we first apply three Hadamard gates on the control qubits and then three CNOT gates targeting the fourth [56].This procedure is similar to the first 3-CNOT sequence of the circuit depicted in Fig. 3(a), with the addition of Hadamard gates on the qubits |p 1 , |p 2 , and |p 3 .
These operations can be performed in parallel on plaquettes belonging to a single row and then repeated L times to cover the whole lattice.Operations on different columns (or rows), however, cannot be parallelized because in each plaquette the CNOT gates must be applied before using one of the control qubits as the target for the neighboring column (row).If we consider Fig. 3(c), the plaquette operations must start from the rightmost column, in such a way that the qubit |p 4 is used as a control before becoming the target of the plaquette on the left.This is an important difference with respect to the application of the gate e −iγ H B : despite the preparation of the ground state of H B on a single plaquette requires a smaller number of gates, for large systems the initialization of |Ω B requires a deeper circuit than the Hamiltonian gate, which, instead, can be run in parallel on all the even or odd columns (rows).This reflects the necessity of having a circuit with depth O(L) to prepare a state with long range entanglement, such as |Ω B , which has been well studied in the literature [33][34][35].In conclusion, when we compare the QAOA results with different choices of the initial state |ψ 0 , we need to take into account the overhead required for preparing |Ω B .
The construction of the QAOA layers we presented so far was restricted to the case of an Abelian Z 2 LGT.We stress, however, that the same procedure can be generalized to pure 2D LGTs with arbitrary discrete gauge groups.In particular, by suitably extending the Hilbert space associated with each link of the square lattice, it is possible to implement the time evolution steps of both the electric and magnetic Hamiltonian based on local unitary operators [52].
In this respect, the simplest generalization is provided by Z n LGTs (see, for instance [54,[57][58][59][60][61][62][63]).In this case, a gauge degree of freedom is encoded into an n-dimensional Hilbert space, as common in quantum clock models with Z n symmetries [64,65].The electric field assumes indeed n different values, which may be represented by suitable qudits (or by embedding each degree of freedom in a set of qubits).H E remains a local Hamiltonian, whose time evolution can be perfomed in a parallel way over all links.
As in the case of the Z 2 theory, also for Z n symmetries there is a suitable unitary transformation mapping the eigenstates of the electric Hamiltonian into the eigenstates of the magnetic operators adopted to build the plaquette terms (the so-called connection operators).Such unitary transformations generalize the Hadamard gates we adopted and correspond to a quantum Fourier transform.Additionally, the plaquette term maintains the same 4-body interaction form through a suitable replacement of σz with quantum clock operators.The implementation of the plaquette operator thus requires generalizing the CNOT to controlled Z n clock gates.Again, the phase diagram of pure 2D Z n LGT models presents a deconfined and topological phase at large h, whose topological order matches the Z n generalization of the toric code [66,67], and a confined phase whose ground state becomes a trivial product for h = 0.
By following the Kogut-Susskind Hamiltonian construction, a further generalization can be implemented to investigate ground states of discrete non-Abelian 2D LGTs (see, for example, Refs.[52,68,69]).In this case, the gauge degrees of freedom can be represented either in an eigenbasis associated with the irreducible representations of the group, which is diagonal in the electric term of the Hamiltonian, or in an eigenbasis associated to the group elements, which is diagonal in the magnetic term of the Hamiltonian.The general structure of a quantum algorithm for implementing the QAOA steps in this case is analogous to the previous one and can be based on the construction in Ref. [52].Given the non-Abelian nature of the group, however, the implementation of the magnetic time evolution requires a further technical generalization.In this case the irreducible representations are not one-dimensional and correspondingly the connection operators acquire a tensor form, therefore the gaugeinvariant plaquette terms must be written in terms of their trace [43], requiring, in turn, to extend the rotation operators U p (γ) to more general single-link gates, which apply phases given by the traces of gauge group matrices.

IV. NUMERICAL RESULTS
In this section we analyze the QAOA performance on our LGT model, showing that the ground state can be prepared through shallow circuits with good fidelity both in the confined and in the topological phase.Unless otherwise stated, the numerical analysis is performed on a lattice with 3×3 plaquettes (18 qubits) and implemented through the python package Qiskit [70], using the circuit sketched in Sec.III B. Simulations of larger systems (L = 4, 5), instead, exploit the mapping onto the 2D-TFIM to reduce the Hilbert space dimension and allow for the exact evaluation of the QAOA Ansatz.

A. Energy landscape
The first feature we are interested in is the structure of the energy landscape, because it determines whether the classical optimization of the variational parameters can be performed efficiently or not [2].Indeed, the presence of rugged energy landscapes is a common problem that severely affects the classical optimization loop of VQAs by making it prone to remaining stuck in local minima, some of which might be far from the ground state energy.To characterize the energy landscape and to quantify the quality of the optimized variational Ansatz |ψ P (γ * , β * ) , we use either the residual energy or the fidelity.Given a target Hamiltonian H targ (h), we denote with |ψ targ its ground state and with E GS the corresponding energy, both obtained with exact diagonalization.The residual energy is simply the difference between the minimized variational energy E(γ * , β * ), defined in Eq. ( 14), and E GS , while the fidelity follows the usual definition In this work, we use the fidelity as a precise estimate of the accuracy of the approximation of the target state, even though, in an actual experiment, it is hardly accessible because it requires an exponential number of measurements.The energy, instead, is easily estimated: the magnetic contribution is diagonal in the computational basis, while the electric contribution is obtained by applying a basis rotation on each qubit, i.e. a set of Hadamard gates, before the measurement.
In the Z 2 LGT model, the energy landscape emerging from the QAOA Ansatz is characterized by many local minima covering a wide energy interval, making random-start local optimization impractical.This effect is particularly severe if the target state is in a phase different from the initial one |ψ 0 .This is illustrated in Fig. 5, where we show the residual energy and the infidelity 1 − F P (γ * , β * ) for 100 different random-start local optimizations with P = 5 QAOA layers and three different values of the magnetic coupling, around or above the topological phase transition.The initial state is the product state |Ω E in the extreme confined limit (h = 0) and the local minimizations are performed with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [71].Although there is a clear concentration of data in the corner corresponding to successful optimizations, where both the infidelity and the residual energy tend to zero, there are many local minima far from the ground state, suggesting that more refined optimization techniques are needed for this model.
However, the clear correlation between the energy and the fidelity is reassuring since the absence of low-energy minima with small projection on the ground state guarantees that any scheme that allows for a reliable minimization of the energy will also lead to a good approximation of the target state.This correlation is intuitively justified by the existence of a gap in the topological phase.Indeed, the construction in Eq. ( 12) cannot mix different topological sectors of the model and, for each topological sector, there is only one ground state.Away from the critical point the ground state is protected by a finite gap, such that the correlation between infidelity and residual energy must hold below this energy scale.Close to the critical point, other orthogonal low-energy states may appear and spoil the correlation.This, however, seems not to be the case.Its resilience is not surprising for small system sizes in which the gap does not close even at h c .However, we observe that the correlation between infidelity and residual energy holds also when we increase the linear dimension L of the lattice, and it actually appear to be even sharper, as shown in App.B 3. The variational energy is therefore a reliable figure of merit for the optimization, it can be efficiently measured in experiments, and the procedure is still effective when the system size is scaled up.All these represent positive factors for the feasibility of the implementation of QAOA on the Z 2 LGT model in near term quantum devices.
B. Heuristic local optimization: two-step QAOA Because of the large number of suboptimal minima present in the energy landscape, it is important to adopt an efficient strategy in order to reliably find a good approximation of the true ground state of the target Hamiltonian.This is, indeed, a crucial task for QAOA and VQAs in general, where the classical optimization outerloop is often the main computational bottleneck and several strategies have been proposed that go beyond a local search from a random start.These strategies range from problem-specific methods to general iterative techniques, based on observed patterns in the optimal schedules [5,[72][73][74].
We adopted here a simple two-step minimization protocol inspired by a digitized Quantum Annealing [4,48] turning-on of one of the two terms of the Hamiltonian.The idea behind the two-step optimization is to leverage on the formal analogy between QAOA and digitized Quantum Annealing (QA) [4]: for depth-P QAOA, we first optimize the total run time of a digitized linear QA [48] of the same depth, and then fine-tune the variational parameters around this schedule.This approach can be used effectively when the system is initialized either in the electric state |Ω E or in the magnetic one |Ω B .For generic applications, the confined electric ground state |Ω E -a uniform superposition of all possible states in the computational (magnetic) basis -is a standard choice for initializing the variational circuit, because it is easy to prepare.However, this choice is non-optimal when we target states with long-range entanglement in the deconfined/topological phase, since a circuit of local unitary gates with bounded depth cannot create states with topological order beyond a certain system size [34].
When the initial state is set to be |Ω E , an adiabatic turning-on in a time τ of the magnetic coupling through H(t) = H E + (t/τ )h H B suggests -after digitization by Trotter decomposition of e −i∆t H(tm) ≈ e −i∆t H E e −i m∆t P h H B , with t m /τ = m/P -setting γ 0 m = m ∆t P h and β 0 m = ∆t in the state in Eq. (12).In our first QAOA step, these linear-schedule parameters are opti-mized by searching for the optimal digitized QA [48] ∆t * -a simple one-dimensional minimization -which leads to setting: The second step in our QAOA procedure is to perform 10 local BFGS optimizations in the 2P-dimensional parameters space, starting from (γ dQA , β dQA ) + , where is a small 2P-dimensional vector with random numbers uniformly distributed in the interval [−0.025, 0.025), keeping the best outcome out of these local optimizations.Schematically: The toric code ground states |Ω B , corresponding to the extreme deconfined limit h → ∞, provide a better initial state |ψ 0 when targeting ground states in the topological phase: they can be exactly prepared with local circuits whose depth scales with the width of the system [19,36].Proceeding once again with an adiabatic turning-on, now of the electric part of the Hamiltonian, through H(t) = h( H B + (t/hτ ) H E ), suggests -after digitization by Trotter decomposition of e −i(∆t/h) H(tm) ≈ e −i∆t H B e −i m∆t hP H E , with t m /τ = m/P -setting γ 0 m = ∆t and β 0 m = m ∆t hP in the state in Eq. (12).Once again, these can be optimized by searching for the optimal digitized QA ∆t * , which leads to: The second step in our QAOA procedure is identical to the previous case, as schematically indicated in Eq. (17).Two noteworthy features of the QAOA minima obtained by applying our two-step QAOA procedure are the smoothness of the schedules (γ , β ), illustrated in Appendix B 2, and the closely related transferability of such smooth schedules from a smaller to a larger L > L sample, discussed in Sec.IV C. We benchmarked our heuristic two-step QAOA approach against a computationally expensive global optimization, finding comparable quality results in terms of ground-state fidelity, both in the confined and in the deconfined phase: we illustrate this in Appendix B 1.
In the following, we will compare the performance of our two-step QAOA for systems prepared either in the electric ground state |Ω E , or in the toric code ground state |Ω B = | + + .For a fair comparison, a remark is in order: while |Ω E is trivially prepared with one layer of single-qubit Hadamard gates, for the preparation of |Ω B one should include an overhead circuit with O(3L 2 ) gates, organized in L layers applied sequentially.As explained in Sec.III, although no optimization is necessary for this preliminary step, it is still required to apply 3 CNOT gates for each plaquette.The QAOA results obtained from the initial product state |Ω E are reported in Fig. 6, where we show the infidelity 1 − F P (γ, β) as a function of the circuit depth P for several values of the magnetic coupling, both in the confined phase (h 3) and in the deconfined one (h 3).As expected, the variational Ansatz converges faster to states in the same phase (e.g.h = 1, 2) but QAOA can reach very good fidelity 1 − F P < 10 −3 , when P ≥ 5, for all the couplings we considered.
With the "reversed" protocol, starting from the toric code ground state |Ω B = | + + , we obtain an overall behavior similar to what observed for |ψ 0 = |Ω E , see Fig. 7(a), with the important difference that now the optimization converges faster when targeting the deconfined phase.Indeed, only P = 3 QAOA layers are now needed to reach 1 − F P < 10 −3 when h > h c , see data for h = 4 or h = 5, while confined states require more QAOA layers to reach comparable accuracy.
For both choices of initial state, we observe that the infidelity decreases exponentially with the circuit depth; the only exceptions for P = 5, 6 can be ascribed to the algorithm remaining stuck in a (high-quality) local minimum, when the target state is very close to the initial one (see Appendix B 2).However, if we focus on the minimal resources to approximate the target state within a certain fidelity threshold, we can further reduce the number of parameters required.Figure 7(b) shows a comparison of QAOA performance with the two possible choices of the initial state, for P = 2 and P = 3, by looking at the best fidelity reached by the two-step optimization as a function of the coupling h.Remarkably, such shallow variational circuits are enough to prepare with high fidelity the ground states in the confined and deconfined phases, provided the initial state is selected in the same phase as the target ground state.Unsurprisingly, the region that requires a larger number of parameters, i.e. a Comparison between two-step QAOA performance by starting from the electric and the magnetic ground states: we plot the fidelity vs magnetic coupling h at fixed values of P.Here P only takes into account the number of QAOA layers with parameterized gates, while it does not include the computational overhead for the preparation of |ΩB , compared to preparing |ΩE .
deeper variational circuit, corresponds to the crossover between the two regimes, where 2 h 3.
We finally observe that the choice of the initial state based on the target value of h plays a role analogous to the choice of the electric or magnetic representation of the LGT Hamiltonians applied in the quantum simulation protocols presented in Refs.[75,76].

C. Schedule transferability
A promising route to reduce the computational cost of the outer-loop classical optimization in VQAs is the transferability of the optimal parameters from small to large instances of the same model.Indeed, as empirically observed or proven in specific applications of VQAs, if you consider two instances of the same model and a fixed variational circuit depth P, the optimal parameters obtained for the small system of size L may serve as a very good warm-start (or educated guess) for a local optimization for the L -size model (L > L) [77][78][79].
Classical numerical simulations soon become unfeasible even for modest sizes, often hindering a more systematic analysis on this issue: for our LGT model, which requires 2L 2 quantum spins, even sizes as small as L = 4, 5 can be extremely challenging to simulate exactly.To partially overcome the size limitation, we exploit the mapping onto the 2D-TFIM, explained in Sec.II, which involves only L 2 spins on a square lattice, by taking advantage of the restrictions imposed by the gauge constraints.This allows us to simulate exactly the variational optimization for L = 4, 5. To study the schedule transferability, we first perform the two-step QAOA on the system with L = 3, as described in Sec.IV B. The optimal angles (γ , β ) found for L = 3 are then used as warm-start points for a local optimization on larger sizes.In particular, we keep the best run out of 10 BFGS optimizations on the larger systems, each of them starting in the neighborhood of (γ , β ) L=3 , similarly to the strategy used in the second part of the two-step QAOA protocol.This procedure is repeated for different values of the coupling.
The results obtained are reported in Fig. 8, where we compare the fidelity F P (γ * , β * ) vs h, for circuit depth P = 6 -which allows us to prepare the ground state for arbitrary h with an error 1 − F P (γ * , β * ) < 10 −3 for L = 3 -and both possible initial states: |ψ 0 = |Ω E (full symbols) and |ψ 0 = |Ω B (empty symbols).The transferability of the parameters is almost perfect when the initial and target states are in the same topological phase, leading to very high fidelities both in the small and large magnetic coupling regimes.Even when we target a ground state in a different phase than the initial one -for instance, |ψ 0 = |Ω E and h > h c ∼ 3.0the final fidelity is still large, allowing us to characterize the topological properties of the final state.The poorer fidelity observed for L = 5, when targeting a state in the opposite phase than |ψ 0 , is due only to the reduced ratio between the circuit depth P and the size L (particularly relevant when crossing a topological phase transition [34]) and not to the size of the system on which the original schedule has been optimized (L = 3).Indeed, transferring the schedule from L = 4 to L = 5, or from L = 3 to L = 5, leads to the same final performance.
Interestingly, the warm-start initialization provided by the L = 3 optimal parameters leads to a successful local minimum search for L = 4, 5, with an accuracy close to what can be achieved with a full global minimization, discussed in Appendix B 1.Moreover, the number of iterations needed for the local optimization is rather small (N iter 50), confirming the benefit of the transferability of optimal solutions: once the L = 3 two-step solution is provided, only a small overhead in computation resources is required to fine-tune the parameters for larger sizes.Hence, transferability provides a speed advantage over starting from scratch a two-step optimization: even though the fidelity reached is comparable, the latter requires more runs of the quantum circuit, making it less efficient when the optimal schedule for a smaller system is already known.
This transferability evidence may be linked to the observation of the smooth schedules we found with the twostep optimization, as shown in Appendix B 2. It is also important to remark that the schedule transferability is not a general property of any minimum in the energy landscape, but it is associated with the smooth solution found with the two-step protocol.For instance, a global optimization yields slightly better results on the L = 3 system, but it often represents a poor choice as an educated guess to initialize a local minimum search on larger sizes, as discussed in Appendix B 1: this phenomenon is similar to overfitting in machine learning [80].In this respect, the two-step scheme appears to outperform an extensive global search.

D. Ground state characterization
In the following, we turn our attention to the properties of the approximate ground states we prepare with QAOA.Despite the finite size limitations of our simulations, the states obtained through QAOA display most of the main features associated with the appearance of topological order and the crossover from a confined to a deconfined phase as h increases.The main observables to distinguish these two regimes are the Wilson loops, as defined in Eq. ( 6).We consider in particular Wilson loops W lx,ly defined over rectangles of size l x × l y .
As explained in Sec.II, it is known that the deconfined phase is characterized by an exponential decay of W lx,ly with the perimeter P of the loop, whereas the confined phase displays a decay dictated by the area A of the loop [42].In particular, the magnetic ground states |Ω B are such that W lx,ly = 1, while in the electric ground state |Ω E Wilson loops always present vanishing expectation values.Recalling Eq. ( 7), the overall behavior of a Wilson loop can be approximated by W ∝ e −χA−δP .Indeed, if χ > 0, the exponential decay with the area dominates for large loops, while if instead χ = 0, the decay is dictated by the perimeter law only.To extract the information about the χ coefficient we estimate the so-called Creutz ratio [44]: This ratio is indeed built to cancel the perimeter contribution to the decay of the observables and approximate the coefficient χ, which is recovered for large l.Fig. 9 displays the Creutz ratio in a system with L = 5 and periodic boundary conditions for states obtained with P = 6 QAOA steps applied either to the state |Ω E (for h < 3) or to |Ω B (for h ≥ 3).The optimization on the L = 5 systems was initialized with the best result obtained with the two-step protocol for L = 3, on top of which we performed 20 local minimum searches, out of which we consider the most successful outcome.
Analogously to other LGT studies on small lattices [31,81,82], the finite size effects in our computation are strong.When considering a Wilson loop of width 3, its opposite sides lay at distance 2. This implies that what we observe in Fig. 9 may provide a quantitative estimate of the behavior in thermodynamic systems only if the correlation length is sufficiently smaller than this distance, thus only sufficiently far from the phase transition.
Despite this limitation, the Creutz ratio χ(3, 3) presents a behavior that clearly distinguishes the confined phase (χ > 0) and the deconfined phase (χ → 0) appearing for h 3, although a quantitative identification of h c is beyond the possibilities of these small systems and loops.
The inset of fig. 9 reports the expectation value of the Wilson loop operators corresponding to the Creutz ratios shown in the main plot.It clearly shows a crossover between the trivial, confined state with W Γ → 0 and the topologically ordered, deconfined limit W Γ → 1.With the chosen scheme, i.e. starting from |Ω E or |Ω B depending on the target state, they perfectly match the results from exact diagonalization (not shown) as expected from the high fidelity reached, see Fig. 8.We emphasize that the possibility of obtaining a reliable estimate of the expectation value of the Wilson loops yields further important implications: Ref. [76] shows indeed that, in a U (1) LGT, even the expectation value of the single plaquette operator can be used to extract the running coupling of the model, which is a fundamental quantity related to its renormalization.
If we chose to start always from the electric ground state, the deviation from exact results would become larger in the deconfined phase, as also expected from the fidelity drop observed in Fig. 8.However, the results obtained in this non-optimal case still provide an acceptable scaling of the Wilson loop for the deconfined regime (empty symbols in Fig. 9): even without perfect reconstruction of the target state, it is still possible to identify the deconfined phase.This is, indeed, useful for experimental investigation, where realistic setups are limited to shallow circuits and noise would decrease the quality of the approximated ground state.
Another observable that marks the onset of topological order is the topological entropy [83,84].Given a connected subsystem A ∪ B ∪ C of the whole lattice, we define its topological entropy as in the equation below Here S X is the von Neumann entanglement entropy of a generic subsystem X, obtained by tracing out all degrees of freedom in the complement of X with respect to the whole system, and {A, B, C} is a tripartition of the region of which we compute the topological entropy.In the toric code state |Ω B , the topological entropy of any subsystem is S topo = − ln 2 and the total entanglement entropy is where N v is the number of vertex operators A v cut by the edge of the bipartition X [19,85].In a product state, such as |Ω E , we expect both quantities to be zero, while for generic values of h the entropy should interpolate between the two limits.To compute the entropy, we choose a subsystem X with 6 qubits, as depicted in Fig. 10(a), and we divide it into three further regions A, B, and C with two qubits each.We compute the entanglement entropy of all the subsets used in Eq. ( 20) by tracing out explicitly their complements and obtain the data plotted in Fig. 10(b).Despite the small dimension of the lattice and its subsystem, our results agree perfectly with the theoretical prediction: in the deconfined phase, the total entanglement entropy is S ent = 4 ln 2, since the partition ABC cuts five vertices and the topological entropy approaches S topo = − ln 2. Finally, we would like to show that it is possible to manipulate the state constructed with QAOA to change its symmetry sector when the system has PBC.Let | + + P denote the approximate ground state constructed with a QAOA circuit of P layers.We then construct approximate candidate ground states in the other topologial sectors by applying non-contractible Wilson loops The subscript label P is here used to distinguish the states obtained via QAOA from the exact eigenstates of the Hamiltonian.)Noncontractible Wilson loops W v/h are immediately implemented via L single-qubit gates σz acting on a vertical or horizontal line.By doing so, however, we introduce an extra error on top of the finite accuracy of the QAOA state: Indeed, choosing a specific vertical or horizontal Wilson loop to change the symmetry sector of the system breaks the translational invariance of the constructed state, producing a small excitation.This effect is visible in Fig. 11, where we show the energies of the state approximated with QAOA, denoted by | + + P , and of the other three states obtained by applying W h and W v on | + + P .For comparison, we plot also the low eigenvalues obtained by exact diagonalization (drawn with solid blue lines).For large h, the four lowest energy levels should be almost degenerate and, indeed, the exact diagonalization results are almost indistinguishable for h ≥ 4. In the same region, the excess energy of the approximate states |τ h , τ v is instead clearly visible, although well below the topological gap with the first proper excited state.
An alternative procedure to explore the different topological sectors in the deconfined regime, is to apply first the relevant Wilson loop on |Ω B and then the QAOA unitaries.In such a way, the initial state W h/v |Ω B is exactly degenerate with |Ω B .We find that the optimal schedule (γ * , β * ) used to prepare the state | + + P minimizes also the expectation value of the energy in the other topological sectors, so no further optimization is required.However, the picture presented in Fig. 11 remains valid and small excitations are created in the other topological sectors.In other words, by inverting the order of application of the operators W h/v and U (γ * , β * ), we observe nearly irrelevant changes on the expectation value of the energy; U (γ * , β * ) is the QAOA evolution operator with optimal parameters for the state | + + P .
The expectation value of the 't Hooft loops τ h and τ v , which distinguish the different topological sectors, is perfectly reconstructed by the algorithm.This last feature is, however, independent from the specific values of h and P, since the QAOA evolution respect the global Z 2 × Z 2 symmetry and τ h/v always anticommutes with W v/h .

V. CONCLUSIONS
In this article, we presented a method to study the ground-state properties of a two-dimensional Z 2 lattice gauge theory using the quantum approximate optimization algorithm.With this method, we are able to get good quality variational approximations while keeping circuits with a small depth.Hence, this allows to prepare the target state with a number of standard single-qubit rotations and CNOT gates comparable with the realistic expectations for near-term quantum technologies.
We focused on the minimal resources needed for an accurate description of the ground state in a quantum circuit setup, to show that interesting physics can be indeed observed despite the small size.In particular, we showed that both the behavior of Wilson loops and the entanglement entropy clearly distinguish the trivial and the topological phase, and they characterize the confinement/deconfinement transition as well.To reliably find good approximations of the ground state, we proposed a two-step protocol for QAOA, which produces regular optimal schedules that can be successfully transferred to larger sizes.In this respect, the two-step protocol outperforms a resource-costly global minimum search, as well as other local optimizations strategies that are prone to remaining trapped in "bad" local minima.However, the role of the noise brought by measurements and gates has been neglected, even though it will inevitably appear in a realistic implementation of our proposal.In general, the effect of noise on VQAs is still an important open question [2] and thus the study of the robustness of our proposal is worthy of further investigation.In the worst-case scenario, where noise prevents an accurate reconstruction of the cost function E P (γ, β) for optimization purposes, one might still use a "simulated" QAOA to infer good variational parameters to be provided to the actual quantum circuit, which will be used mainly for measuring physical properties.
We emphasize that the QAOA technique we propose can be easily combined to extend several proposals for the study of 2D LGTs through digital quantum simulations [52-54, 69, 76, 86-89].Digital quantum simulations of LGTs on small systems have already been implemented in trapped ion experiments [90,91] and superconducting qubit platforms [92][93][94][95].These experiments inspired several theoretical studies aimed at investigating the dynamics of the most important LGT excitations [82,[96][97][98][99].Our results provide a tool to efficiently initialize the ground states of LGTs, which, in turn, make it possible to engineer in a controlled way several of the excited states studied to explore the dynamical and topological properties of LGTs, including, for example, flux excitations and mesons.The system we considered can be regarded as a surface code perturbed by onsite interactions that provide a kinetic energy to its plaquette excitations [100,101].Hence, the study of its dynamics delivers information on the resilience of topological quantum memories in which anyons acquire a non-trivial dispersion.Furthermore, the topological order of the Z 2 LGT is the same of the most common topological quantum spin liquids, and our QAOA approach can be extended, for instance, to the study of quantum dimer models based on plaquette interactions, such as the Rokhsar Kivelson model [102], which displays this kind of topological phases and transitions on suitable lattices [46].
More in general, our variational quantum optimization successfully enables to explore the properties of Hamiltonians with non-trivial four-body interactions, which represent not only an essential element for designing topological phases but also a useful tool for encoding classical optimization problems [50,51,103,104].Such interactions are compatible with the native geometry and qubit gate connectivity of several recently developed quantum computation platforms, encompassing both twodimensional superconducting architectures, such as the Google Sycamore array [16,19], and programmable arrays of Rydberg atoms [105][106][107][108].In these systems, no additional overhead would be needed to map logical into physical qubits and measurements would give direct information on the addressed models, as in the case of the Z 2 LGT.
In conclusion, the combination of QAOA, initialization of the excitations and digital quantum simulation of their time evolution opens the path to study many aspects of the dynamics of the confined and deconfined phases in LGTs as well as the anyonic excitations appearing in topologically ordered phases.Quantum Technologies (TQT).GM acknowledges support from Austrian Science Fund through the SFB Be-yondC Project No. F7108-N38.Comparison of the accuracy with respect to the magnetic coupling h between global optimization (empty squares) and the two-step approach (full squares) described in the text.The data refer to P = 6 and |ψ0 = |ΩE .under investigation.In particular, we focus on a benchmark of our heuristic approach against a global minimum search, which, remarkably, yields similar-quality results for both phases, in terms of ground state fidelity, offering a good numerical validation of our scheme.
In addition, we comment on the transferability of the optimal schedules, obtained by either two-step or basin hopping, to larger system sizes, a strategy that could provide an educated guess to lower the computational cost for a new optimization.
Finally, we observe some patterns for optimal QAOA variational parameters obtained with the two-step scheme, in particular their smoothness as a function of m = 1 • • • P, similarly to other results for different QAOA applications [4,5,13].

Global optimization vs two-step scheme
In order to prove the effectiveness of the two-step optimization protocol, we compare it with a global minimum search, based on the basin hopping method [109] from the scipy.optimizePython library.In the latter algorithm, we allow up to 500 local minimizations, each of them initialized in the proximity of a previously found local minimum.The parameter space is explored with an effective temperature chosen to allow jumps between typical low energy minima.To reliably find the absolute minimum we run the basin hopping optimization 100 times and take the best result.
Figure 13 shows a comparison between the fidelity obtained with the global and the two-step optimizations, as a function of h for fixed P = 6.The initial state is |ψ 0 = |Ω E .For h < 3 the two-step approach can match the global optimization performance and it yields the same results, while for h ≥ 3 it finds a sub-optimal local minimum.However, we stress that, even in this case, the final fidelity is almost one: the difference in the accuracy between the two methods is much lower than a realistic experimental resolution.Moreover, the two-step protocol has the clear advantage of requiring only a single local optimization -on top of a modest computational overhead for the one-dimensional optimal ∆t grid search -to be compared with 500 × 100 local optimizations for the basin hopping method.Consequently, the twostep heuristics certainly requires drastically fewer function evaluations and it is, therefore, a better candidate to be implemented on a realistic quantum device and also much faster to simulate on a classical computer.
Regarding the transferability of the optimal schedules to larger system sizes, we use the optimal angles obtained for L = 3, either with basin hopping or with a two-step optimization, as initial guess for local optimizations of larger system sizes L = 4, 5. Specifically, for each value of h, we compare the best fidelity out of 10 BFGS local search runs, each of them initialized with the optimal 2P parameters previously found for L = 3, plus a small noise to facilitate the exploration of the energy landscape.
Our results are reported in Fig. 14(a), where we compare local minimizations starting from the L = 3 two-step optimal schedules (star symbols) local searches starting from to the L = 3 global minimum (circles), and the two-step process applied directly on the lerger system (squares).We find that the optimal angles returned by the two-step algorithm provide a better guess for larger systems, resulting in higher fidelity than a local search initialized with the global minimum for L = 3.This may be linked to the existence of some patterns in the optimal parameters found with the two-step scheme, in particular their smoothness as a function of m = 1 • • • P, as summarized in Appendix B 2. The performance of the two-step optimization applied directly on the target system L = 4 or L = 5 is instead comparable with the transfer of the schedule from L = 3, although the latter is slightly better at large magnetic fields.
However, once the optimal schedule for a given system size is known, it is convenient to leverage on that result to initialize the QAOA search for larger systems, instead of running a new two-step optimization from scratch.In fact, although the performances in terms of final fidelity are similar, the schedule transfer requires less iterations than the two-steps optimization.This is shown in Fig. 14(b), where we compare the number of BFGS iterations required in the final local minimum search on a system with L = 5 for different optimization strategies: schedule transferring from L = 3 to L = 5, from L = 4 to L = 5, and the two-steps protocol directly on L = 5.The latter requires in general a larger number of iterations and its overall cost must be added to the resources required for the optimization of the time step ∆t.

Smooth schedules
In the Z 2 LGT we studied in this paper, we found that the two-step optimization scheme produces smooth protocols for the optimal variational parameters more easily than other heuristic methods present in the literature, such as the standard application of iterative schemes based on parameter interpolation or Fourier component optimization [4,5].Moreover, it provides a minimum (γ * , β * ) for a chosen circuit depth P without requiring the solution for shallower circuits with P < P, contrarily to both the iterative methods just mentioned.
The presence of regular patterns in the optimal parameters suggests a comparison with a digitized quantum annealing scheme, such as the ones adopted to initialize the  16), with the corresponding value of h and P = 6 steps.The quantum circuit is initially prepared in the electric ground state |ΩE .two-step QAOA.In Sec.IV B we defined two possible annealing protocols, depending on the choice of the initial state.If |ψ 0 = |Ω E , we construct the time dependent Hamiltonian while if |ψ 0 = |Ω B we use In both cases, t ∈ [0, τ ] and at the end of the protocol H(t = τ ) = H targ .The corresponding parameters γ m and β m of a digitized quantum annealing are reported in Eq. ( 16) and Eq. ( 18), respectively.For a graphical representation of smooth optimal twostep schedules and a direct comparison with digitized Quantum Annealing, it is useful to consider the following more general protocol, as customary in Adiabatic Quantum Computation [110]: where s(t) ∈ [0, 1] is a monotone time-dependent parameter that interpolates between H E and H B .With this notation, we can identify Eq. (B1) with a process starting from s(0) = 0 and ending in s(τ ) = s f , with s f = h h+1 ; Eq. (B2), instead, corresponds to a process with s(0) = 1, ending again in s(τ ) = s f (both identifications are valid modulo an overall multiplicative factor).
A digitized Quantum Annealing process [4] consists in choosing a discretization of the time interval [0, τ ] into P small time steps ∆t m , such that m ∆t m = τ .Correspondingly, the continuous schedule s(t) is discretized into a sequence of short-time evolutions generated by H(s m ), where (B5) Once we have found optimal smooth QAOA parameters γ , β with our two-step QAOA scheme discussed in Sec.IV B, we can extract the corresponding digitized schedule s * m and compare it with the linear digitized quantum annealing protocol s dQA m that we used as an educated guess for the local minimization.
As examples of typical smooth QAOA optimal parameters, we report in Fig. 15(a)-(d) the schedules s * m corresponding to four different values of the coupling h, both below and above the "topological transition", with P ≥ 3 and initial state |ψ 0 = |Ω E .The dashed black lines correspond to the linear annealing schedule of Eq.( 16) we used as starting point for the local minimizations, with P = 6.In all four cases, it appears evident that as P increases, the parameters gradually approach a smooth continuous behavior, with the possible exception of a single localized irregularity, which seems to appear in Fig. 15(a) for P = 5, 6.This is not surprising, however, since we are preparing a state very close to the initial one.Thus, a large value of P could "overfit" the target state and many different parameters choices, usually nonsmooth, could yield similar accuracy.A comparison with Fig. 6 for the case h = 1, clearly shows a degradation of performance (almost-flat curve) of the infidelity vs P, exactly for P = 5, 6: this irregularity can thus be interpreted as a local lower-quality minimum or a saturation of the numerical precision of the algorithm.
For larger values of the coupling h, instead, we observe a clear continuity in the optimal schedule s * m , as we change both P and h.This leads to the interesting consequence that the optimal schedule for a given H targ and circuit depth P could be used as a seed to initialize the optimization for different values of h, requiring only a small fine-tuning of the parameters to adapt the schedule to the new target ground state.
A similarly smooth pattern is observed when we initialize the system in the magnetic ground state |Ω B , as reported in Fig. 16(a)-(d).The dashed black lines correspond here to the schedules s dQA m extracted from Eq.( 18), with P = 6.The main difference is that the smoothness now is more easily lost when targeting the deconfined phase, see panel(d), which is closer to the initial state.Similar comments on this irregularity apply as for the previous case, by comparing with Fig. 7(a).On a side note, we notice that the evident irregularity in panel(d) for P = 5, 6 involves a single point with a numerical value smaller than 0.4: this is not a significant feature, and it could easily be eliminated by an appropriate smoothing procedure with a likely improve in performance.
Unsurprisingly, the two-step optimization might get trapped in a (high-quality) local minimum even when we target the opposite phase: this is seen, e.g., for the outlier set of h = 1 and P = 5 in Fig. 16(a), which might be associated to a suboptimal minimum.This observation is once again consistent with the corresponding data in Fig. 7(a), where the curve for h = 1 shows a small spike in correspondence to P = 5.
Regarding the comparison with the linear dQA protocol s dQA m (dashed black lines), in both Figs. 15 and 16, the overall monotonicity of optimized s * m is the same of the original schedule, i.e., an increasing function of m when |ψ 0 = |Ω E , and a decreasing function when |ψ 0 = |Ω B .However, when targeting states in a phase that differs from the initial one, the optimal schedule deviates more and more from the original Ansatz, highlighting the importance of the local optimization of the parameters.

Energy landscape
In the Z 2 LGT model, the energy landscape associated to the QAOA Ansatz is characterized by the presence of many local minima, as discussed in Sec.IV A. This makes the employment of a clever optimization strategy, such as the two-step protocol or schedule transferability, a necessity to target reliably low energy minima.However, for general variational problems, it might happen that it exists a deep minimum in the energy landscape, associated to a state with small or no overlap with the target one.This is not the case for the problem under investigation, where there is a clear correlation between the fidelity and the variational energy for the minima in the energy landscape, see Fig. 5 in the main text.
Here, we show that this correlation holds also for the larger systems considered in this paper, L = 4 and L = 5, corresponding to 32 and 50 qubits respectively.We repeat the analysis of Sec.IV A: focusing on |ψ 0 = |Ω E , we perform 100 QAOA runs with random initial parameters, targeting states in the deconfined phase h ≥ 3.In Fig. 17 we plot the infidelity 1−F P vs the residual energy of the minima found with this procedure, for both L = 4 (a.panel) and L = 5 (b.panel) and a circuit depth of P = 6.Interestingly, it appears that increasing the system size leads to a sharper correlation between energy and fidelity, compared to the data presented in Fig. 5.

Figure 1 .
Figure 1.Representation of a star operator Av (in red) and a plaquette operator Bp (in blue), with the corresponding qubits on the links (solid circles).

Figure 2 .
Figure 2. The two non-contractible 't Hooft loops τ h , τv and a simple example of how a contractible 't Hooft loop operator τC is written as a product of star operators.The dotted link indicates a cancellation due to (σ x l ) 2 = 1.

Figure 3 .
Figure 3. (a): Circuit implementation of the operator e −iγ H B acting on a single plaquette, with a target and 3 control qubits.The states |pi are expressed in the σz eigenbasis and Up(γ) = e iγ σz is a single-qubit rotation around the z axis.(b): implementation of the single-qubit operations that describe e −iβ H E .(c): example of the sequential implementation of the operator e −iγ H B on two neighboring plaquettes.The qubit |p4 is used first as the target qubit for the first plaquette and aftwerwards as a control qubit for the second.

Figure 5 .
Figure 5. Infidelity vs residual energy rescaled over the magnetic coupling h, for h = 5, 4, 3. Data refer to 100 local optimizations on a system with linear dimension L = 3, initial random guess of the QAOA parameters, and circuit depth P = 5.The initial state is |ψ0 = |ΩE .

Figure 6 .
Figure 6.Infidelity vs the number of QAOA layers P for L = 3 and several values of the magnetic coupling h.The initial state is the electric ground state |ΩE , thus the convergence towards the exact GS is faster for smaller values of h.We show the best result out of ten local BFGS minimizations following the heuristic two-step optimization method.

Figure 7 .
Figure 7. (a): Infidelity vs the number of QAOA layers P for L = 3 and several values of the magnetic coupling h.Data correspond to the best out of 10 results obtained in the twostep optimization, performed on a state initially prepared in the toric code ground state |ΩB , hence the convergence is now faster for larger couplings h.(b):Comparison between two-step QAOA performance by starting from the electric and the magnetic ground states: we plot the fidelity vs magnetic coupling h at fixed values of P.Here P only takes into account the number of QAOA layers with parameterized gates, while it does not include the computational overhead for the preparation of |ΩB , compared to preparing |ΩE .

F 5 Figure 8 .
Figure 8. Fidelity vs the magnetic coupling h.The data are obtained by using the two-step optimal schedules for L = 3, as initial guess for 10 local BFGS optimizations performed on the larger systems L = 4, 5 (we fixed P = 6).Here, we report the best result out of the 10 runs.The full and empty symbols respectively correspond to |ψ0 = |ΩE and |ψ0 = |ΩB .

Figure 9 .
Figure 9. Creutz ratio χ(l, l), defined in Eq.(19), for two different loops in a system with L = 5.Inset: expectation value of Wilson operators W l,l , corresponding to the data in the main plot, vs the coupling strength h.The vertical dashed line indicates the critical value of the coupling hc.All data refer to the best energy out of 20 BFGS local optimizations, with P = 6 and |ψ0 = |ΩE , |ΩB , performed on a system of linear size L = 5 and initialized with the optimal parameters found for L = 3. Empty symbols in the inset are data obtained with |ψ0 = |ΩE .

Figure 10 .
Figure 10.(a): Graphical representation of the subsystem used to compute the topological and entanglement entropy, with the tripartition A, B, C highlighted.Notice that a total of Nv = 5 vertex are cut by the outer edge of A ∪ B ∪ C. Empty dots indicate the presence of PBC in the lattice, thus identifying the upper edge with the lower one and the right edge with the left one.(b): Entanglement and topological entropy as a function of the coupling h.Notice that we plot −Stopo to make it positive.

Figure 11 .
Figure 11.Four lowest energy states in different symmetry sectors of the non-contractible 't Hooft loop operators vs coupling strength h.All data refer to the best solution found for L = 3, P = 6, and |ψ0 = |ΩE .States different from | + + P are obtained by acting with non-contractible Wilson loops after the unitary evolution, while the subscript P indicates that they are obtained by QAOA and thus are not exact eigenstates.The solid lines correspond to the 5 lowest eigenvalues obtained with exact diagonalization: The 5th eigenvalue, corresponding to the first excited state, is shown to highlight the topological gap in the deconfined phase.Only four lines are visible because | + − and | − + are exactly degenerate, as their approximations obtained with QAOA.

Figure 12 .
Figure 12.(a) A stripe of the 3 × 3 lattice with periodic boundary conditions (b) Quantum circuit to implement the plaquette operator e iγB on all three plaquettes.The labelling of the qubit lines emphasizes that all the boundary qubits are shared with the plaquette stripes above and below.The partially depicted green CNOTS are related to the simultaneous implementation of the same algorithm on the neighboring stripes: qubits 2,3,5,6,8,9 act as controls also for the circuit in neighboring plaquettes.

Figure 13 .
Figure13.Comparison of the accuracy with respect to the magnetic coupling h between global optimization (empty squares) and the two-step approach (full squares) described in the text.The data refer to P = 6 and |ψ0 = |ΩE .

Figure 14 .
Figure 14.(a) Comparison of the two-step vs global schedules for L = 3 as an educated guess for optimization on larger system sizes L = 4, 5.The data refer to the best out of 10 local BFGS minimizations, run by starting close to the optimal two-step schedule (blue stars) or optimal global-schedule (green circles).The plot also shows the fidelity of optimal two-step schedule (orange squares) for L = 4, 5. (b) Number of iterations required for the convergence of the final BFGS optimization; comparison between transferring the schedule from L = 3 to L = 5 (blue stars), transferring from L = 4 to L = 5 (red triangles), and two-steps optimization directly on L = 5 (orange squates).All data refer to P = 6 and |ψ0 = |ΩE .

Figure 15 .
Figure 15.Variational parameters associated to the best result out of 10 local minimum searches, following the two-step optimization procedure.In each panel, the black dashed line corresponds, through Eq. (B5), to the linear annealing schedule defined by Eq.(16), with the corresponding value of h and P = 6 steps.The quantum circuit is initially prepared in the electric ground state |ΩE .
B4)with m = 1 • • • P. The resulting expression can be further simplified with a first order Trotter split-up, neglecting quadratic terms in ∆t m .Thus, the final state is written as the variational Ansatz in Eq. (12), with fixed param-γ m γ m + β m , ∆t m = γ m + β m .

Figure 16 .
Figure 16.Variational parameters associated to the best result out of 10 local minimum searches, following the two-step optimization procedure.In each panel, the black dashed line corresponds, through Eq. (B5), to the linear annealing schedule defined by Eqs.(18), with the corresponding value of h, and P = 6 steps.The quantum circuit is initially prepared in the magnetic ground state |ΩB .

Figure 17 .
Figure 17.Infidelity vs residual energy for 100 minima in the energy landscape found with random initialization of local BFGS searches, with |ψ0 = |ΩE and P = 6.Panel a. corresponds to the lattice with 4 × 4 plaquettes, panel b. to 5 × 5 plaquettes.We only show the data that falls in the interval (EP − EGS)/h ∈ [0, 4.2].