Circuit-based digital adiabatic quantum simulation and pseudoquantum simulation as new approaches to lattice gauge theory
aa r X i v : . [ qu a n t - ph ] S e p Prepared for submission to JHEP
Circuit-based digital adiabatic quantum simulationand pseudoquantum simulation as new approachesto lattice gauge theory
Xiaopeng Cui a Yu Shi a, Ji-Chong Yang a,b a Department of Physics & State Key Laboratory of Surface Physics, Fudan University, Shanghai200433, China b Department of Physics, Liaoning Normal University, Dalian 116029, China
E-mail: , [email protected] , [email protected] Abstract:
Gauge theory is the framework of the Standard Model of particle physicsand is also important in condensed matter physics. As its major non-perturbative ap-proach, lattice gauge theory is traditionally implemented using Monte Carlo simulation,consequently it usually suffers such problems as the Fermion sign problem and the lack ofreal-time dynamics. Hopefully they can be avoided by using quantum simulation, whichsimulates quantum systems by using controllable true quantum processes. The field ofquantum simulation is under rapid development. Here we present a circuit-based digitalscheme of quantum simulation of quantum Z lattice gauge theory in and di-mensions, using quantum adiabatic algorithms implemented in terms of universal quantumgates. Our algorithm generalizes the Trotter and symmetric decompositions to the case thatthe Hamiltonian varies at each step in the decomposition. Furthermore, we carry through acomplete demonstration of this scheme in classical GPU simulator, and obtain key featuresof quantum Z lattice gauge theory, including quantum phase transitions, topological prop-erties, gauge invariance and duality. Hereby dubbed pseudoquantum simulation, classicaldemonstration of quantum simulation in state-of-art fast computers not only facilitates thedevelopment of schemes and algorithms of real quantum simulation, but also represents anew approach of practical computation. J. High Energ. Phys. , 160 (2020). https://doi.org/10.1007/JHEP08(2020)160 Corresponding author. ontents Z lattice gauge theory 52.2 Circuit-based quantum adiabatic algorithm 72.3 Preparation of the initial ground state 112.4 Measurement of physical quantities 132.5 Computational hardware 142.6 Adiabaticity and parameter values 14 Quantum simulation and quantum computation can efficiently solve some problems thatcannot be efficiently solved in classical computers [1–4], and is under extensive studiesworldwide, thanks to the rapid development of quantum science and technology. A control-lable quantum system, which may even be universal or programmable, simulates variousquantum systems, whose physical properties can be conveniently investigated with variousparameter values. Even with only tens of qubits, far less than those in full fault-tolerantquantum computing, and even in the presence of some noises, a quantum machine canperform some tasks surpassing its classical counterparts, exhibiting the so-called quantumsupremacy [5, 6]. Many quantum simulations are such tasks. Hopefully they not only cansolve specific problems, but also represent new scientific methods [3].Quantum simulations of important models in theoretical physics are enabled by somequantum algorithms, including the Trotter decomposition [2], the sparse Hamiltonian quan-tum walk [7], the dense Hamiltonian density matrix exponentiation [8–10], the adiabaticalgorithm [11, 12], and so on. It is timely even to develop various quantum softwares [13].An important battlefield of quantum simulation appears to be lattice gauge theory(LGT) [14], the major non-perpurbative approach to gauge theory. In particle physics,– 1 –auge theory is the framework of the Standard Model, describing both electroweak andstrong interactions among elementary particles, and is also a guide beyond the StandardModel. A LGT is a gauge theory defined on a spacetime lattice in path integral formalismor on a space lattice in Hamiltonian formalism, making the degrees of freedom countableand convenient for numerical calculations. Gauge theory is also important in condensedmatter physics, where the lattice can be a real structure. It is often an effective descriptionof constraints in strongly correlated systems, and uses emergent gauge fields to characterizetopological orders, which exist in fractional quantum Hall effect, spin liquids and possiblyin high temperature superconductivity, as well as in topological quantum computing, etc.Topological order represents a new and active paradigm beyond the traditional frameworksof symmetry breaking and Fermi liquid theory.Using Monte Carlo (MC) simulation, LGT has made great achievements [15–17]. Butthere are also difficulties, including the lack of real-time dynamics because of the use ofEuclidean spacetime, and the notorious Fermion sign problem [18], which was shown tobe NP hard [19]. The cause of the Fermion sign problem is that in presence of Fermions,the Boltzmann weight of an appropriately defined configuration, to which the quantumproblem is mapped, may become negative, therefore the partition function oscillates vio-lently, consequently the importance sampling in MC becomes invalid. These difficulties arerelated to the unsolved problems in quantum chromodynamics, such as color confinementand phase diagram of quark-gluon plasma. The Fermion sign problem exists in most ofthe MC-based methods such as quantum MC (QMC), except special algorithms for somespecific models [20], and in some special issues to be told below.As a new approach avoiding Fermion sign problem and an ideal avenue to study quan-tum phase transition (QPT) and real-time quantum dynamics, quantum simulations ofLGTs 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 [21–26]. Otherproposals were also made [27, 28]. Among these schemes, some are analog, based on phys-ical mapping between the Hamiltonians of simulated and simulating systems, while othersare digital, based on Trotter decomposition of the finite-time evolution to many small stepsthat are much easier to be implemented. A quantum-classical algorithm was developed fortwo-site Schwinger model [29].Experimentally, quantum simulations of (1+1)-dimensional quantum electrodynamics(QED) or U(1) theory were performed using trapped ions [30] and cold atoms [31]. Initialattempts were also made in quantum simulation of quantum Z LGT using cold atoms [32].A practical proposal using trapped ions was also made on quantum simulations of QED,Chern-Simons theory and Z theory [33].Quantum Z LGT is the simplest quantum LGT [16, 34–36]. On one hand, Z n theorycan be obtained as the discretization of U(1) theory, suitable for quantum simulations [37].On the other hand, quantum Z LGT is also important in condensed matter physics [36, 38]. Z toric code model [39], which is important in topological quantum computing, can beregarded as a variant of quantum Z LGT, and has been experimentally demonstrated ina quantum simulation using nuclear magnetic resonance [40, 41].In two spatial dimensions, quantum Z LGT is dual to quantum Ising model in a– 2 –ransverse field [16, 34–36], which is thus often invoked for QPT properties, especiallythe critical point [42–45]. Besides, Z n theory in one spatial dimension was studied byusing density matrix renormalization group (DMRG) [37], and was studied in two spatialdimensions by using tensor network techniques, as the low energy limit of toric modelin a magnetic field under the constraint of gauge invariance [46]. Direct DMRG studyof quantum Z LGT is difficult as realizing plaquette interactions consistently with thegauge 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 [47].Interestingly, the duality between such a spin gauge theory and a generalized Ising modelallows for scalable quantum simulation with Rydberg atoms [48].Coupling of Z gauge field with various kinds of matter has also been studied, start-ing with Ising matter [35]. Fractionalization of electrons is obtained in theories of strongcorrelations with Z gauge fields, giving rise to the so-called orthogonal metals [49]. In-terestingly, the issue of Fermions coupled with Z gauge field is exactly one of the specialissues free of sign problem in MC-based methods [50], another issue being Fermions withan even number of flavors [51]. The cases with both of these characteristics were studied byusing QMC [52]. Coupling of Fermions with Z gauge field with Gauss law not imposed butemerged was also studied by using QMC [53]. These QMC studies were all in two spatialdimensions. A modified Z gauge theory coupled with Fermions was studied analyticallyin one and two spatial dimensions [54]. One-dimensional Bosons coupled with Z gaugefield with Gauss law not imposed but emerged was studied also by using DMRG [55]. Aquantum link model of QED was approached by using tensor network method [56].In quantum Z LGT, defined on a square lattice, each link is occupied by one qubit. Inthe × × lattice, which is the smallest three-dimensional lattice, there are 24 links (Fig. 1).In the × lattice, the second smallest two-dimensional lattice, there are 18 links (Fig. 1).In a quantum algorithm, additional qubits may also be needed as ancillas in simulatingthe 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 nonlocalinteractions, for example. Regarding the quantum adiabatic algorithm, tens of thousandssteps are needed in Trotter decomposition. Therefore, at present time, it seems difficult forthe experimental platforms to fully meet the requirements on the qubit number, error rateand coherence time.Therefore, it appears interesting to use classical high-performance computing platformsto demonstrate quantum simulation in general, and that of quantum Z LGT in particu-lar. We call it pseudoquantum simulation, which serves as a benchmark for real quantumsimulation, and facilitates the development of quantum algorithm and quantum softwares.Meanwhile, it is also a new method of computation and simulation, providing useful resultson the computed problems.Such pseudoquantum simulation is realistic if it is run on a fast enough classical comput-ing platform. An example is the latest graphics processing unit (GPU) parallel computingarchitecture, 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 [57]. It is based on the software– 3 – a) (b)
Figure 1 . (a) Three-dimensional × × lattice with periodic boundary condition. (b) two-dimensional × lattice with periodic boundary condition. There is one qubit on each link, withan ordering numbers indicated. In (a), for example, qubits numbered 1, 3, 4 and 15 form anelementary plaquette, while qubits numbered 4, 6, 1, and 18 form another, and there are qubits.In (b), for example, qubits numbered 3, 4, 5 and 10 form an elementary plaquette, while qubitsnumbered 5, 6, 1, and 12 form another, and there are 18 qubits. platform and application programming interface CUDA created by Nvidia, which allowsthe development of parallel program using a CUDA-enabled GPU. QuEST is designed asa 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 de-tailed below. Therefore it fits well the need of our pseudoquantum simulation of quantum Z LGT.In this paper, we present a scheme of circuit-based digital quantum simulation of (2+1)-dimensional and (3+1)-dimensional quantum Z LGT. We use universal quantum circuit toimplement the quantum adiabatic algorithm, and we generalize the Trotter decompositionand symmetrized Trotter decomposition to the case of adiabatic varying the Hamiltonian ineach step of decomposition. As the first work and proof of principle, here we consider puregauge theory, without coupling to Fermions yet. Moreover, we perform a complete pseudo-quantum 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.682TFLOPSand 7.834TFLOPS, with double precisions, respectively. Our work demonstrates the ad-vantages of quantum simulation as well as the usefulness of pseudoquantum simulation.Meanwhile we obtain useful results regarding various properties of quantum Z LGT. Itseems to contain the first numerical result on quantum Z LGT in three spatial dimensions,including the features of first-order QPT, as hinted by the fact that the thermal phasetransition in the classical Z LGT is first-order in 4 spatial dimensions [58].The rest of the paper is organized as the following. In Sec. 2, we introduce the quantum Z LGT, the lattices we consider, the quantum adiabatic algorithms and their realizationin terms of the quantum circuits, our generalization of the Trotter and symmetrized Trot-– 4 –er decompositions, the simulation of measurement, as well as the hardware we use todemonstrate the scheme of quantum simulation. In Sec. 3, we present the results fromour GPU computation simulating the quantum simulation, including the expectations ofWegner-Wilson loop operators and of the Hamiltonian, which are used to determine thecritical points and orders of the QPTs, and confirm the self-duality in 3+1 dimensions. Wealso calculate the densities of states, which confirm the gauge invariance, the self-duality in3+1 dimensions and its absence in 2+1 dimensions. We will also present the evidences oftopological nature of QPT. Finally, a summary is made in Sec. 4. Z lattice gauge theory Z 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 spontaneousmagnetization at the phase transition, which is characterized by onset of topological orderrather than symmetry breaking. So there is no local order parameter [34]. On the otherhand, 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 gaugepotential and thus electric field on the links between sites.Here we focus on the quantum Z LGT with the Hamiltonian [36, 38] H = Z + gX, (2.1)where g is the coupling constant, X ≡ − X l σ xl , Z ≡ X (cid:3) Z (cid:3) , Z (cid:3) ≡ − Y l ∈ (cid:3) σ zl , (2.2)with Z (cid:3) defined for each elementary plaquette (the smallest square) (cid:3) (Fig. 1). For qubit l , σ zl | i l = | i l , σ zl | i l = −| i l . The Hamiltonian of the classical Z theory is Z only, witheach operator σ zl reduced to a classical variable. The quantum nature of H is due to thenoncommutativity between σ xl and σ zl , resulting the competition between Z and X , in away like that between energy and entropy in a thermal phase transition. The couplingconstant g is a control parameter, playing a role in QPT similar to the role of temperaturein thermal phase transition.For convenience, in our numerics, the links and thus the qubits are numbered as inFig. 1. For three spatial dimensions d=3, we use × × lattice with periodic boundarycondition, where there are 24 links and 24 elementary plaquettes. We use 25 qubits, oneof which is the ancilla. For two spatial dimensions d=2, we use × lattice with periodicboundary condition, where there are 18 links and 9 elementary plaquettes. We use 19qubits, one of which is the ancilla. We assume periodic boundary condition, that is, eachlattice is a torus. Although the lattice sizes are very small, the key features of the quantum Z LGT do appear. In general, for a d -dimensional square lattice with linear size L , thenumber of links is N l = dL d , the number of plaquettes is N p = N l ( d − / .– 5 –his theory possesses Z gauge invariance, similar to Gauss law, dictating that eacheigenstate | ψ i of H must satisfy G i | ψ i = | ψ i (2.3)where G i ≡ Q l ∋ i σ xl is the product of the σ xl ’s on all the links ending at each lattice point i . One gauge invariant operator is Wegner-Wilson loop operator, which is defined as W C = Y l ∈ C σ zl (2.4)along a closed loop C on the direct lattice. G i W C G − i = W C . A Wegner-Wilson loopoperator is not necessarily along a non-contractible loop, for example, Z (cid:3) is a Wegner-Wilson loop. But those along non-contractible loops play special roles. Generalizing thewell known case in d=2, here we define the special Wegner-Wilson loop operators W µ = Y l ∈ C µ σ zl , (2.5)along non-contractible loop C µ , µ = 1 , · · · , d (Fig. 2). Periodic boundary conditions meansthat topologically they encircle a torus. As far as it encircles the lattice in µ direction, thedetails of C µ does not matter. σ x and any product of σ x ’s are also gauge invariant. G i σ x G − i = σ x . Generalizing thewell known case in d=2, here we make a definition of the so-called ’t Hooft loop operatorin any dimention d, V µ ≡ Y l ∈ C µ σ xl , (2.6)which is a product of σ xl pierced by a non-contractible ( d − -dimensional surface C µ onthe dual lattice, µ = 1 , · · · , d (Fig. 2). A non-contractible ( d − -dimensional surface canbe deformed by using G i operators, because of gauge invariance. C µ is subscripted as µ asit can be deformed, if needed, to be parallel to C µ . Periodic boundary conditions meansthat topologically they encircle a torus. In d=2, C µ is a loop. V µ commutes with H , while W µ does not unless g = 0 . For µ = ν , W µ V ν = − V ν W µ . (2.7)Consequently W µ acting on an eigenstate of V ν ( ν = µ ) yields another eigenstate of V ν , withthe eigenvalue of the opposite sign. Therefore, from the state with | V ν = 1 , ν = 1 , · · · , d i ,using consecutive actions of W µ ’s, one can generate | V ν = − , V d − ν = 1 , ν = 1 , · · · m i = Q mν =1 W ν | V ν = 1 , ν = 1 , · · · , d i , m = 1 , · · · , d , hence P dm =1 C mn = 2 d − states in totalcan be generated. These d common eigenstates | V , · · · , V d i of V µ with eigenvalues ± ,( µ = 1 , · · · , d ), are also eigenstates of H . They are degenerate ground states when g = 0 .The degeneracy on a d-dimensional lattice with genus G is d G . On the lattice consideredhere with period boundary condition, the degeneracy is d . The degeneracy changes when g = 0 , with only | V ν = 1 , ν = 1 , · · · , d i remaining as the ground state. But they represent– 6 – a) (b) C x C y C x C y Figure 2 . (a) d=3 Lattice, with the indication of non-contractible loops C x , C y and C z on thedirect lattice, and non-contractible C y on the dual lattice. C x , C y and C z are used in definingspecific Wegner-Wilson operators W x , W y and W z . C x and C z are similar to C y . They are usedin 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 asnon-contractible loops C x and C y in defining ’t Hooft operators V x and V y . different topological sectors, because V µ is conserved while each non-contractible ( d − -dimensional surface can deform.Depending on g , there are two different phases, deconfined phase at small g and confinedphase at large g , separated by a QPT at the critical point g c , where there is a phasetransition, which however cannot be characterized by a change of symmetry. There is Z topological order in the deconfined phase, while the confined phase is trivial, and QPT inthis theory is topological [34, 36, 38]. Our digital scheme of quantum simulation of the quantum Z LGT is based on the quantumadiabatic algorithm, and we use quantum circuits consisting of one-qubit and two-qubitquantum gates to implement it. In the adiabatic evolution, g varies from to a large enoughvalue g f , passing g c . The slow variation allows the system to adapt to the instantaneousground state. According to the adiabatic theorem, the state of the system, starting as aground 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 ) .In our simulation, each time g is updated, the state evolves for a very short time, asrealized by the quantum circuit under the present g value, then g is updated again, and the– 7 –tate 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. Thisso-called “substep” really corresponds to the “step” in Trotter decomposition. The reasonof referring to the decomposition steps as substeps is that the evolution is paused after aperiod called step, and calculations, or called pseudo-measurements, are done, afterwardsthe evolution is resumed. In real quantum simulation, the measurement depends on theactual situation.Therefore, the stepwise Hamiltonian, for the m -th substep within k -th step, is H k,m = Z + g k,m X, (2.8)where g k,m = ( k − g s + mδ, (2.9) ( m = 1 , · · · , n ) , n is the total number of substeps in each step, which is freely set. g s is theincrease 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 Y m =1 e − iH k,m tsn , (2.10)while the total evolution is N s Y k =1 n Y m =1 e − iH k,m tsn . (2.11)The evolution in each substep of t s /n , under H k,m , consisting of two noncommutativeparts Z and g k,m X , can be decomposed, in an asymmetric way, into consecutive evolutionof Z and g k,m X , e − iH k,m tsn ≈ e − iZ tsn e − ig k,m X tsn . (2.12)Therefore the evolution in k -th step is n Y m =1 (cid:16) e − iZ tsn e − ig k,m X tsn (cid:17) . (2.13)which generalizes Trotter decomposition. It reduces to Trotter decomposition if H k,m isindependent of m .From the identity e A + B = e A e B e − [ A,B ]+ ··· for two operators A and B , we have e − i ( H + H ) τ = e − iH τ e − iH τ e τ [ H ,H ]+ ··· , thus || e − i ( H + H ) τ − e − iH τ e − iH τ || ≈ τ || [ H , H ] || . (2.14)Therefore, the error for one substep, in the asymmetric decomposition (2.12), is ∆ asyk,m ≡ || e − iH k,m tsn − e − iZ tsn e − ig k,m X tsn || ≈ g k,m N p t s n , (2.15)– 8 –here we have considered that each Z (cid:3) is noncommutative with 4 σ x ’s. It is calculatedthat N s X k =1 n X m =1 g k,m = 12 ( N s n + N s ) g s ≈ N s ng s / . (2.16)Therefore the total error of the asymmetric decomposition is ∆ asy = N s X k =1 n X m =1 ∆ asyk,m = ( N s + N s n ) N p g s t s n ≈ N s N p g s t s n . (2.17)We have also used the symmetric decomposition e − iH k,m tsn ≈ e − iZ ts n e − ig k,m X tsn e − iZ ts n . (2.18)Therefore the evolution in k -th step is n Y m =1 (cid:16) e − iZ ts n e − ig k,m X tsn e − iZ ts n (cid:17) , (2.19)which generalizes the symmetrized Trotter decompositon. It can be rewritten as e − iZ ts n n − Y m =1 (cid:16) e − ig k,m X tsn e − iZ tsn (cid:17) e − ig k,n X tsn e − iZ ts n , (2.20)which indicates a more convenient way of execution in our computation.Using ln( e A/ e B e A/ ) = A + B − ([ A, [ A, B ]] + 2[ B, [ A, B ]]) /
24 + · · · [59], we find e − iH k,m tsn − e − iZ ts n e − ig k,m X tsn e − iZ ts n ≈ i t s n (cid:0) g k,m [ Z, [ Z, X ]] + 2 g k,m [ X, [ Z, X ]] (cid:1) , (2.21)where [ Z, [ Z, X ]] is of the order of d − N p , while [ X, [ Z, X ]] is of the order of N p ,for the following reason. Each Z (cid:3) is noncommutative with 4 σ x ’s, hence [ Z (cid:3) , X ] is thesum of 4 products of one σ y and 3 σ z ’s. Each σ y is noncommutative with the d − Z (cid:3) ’s of the plaquettes sharing with the link l . On the other hand, each product of one σ y and 3 σ z ’s is noncommutative with 4 σ x ’s. Therefore [ Z, [ Z, X ]] = O [8( d − N p ] , while [ X, [ Z, X ]] = O (16 N p ) . Another way of reasoning is the following. Each σ x is shared by d − plaquettes, thus [ Z, σ xl ] yields d − products of σ y and 3 σ z ’s. Each productis noncommutative with Z (cid:3) ’s of the d − plaquettes, and with the 4 σ x ’s on the sameplaquette. Consequently, [ Z, [ Z, X ]] = O [4( d − N l ] , [ X, [ Z, X ]] = O [8( d − N l ] . With N p = N l ( d − / , this is the same as above.Therefore the error, in the symmetric decomposition (2.18), is ∆ symk,m ≡ || e − iH k,m tsn − e − iZ ts n e − ig k,m X tsn e − iZ ts n || ≈ ( d − g k,m + 4 g k,m N p t s n . (2.22)It is calculated that N s X k =1 n X m =1 g k,m = ( N s − N s (2 N s − ng s N s n ( n + 1)(2 n + 1) g s n + N s ( N s − n + 1) g s ≈ N s ng s . (2.23)– 9 –herefore the total error of the symmetric decomposition is ∆ sym = N s X k =1 n X m =1 ∆ symk,m = [ ( d − N s g s + 49 N s g s ] N p t s n ≈ N s N p g s t s n . (2.24)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 oneby a factor of t s /n . We have done our simulations using both decompositions. The resultsfrom the symmetrized decomposition is clearly better and thus presented below.Note that our decompositions are different from, and generalize, the usual Trotterdecomposition and symmetrized Trotter decomposition, each of which repeats a constantevolution for a number of times.We now discuss how to implement the evolution in each substep, or called each decom-position step. First, e − iZt s /n = Q (cid:3) e − iZ (cid:3) t s /n , hence the evolution of Z can be realized bythe consecutive evolution of all the plaquettes.Evolution of e − iZ (cid:3) t s /n , for each palquette (cid:3) , is realized in terms of a quantum circuit,where there is also an ancilla, shown in Fig. 3, e − iZ (cid:3) tsn = A − R az ( − t s n ) A, (2.25)with A = Y l ∈ (cid:3) CN OT l,a , (2.26)where CN OT l,a is a controlled-NOT gate controlled by the qubit l and targeting on theancilla a , A − is the product of these CNOT gates in reversed order, R az ( φ ) ≡ e − iσ az φ/ issingle-qubit gate on ancilla representing rotation of angle φ around z-axis. Initially, theancilla is set to be | r i = | i . Preceding R az ( − t s /n ) , each CN OT l,a flips | r i if and only ifthe control qubit l is | i . Therefore R az ( − t s /n ) acts as e it s /n when there are even numberof | i ’s on the plaquette, and acts as e − it s /n when there are odd number of | i ’s on theplaquette. This is precisely the effect of e − iZ (cid:3) t s /n . Afterwards, the four CNOT gates after R az ( − t s /n ) return | r i to | i , which can be used for the next plaquette. So we only needone ancilla. If we set | r i = | i initially, this circuit can be used to simulate the reversedevolution e iZ (cid:3) t s /n , in other words e − iZ (cid:3) t for reversed time t = − t s /n .It is straightforward to realize the evolution under the other part gX in the Hamiltonian, e − igX tsn = Y l e − i ( − σ xl ) g tsn = Y l R lx ( − g t s n ) , (2.27)where R lx ( φ ) on qubit l represents rotation of angle φ around the x-axis. For m -th substepof k -th step, g = g k,m .Our realization of e − iZ (cid:3) t s /n , as given in (2.25) and Fig. 3, is a direct application ofthe standard strategy for the evolution under an interaction that is a tensor product of σ z operators [60]. In some previous schemes of quantum simulation of LGTs, the four-bodyinteractions are obtained stroboscopically through a sequence of two-body interactions with– 10 – igure 3 . The quantum circuit realizing e − iZ (cid:3) t s /n or e iZ (cid:3) t s /n , depending on whether the initialstate of the ancila | r i is set to be | i or | i . ancillary degrees of freedom, and gauge invariance in each step of Trotter decomposition isemphasized [23, 25].Here, as Z (cid:3) and σ xl are gauge invariant operators, the evolution operators are gaugeinvariant in every substep of the digital decompositions. The adiabatic quantum simulation starts with g = 0 , i.e. H = Z , of which there aredegenerate ground states satisfying Z (cid:3) = − , ∀ (cid:3) , (2.28)which can be prepared by using the quantum circuit in Fig. 4.Each qubit is initially in | i and is then transformed to be ( | i + | i ) / √ by using aHadamard gate H . Therefore the state of the system is in the equal superposition of allbasis states, | ψ i = 1 √ N l X i =0 · · · X i Nl =0 | i · · · i N l i . (2.29)Each qubit is not entangled any other qubit, each plaquette is also in the equal superpositionof all its basis states.Then the four CNOT gates between the four qubits of one plaquette and the ancillainitially in | i a produce the state √ ( | r = 0 i a | Z (cid:3) = − i + | r = 1 i a | Z (cid:3) = 1 i ) , where | Z (cid:3) = ± i are states of all the qubits on the lattice satisfying Z (cid:3) | Z (cid:3) = ± i = ±| Z (cid:3) = ± i .The ancilla is entangled with the qubits on the lattice, with σ az ≡ − r = − Z (cid:3) in eachbranch. Then the measurement operation M on the ancilla projects it to | r = 0 i , therebyselects the state of qubits on the lattice to be | Z (cid:3) = − i .Afterwards, the ancilla is returned to | i a by other four CNOT gates (Fig. 4). Thesame ancilla is ready to work on another plaquette, which may or may not share a qubitwith a plaquette already worked on. Similar procedure goes on, till all plaquettes have been– 11 – igure 4 . The quantum circuit preparing the initial ground state of Z (cid:3) for one plaquette. worked on. The state of the system is then an eigenstate of Z (cid:3) for each (cid:3) , with eigenvalue − , and is an equal superposition of all configurations satisfying Z (cid:3) = − for each (cid:3) .To summarize, we perform consecutive projections Q (cid:3) P (cid:3) , where P (cid:3) = | Z (cid:3) = − ih Z (cid:3) = − | . Then the state Q (cid:3) P (cid:3) | ψ i is exactly a ground state of Z , and is the eigenstate of the’t Hooft operators V µ with eigenvalue , for all µ ’s. It satisfies the gauge invariance. Thisground state can adiabatically evolve to the ground state for g = 0 .The evolution of each decomposition substep is also gauge invariant, so the state pre-serves gauge invariance during the evolution. For a gauge invariant state | ψ i , satisfying G i | ψ i = | ψ i , and the evolution of a time period τ under a gauge invariant operator O satisfying G i OG − i = O , one has G i ( e − iOτ | ψ i ) = G i e − iOτ G − i G i | ψ i = e − iOτ | ψ i . Therefore,with the initial state gauge invariant, gauge invariance is always preserved in each evolutionunder Z (cid:3) or σ xl in the digital decompositions.The other degenerate ground states at g = 0 can be obtained by using W µ ’s as describedabove. This action changes the topological sector. When g = 0 , the state adiabaticallyevolves 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 redefinethe Hamiltonian as X + KZ , with the coupling constant K corresponding to /g . Theground state of X is with σ x = 1 for all qubits. With the increase of K , the groundstate evolves from confined phase to deconfined phase. During the adiabatic evolution, theground state remains in the topological sector of V µ = 1 for all µ ’s. To enter other sectorsand adiabatically approaches the other degenerate ground states of H ( g = 0) , one can alsouse 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 thecollapse 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, makethe success rate only − N p . Hence the second approach is preferred. We have actually alsotried it in our demonstration, which yields result consistent with the first approach, so isomitted here, as the emphasis is on the deconfined phase.– 12 – .4 Measurement of physical quantities The energy E in state | ψ i is the expectation value of H , h H i = h Z i + g h X i , where h Z i = h ψ | Z | ψ i = P i z i P ( z i ) , h X i = h ψ | X | ψ i = P i x i P ( x i ) , where { z i } and { x i } represent theeigenvalues or measurement results of Z and X , respectively, P ( z i ) and P ( x i ) are thecorresponding probability distributions, called densities of states hereby. In our simulator,we obtain E by summing up h Z i and g h X i . We also calculate the expectation values ofWegner-Wilson loop operators.We use CUDA parallel acceleration method to count the statistical summation of allbasis vectors on GPU. We write our own codes for the calculation of the measurementresults, 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 ) } inthe representation { σ zl } , as h W C i = h ψ | Q l ∈ C σ zl | ψ i , h Z i = − P (cid:3) h ψ | Q l ∈ (cid:3) σ zl | ψ i , h X i = − P l h ψ | H σ zl H | ψ i , where H is Hadamard gate.We mention that in real quantum simulation, E can also be measured by using quantumphase estimation. For an eigenstate | u i of a time-independent H , the evolution for a timeduration t is e − iHt | u i = e − iφ | u i , where φ ≡ Et . e − iHt can be realized in a way similar to the decompositions described in Sec. 2.2, butnow g is fixed. That is, e − iHt ≈ ( e − iZ tn e − igX tn ) n , (2.30)or e − iHt ≈ ( e − iZ t n e − igX tn e − iZ t n ) n = e − iZ t n ( e − igX tn e − iZ tn ) n − e − igX tn e − iZ t n , (2.31)where n is the number of decomposition steps here. These are Trotter decomposition andsymmetrized Trotter decomposition. The errors are just n times those for t/n , given inEq. (2.15) and (2.21) with g k,m replaced as g , that is, gN p t n and || i t n (cid:0) g [ Z, [ Z, X ]] + 2 g [ X, [ Z, X ]] (cid:1) || ≈ ( d − g + 4 g N p t n , respectively.For the purpose of quantum phase estimation, we need the conditional evolution con-trolled by an ancilla, U ( t ) = | i h | ⊗ e − iHt + | i h | ⊗ e iHt , (2.32)which can be realized in terms of controlled gates, and ancilla controlling the time direction.For Z part, the method is as shown in Fig. 3. For X part, the method is as shown in Fig. 5.Then U ( t ) √ ( | i + | i ) | u i = √ ( e − iφ | i + | i e iφ ) | u i . Therefore, one can prepare U (2 k t ) √ ( | i + | i ) | u i = e − i k φ √ ( | i + e i k · φ | i ) | u i , ( k = 0 , , , ... ) , as required by thealgorithm of quantum phase estimation. Subsequently, the probability distribution of φ ,and thus E , can be obtained by using the standard procedure of quantum phase estimation.– 13 – igure 5 . The quantum circuit for Q l ∈ (cid:3) e iσ xl gt s /n or Q l ∈ (cid:3) e − iσ xl gt s /n , depending on whether thestate of the ancila is | i or | i . GPU card Simulation Precision capacity (TFLOPS) N f N qubit M K40m float 5.046 4 30 8.2GBK40m double 1.682 8 29 8.2GBV100 float 15.67 4 31 16.2GBV100 double 7.834 8 30 16.2GB
Table 1 . Maximal scales that QuEST can perform with one Nvidia Tesla K40m-12GB or NvidiaTesla V100-SXM2-32GB GPU card. N f is the byte number of a floating point number, N qubit isthe number of qubits, M is the memory requirement. We use a Nvidia Tesla K40m GPU card, which was used by QuEST team in their simulationof 29 qubits with float decision [57], as well as a Nvidia Tesla V100-SXM2-32GB GPU card.We estimated the maximal scales of quantum simulations that QuEST can simulate underdifferent precisions, as listed in Table.1. We use double precisions.
We now estimate the total errors in the adiabatic process of the quantum simulation, using ∆ asy ≈ N s N p g s t s n , ∆ sym ≈ N s N p g s t s n , given in Eqs. (2.17) and (2.24). We choose thefinal value of g to be g f = 2 . For each step, the increase of g is set to be g s = 0 . .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 n . n is different in different casesas described in the following.First consider the asymmetric decomposition. For d=3 lattice, the number of plaquettesis N p = 24 . The number of substeps in the decomposition is set to be n = 200 . The timefor each step is chosen to be t s = 0 . . The total evolution time is then t f = N s t s = 200 . t s /n = 0 . × − . The total error is about . × − . For d=2, the number of plaquettesis N p = 9 . As the number of qubits are less than in d=3, we set the number of substepsto be n = 5000 . We choose the time step to be t s = 0 . . The total evolution time is then t f = N s t s = 400 . Now t s /n ≈ . × − . The total error is about . × − .– 14 –ow consider the asymmetric decomposition. As the error is larger than the symmetricone by one factor of t s /n , we use larger value of n = 500 . We use smaller value of t s = 0 . for d=3, and t s = 0 . for d=2. Therefore, t s /n = 0 . × − for d=3, and t s /n = 0 . × − for d=2. Therefore the total error is . × − for d=3, . × − for d=2.These parameter values are chosen to allow the computation to be completed in anacceptable time under adiabatic condition. The computation is proportional to N q , where N q is the number of qubits, which is equal to N l + 1 in the present model. It is also pro-portional to N s and n . For the parameter values given above, the time for the computationbased on asymmetric decomposition, run on a K40m server, is about 7 days for d = 3 and16 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 accuraciesare all acceptable. The difference between the computation times is due to the hardwaredifference rather than the decomposition methods.Our simulation satisfies the adiabatic condition, which says that the variation of theHamiltonian should be slower than the dynamical time scale, in other words, the matrixelement of ∂H/∂t should be smaller than the square of the energy gap [12, 61]. In oursimulation, the matrix element of ∂H/∂t is of the order of ∂g/∂t = g s /t s , which varies from . to . .In the weak coupling limit g → , the ground state is with all Z (cid:3) = − , while the firstexcited state is one with a pair of visons, that is, a pair of plaquettes with Z (cid:3) = 1 , whichcan be created by flipping the qubits on the links pierced by a string on the dual lattice.Thus the gap is . Its square is a lot larger than g s /t s . At g = 0 , the ground state is d -fold degenerate, with different eigenvalues ± of the d ’t Hooft operators. They becomenondegenerate when g = 0 , and the one with eigenvalues of all V µ ’s being 1 becomes theunique ground state. Hence there are small energy splittings between the ground stateand 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 isprepared 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. Thegauge invariance dictates that in the first excited state, there is a plaquette with σ xl = 1 on its four links. Hence the energy gap is g , which is very large. It is known that groundstates and spectra in weak and strong coupling limits are all stable up the QPT point [36].For a finite size L , the adiabatic condition is satisfied for general values of g , includingthe critical point as the energy gap is of the order of /L , L being the linear size of thesystem [12]. Hence the gap is about . for a general value of g . Now we turn to the results of our actual pseudoquantum simulations of quantum Z LGTon 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 to in substeps of t s /n . Aftereach step of t s , several quantities are calculated.– 15 – .1 Wegner-Wilson loops As g increases, the ground state evolves, from the equal superposition of the configurationswith all plaquettes in Z (cid:3) = − , to be near the state with all qubits in σ x = 1 . During thisprocess, it undergoes a QPT, which can be characterized in terms of the Wegner-Wilsonloop operator W C defined along a loop C on the direct lattice [36, 38]. In the confinedphase g > g c , W C obeys the area law h W C i ∝ g − A C = exp( − A C ln g ) , where A C is the areaenclosed by the contour C . In the deconfined phase at g < g c , W C obeys the perimeterlaw h W C i = exp( − B ( g ) P C ) , where P C is the perimeter of the contour C , B ( g ) is a smoothfunction of g and vanishes as g → . At g = 0 , all Z (cid:3) = − , the Z flux is expelled, h W C i = 1 . At small values of g , the fluctuations lead to the perimeter law. Our simulationsconfirm this picture.In our simulation, for d=3 and d=2 lattices respectively, we choose three contoursdenoted as c1, c2, c3, as shown in Fig. 6. On each lattice, the perimeter ratios are 1:1.5:2while the area ratios are 1:2:3. (a) (b) (c)(d) (e) (f) Figure 6 . (a-c) Wegner-Wilson loops c1, c2 and c3 on d=3 lattice. (d-f) Wegner-Wilson loops c1,c2 and c3 on d=2 lattice.
Features on Wegner-Wilson loops are shown in Fig. 7. As g → , h W C i → . With theincrease of g , h W C i decrease relatively slowly when g is small, as indicated in subfigures(a) and (c), and the ratios between log h W C i ’s of different loops equal the ratios betweenthe perimeters, as indicated in subfigures (b) and (d). When g is relatively large, h W C i ’sdecrease as some powers of g , and the powers are shown to be areas, since the ratios between log h W C i ’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 h W C i ratio, which will be discussed below.Previous tensor network calculation for d=2 gave h W C i as functions of perimeter andarea for a small and a large values of coupling constant, respectively [46]. Complementarily,– 16 – a) < W c > g
044 =0 . [42]. A cluster Monte Carlo calculation implies g c = 1 / . . , while therewere other calculations implying g c between / .
046 = 0 . and / .
742 = 0 . [44].An exact diagonalization calculation of TIM on L = 6 lattice implies g c = 0 . [43]. Anentanglement renormalization calculation of TIM on L = 54 lattice implies g c = 1 / .
075 =0 . [45].In the tensor network calculation on lattices with L ≤ [46], g c is . from theenergy gap in the same topological sector of the ground state, and is . from the stringoperator expectation, and is . from the overlap between the ground state and thelowest energy states in other topological sectors acted by Wegner-Wilson operators on thenon-contractible loops. – 18 – a) g
24 + 4 n , which are not considered in the plot. The prominent feature is that it also vanishesat x = ± . (c) DOS of Z eigenstates for d=2 × lattice. Note that it identically vanishes when z = − n , which are not considered in the plot. (d) DOS of X eigenstates for d=2 × lattice.It vanishes when x = −
18 + 4 n , which are not considered in the plot. The prominent feature is thatit also vanishes at x = ± . (a) and (b) satisfy D Z ( z, g ) = D X ( x = z, /g ) , confirming self-dualityin d=3. (c) and (d) confirm the absence of self-duality in d=2. − and towards , while the expectation value of X decreases from towards − . A first-order phase transition is one where there is a discontinuity of a first derivative ofthe free energy or energy, at the critical value of the control parameter, in contrast to asecond-order phase transition, where the first derivatives are continuous while there is adiscontinuity of its first derivatives at the critical point.The thermal phase transition in the classical Z LGT is first-order in 4 spatial dimen-sions, and is second-order in 3 spatial dimensions [58]. This suggests that for quantum Z LGT, the QPT is first-order in d=3 spatial dimensions, and is second-order in thed=2 spatial dimensions, because the classical theory in D spatial dimensions corresponds– 22 – a) -24 -20 -16 -12 -8 -4 0 4 8 12 16 20 24-24 -20 -16 -12 -8 -4 0 4 8 12 16 20 240.00.10.20.30.40.50.60.7 0.00.10.20.30.40.50.60.7 D O S z g 0.9 1 1.1 (b) -24 -20 -16 -12 -8 -4 0 4 8 12 16 20 24-24 -20 -16 -12 -8 -4 0 4 8 12 16 20 240.00.10.20.30.40.50.60.7 0.00.10.20.30.40.50.60.7 D O S x g 1.111 1.0 0.909 (c) -9 -5 -1 3 7-9 -5 -1 3 70.00.10.20.30.40.50.60.70.80.91.0 0.00.10.20.30.40.50.60.70.80.91.0 D O S z g 0.3 0.4 0.5 (d) -18 -14 -10 -6 -2 2 6 10 14 18-18 -14 -10 -6 -2 2 6 10 14 180.00.10.20.30.40.50.60.70.80.91.0 0.00.10.20.30.40.50.60.70.80.91.0 D O S x g 3.333 2.5 2.0 Figure 13 . (a) Ground-state DOS of Z eigenstates as a function of z , with g = 0 . , . , . , ond=3 × × lattice. (b) Ground-state DOS of X eigenstates as a function of x , with g = 1 . , . , . , which are the inverses of the values of g in (a), on d=3 × × lattice. (c) Ground-stateDOS of Z as a function of z , with g = 0 . , . , . , on d=2 × lattice. (d) Ground-state DOS of X as a function of x , with g = 2 . , , . , which are the inverses of the values of g in (c), on d=2 × lattice. To reach these values, we specifically extend the range of g in the adiabatic evolution.(a) and (b) satisfy D Z ( z, g ) = D X ( x = z, /g ) , confirming self-duality in d=3. (c) and (d) confirmthe absence of self-duality in d=2. to the quantum theory in d=D-1 spatial dimensions, in other words, D=d+1 spacetimedimensions, where 1 represents the time dimension.Limited by the smallness of the lattice size, one cannot make conclusion about whetherthere exists discontinuity in the slope of h H i in Fig. 8 or ∂ h H i /∂g in Fig. 9.Therefore, we make a comparison between d=3 and d=2, by putting together h Z i andsecond order derivatives ∂ h H i /∂g for d=2 and d=3, as functions of ( g − g c ) /g c , whichis dimensionless (Fig. 14). Then it can be seen clearly that in d=3, the change of h Z i atQPT is much steeper and the valley is much sharper. This would be consistent with the– 23 – a) -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0-1.0-0.8-0.6-0.4-0.20.0 -1.0-0.8-0.6-0.4-0.20.0 < Z > / | < Z ( g = ) > | (g-g c )/g c d=3 d=2 (b) -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0-50-40-30-20-100 -50-40-30-20-100 < H > / g (g-g c )/g c d=3 d=2 Figure 14 . Comparison of results for d=3 × × and d=2 × lattices, as functions of ( g − g c ) /g c , where g c is the respective critical value in each dimensionality. (a) h Z i , relative to theabsolute value at g = 0 . (b) Second derivative of the energy with respect to g . It can be seen clearlythat in d=3, the change of h Z i at QPT is much steeper, and the valley of second derivative of theenergy is much sharper. claim that in d=3, QPT is first-order, even though the discontinuity and the singularityare rounded out, while in d=2, QPT is second-order.Now we examine the properties of the ground states. As can be seen in Fig. 13, in d=3,distributions of DOS of Z appear significantly different before and after QPT. For example,when g = 0 . < g c , the peaks are at z = − and z = − ; when g = 1 . > g c , the peak isat z = − . When g = g c = 1 , there are peaks at z = − and z = − , i.e. the locationsof peaks when g < g c , and at z = − , i.e. the location of peak when g > g c . Same featureexists in DOS of X, by replacing g as /g because of self-duality.This feature indicates a quantum version of phase coexistence at the critical point, ahallmark 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 alsoslowly varies with g , as a same phase within this regime. When g = g c , the ground state isroughly a nearly-equal superposition of the two states. This feature is absent in d=2, alsoseen in Fig. 13.Therefore, our simulations support the implication from 4-dimensional classical Z LGTthat QPT in quantum Z 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(Fig 7). In d=3, there is a dip in the ratio of the logarithms, which is absent in d=2. Wehave examined that the dip is right at g c determined from the lowest point of the ∂ h H i /∂g in Fig. 14. It is not excluded that the dip is a finite-size effect. It is also possible that thedip is a signature of first-order QPT at g = g c , as a consequence of the quantum phasecoexistence. The matrix element of the Wegner-Wilson loop operator between the two– 24 –tates representing the two phases possibly decrease the expectation in their superpositionrepresenting the phase coexistence. The dip in the ratio between c3 and c1 is deeper thanthat 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 studiedin a quantum adiabatic algorithm for the exact cover problem, which was regarded as acriterion for first-order QPT [62]. Nevertheless, no definite conclusion can be drawn yet. The phase transition in the classical Z LGT is between states that cannot be distinguishedin symmetry [16, 34]. It is now called topological phase transition [38]. Our simulationverifies that QPT in quantum Z LGT is also topological.First, the absence of symmetry breaking is directly verified in h σ zl i of qubit l , as afunction of g . As can be seen in Fig. 15, h σ zl i remains consistent with for all values of g .For this qubit l , we have also studied the expectation values of x l ≡ − σ xl , z l = P (cid:3) ∋ l Z (cid:3) ,which is the average of Z (cid:3) ’s of the plaquettes sharing the qubit l normalized by the numberof plaquett per link, and h l = z l + x l . (a) g
1) = E (1 , − ,
1) = E (1 , , − , E = E ( − , − ,
1) = E (1 , − , −
1) = E ( − , , − , E = E ( − , − , − . The remaining degeneracy is due to rotational symme-try. At g = 0 , we prepare the ground state ( V x , V y , V z ) = (1 , , and other three ones (1 , − , , ( − , − , and ( − , − , − , using the corresponding W µ operators, definedalong 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 withone of the states studied. Then we execute adiabatic algorithm to obtain the dependenceof the energies and the splittings on g before and after the QPT, as shown in Fig. 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 = ± and V y = ± . Thequantum adiabatic algorithm is executed on each of them. They become nondegeneratewhen g = 0 . In our simulation, we calculate the four energies E ( V x , V y ) , which are E = E (1 , , , E = E ( − , and E = E (1 , − , E = E ( − , − . The simulation confirmsthat E is the lowest while E = E . As shown in the inset of in (d) of Fig. 16, the splitting E i ≡ E i − E can be fitted by least-square method as E L +1 = E L +1 = 1 . g + 0 . ,E L +1 = 2 . g + 0 . , (3.8)with L = 3 , which verify E i ∝ g L +1 , (3.9)with L = 3 fixed.Previous tensor network calculation found that the splittings between the four eigen-states of V x and V y exponentially decay with L for given small g [46]. This is consistent withand complements our result, as g L +1 = ge − L/ξ , with ξ ≡ − / ln g . Our result reveals the g dependence but does not give the L dependence, while their result gives the L dependencebut does not give the g dependence. – 26 – a) E i g E E E E E i g (b) E E E E i g E i g (c) E i g E E E E E i g (d) E Li g E E E E Li g Figure 16 . (a) E ( V x , V y , V z ) on d = 3 2 × × lattice. E = E (1 , , , E = E ( − , ,
1) = E (1 , − ,
1) = E (1 , , − , E = E ( − , − ,
1) = E (1 , − , −
1) = E ( − , , − , E = E ( − , − , − ,as functions of g . (b) On d = 3 2 × × lattice, E i ≡ E i − E as a function of g , ( i = 2 , , , (c) E ( V x , V y ) on d = 2 3 × lattice. E = E (1 , , , E = E ( − , , E = E (1 , − , E = E ( − , − .(d) On d = 2 3 × lattice, E i L +1 , ( i = 2 , , , as functions of g , where E i ≡ E i − E . Implementing quantum adiabatic algorithm in terms of quantum circuit consisting of uni-versal quantum gates of one or two qubits, we present a digital scheme of quantum adiabaticsimulation of quantum Z LGT, which is important in both high energy physics and con-densed matter physics. Furthermore, we classically demonstrate this quantum simulationscheme by running the GPU simulator QuEST on a Nvidia GPU server, and obtain resultsuseful in this field of physics.In our algorithm, we have generalized Trotter decomposition and symmetrized Trotterdecomposition such that the Hamiltonian adiabatically varies during the decomposition.Hence we have proposed a scheme of digital adiabatic quantum simulation. Hopefully, this– 27 –ethod can be used for various problems.We have studied quantum Z 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 evo-lution. Since the lattices are very small, singularities in QPT are rounded out. But keyfeatures are observed. It is indicated that QPT is first-order in d=3 and is second-order ind=2, with critical point consistent with previously known results. In d=3 and d=2 respec-tively, we also clearly observe topological characteristics of QPT, including the vanishingof the “magnetization” for all values of the coupling constant, the excitation properties, thelowest energies in different topological sectors and the their splittings, which are propor-tional 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 presentedthe unique numerical result on quantum Z LGT in 3 spatial dimensions.This work seems to be the first complete demonstration as a proof of principle, albeitin a classical simulator, carrying through quantum simulation of a LGT and observe thekey features of physics. Thereby it demonstrates that real quantum simulations of LGTs infuture can be done in the proposed way.On the other hand, this work also shows that high-performance classical demonstrationof quantum simulation, which may be dubbed pseudoquantum simulation, represents a newway of computation, in addition to facilitating the development of quantum software forexperimental 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 isvaried, so it is very natural and convenient to obtain various quantities as functions ofthis parameter, which may not be convenient in other methods, hence adiabatic quantumsimulation and pseudoquantum simulation are convenient tools of QPT study.The classical demonstration proceeds according to rules of quantum mechanics, there-fore pseudoquantum simulation is legitimately a reliable approach, as far as computationalresources allow. Compared with other computational approaches, it can be used withoutintricate algorithmic design depending on the details of the computed problem, as in manyother numerical methods, which are often applicable only to one or two dimensions.As the next step regarding quantum and pseudoquantum simulations of LGT, we shallstudy Fermions coupled with the Z gauge field using our method, on which the results canbe compared with the existing results of QMC calculations, which are free of Fermion signproblem. Thereby the method can be further benchmarked, which can then be applied toother LGTs, for which Fermion sign problem exists in MC-based methods, as well as theLGTs that recently have been studied by using tensor network methods [47, 48].After the release of the present work as a preprint (arXiv:1910.08020), there appearedmore recent progress in the related fields, including DMRG study of the one-dimensionalspinless Fermions coupled with Z gauge field [63], schemes of measuring nonlocal observ-ables in the quantum simulation of a general LGT [64], experimental observation of gaugeinvariance in a 71-site quantum simulator of an extended U(1) LGT [65], and a sign-freeQMC study of Z gauge field coupled with both Fermions and Bosons, demonstrating thetransition from conventional metal to orthogonal metal [66].– 28 – cknowledgments This work was supported by National Science Foundation of China (Grant No. 11574054).
References [1] Feynman, R. P.,
Simulating Physics with Computers , Int. J. Theor. Phys. , 467 (1982).[2] Lloyd, S., Universal quantum simulators Science , , 1073 (1996).[3] Cirac J. I. and Zoller P., Goals and opportunities in quantum simulation , Nature Physics ,264 (2012).[4] Georgescu I. M., Ashhab S. and Rev. Mod. Phys. , 153 (2014).[5] Preskill J., Quantum computing in the NISQ era and beyond , Quantum , 79 (2018).[6] Arute F et al. , Quantum supremacy using a programmable superconducting processor , Nature,574, 505 (2019); Wang H., et al. , Boson Sampling with 20 Input Photons and a 60-ModeInterferometer in a -Dimensional Hilbert Space , Phys. Rev. Lett. 123, 250503 (2019).[7] Childs, A. M. et al. Exponential algorithmic speedup by quantum walk ,arXiv:quant-ph/0209131.[8] Lloyd, S., Mohseni, M. and Rebentrost, P., Quantum principal component analysis , NaturePhysics , 631–633 (2014).[9] Rebentrost, P., Steffens, A., Marvian, I. and Lloyd, S., Quantum singular-valuedecomposition of nonsparse low-rank matrices , Phys. Rev. A , 6 (2018).[10] Wossnig, L., Zhao, Z. K. and Prakash, A., Quantum linear system algorithm for densematrices , Phys. Rev. Lett. , 5 (2018).[11] Farhi, E., Goldstone, J., Gutmann, S., Sipser, M.,
Quantum computation by adiabaticevolution , eprint arXiv:quant-ph/0001106 (2000); Farhi, E., Goldstone, J., Gutmann, S.,Sipser, M.,
A Quantum adiabatic evolution algorithm applied to random instances of anNP-complete problem , Science , 472 (2001).[12] Hamma, A. and Lidar, D. A.,
Adiabatic preparation of topological order , Phys. Rev. Lett. , 4 (2008).[13] Mueck,L.,
Quantum Software , Nature , 171 (2017); Chong, F. T., Franklin, D. andMartonosi, M.,
Programming languages and compiler design for realistic quantum hardware ,Nature , 180 (2017).[14] For reviews, see Uwe-Jens Wiese,
Ultracold quantum gases and lattice systems: quantumsimulation of lattice gauge theories , Annalen der Physik 525, 777 (2013); Zohar E., Cirac J.I., Reznik B.,
Quantum simulations of lattice gauge theories using ultracold atoms in opticallattices , Rep. Prog. Phys. , 014401 (2016); Dalmonte M. and Montangero S., Lattice gaugetheory simulations in the quantum information era , Contemporary Phys. , 388 (2016).[15] Wilson K. G., Confinement of quarks , Phys. Rev. D , 2445 (1974).[16] Kogut J. B. An introduction to lattice gauge theory and spin systems , Rev. Mod. Phys. ,659 (1979); Kogut J. B., The lattice gauge theory approach to quantum chromodynamics ,Rev. Mod. Phys. , 775 (1983). – 29 –
17] Creutz M.,
Quarks, Gluons, and Lattices , Cambridge University Press, Cambridge (1985);Gattringer C. and Lang C.,
Quantum Chromodynamics on the Lattice: An IntroductoryPresentation , Lecture Notes in Physics, Springer, Berlin (2010).[18] Hirsch, J. E. et al. , Monte Carlo Calculations of One-dimensional Fermion systems , Phys.Rev. B , 5033 (1982); Blankenbecler R., Scalapino, D. J. and Sugar R. L., Monte CarloCalculations of Coupled Boson-Fermion Systems,
Phys. Rev. D , 2278 (1981).[19] Troyer M. and Wiese U. J., Computational complexity and fundamental limitations toFermionic quantum Monte Carlo simulations , Phys. Rev. Lett. , 170201 (2005).[20] Li, Z. and Yao, H., Sign-problem-free Fermionic quantum Monte Carlo: developments andapplications , Ann. Rev. Condens. Matter Phys. , 337 (2019).[21] Osterloh, K., Baig, M., Santos, L., Zoller, P. and Lewenstein, M., Cold atoms in non-abeliangauge potentials: From the Hofstadter “moth” to lattice gauge theory , Phys. Rev. Lett. ,010403 (2005); Büchler, H. P., Hermele, M., Huber, S. D. , Fisher, M. P. A. and Zoller, P., Atomic quantum simulator for lattice gauge theories and ring exchange models , Phys. Rev.Lett. , 040402 (2005); Weimer H., Müller M., Lesanovsky I, Zoller P and Büchler, ARydberg quantum simulator , Nature Physics , 382 (2010); Banerjee D. et al. , AtomicQuantum Simulation of Dynamical Gauge Fields Coupled to Fermionic Matter: From StringBreaking to Evolution after a Quench , Phy. Rev. Lett. , 175302 (2011); Banerjee D. etal. , Atomic Quantum Simulation of U(N) and SU(N) Non-Abelian Lattice Gauge Theories ,Phys. Rev. Lett. , 125303 (2013).[22] Zohar E. and Reznik B.,
Confinement and Lattice Quantum-Electrodynamic Electric FluxTubes Simulated with Ultracold Atoms , Phy. Rev. Lett. , 275301 (2011); Zohar E., CiracJ. I., Reznik B.,
Simulating Compact Quantum Electrodynamics with Ultracold Atoms:Probing Confinement and Nonperturbative Effects , Phy. Rev. Lett. , 125302 (2012);Zohar E., Cirac J. I., Reznik B.,
Simulating (2+1)-Dimensional Lattice QED with DynamicalMatter Using Ultracold Atoms , Phy. Rev. Lett. , 055302 (2013); Zohar E., Cirac J. I.,Reznik B.,
Cold-Atom Quantum Simulator for SU(2) Yang-Mills Lattice Gauge Theory ,Phys. Rev. Lett. , 125304 (2013).[23] Tagliacozzo, L., Celi, A., Orland, P., Mitchell, M. W. and Lewenstein, M.,
Simulation ofnon-abelian gauge theories with optical lattices , Nature Communications , 8 (2013).[24] Zohar E., Cirac J. I., Reznik B., Quantum simulations of gauge theories with ultracold atoms:Local gauge invariance from angular-momentum conservation , Phys. Rev. A , 023617(2013).[25] Zohar, E., Farace, A., Reznik, B. and Cirac, J. I. Digital quantum simulation of Z(2) latticegauge theories with dynamical fermionic matter , Phys. Rev. Lett. , 070501 (2017); ZoharE., Farace A., Reznik B., Cirac J. I.,
Digital lattice gauge theories , Phys. Rev. A , 023604(2017); Bender J., Zohar E., Farace A., Cirac J., Digital quantum simulation of lattice gaugetheories in three spatial dimensions , New J. Phys. , 093001 (2018).[26] L. Tagliacozzo, Celi A., Zamoraa A., Lewenstein M., Optical Abelian lattice gauge theories ,Ann. Phys. 330, 160 (2013).[27] Byrnes, T. and Yamamoto, Y.
Simulating lattice gauge theories on a quantum computer ,Phys. Rev. A , 16 (2006).[28] Lamm, H., Lawrence, S., Yamauchi, Y., General Methods for Digital Quantum Simulation ofGauge Theories , Phys. Rev. D , 034518 (2019). – 30 –
29] Kolco, N.,
Quantum-classical computation of Schwinger model dynamics using quantumcomputers , Phys. Rev. A , 032331 (2018).[30] Martinez E. A. et al. , Real-time dynamics of lattice gauge theories with a few-qubit quantumcomputer , Nature 534, 516 (2016); Kokail C.,
Self-verifying variational quantum simulationof lattice models , Nature , 355 (2019).[31] Kasper V. et al. , Implementing quantum electrodynamics with ultracold atomic systems , NewJ. Phys. 19, 023030 (2017); Mil A. et al. , Realizing a scalable building block of a U(1) gaugetheory with cold atomic mixtures , arXiv:1909.07641.[32] Görg, F. et al. , Realization of density-dependent Peierls phases to engineer quantized gaugefields coupled to ultracold matter , Nat. Phys. , 1161 (2019); Schweizer C. et al. , Floquetapproach to Z lattice gauge theories with ultracold atoms in optical lattices , Nature Physics , 1168 (2019); Barbiero, L., Coupling ultracold matter to dynamical gauge fields in opticallattices: From flux-attachment to Z2 lattice gauge theorie , arXiv:1810.02777.[33] Davoudi Z. et al. , Towards analog quantum simulations of lattice gauge theories with trappedions , arXiv:1908.03210.[34] Wegner, F. J.,
Duality in generalized Ising models and phase transitions without local orderparameters , Journal of Mathematical Physics , 2259–and (1971).[35] Fradkin, E., Susskind, L., Order and disorder in gauge systems and magnets , Phys. Rev. D , 2637-2658 (1978).[36] Fradkin, E., Field Theories of Condensed Matter Physics , Cambridge University Press,Cambridge (2013).[37] Ercolessi E., Facchi P., Magnifico G., Pascazio S. and Pepe F. V.,
Phase transitions in Z n gauge models: towards quantum simulations of the Schwinger-Weyl QED , Phys. Rev. D ,074503 (2018); Magnifico G., Dalmonte M., Facchi P., Pascazio S., Pepe F. V. and ErcolessiE., Real time dynamics and confinement in the Z n Schwinger-Weyl model for 1+1 QED ,arXiv:1909.04821.[38] Sachdev, S.,
Topological order, emergent gauge fields, and fermi surface reconstruction ,Reports on Progress in Physics , 014001 (2019).[39] Kitaev, A. Y., Fault-tolerant quantum computation by anyons , Annals of Physics , 2–30(2003).[40] Li, K. R. et al. , Experimental identification of non-abelian topological orders on a quantumsimulator , Phys. Rev. Lett. , 5 (2017).[41] Luo, Z. H. et al. , Experimentally probing topological order and its breakdown through modularmatrices , Nature Physics , 160 (2018).[42] Rieger, H. and Kawashima, N., Application of a continuous time cluster algorithm to thetwo-dimensional random quantum Ising ferromagnet , Eur. Phys. J. B. , 233 (1999).[43] Hamer, C. J., Finite-size scaling in the transverse Ising model on a square lattice , J. Phys. A:Math. Gen. , 6683 (2000).[44] Blote, H. W. J. and Deng, Y. J., Cluster Monte Carlo simulation of the transverse Isingmodel , Phys. Rev. E , 8 (2002).[45] Evenbly G., and Vidal G., Entanglement Renormalization in Two Spatial Dimensions , Phys.Rev. B , 180406 (2009). – 31 –
46] Tagliacozzo, L. and Vidal G.,
Entanglement renormalization and gauge symmetry , Phys.Rev. B , 115127 (2011).[47] Tagliacozzo L., Celi A., and Lewenstein M., Tensor Networks for Lattice Gauge Theories withContinuous Groups , Phys. Rev. X (4), 041024 (2014); Tschirsich F., et al. , Phase Diagramand Conformal String Excitations of Square Ice using Gauge Invariant Matrix ProductStates , SciPost Phys. 6, 028 (2019); J. Haegeman J., et al. , Gauging Quantum States: FromGlobal to Local Symmetries in Many-Body Systems , Phys. Rev. X (1), 011024 (2015).[48] Celi, A., et al. , Emerging 2D Gauge theories in Rydberg configurable arrays ,arXiv:1907.03311.[49] Senthil, T. and Fisher, M. P. A., Z gauge theory of electron fractionalization in stronglycorrelated systems , Phys. Rev. B , 7850 (2000); Rĺźegg A., Huber S. D. and Sigrist M,Phys. Rev. B , 155118 (2010); Rahul Nandkishore, R., Metlitski, M. A., and Senthil T., Orthogonal metals: The simplest non-Fermi liquids , Phys. Rev. B , 045128 (2012).[50] Trebst, S. et al. , Breakdown of a Topological Phase: Quantum Phase Transition in a LoopGas Model with Tension , Phys. Rev. Lett. , 070602 (2007).[51] Dagotto, E., Kogut, J. B. and Kocić, A., Computer simulation of chiral-symmetry breaking in(2+1)-dimensional QED with N flavors , Phys. Rev. Lett. , 1083 (1989).[52] Gazit, S., Randeria, M. and Vishwanath, A., Emergent dirac fermions and broken symmetriesin confined and deconfined phases of Z gauge theories , Nature Physics , 484 (2017);Gazit, S., Assaad, F. F., Sachdev, S., Vishwanath, A. and Wang, C., Confinement transitionof Z gauge theories coupled to massless fermions: Emergent quantum chromodynamics andSO(5) symmetry , Proc. Nat. Acad. Sci. (USA) , E6987 (2018); Gazit, S., Assaad, F. F.and Sachdev, S. , Fermi-surface reconstruction without symmetry breaking , arXiv:1906.11250.[53] Assaad, F. F. and Grover, T.
Simple fermionic model of deconfined phases and phasetransitions , Phys. Rev. X , 041049 (2016); Frank, J., Huffman, E. and Chandrasekharan, S. Emergence of Gauss’ Law in a Z Lattice Gauge Theory , arXiv:1904.05414.[54] Prosko, C., Lee, S. and Maciejko J.,
Simple Z lattice gauge theories at finite fermiondensity , Phys. Rev. B , 205104 (2017).[55] González-Cuadra, D., et al. , Intertwined topological phases induced by emergent symmetryprotection , Nat. Commun. , 2694 (2019); González-Cuadra, D., et al. , Symmetry-breakingtopological insulators in the Z Bose-Hubbard model , Phys. Rev. B , 045139 (2019).[56] Felser, T. et al. , Two-dimensional quantum-link lattice Quantum Electrodynamics at finitedensity , arXiv:1911.09693.[57] Jones, T., Brown, A., Bush, I. and Benjamin, S.,
Quest and high performance simulation ofquantum computers , Scientific Report , 10736 (2019).[58] Creutz, M., Jacobs, L. and Rebbi, C., Experiments with a gauge-invariant ising system ,Phys. Rev. Lett. , 1390 (1979); Creutz, M., Jacobs, L. and Rebbi, C., Monte Carlo studyof Abelian lattice gauge theories , Phys. Rev. D , 1915 (1979).[59] Kennedy A. D., Clark M. A. and Silva P. J., Force Gradient Integrators , arXiv:0910.2950.[60] Nielsen, M. A. and Chuang I. L.,
Quantum computation and quantum information ,Cambridge University Press, Cambridge (2000), P. 210.[61] Shi, Y. and Wu, Y.-S.,
Perturbative formulation and nonadiabatic corrections in adiabaticquantum-computing schemes , Phys. Rev. A , 024301 (2004). – 32 –
62] Young, A. P., Knysh, S. and Smelyanskiy, V. N.,
First-order phase transition in the quantumadiabatic algorithm , Phys. Rev. Lett. , 4 (2010).[63] Borla U. et al. , Confined phases of one-dimensional spinless Fermions couple to Z gaugetheory, Phys. Rev. Lett. , 120503 (2020).[64] Zohar E.,
Local Manipulation and Measurement of Nonlocal Many-Body Operators in LatticeGauge Theory Quantum Simulators , arXiv:1911.11156.[65] Yang B. et al. , Observation of gauge invariance in a 71-site quantum simulator ,arXiv:2003.08945.[66] Chen C., Xu X. Y., Qi Y. and Meng Z. Y.,
Metal to Orthogonal Metal Transition , Chin.Phys. Lett. (4), 047103 (2020).(4), 047103 (2020).