Circuit-based digital adiabatic quantum simulation and pseudoquantum simulation as new approaches to lattice gauge theory

Gauge theory is the framework of the Standard Model of particle physics and is also important in condensed matter physics. As its major non-perturbative approach, lattice gauge theory is traditionally implemented using Monte Carlo simulation, consequently it usually suffers such problems as the Fermion sign problem and the lack of real-time dynamics. Hopefully they can be avoided by using quantum simulation, which simulates quantum systems by using controllable true quantum processes. The field of quantum simulation is under rapid development. Here we present a circuit-based digital scheme of quantum simulation of quantum ℤ2 lattice gauge theory in 2 + 1 and 3 + 1 dimensions, using quantum adiabatic algorithms implemented in terms of universal quantum gates. Our algorithm generalizes the Trotter and symmetric decompositions to the case that the Hamiltonian varies at each step in the decomposition. Furthermore, we carry through a complete demonstration of this scheme in classical GPU simulator, and obtain key features of quantum ℤ2 lattice gauge theory, including quantum phase transitions, topological properties, gauge invariance and duality. Hereby dubbed pseudoquantum simulation, classical demonstration of quantum simulation in state-of-art fast computers not only facilitates the development of schemes and algorithms of real quantum simulation, but also represents a new approach of practical computation.


Introduction
Quantum simulation and quantum computation can efficiently solve some problems that cannot be efficiently solved in classical computers [1][2][3][4], and is under extensive studies worldwide, thanks to the rapid development of quantum science and technology. A controllable quantum system, which may even be universal or programmable, simulates various quantum systems, whose physical properties can be conveniently investigated with various parameter values. Even with only tens of qubits, far less than those in full fault-tolerant quantum computing, and even in the presence of some noises, a quantum machine can perform some tasks surpassing its classical counterparts, exhibiting the so-called quantum supremacy [5][6][7]. Many quantum simulations are such tasks. Hopefully they not only can solve specific problems, but also represent new scientific methods [3]. Quantum simulations of important models in theoretical physics are enabled by some quantum algorithms, including the Trotter decomposition [2], the sparse Hamiltonian quantum walk [8], the dense Hamiltonian density matrix exponentiation [9][10][11], the adiabatic algorithm [12][13][14], and so on. It is timely even to develop various quantum softwares [15,16].
An important battlefield of quantum simulation appears to be lattice gauge theory (LGT) [17][18][19], the major non-perpurbative approach to gauge theory. In particle physics, JHEP08(2020)160 gauge theory is the framework of the Standard Model, describing both electroweak and strong interactions among elementary particles, and is also a guide beyond the Standard Model. A LGT is a gauge theory defined on a spacetime lattice in path integral formalism or on a space lattice in Hamiltonian formalism, making the degrees of freedom countable and convenient for numerical calculations. Gauge theory is also important in condensed matter physics, where the lattice can be a real structure. It is often an effective description of constraints in strongly correlated systems, and uses emergent gauge fields to characterize topological orders, which exist in fractional quantum Hall effect, spin liquids and possibly in high temperature superconductivity, as well as in topological quantum computing, etc. Topological order represents a new and active paradigm beyond the traditional frameworks of symmetry breaking and Fermi liquid theory.
Using Monte Carlo (MC) simulation, LGT has made great achievements [20][21][22][23][24]. But there are also difficulties, including the lack of real-time dynamics because of the use of Euclidean spacetime, and the notorious Fermion sign problem [25,26], which was shown to be NP hard [27]. The cause of the Fermion sign problem is that in presence of Fermions, the Boltzmann weight of an appropriately defined configuration, to which the quantum problem is mapped, may become negative, therefore the partition function oscillates violently, consequently the importance sampling in MC becomes invalid. These difficulties are related to the unsolved problems in quantum chromodynamics, such as color confinement and phase diagram of quark-gluon plasma. The Fermion sign problem exists in most of the MC-based methods such as quantum MC (QMC), except special algorithms for some specific models [28], and in some special issues to be told below.
As a new approach avoiding Fermion sign problem and an ideal avenue to study quantum phase transition (QPT) and real-time quantum dynamics, quantum simulations of LGTs are under study. They were first theoretically explored, mostly but not exclusively, in cold atom platforms for the simulations of U(N), SU(N) and Z n LGTs [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43]. Other proposals were also made [44,45]. Among these schemes, some are analog, based on physical mapping between the Hamiltonians of simulated and simulating systems, while others are digital, based on Trotter decomposition of the finite-time evolution to many small steps that are much easier to be implemented. A quantum-classical algorithm was developed for two-site Schwinger model [46].
Quantum Z 2 LGT is the simplest quantum LGT [21,22,[55][56][57]. On one hand, Z n theory can be obtained as the discretization of U(1) theory, suitable for quantum simulations [58,59]. On the other hand, quantum Z 2 LGT is also important in condensed matter physics [57,60]. Z 2 toric code model [61], which is important in topological quantum computing, can be regarded as a variant of quantum Z 2 LGT, and has been experimentally demonstrated in a quantum simulation using nuclear magnetic resonance [62,63].

JHEP08(2020)160
In two spatial dimensions, quantum Z 2 LGT is dual to quantum Ising model in a transverse field [21,22,[55][56][57], which is thus often invoked for QPT properties, especially the critical point [64][65][66][67]. Besides, Z n theory in one spatial dimension was studied by using density matrix renormalization group (DMRG) [58,59], and was studied in two spatial dimensions by using tensor network techniques, as the low energy limit of toric model in a magnetic field under the constraint of gauge invariance [68]. Direct DMRG study of quantum Z 2 LGT is difficult as realizing plaquette interactions consistently with the gauge symmetry is challenging. But such study is technically possible for Abelian and non-Abelian gauge theory in 2D by using symmetry-preserving tensor network techniques [69][70][71]. Interestingly, the duality between such a spin gauge theory and a generalized Ising model allows for scalable quantum simulation with Rydberg atoms [72].
Coupling of Z 2 gauge field with various kinds of matter has also been studied, starting with Ising matter [56]. Fractionalization of electrons is obtained in theories of strong correlations with Z 2 gauge fields, giving rise to the so-called orthogonal metals [73][74][75]. Interestingly, the issue of Fermions coupled with Z 2 gauge field is exactly one of the special issues free of sign problem in MC-based methods [76], another issue being Fermions with an even number of flavors [77]. The cases with both of these characteristics were studied by using QMC [78][79][80]. Coupling of Fermions with Z 2 gauge field with Gauss law not imposed but emerged was also studied by using QMC [81,82]. These QMC studies were all in two spatial dimensions. A modified Z 2 gauge theory coupled with Fermions was studied analytically in one and two spatial dimensions [83]. One-dimensional Bosons coupled with Z 2 gauge field with Gauss law not imposed but emerged was studied also by using DMRG [84,85]. A quantum link model of QED was approached by using tensor network method [86].
In quantum Z 2 LGT, defined on a square lattice, each link is occupied by one qubit. In the 2 × 2 × 2 lattice, which is the smallest three-dimensional lattice, there are 24 links (figure 1). In the 3 × 3 lattice, the second smallest two-dimensional lattice, there are 18 links (figure 1). In a quantum algorithm, additional qubits may also be needed as ancillas in simulating the adiabatic evolution, and in phase estimation simulating the measurement, and so on. Symmetries may reduce some degrees of freedom, but at the price of introducing nonlocal interactions, for example. Regarding the quantum adiabatic algorithm, tens of thousands steps are needed in Trotter decomposition. Therefore, at present time, it seems difficult for the experimental platforms to fully meet the requirements on the qubit number, error rate and coherence time.
Therefore, it appears interesting to use classical high-performance computing platforms to demonstrate quantum simulation in general, and that of quantum Z 2 LGT in particular. We call it pseudoquantum simulation, which serves as a benchmark for real quantum simulation, and facilitates the development of quantum algorithm and quantum softwares. Meanwhile, it is also a new method of computation and simulation, providing useful results on the computed problems. JHEP08(2020)160 Such pseudoquantum simulation is realistic if it is run on a fast enough classical computing platform. An example is the latest graphics processing unit (GPU) parallel computing architecture, which greatly accelerates large-scale complex scientific computation. Indeed, a GPU simulator called Quantum Exact Simulation Toolkit (QuEST) has been developed, as a new software platform simulating the quantum circuit [87]. It is based on the software platform and application programming interface CUDA created by Nvidia, which allows the development of parallel program using a CUDA-enabled GPU. QuEST is designed as a C library, and allows quantum codes to be deployed in a variety of computing platforms. With high precision, it can simulate 29 to 31 qubits on a single Nvidia GPU card, as detailed below. Therefore it fits well the need of our pseudoquantum simulation of quantum Z 2 LGT.
In this paper, we present a scheme of circuit-based digital quantum simulation of (2+1)dimensional and (3+1)-dimensional quantum Z 2 LGT. We use universal quantum circuit to implement the quantum adiabatic algorithm, and we generalize the Trotter decomposition and symmetrized Trotter decomposition to the case of adiabatic varying the Hamiltonian in each step of decomposition. As the first work and proof of principle, here we consider pure gauge theory, without coupling to Fermions yet. Moreover, we perform a complete pseudoquantum simulation, that is, we classically demonstrate the scheme of quantum simulation, by using QuEST in Nvidia Tesla K40m and V100 GPU cards, which operate 1.682TFLOPS and 7.834TFLOPS, with double precisions, respectively. Our work demonstrates the advantages of quantum simulation as well as the usefulness of pseudoquantum simulation. Meanwhile we obtain useful results regarding various properties of quantum Z 2 LGT. It seems to contain the first numerical result on quantum Z 2 LGT in three spatial dimensions, including the features of first-order QPT, as hinted by the fact that the thermal phase transition in the classical Z 2 LGT is first-order in 4 spatial dimensions [88,89].

JHEP08(2020)160
The rest of the paper is organized as the following. In section 2, we introduce the quantum Z 2 LGT, the lattices we consider, the quantum adiabatic algorithms and their realization in terms of the quantum circuits, our generalization of the Trotter and symmetrized Trotter decompositions, the simulation of measurement, as well as the hardware we use to demonstrate the scheme of quantum simulation. In section 3, we present the results from our GPU computation simulating the quantum simulation, including the expectations of Wegner-Wilson loop operators and of the Hamiltonian, which are used to determine the critical points and orders of the QPTs, and confirm the self-duality in 3+1 dimensions. We also calculate the densities of states, which confirm the gauge invariance, the self-duality in 3+1 dimensions and its absence in 2+1 dimensions. We will also present the evidences of topological nature of QPT. Finally, a summary is made in section 4.
2 Schemes and algorithms of quantum simulation and pseudoquantum simulation 2.1 Quantum Z 2 lattice gauge theory LGT is defined on a lattice. It was first proposed as a generalization of Ising model, elevating the global up-down symmetry to a local symmetry, but without spontaneous magnetization at the phase transition, which is characterized by onset of topological order rather than symmetry breaking. So there is no local order parameter [55]. On the other hand, with some differences in details, Z n LGT can also be obtained by discretizing U(1) gauge theory, by defining the matter field on the sites of a lattice, and the n-valued gauge potential and thus electric field on the links between sites.
Here we focus on the quantum Z 2 LGT with the Hamiltonian [57,60] H = Z + gX, where g is the coupling constant, with Z defined for each elementary plaquette (the smallest square) (figure 1). For qubit l, σ z l |0 l = |0 l , σ z l |1 l = −|1 l . The Hamiltonian of the classical Z 2 theory is Z only, with each operator σ z l reduced to a classical variable. The quantum nature of H is due to the noncommutativity between σ x l and σ z l , resulting the competition between Z and X, in a way like that between energy and entropy in a thermal phase transition. The coupling constant g is a control parameter, playing a role in QPT similar to the role of temperature in thermal phase transition.
For convenience, in our numerics, the links and thus the qubits are numbered as in figure 1. For three spatial dimensions d=3, we use 2 × 2 × 2 lattice with periodic boundary condition, where there are 24 links and 24 elementary plaquettes. We use 25 qubits, one of which is the ancilla. For two spatial dimensions d=2, we use 3 × 3 lattice with periodic boundary condition, where there are 18 links and 9 elementary plaquettes. We use 19 JHEP08(2020)160 qubits, one of which is the ancilla. We assume periodic boundary condition, that is, each lattice is a torus. Although the lattice sizes are very small, the key features of the quantum Z 2 LGT do appear. In general, for a d-dimensional square lattice with linear size L, the number of links is N l = dL d , the number of plaquettes is N p = N l (d − 1)/2. This theory possesses Z 2 gauge invariance, similar to Gauss law, dictating that each eigenstate |ψ of H must satisfy where G i ≡ l∋i σ x l is the product of the σ x l 's on all the links ending at each lattice point i. One gauge invariant operator is Wegner-Wilson loop operator, which is defined as operator is not necessarily along a non-contractible loop, for example, Z is a Wegner-Wilson loop. But those along non-contractible loops play special roles. Generalizing the well known case in d=2, here we define the special Wegner-Wilson loop operators along non-contractible loop C µ , µ = 1, · · · , d (figure 2). Periodic boundary conditions means that topologically they encircle a torus. As far as it encircles the lattice in µ direction, the details of C µ does not matter. σ x and any product of σ x 's are also gauge invariant. G i σ x G −1 i = σ x . Generalizing the well known case in d=2, here we make a definition of the so-called 't Hooft loop operator in any dimention d, which is a product of σ x l pierced by a non-contractible (d − 1)-dimensional surface C µ on the dual lattice, µ = 1, · · · , d (figure 2). A non-contractible (d − 1)-dimensional surface can be deformed by using G i operators, because of gauge invariance. C µ is subscripted as µ as it can be deformed, if needed, to be parallel to C µ . Periodic boundary conditions means that topologically they encircle a torus. In d=2, C µ is a loop.
JHEP08(2020)160 Figure 2. (a) d=3 Lattice, with the indication of non-contractible loops C x , C y and C z on the direct lattice, and non-contractible C y on the dual lattice. C x , C y and C z are used in defining specific Wegner-Wilson operators W x , W y and W z . C x and C z are similar to C y . They are used in defining 't Hooft operators V x , V y , V z . (b) d=2 Lattice, with non-contractible loops C x and C y on the direct lattice, used in defining specific Wegner-Wilson loop operators W x and W y , as well as non-contractible loops C x and C y in defining 't Hooft operators V x and V y .
The degeneracy on a d-dimensional lattice with genus G is 2 dG . On the lattice considered here with period boundary condition, the degeneracy is 2 d . The degeneracy changes when g = 0, with only |V ν = 1, ν = 1, · · · , d remaining as the ground state. But they represent different topological sectors, because V µ is conserved while each non-contractible (d − 1)dimensional surface can deform.
Depending on g, there are two different phases, deconfined phase at small g and confined phase at large g, separated by a QPT at the critical point g c , where there is a phase transition, which however cannot be characterized by a change of symmetry. There is Z 2 topological order in the deconfined phase, while the confined phase is trivial, and QPT in this theory is topological [55,57,60].

Circuit-based quantum adiabatic algorithm
Our digital scheme of quantum simulation of the quantum Z 2 LGT is based on the quantum adiabatic algorithm, and we use quantum circuits consisting of one-qubit and two-qubit quantum gates to implement it. In the adiabatic evolution, g varies from 0 to a large enough value g f , passing g c . The slow variation allows the system to adapt to the instantaneous ground state. According to the adiabatic theorem, the state of the system, starting as a ground state of the initial Hamiltonian H(g = 0) = Z, evolves as ground state of H(g), and ends up as the ground state of the final Hamiltonian H(g f ).

JHEP08(2020)160
In our simulation, each time g is updated, the state evolves for a very short time, as realized by the quantum circuit under the present g value, then g is updated again, and the state evolves under the new value of g. The iteration continues until a final value g f of g. Therefore, the adiabatically varying Hamiltonian is stepwise. We divide the evolution to N s steps, and each step is further divided to n substeps, and g varies at each substep. This so-called "substep" really corresponds to the "step" in Trotter decomposition. The reason of referring to the decomposition steps as substeps is that the evolution is paused after a period called step, and calculations, or called pseudo-measurements, are done, afterwards the evolution is resumed. In real quantum simulation, the measurement depends on the actual situation.
Therefore, the stepwise Hamiltonian, for the m-th substep within k-th step, is (m = 1, · · · , n), n is the total number of substeps in each step, which is freely set. g s is the increase of g in each step, which lasts time t s , δ = g s /n is the increase of g in each substep. The total number of steps is N s = g f /g s , and the total number of substeps is nN s . The evolution in k-th step is n m=1 e −iH k,m ts n , (2.10) while the total evolution is Ns k=1 n m=1 e −iH k,m ts n . (2.11) The evolution in each substep of t s /n, under H k,m , consisting of two noncommutative parts Z and g k,m X, can be decomposed, in an asymmetric way, into consecutive evolution of Z and g k,m X, e −iH k,m ts n ≈ e −iZ ts n e −ig k,m X ts n . (2.12) Therefore the evolution in k-th step is n m=1 e −iZ ts n e −ig k,m X ts n . ( (2.14) Therefore, the error for one substep, in the asymmetric decomposition (2.12), is

JHEP08(2020)160
where we have considered that each Z is noncommutative with 4 σ x 's. It is calculated that Therefore the total error of the asymmetric decomposition is We have also used the symmetric decomposition Therefore the evolution in k-th step is which generalizes the symmetrized Trotter decompositon. It can be rewritten as which indicates a more convenient way of execution in our computation.
is of the order of 16N p , for the following reason. Each Z is noncommutative with 4 σ x 's, hence [Z , X] is the sum of 4 products of one σ y and 3 σ z 's. Each σ y is noncommutative with the 2(d − 1) Z 's of the plaquettes sharing with the link l. On the other hand, each product of one σ y and 3 σ z 's is noncommutative Another way of reasoning is the following. Each σ x is shared by 2(d − 1) plaquettes, thus [Z, σ x l ] yields 2(d − 1) products of σ y and 3 σ z 's. Each product is noncommutative with Z 's of the 2(d − 1) plaquettes, and with the 4 σ x 's on the same plaquette. Consequently, this is the same as above.
Therefore the error, in the symmetric decomposition (2.18), is It is calculated that

JHEP08(2020)160
Therefore the total error of the symmetric decomposition is The ratio of the errors of the symmetric and asymmetric decompositions is ∼ g f ts/n. With g f = O(1), the total error of the symmetric decomposition is less than the asymmetric one by a factor of t s /n. We have done our simulations using both decompositions. The results from the symmetrized decomposition is clearly better and thus presented below. Note that our decompositions are different from, and generalize, the usual Trotter decomposition and symmetrized Trotter decomposition, each of which repeats a constant evolution for a number of times.
We now discuss how to implement the evolution in each substep, or called each decomposition step. First, e −iZts/n = e −iZ ts/n , hence the evolution of Z can be realized by the consecutive evolution of all the plaquettes.
Evolution of e −iZ ts/n , for each palquette , is realized in terms of a quantum circuit, where there is also an ancilla, shown in figure 3, where CN OT l,a is a controlled-NOT gate controlled by the qubit l and targeting on the ancilla a, A −1 is the product of these CNOT gates in reversed order, R a z (φ) ≡ e −iσ a z φ/2 is single-qubit gate on ancilla representing rotation of angle φ around z-axis. Initially, the ancilla is set to be |r = |0 . Preceding R a z (−2t s /n), each CN OT l,a flips |r if and only if the control qubit l is |1 . Therefore R a z (−2t s /n) acts as e its/n when there are even number of |1 's on the plaquette, and acts as e −its/n when there are odd number of |1 's on the plaquette. This is precisely the effect of e −iZ ts/n . Afterwards, the four CNOT gates after R a z (−2t s /n) return |r to |0 , which can be used for the next plaquette. So we only need one ancilla. If we set |r = |1 initially, this circuit can be used to simulate the reversed evolution e iZ ts/n , in other words e −iZ t for reversed time t = −t s /n.
It is straightforward to realize the evolution under the other part gX in the Hamiltonian, where R l x (φ) on qubit l represents rotation of angle φ around the x-axis. For m-th substep of k-th step, g = g k,m .
Our realization of e −iZ ts/n , as given in (2.25) and figure 3, is a direct application of the standard strategy for the evolution under an interaction that is a tensor product of σ z operators [91]. In some previous schemes of quantum simulation of LGTs, the four-body interactions are obtained stroboscopically through a sequence of two-body interactions with JHEP08(2020)160 Figure 3. The quantum circuit realizing e −iZ ts/n or e iZ ts/n , depending on whether the initial state of the ancila |r is set to be |0 or |1 . ancillary degrees of freedom, and gauge invariance in each step of Trotter decomposition is emphasized [38,[40][41][42].
Here, as Z and σ x l are gauge invariant operators, the evolution operators are gauge invariant in every substep of the digital decompositions.

Preparation of the initial ground state
The adiabatic quantum simulation starts with g = 0, i.e. H = Z, of which there are degenerate ground states satisfying which can be prepared by using the quantum circuit in figure 4. Each qubit is initially in |0 and is then transformed to be (|0 + |1 )/ √ 2 by using a Hadamard gate H. Therefore the state of the system is in the equal superposition of all basis states, Each qubit is not entangled any other qubit, each plaquette is also in the equal superposition of all its basis states. Then the four CNOT gates between the four qubits of one plaquette and the ancilla initially in |0 a produce the state 1 √ 2 (|r = 0 a |Z = −1 + |r = 1 a |Z = 1 ), where |Z = ±1 are states of all the qubits on the lattice satisfying Z |Z = ±1 = ±|Z = ±1 . The ancilla is entangled with the qubits on the lattice, with σ a z ≡ 1 − 2r = −Z in each branch. Then the measurement operation M on the ancilla projects it to |r = 0 , thereby selects the state of qubits on the lattice to be |Z = −1 .
Afterwards, the ancilla is returned to |0 a by other four CNOT gates (figure 4). The same ancilla is ready to work on another plaquette, which may or may not share a qubit with a plaquette already worked on. Similar procedure goes on, till all plaquettes have been JHEP08(2020)160 worked on. The state of the system is then an eigenstate of Z for each , with eigenvalue −1, and is an equal superposition of all configurations satisfying Z = −1 for each .
To summarize, we perform consecutive projections P , where P = |Z = −1 Z = −1|. Then the state P |ψ 0 is exactly a ground state of Z, and is the eigenstate of the 't Hooft operators V µ with eigenvalue 1, for all µ's. It satisfies the gauge invariance. This ground state can adiabatically evolve to the ground state for g = 0.
The evolution of each decomposition substep is also gauge invariant, so the state preserves gauge invariance during the evolution. For a gauge invariant state |ψ , satisfying G i |ψ = |ψ , and the evolution of a time period τ under a gauge invariant operator O satisfying Therefore, with the initial state gauge invariant, gauge invariance is always preserved in each evolution under Z or σ x l in the digital decompositions. The other degenerate ground states at g = 0 can be obtained by using W µ 's as described above. This action changes the topological sector. When g = 0, the state adiabatically evolves to the lowest energy state in this sector, which is not the ground state.
Alternatively, for the adiabatic preparation of the ground states for various values of g, one can also start with the ground state of X. For this approach, one had better redefine the Hamiltonian as X + KZ, with the coupling constant K corresponding to 1/g. The ground state of X is with σ x = 1 for all qubits. With the increase of K, the ground state evolves from confined phase to deconfined phase. During the adiabatic evolution, the ground state remains in the topological sector of V µ = 1 for all µ's. To enter other sectors and adiabatically approaches the other degenerate ground states of H(g = 0), one can also use the actions of W µ 's.
In our demonstration, we use the first method, as it is easily simulated in our computing. Moreover, QuEST provides a method function of controlling the ancilla to simulate the collapse to the destined state in the GPU simulator. In the real quantum simulation, ancilla measurements of N p times, each conditioned on the result of the previous one, make the success rate only 2 −Np . Hence the second approach is preferred. We have actually also tried it in our demonstration, which yields result consistent with the first approach, so is omitted here, as the emphasis is on the deconfined phase. JHEP08(2020)160

Measurement of physical quantities
The energy E in state |ψ is the expectation value of H, H = Z + g X , where Z = ψ|Z|ψ = i z i P (z i ), X = ψ|X|ψ = i x i P (x i ), where {z i } and {x i } represent the eigenvalues or measurement results of Z and X, respectively, P (z i ) and P (x i ) are the corresponding probability distributions, called densities of states hereby. In our simulator, we obtain E by summing up Z and g X . We also calculate the expectation values of Wegner-Wilson loop operators.
We use CUDA parallel acceleration method to count the statistical summation of all basis vectors on GPU. We write our own codes for the calculation of the measurement results, which are not included in QuEST.
The measurement is simulated at the end of each step, i.e. when g = kg s , (k = 1, · · · , N s ). All these quantities can be calculated in terms of the distribution {P (z i )} in the representation {σ z l }, as W C = ψ| l∈C σ z l |ψ , Z = − ψ| l∈ σ z l |ψ , X = − l ψ|Hσ z l H|ψ , where H is Hadamard gate. We mention that in real quantum simulation, E can also be measured by using quantum phase estimation. For an eigenstate |u of a time-independent H, the evolution for a time duration t is e −iHt |u = e −iφ |u , where φ ≡ Et. e −iHt can be realized in a way similar to the decompositions described in section 2.2, but now g is fixed. That is, where n is the number of decomposition steps here. These are Trotter decomposition and symmetrized Trotter decomposition. The errors are just n times those for t/n, given in eq. (2.15) and (2.21) with g k,m replaced as g, that is, 2gN p t 2 n 2 and || i 24 For the purpose of quantum phase estimation, we need the conditional evolution controlled by an ancilla, which can be realized in terms of controlled gates, and ancilla controlling the time direction. For Z part, the method is as shown in figure 3. For X part, the method is as shown in figure 5.

Computational hardware
We use a Nvidia Tesla K40m GPU card, which was used by QuEST team in their simulation of 29 qubits with float decision [87], as well as a Nvidia Tesla V100-SXM2-32GB GPU card. We estimated the maximal scales of quantum simulations that QuEST can simulate under different precisions, as listed in table. 1. We use double precisions.

Adiabaticity and parameter values
We now estimate the total errors in the adiabatic process of the quantum simulation, using ∆ asy ≈ N 2 s N p g s t 2 s n , ∆ sym ≈ 4 9 N 3 s N p g 2 s t 3 s n 2 , given in eqs. (2.17) and (2.24). We choose the final value of g to be g f = 2. For each step, the increase of g is set to be g s = 0.001. For each substep in the decomposition, the increase is g s /n. Thus the number of steps is N s = g f /g s = 2000, while the number of substeps is 2000n. n is different in different cases as described in the following.
First consider the asymmetric decomposition. For d=3 lattice, the number of plaquettes is N p = 24. The number of substeps in the decomposition is set to be n = 200. The time for each step is chosen to be t s = 0.1. The total evolution time is then t f = N s t s = 200. t 3 s /n 2 = 0.25×10 −7 . The total error is about 2.1×10 −3 . For d=2, the number of plaquettes is N p = 9. As the number of qubits are less than in d=3, we set the number of substeps to be n = 5000. We choose the time step to be t s = 0.2. The total evolution time is then t f = N s t s = 400. Now t 3 s /n 2 ≈ 0.32 × 10 −9 . The total error is about 1.1 × 10 −5 .

JHEP08(2020)160
Now consider the asymmetric decomposition. As the error is larger than the symmetric one by one factor of t s /n, we use larger value of n = 500. We use smaller value of t s = 0.01 for d=3, and t s = 0.02 for d=2. Therefore, t 2 s /n = 0.2×10 −6 for d=3, and t 2 s /n = 0.8×10 −7 for d=2. Therefore the total error is 1.92 × 10 −3 for d=3, 2.9 × 10 −3 for d=2.
These parameter values are chosen to allow the computation to be completed in an acceptable time under adiabatic condition. The computation is proportional to 2 Nq , where N q is the number of qubits, which is equal to N l + 1 in the present model. It is also proportional to N s and n. For the parameter values given above, the time for the computation based on asymmetric decomposition, run on a K40m server, is about 7 days for d = 3 and 16 hours for d = 2, while the time for the computation based on symmetric decomposition, run on a V100 server, is about 28 hours for d = 3 and 4 hours for d = 2. The accuracies are all acceptable. The difference between the computation times is due to the hardware difference rather than the decomposition methods.
Our simulation satisfies the adiabatic condition, which says that the variation of the Hamiltonian should be slower than the dynamical time scale, in other words, the matrix element of ∂H/∂t should be smaller than the square of the energy gap [14,92]. In our simulation, the matrix element of ∂H/∂t is of the order of ∂g/∂t = g s /t s , which varies from 0.005 to 0.1.
In the weak coupling limit g → 0, the ground state is with all Z = −1, while the first excited state is one with a pair of visons, that is, a pair of plaquettes with Z = 1, which can be created by flipping the qubits on the links pierced by a string on the dual lattice. Thus the gap is 4. Its square is a lot larger than g s /t s . At g = 0, the ground state is 2 d -fold degenerate, with different eigenvalues ±1 of the d 't Hooft operators. They become nondegenerate when g = 0, and the one with eigenvalues of all V µ 's being 1 becomes the unique ground state. Hence there are small energy splittings between the ground state and other eigenstates of V µ . However, V µ 's are conserved because of gauge invariance, consequently these states belong to different topological sectors. Once the initial state is prepared in one of the topological sectors, it remains there when g is varied. Consequently, the small splittings between the ground state and other V µ eigenstates do not matter.
In the strong-coupling limit g → ∞, the ground state is with σ x = 1 for all qubits. The gauge invariance dictates that in the first excited state, there is a plaquette with σ x l = 1 on its four links. Hence the energy gap is 8g, which is very large. It is known that ground states and spectra in weak and strong coupling limits are all stable up the QPT point [57].
For a finite size L, the adiabatic condition is satisfied for general values of g, including the critical point as the energy gap is of the order of 1/L, L being the linear size of the system [14]. Hence the gap is about 0.5 for a general value of g.

Results of pseudoquantum simulation
Now we turn to the results of our actual pseudoquantum simulations of quantum Z 2 LGT on d=3 and on d=2 lattices. In each case, we first prepare the initial ground state at g = 0, then execute the adiabatic algorithm by varying g from 0 to 2 in substeps of t s /n. After each step of t s , several quantities are calculated. JHEP08(2020)160

Wegner-Wilson loops
As g increases, the ground state evolves, from the equal superposition of the configurations with all plaquettes in Z = −1, to be near the state with all qubits in σ x = 1. During this process, it undergoes a QPT, which can be characterized in terms of the Wegner-Wilson loop operator W C defined along a loop C on the direct lattice [57,60]. In the confined phase g > g c , W C obeys the area law W C ∝ g −A C = exp(−A C ln g), where A C is the area enclosed by the contour C. In the deconfined phase at g < g c , W C obeys the perimeter law W C = exp(−B(g)P C ), where P C is the perimeter of the contour C, B(g) is a smooth function of g and vanishes as g → 0. At g = 0, all Z = −1, the Z 2 flux is expelled, W C = 1. At small values of g, the fluctuations lead to the perimeter law. Our simulations confirm this picture. In our simulation, for d=3 and d=2 lattices respectively, we choose three contours denoted as c1, c2, c3, as shown in figure 6. On each lattice, the perimeter ratios are 1:1.5:2 while the area ratios are 1:2:3.
Features on Wegner-Wilson loops are shown in figure 7. As g → 0, W C → 1. With the increase of g, W C decrease relatively slowly when g is small, as indicated in subfigures (a) and (c), and the ratios between log W C 's of different loops equal the ratios between the perimeters, as indicated in subfigures (b) and (d). When g is relatively large, W C 's decrease as some powers of g, and the powers are shown to be areas, since the ratios between log W C 's for different loops equal the ratios between areas, as can be seen in subfigures (b) and (d). Features in d=3 and d=2 are similar, except that in d=3, there is a dip in the g-dependence of the log W C ratio, which will be discussed below.
Previous tensor network calculation for d=2 gave W C as functions of perimeter and area for a small and a large values of coupling constant, respectively [68]. Complementarily, here we give W C as a function of g, satisfying the area and perimeter laws for three coutours in d=2 and d=3 respectively.

Critical points and duality
In the adiabatic evolution, in steps of t s , we pause the evolution, and calculate Z and X , which are summed to give H , as shown in figure 8, from which we also obtain the first and second derivatives with respect to g, as shown in figure 9.
Our simulations are on very small lattices, therefore, the singularities at the critical points of QPT are rounded out. Nevertheless one can observe the critical points from the energy properties, where the transition is clearer than the transition of the Wegner-Wilson loops from the perimeter law to the area law.
We have determined the critical points g c 's using two methods. observed that which are not precise enough without finite-size scaling.
In the tensor network calculation on lattices with L ≤ 8 [68], g c is 0.3267 from the energy gap in the same topological sector of the ground state, and is 0.327 from the string operator expectation, and is 0.3285 from the overlap between the ground state and the lowest energy states in other topological sectors acted by Wegner-Wilson operators on the non-contractible loops.
(3.5) Therefore ψ(g)|Z|ψ(g) = ψ(1/g)|X|ψ(1/g) . In our simulation in d = 3, as the second approach determining g c , we find that the crossing point of Z and X is right at g c = 1.0 (figure 10), the same as the theoretical result.
As shown in figure 10, our simulation also demonstrates (3.6) for various values of g, further confirming self-duality in d=3. It also clearly shows the absence of self-duality in d=2.

Densities of states
In the ground state |ψ(g) for g equal to multiplies of g s , we numerically calculate the density of eigenstates {|z i } of Z and density of eigenstates {|x j } of and X, respectively. Z|z i = z i |z i , X|x j = x j |x j . From the decompositions |ψ(g) = i α i (g)|z i = j β j (g)|x j , it is known that the densities of states are D(g, z i ) = |α i (g)| 2 and D(g, x j ) = |β j (g)| 2 .
First consider d=3 2 × 2 × 2 lattice, with periodic boundary condition. As shown in figure 11, because of the geometric constraint, flipping one qubit between σ z eigenstates reverses the signs of the eigenvalues of Z 's of 4 plaquettes sharing this qubit, thus reverses the eigenvalue of Z by 8. Flipping the qubits on two crossing or neighboring parallel links reverses the signs of the eigenvalues of Z 's of 6 plaquettes, thus changes the eigenvalue of Z by 12. In general, 4 + 2n plaquettes can be flipped (n = 0, 1, 2 . . .). Consequently, the possible eigenvalues of Z are ±24, ±16, ±12, ±8, ±4, 0. The prominent feature here is that ±20 are forbidden, because at least 4 plaquettes are reversed.
As shown in figure 12 and figure 13, the above deduction is fully verified by our simulations, in which the possible eigenvalues z i 's of Z are indeed −24, −16, −12, · · · , 12, 16, 24. Also note that in the ground state, those Z eigenstates with large positive eigenvalues are difficult to occupy, as can be seen in figure 12 and figure 13. Interestingly, as can be seen in figure 12 and figure 13, the possible eigenvalues of X are the same as those of Z, that is, −24, −16, −12, · · · , 12, 16, 24, despite flipping the σ x eigenstate of one qubit changes the eigenvalue of X only by 2. This is a consequence of gauge invariance, which dictates that the qubits flipped in σ x eigenstates must be in closed loops, whose possible perimeters are 4, 6, 8, · · · , or 4 + 2n, (n = 0, 1, · · · ), thus X can only be changed by 8 + 4n. This is valid in any d ≥ 2. In d=3, the fact that Z 's of at least 4 plaquettes are reversed corresponds to the fact that the perimeter of a loop is at least 4.
This also confirms the self-duality in d=3, which implies that the possible eigenvalues of X must be the same as those of Z. Moreover, self-duality implies D(g, z) = D(1/g, x = z), which is also clearly confirmed in our simulation ( figure 12 and figure 13).
From DOS' of Z and X, it is calculated that with the increase of g from 0, the expectation value of Z increases from −24 towards 0, and while the expectation value of X decreases from 0 to −24 (figure 8).
In d=2, Z 's of two plaquettes are reversed by flipping the σ z eigenstate of one qubit. More generally, it is possible to create 2 visons by flipping σ z l 's of a string of qubits. Therefore the number of flipped plaquettes is 2n, (n = 1, 2, · · · ), as shown in figure 11. In d=2 3 × 3 lattice with periodic boundary condition, the possible eigenvalues of Z are −9, −5, · · · , 3, 7.
In d = 2, the possible eigenvalues of X are ±18, ±10, ±6, ±2, because the number of qubits flipped in σ x eigenstates must be 4 + 2n, that is, X is changed by 8 + 4n, (n = 0, 1, · · · ), as said above for any d ≥ 2. In consistency with the absence of self-duality, there is no identity D Z (z, g) = D X (x = z, 1/g), as can be seen in figure 12 and figure 13. As already shown in figure 10, with the increase of g from 0, the expectation value of JHEP08(2020)160 Z increases from −9 and towards 0, while the expectation value of X decreases from 0 towards −18.

Orders of quantum phase transitions
A first-order phase transition is one where there is a discontinuity of a first derivative of the free energy or energy, at the critical value of the control parameter, in contrast to a second-order phase transition, where the first derivatives are continuous while there is a discontinuity of its first derivatives at the critical point. The thermal phase transition in the classical Z 2 LGT is first-order in 4 spatial dimensions, and is second-order in 3 spatial dimensions [88,89]. This suggests that for quantum JHEP08(2020)160 LGT, the QPT is first-order in d=3 spatial dimensions, and is second-order in the d=2 spatial dimensions, because the classical theory in D spatial dimensions corresponds to the quantum theory in d=D-1 spatial dimensions, in other words, D=d+1 spacetime dimensions, where 1 represents the time dimension.
Limited by the smallness of the lattice size, one cannot make conclusion about whether there exists discontinuity in the slope of H in figure 8 or ∂ H /∂g in figure 9.  Figure 14. Comparison of results for d=3 2×2×2 and d=2 3×3 lattices, as functions of (g−g c )/g c , where g c is the respective critical value in each dimensionality. (a) Z , relative to the absolute value at g = 0. (b) Second derivative of the energy with respect to g. It can be seen clearly that in d=3, the change of Z at QPT is much steeper, and the valley of second derivative of the energy is much sharper. Therefore, we make a comparison between d=3 and d=2, by putting together Z and second order derivatives ∂ 2 H /∂g 2 for d=2 and d=3, as functions of (g − g c )/g c , which is dimensionless ( figure 14). Then it can be seen clearly that in d=3, the change of Z at QPT is much steeper and the valley is much sharper. This would be consistent with the claim that in d=3, QPT is first-order, even though the discontinuity and the singularity are rounded out, while in d=2, QPT is second-order. Now we examine the properties of the ground states. As can be seen in figure 13, in d=3, distributions of DOS of Z appear significantly different before and after QPT. For example, when g = 0.9 < g c , the peaks are at z = −24 and z = −16; when g = 1.1 > g c , the peak is at z = −8. When g = g c = 1, there are peaks at z = −24 and z = −16, i.e. the locations of peaks when g < g c , and at z = −8, i.e. the location of peak when g > g c . Same feature exists in DOS of X, by replacing g as 1/g because of self-duality.
This feature indicates a quantum version of phase coexistence at the critical point, a hallmark of first-order phase transition. When g < g c , the ground state slowly varies with g, as a same phase. When g > g c , the ground state is dramatically different from that for g < g c , as indicated by their distinct DOS distributions. The ground state for g > g c also slowly varies with g, as a same phase within this regime. When g = g c , the ground state is roughly a nearly-equal superposition of the two states. This feature is absent in d=2, also seen in figure 13.
Therefore, our simulations support the implication from 4-dimensional classical Z 2 LGT that QPT in quantum Z 2 LGT is first-order in d=3, while it is second-order in d=2. Now we go back to the expectations of Wegner-Wilson operators of different loops (figure 7). In d=3, there is a dip in the ratio of the logarithms, which is absent in d=2. We have examined that the dip is right at g c determined from the lowest point of the ∂ 2 H /∂g 2 in figure 14. It is not excluded that the dip is a finite-size effect. It is also possible that the dip is a signature of first-order QPT at g = g c , as a consequence of the quantum phase coexistence. The matrix element of the Wegner-Wilson loop operator between the two states representing the two phases possibly decrease the expectation in their superposition representing the phase coexistence. The dip in the ratio between c3 and c1 is deeper than that in the ratio between c2 and c1, in consistency with property that the larger the loop, the stronger the effect. This conjecture is supported by the absence of such a dip in d=2, where QPT is second-order, and also by the existence of a similar dip in a quantity studied in a quantum adiabatic algorithm for the exact cover problem, which was regarded as a criterion for first-order QPT [93]. Nevertheless, no definite conclusion can be drawn yet.

Topology
The phase transition in the classical Z 2 LGT is between states that cannot be distinguished in symmetry [21,22,55]. It is now called topological phase transition [60]. Our simulation verifies that QPT in quantum Z 2 LGT is also topological.
First, the absence of symmetry breaking is directly verified in σ z l of qubit l, as a function of g. As can be seen in figure 15, σ z l remains consistent with 0 for all values of g. For this qubit l, we have also studied the expectation values of x l ≡ −σ x l , z l = 1 4 ∋l Z , which is the average of Z 's of the plaquettes sharing the qubit l normalized by the number of plaquett per link, and h l = z l + x l .
Second, the above results on the DOS of X indirectly verifies the existence of visons ( figure 12). Two separated visons can only be annihilated by a nonlocal operator or by contacting each other.
Third, we have studied g-dependent energy splittings between the ground state and the other 2 d − 1 common eigenstates of the d 't Hooft loop operators V µ 's, µ = 1, · · · , d, JHEP08(2020)160 (figure 2). The presence of these d lowest energy states with exponentially small energy splittings, scaled as g L+1 and vanishing with L, is a defining characteristics of Z 2 topological order [60].
In d=3, at g = 0, there are 8 degenerate ground states with V x = ±1, V y = ±1 and V z = ±1. For a generic value of g, there are 4 energy levels −1, −1). The remaining degeneracy is due to rotational symmetry.
At g = 0, we prepare the ground state (V x , V y , V z ) = (1, 1, 1) and other three ones (1, −1, 1), (−1, −1, 1) and (−1, −1, −1), using the corresponding W µ operators, defined along non-contractible loops on the direct lattice. Other four ground states are not studied, because theoretically it is known that even when g = 0, each of them remain degenerate with one of the states studied. Then we execute adiabatic algorithm to obtain the dependence of the energies and the splittings on g before and after the QPT, as shown in figure 16. In d=2, at g = 0, we first prepare the ground state with V x = V y = 1, and then use W µ operators to obtain the other 3 degenerate ground state with V x = ±1 and V y = ±1. The quantum adiabatic algorithm is executed on each of them. They become nondegenerate when g = 0. In our simulation, we calculate the four energies E(V x , V y ), which are E 1 = E(1, , 1), E 2 = E(−1, 1) and E 3 = E(1, −1), E 4 = E(−1, −1). The simulation confirms that E 1 is the lowest while E 2 = E 3 . As shown in the inset of in (d) of figure 16, the splitting E i1 ≡ E i − E 1 can be fitted by least-square method as with L = 3 fixed. Previous tensor network calculation found that the splittings between the four eigenstates of V x and V y exponentially decay with L for given small g [68]. This is consistent with and complements our result, as g L+1 = ge −L/ξ , with ξ ≡ −1/ ln g. Our result reveals the g dependence but does not give the L dependence, while their result gives the L dependence but does not give the g dependence. In our algorithm, we have generalized Trotter decomposition and symmetrized Trotter decomposition such that the Hamiltonian adiabatically varies during the decomposition. Hence we have proposed a scheme of digital adiabatic quantum simulation. Hopefully, this method can be used for various problems.
We have studied quantum Z 2 LGT in lattices of spatial dimensions d=2 and d=3. Gauge invariance is realized in the initial state and is preserved during the adiabatic evolution. Since the lattices are very small, singularities in QPT are rounded out. But key features are observed. It is indicated that QPT is first-order in d=3 and is second-order in d=2, with critical point consistent with previously known results. In d=3 and d=2 respectively, we also clearly observe topological characteristics of QPT, including the vanishing of JHEP08(2020)160 the "magnetization" for all values of the coupling constant, the excitation properties, the lowest energies in different topological sectors and the their splittings, which are proportional to g L+1 in d=2. We have observed the change of the degeneracy caused by nonzero g, which is also an indication of topological order. To our knowledge, we have presented the unique numerical result on quantum Z 2 LGT in 3 spatial dimensions.
This work seems to be the first complete demonstration as a proof of principle, albeit in a classical simulator, carrying through quantum simulation of a LGT and observe the key features of physics. Thereby it demonstrates that real quantum simulations of LGTs in future can be done in the proposed way.
On the other hand, this work also shows that high-performance classical demonstration of quantum simulation, which may be dubbed pseudoquantum simulation, represents a new way of computation, in addition to facilitating the development of quantum software for experimental quantum simulation.
As shown in the comparison with previous results from tensor network calculation, it is a basic element of our adiabatic approach that the parameter in the Hamiltonian is varied, so it is very natural and convenient to obtain various quantities as functions of this parameter, which may not be convenient in other methods, hence adiabatic quantum simulation and pseudoquantum simulation are convenient tools of QPT study.
The classical demonstration proceeds according to rules of quantum mechanics, therefore pseudoquantum simulation is legitimately a reliable approach, as far as computational resources allow. Compared with other computational approaches, it can be used without intricate algorithmic design depending on the details of the computed problem, as in many other numerical methods, which are often applicable only to one or two dimensions.
As the next step regarding quantum and pseudoquantum simulations of LGT, we shall study Fermions coupled with the Z 2 gauge field using our method, on which the results can be compared with the existing results of QMC calculations, which are free of Fermion sign problem. Thereby the method can be further benchmarked, which can then be applied to other LGTs, for which Fermion sign problem exists in MC-based methods, as well as the LGTs that recently have been studied by using tensor network methods [69][70][71][72].
After the release of the present work as a preprint (arXiv:1910.08020), there appeared more recent progress in the related fields, including DMRG study of the one-dimensional spinless Fermions coupled with Z 2 gauge field [94], schemes of measuring nonlocal observables in the quantum simulation of a general LGT [95], experimental observation of gauge invariance in a 71-site quantum simulator of an extended U (1) LGT [96], and a sign-free QMC study of Z 2 gauge field coupled with both Fermions and Bosons, demonstrating the transition from conventional metal to orthogonal metal [97].