Toward Quantum Simulations of Z2 Gauge Theory Without State Preparation

Quantum computers hold the promise of solving problems beyond the capability of classical computers in many aspects of high energy physics, in particular calculations involving finite density and real time evolution [1]. Achieving this promise requires the development not only of hardware with reduced noise but efficient algorithms as well. At present, a major roadblock to such simulations is the creation of the initial, strongly-coupled quantum states on the quantum computer. Current methods for this state preparation are expensive to implement (usually dominating the circuit depth of simulations) and often difficult to generalize beyond ground states. In particular, the preparation of scattering states requires substantial complexity [2–10]. Consequently, attempts to develop more efficient procedures have been explored [11–15]. One proposal [12, 15] avoids the issue of initial state preparation entirely by stochastically sampling the density matrix, ρ, classically, and then passing simpler basis states to the quantum computer with proper weights. This method can be understood as a Schwinger-Keldysh contour [16, 17], where the classically sampled Euclidean path integral is matched at its boundaries to a Minkowski path integral computed on the quantum computer. For quantum field theories, the natural way to sample ρ is through a lattice field theory (LFT) calculation with open boundary conditions. In this picture, the quantum computer acts as an operator insertion into a standard LFT calculation, except it returns time-dependent quantities. For thermal states, one need only sample from the Euclidean path integral at the desired inverse temperature β [12]. For other states (e.g. the pion or scattering protons), it was proposed to use projections of the configurations onto the quantum numbers of the desired states [15]. In this paper, we demonstrate the nonthermal state preparation using projection operators suggested in [15] by computing results for oneand two-“particle” plane waves in the 2+1d Z2 lattice gauge theory. To do this, we begin begin with a discussion of the action and Hamiltonian formulations of the Z2 lattice gauge theory in Sec. II. This is followed in Sec. III, with a review of the path integral


I. INTRODUCTION
Quantum computers hold the promise of solving problems beyond the capability of classical computers in many aspects of high energy physics, in particular calculations involving finite density and real time evolution [1]. Achieving this promise requires the development not only of hardware with reduced noise but efficient algorithms as well. At present, a major roadblock to such simulations is the creation of the initial, strongly-coupled quantum states on the quantum computer. Current methods for this state preparation are expensive to implement (usually dominating the circuit depth of simulations) and often difficult to generalize beyond ground states. In particular, the preparation of scattering states requires substantial complexity [2][3][4][5][6][7][8][9][10]. Consequently, attempts to develop more efficient procedures have been explored [11][12][13][14][15].
One proposal [12,15] avoids the issue of initial state preparation entirely by stochastically sampling the density matrix, ρ, classically, and then passing simpler basis states to the quantum computer with proper weights. This method can be understood as a Schwinger-Keldysh contour [16,17], where the classically sampled Euclidean path integral is matched at its boundaries to a Minkowski path integral computed on the quantum computer. For quantum field theories, the natural way to sample ρ is through a lattice field theory (LFT) calculation with open boundary conditions. In this picture, the quantum computer acts as an operator insertion into a standard LFT calculation, except it returns time-dependent quantities. For thermal states, one need only sample from the Euclidean path integral at the desired inverse temperature β [12]. For other states (e.g. the pion or scattering protons), it was proposed to use projections of the configurations onto the quantum numbers of the desired states [15].
In this paper, we demonstrate the nonthermal state preparation using projection operators suggested in [15] by computing results for one-and two-"particle" plane waves in the 2+1d Z 2 lattice gauge theory. To do this, we begin begin with a discussion of the action and Hamiltonian formulations of the Z 2 lattice gauge theory in Sec. II. This is followed in Sec. III, with a review of the path integral matching algorithm. Estimates of the required quantum resources are in Sec. IV. Within Sec. V are numerical results obtained on a quantum simulator for 3 2 and 4 2 lattices, and we conclude and consider future work in Sec. VI.

II. MODEL
Simulations on quantum computers are naturally formulated in the language of Hamiltonians, while classical, Euclidean lattice field theory uses actions, therefore we must specify both in order to perform calculations. We begin our discussion with the anisotropic Wilson action for a general gauge theory which can be written as: where t, s label the temporal and spatial directions and the plaquettes U i are formed from gauge links given by elements of the group. After gauge fixing the time-like links, we will classical sample with this action. Using β t = β s introduces an anisotropy ξ = a/a 0 where the physical lattice spacing in the temporal direction a 0 is not equal to the spatial one a. From Eq. (1), we can derive the Kogut-Susskind Hamiltonian [18] via the transfer matrix in the limit of ξ → ∞ (see [19,20]). This Hamiltonian is given by: where l ij are the conjugate variables of U s . We have introduced the Hamiltonian coupling β H = √ β s β t and the bare speed of light c = βs βt . While efficient digitization of gauge groups is a field of active research , specializing to the case of Z 2 gauge theory is straight forward because we map each link to a single qubit. It then follows that Eq. (2) becomes: Although this mapping is relatively simple, the plaquetteterm requires a four-qubit operation. This can be avoided by reformulating this theory in its dual representationthe transverse Ising model in 2+1d [46][47][48]. The relation between the theories maps the flux on the plaquette to a spin on the dual lattice [46][47][48] and is graphically depicted in Fig. 1.
In this representation, the plaquette-term becomes a single-qubit operation, while the l ij term becomes a twoqubit one: where the summations over n correspond to the centers of the plaquettes in the gauge representation and the sum overμ is the unit vectors in thex andŷ directions. The relation between the two representations is seen in Fig. 1. This gauge theory has a confined and deconfined phase. When J Γ (β H 1), the model is deconfined and excitations correspond approximately to plaquettes with flux pointing in the opposite direction. On the other hand when Γ J (β H 1), the system is confined and excitations correspond to droplets with domain walls in theσ x basis. In this work, we use one set of couplings: J = 0.3 and Γ = 1.
The final object we need to define for our model is the projection operators for particle-like states. For this paper, we will investigate quantum states of fixed parity excited by the operator where k = 2π[ kx nx , ky ny ], r = [r x , r y ]. Due to the small lattices, we only consider plane wave excitations i.e. f (r) = 1. As larger lattices become available, it would be interesting to study how wave packets with nontrivial envelopes f (r) evolve since these should have superior overlap with physical particles.

III. ALGORITHM
We are interested in the matrix elements ψ i |O(t)|ψ j of operators O(t) = e iHt Oe −iHt between two states, ψ i and ψ j . The difficulty in preparing these strongly-coupled states on the quantum computer can be avoided by instead computing the matrix elements Ψ i |O(t)|Ψ j between basis states Ψ i which are cheaper to prepare, and then weighting various matrix elements properly. This forms the basis of our hybrid algorithm-the classical portion is used to obtain the weights, and the quantum potion computes matrix elements between easily prepared basis states.
To do this, we consider a thermal state given by a den- Provided this state has overlap with ψ i , ψ j , then there are two operators P, Q respectively that would project out the component of the thermal state. Thus, if one can properly sample from ρ, then the desired matrix elements can be written as where O(t) ji = Ψ j | O(t) |Ψ i . The notation · ρ denotes expectation values sampled from the distribution ρ ij . The overall normalization δ ij ρ measures the weight of i ρ ii which is often unneeded, but can be computed if desired [15]. Efficient classical sampling of the distribution ρ ij can be obtained from standard Euclidean path integral methods provided open boundary conditions (OBC) in time are used [49]. The classical side of the algorithm thus involves Monte Carlo simulations to sample the ρ in the gauge representation. The configurations are then transformed to the dual representation and using Eq. (5) for P, Q we obtain the initial and final states. These operators P, Q naturally belong to the classical portion of the algorithm because they are non-unitary projections. However it is possible to include them in the quantum portion using projective measurements such as [50].
The quantum portion of the algorithm implements a measurement of Ψ i |O(t)|Ψ j . Matrix elements of the form Ψ i |O(t)|Ψ i may be efficiently computed on a quantum processor [51][52][53][54][55]; thus we recast our matrix elements in terms of diagonal ones |Ψ u , |Ψ v = |Ψ i ±|Ψ j . If these are instead used as the initial states on the quantum computer, the desired matrix elements can be obtained via

IV. QUANTUM RESOURCES
Key to the effectiveness of this method is the cost of preparing |Ψ u , |Ψ v on given quantum hardware is reduced compared to the full |ψ i , |ψ j . To study this, we use CNOTs as the universal two-qubit gate and compute how many are required to implement these states.
We consider two ways of preparing the initial superposition states |Ψ u and |Ψ v which each have different costs and benefits. One method which reduces the quantum cost is to perform a quantum simulation for each summand f (r)e −ik·r σ + r of Eq. (5) and then classically perform the sum over r. The benefit of this is the cost of preparing the state is O(N 2 s ) two-qubit gates where N s is the lattice size in one dimension. A downside to this method is the number of circuits needed is (N 4 s ) Np where N p is the number of excitation operators a † k used. While this is tractable for few particles, it clearly scales poorly asymptotically. While this method is O(N 2 s ), the overall coefficient is always less than one. This is because in the case of Z 2 , the dual representation has only two spin states per site. Thus, |ψ i and |ψ j can differ by at most N 2 s /2 sites, depending on correlations between sites. Only in the case of differences are two-qubit operations required; therefore, the number of CNOT gates is ≈ c(J, Γ)N 2 s /2 where 0 ≤ c(J, Γ) ≤ 1 is a coupling-dependent number related to the fraction of sites that actually differ.
Averaging over the different k which are allowed, the number of CNOT gates required for the set of couplings considered for small lattices are listed in Table  I. From these results, we have confirmation of the O(N 2 s ) from a direct implementation for the split summation method. Using this data, we can estimate the asymptotic c(0.3, 1) ≈ 0.7. TABLE I. The average number of CNOT gates for preparing one (1p) and 2 (2p) plane waves on N 2 s lattice. 100 configurations were used in total to extract this estimate.
The second method uses an ancilla qubit. Suppose we have a pair of unitary operations,Û i andÛ j , that prepare the initial states |ψ i and |ψ j . These unitaries are expected to scale like O(N 2Np s ) if we want to prepare a complete plane wave but it is possible this would be O(N 2 s log(N s ) + N p ) if a Fourier transform can be used [56,57]. If the ancilla is prepared in the |± = |0 ± |1 state then controlled implementations of the unitaries can be applied and then transforming the ancilla back to the computational basis will efficiently add or subtract the two states. When the spatial dimensions have an even number of plaquettes, we can write the trotterization of the time evolution operator U (t) up to O N t δt 2 as a stroboscopic set of operators, alternating between even and odd site spins [58]: U V corresponds to the potential energy H V in Eq. (4) and U N,i correspond to kinetic energy in the n = x, y direction for i = e, o even or odd sites: x nσ x n+x (9) where e, o indicate the even or odd spatial sites used to decompose the lattice.Û V (δt) is a product of single qubit rotations, and thus may be applied in one step. The oper-atorsÛ N,i (δt) required two-qubit gates and thus depend on gates available. Together, the four U N,i require 2N 2 s two-qubit gates per Trotter-step if a CNOT or CZ gate are native. A reduction in cost to N 2 s can be obtained if the Moeller-Sorenson XX is available. While the gate count is roughly fixed, the circuit depth is dependent upon the observable investigated. For Hermitian observables (e.g. magnetization) then the two-qubit operations can be parallelized and the circuit depth will be approximately 4−10 two-qubit gates deep depending on the native quantum gates and spatial dimensions of the lattice. If a unitary such as Ψ i |U (t)|Ψ j is desired, then the two-qubit operations cannot be parallelized and ∼ N qubits Toffoli gates will be required instead because controlled time evolution operators will be necessary to measure this operator. These controlled evolution operators arise from needing to use an ancillary to measure the expectation value of this unitary non-Hermitian operator. Since Ψ u |Û (t)|Ψ u is not Hermitian, it is not directly accessible from a quantum computer. However if we prepare our system in the state and then apply a controlled version ofÛ (t), CU (t) and then measure eitherσ x orσ y on the ancilla qubit we will get the real and imaginary parts respectively [20]. Turning our regular evolution operator into a controlled evolution operator simply involves two transformaitons. The first makes all single qubit unitaries into controlled unitaries, which can easily be done with two CNOT gates, and a few single qubit rotations. Turning the two qubit operations into controlled unitaries simply involves turning all CNOT gates into Toffoli gates and all single qubit operations into controlled ones. As a simple metric for comparing the path integral matching procedure, we use the cost of a Trotterization step as a proxy for the gate cost of adiabatic state preparation. Clearly, the c(J, γ)N 2 s /2 two-qubit gate cost of preparing the one-or two-particle states is cheaper than the 2N 2 s required single Trotter step. Since adiabatic state preparation typically requires multiple Trotter steps, this suggests for these states, the path integral matching algorithm yields shallower circuits, albeit at the cost of increased classical resources and the total number of quantum circuits.

V. RESULTS
As a demonstration of how particle states can be studied with the path integral matching algorithm we study the one-and two-particle plane wave excitations. In the simulations we examine the system in the deconfined phase with 100 configurations generated with coupling parameters J = 0.3 and Γ = 1. Additonal information regarding the various choices of quantum simulation are listed in Table II. Different choices of states ψ i (e.g. a pion, two protons) and operators O(t) (e.g. electromagnetic or axial currents J µ ) the matrix element ψ i |O(t)|ψ j will correspond to transition form factors and cross-sections. Here, we consider the unitary operator In order to properly normalize the energies, we must subtract the energy of the vacuum state E V . This is done by modifying Eq. (10) to where N s -dependent E V can accurately be measured using classical methods. In this work, we will focus on the extracting the one-and two-particle energies. For this, our projection operators P = Q and for the case of exact time evolution, the matrix elements for Eq. (10) are where φ l are the eigenbasis ofĤ, and |β is the thermal state produced by the stochastic sampling. For the full ensemble, we compute Eq. (12) for multiple values of t and perform a Fourier transformation so that we can extract the spectral function f (aE). Using this spectral function we can extract the energy of given particle states. While the low-momentum particle states could be obtained from Euclidean calculations or exact diagonalization, we use this computation as a nontrivial test of our method to reproduce these results. Further, extracting these quantities may provide an efficient way toward setting the physical scale of δt in Minkowski lattice field theory similar to how the Euclidean a is determined by measuring the Sommer parameter, string tension, or a mass scale. While Eq. (12) is correct in the limit of δt → 0, Trotterization process introduces additional interactions, which allow for the mixing of states. This can be seen by comparing the generic leading-order Trotter operator to the exact evolution: By neglecting the commutator terms like [H A , H B ], we are effectively simulating a different Hamiltonian which may have reduced symmetries. This means that even in the case where the projection operators exactly pick out a single state, finite δt calculations of Eq. (12) may contain contamination from other states -potentially even ones with the incorrect quantum numbers.

V.1. Dispersion Relation: single spin
First, we computed the matrix element k|O(t)|k of the states excited by using Eq. (5) for fixed momenta k as a projection operator on N 2 s = 3 2 , 4 2 lattices. These states were then evolved for an approximately fixed T = N t δt for δt = 0.15, 0.2, 0.25, 0.3, 0.4. For each δt, we perform a discrete Fourier transform with fixed δE to obtain an approximation of the spectral function. For each δt, we take the peak with largest spectral weight to be the single particle energy aE k (δt). An example of these results is shown in Fig. 2 for the k = 0 state. The dominant error is from the finite δE below which we cannot resolve the peaks of the spectral function. This error could be reduced by taking a longer T . In order to extract the δt → 0 results, we perform an extrapolation. Since the error from the Hamiltonian and the Trotterization should be O(δt 2 ), we fit the data to the function To investigate the potential for circuit depth reduction by using larger δt, we consider two fits. The first is performed only for the points δt ≥ 0.25 with the result of aE 0 = 0.62 (6), while the second is performed including all the data and finds aE 0 = 0.62 (3). Both of these results are in good agreement with the exact 3 2 value of aE 0 = 0.6204. While obviously using the smaller values of δt provide for reduced uncertainty, they come at the cost of larger circuits. For the same fixed T , we increased our longest circuit depth by a factor of 1.6.
In Fig. 3 where s sums over the spatial directions. We also plot a continuum dispersion relation aE k = (am) 2 + (ak) 2 where am = 0.40(2) is given by extrapolating the N s = 2, 3, 4 results of aE 0 = 0.8701, 0.62(2), 0.50(3) assuming O(a 2 ) corrections dominate, as suggested by [59]. The N s = 2 value was computed by exact diagonalization for simplicity. Comparing am to these values, we can see that for N s ≤ 4 there are substantial finite volume effects at k = 0, which increase for k > 0. Despite this, we find qualitative agreement between our results for the dispersion relation at finite volume and the continuum extrapolated one.

V.2. Two-particle states
Next, we considered the case of two plane-waves scattering by computing the matrix element k, p|O(t)|k, p excited byP =â † kâ † p for N 2 s = 4 2 and δt = 0.05. Our final results are the spectral function shown in Fig. 4 obtained from a discrete Fourier transform.
In order to demonstrate time dependence, we need to compare these spectral functions to the initial state of the scattering plane waves. In the inset of Fig. 4, we present the eigenstate decomposition for the initial state of ak = (0, 1), ap = (1, 0), which we can compare to the final result in the larger figure. We find a change in the relative weight of various eigenstates from t = 0. This indicates that the wavepackets are interacting and that scattering processes can occur. Alas, we do not observe a clean spectrum at t = 0 of a single eigenstate. This, in turn, suggests the two plane wave ansatz for the source and sink does not have as strong overlap with two single particle eigenstates on this lattice. This suggests a need both for larger lattices such that the particle states can be physically separated, and that wavepacket-like excitations should be investigated.

VI. CONCLUSIONS
This work has extended the general methods developed in [15] to extract the matrix elements of particle excitations while still reducing the circuit depth compared to adiabatic state preparation. This has been explicitly shown for the dispersion relation and a two-particle spectral function in the Z 2 gauge theory on small lattices. Similar calculations for this model may be tractable to simulate on quantum computers in the near future, with N s = 2 2 lattices potentially feasible already. We have also observed evidence of scattering and particle interactions. A crucial direction of future work would be to study how other matrix elements, e.g. form factors, can be extracted using these techniques. While the plane wave ansatzes used here was capable of extracting meaningful results, they are likely not an ideal choice for finite size particle states because their overlap with multiple states will lead to signal to noise problems which, unlike Euclidean calculations, are not suppressed at longer times. More optimal choices for projection operators will be required, as well as signal-to-noise mitigation techniques. In the case of scattering states, improved operators are currently being developed in the Euclidean lattice field theory community [60,61] and will be an important avenue of study in the future of the method studied here. Novel techniques in state preparation on quantum devices such as projected cooling [62,63] may also prove useful in reducing excited state contamination.