Simulating $\mathbb{Z}_2$ Lattice Gauge Theory with the Variational Quantum Thermalizer

The properties of strongly-coupled lattice gauge theories at finite density as well as in real time have largely eluded first-principles studies on the lattice. This is due to the failure of importance sampling for systems with a complex action. An alternative to evade the sign problem is quantum simulation. Although still in its infancy, a lot of progress has been made in devising algorithms to address these problems. In particular, recent efforts have addressed the question of how to produce thermal Gibbs states on a quantum computer. In this study, we apply a variational quantum algorithm to a low-dimensional model which has a local abelian gauge symmetry. We demonstrate how this approach can be applied to obtain information regarding the phase diagram as well as unequal-time correlation functions at non-zero temperature.


I. INTRODUCTION
Much interest has been generated recently regarding the potential of quantum computing and quantum simulation to address long-standing problems in high-energy physics [1,2].Traditionally, lattice gauge theory has been applied in the Euclidean formulation, where advances in computing power have led to great success in studying properties of strongly-interacting matter in thermal equilibrium.However, due to the notorious sign problem, there has been a growing interest in the Hamiltonian formulation and its application to address non-zero chemical potential and real-time dynamics.For digital quantum computers in the medium term, one problem that has arisen is how to produce thermal states.
In recent years, there has already been an explosion of progress in tackling this problem.Each approach to producing a thermal state has had to address the issue of how to produce a mixed state.This can typically be done by enlarging the system size and then entangling the ancillary system with the original system.One idea which takes this approach to create a Gibbs state via purification is the so-called thermofield double states method [3].Another interesting method is known as active cooling [4].This takes a general initial state, couples it to a thermal heat bath, and uses the concept of the "Maxwell demon" to obtain a thermal state at the desired temperature.Alternatively, as the Boltzmann operator is nonunitary, one can also attempt to directly approximate its action on quantum states.Using judiciously chosen initial states produced by Haar-random circuits, one can use this approximation of the Boltzmann operator to obtain Gibbs states for spin models [5] as well as gauge theories ˚mfromm@itp.uni-frankfurt.de: philipsen@itp.uni-frankfurt.de; michael.spannowsky@durham.ac.uk § winterowd@itp.uni-frankfurt.de[6].
One can also study the problem of producing thermal states using variational methods.One such method, the variational quantum thermalizer (VQT), introduced in [7], serves as a natural extension of the variational quantum eigensolver (VQE) [8,9] to the case of non-zero temperature.The latter uses the variational principle at T " 0 to approximate the ground-state of a Hamiltonian Ĥ by a parametrized ansatz |ψ ξ ⟩, subject to Similarily, for the case of T ‰ 0, VQT addresses the problem of finding a variational approximation ρvar to the thermal state, ρGibbs " e ´β Ĥ for inverse temperature β " 1{T .Similar to its zero temperature analogue in Eq. (1), this problem can also be reformulated as an optimization problem for the free energy with respect to variational parameters ξ, such that where the second term on the right-hand side involves the von Neumann entropy S " ´trpρ log ρq.As the density operator appears non-linearly in the expression for the entropy and thus is not simply determined by a quantum-mechanical expectation value, several avenues have been taken in the literature: The original formulation of VQT [7] employs classical sampling of the distribution corresponding to ρ ξ which necessitates a model ansatz.Through this ansatz, the entropy is then given classically by a closed-form expression.To prepare the latent distribution corresponding to ρ var in a quantum circuit, the noise-assisted variational quantum thermalizer [10] (NAVQT) uses (simulated) parameterized depolarization gates that allow for the preparation of a mixed state, again with a closed-form, approximate expression for the entropy in terms of the noise level λ, which thus becomes a variational parameter.Of particular interest to us, due to both its simplicity and versatility, is a variant of the VQT introduced in [11], see Fig. 1.
. Schematic depiction of the hybrid approach used in the VQT algorithm.The action of the variational circuits V QC1 and V QC2 and the measurements are performed on the quantum device.The results of the measurements are then sent to the classical device, which is depicted in grey.This information allows one to construct the cost function and propose a new value for the variational parameters using one's optimizer of choice.
While the first two examples used a classical sampling of input states [7] or stochastic mixtures of unitary circuits [10], respectively, to create mixed quantum states, the approach in [11] utilizes intermediate projective measurements.As measurements are non-unitary operations, they act as a source of entropy, ultimately leading to a variational approximation of the thermal state 1 .
The subject of our study is the Hamiltonian formulation of the Z 2 lattice gauge theory (LGT) in 1 `1d ), in which we investigate the application of the VQT to this theory.This is performed in its formulation with a gauge-redundant Hilbert space [13] and its recent, resource-efficient form derived systematically in [14,15].In addition, we investigate several observables of interest at both non-zero temperature T and non-zero chemical potential µ.

II. BACKGROUND A. Variational Quantum Thermalizer
As displayed in Fig. 1, starting from an n q -qubit initial state |ψ 0 ⟩ with density matrix ρ 0 " |ψ 0 ⟩⟨ψ 0 |, we set up a concatenation of two separate variational circuits VQC 1 and VQC 2 .The first variational circuit, VQC 1 , creates the state |ψ VQC 1 ⟩, which is given by An intermediate measurement in the computational basis collapses this state and statistically yields the latent density matrix ρ 1 pθq " 1 For a version of this algorithm using entanglement of the system with nq ancillary qubits followed by a projective measurement on the latter, see [12].
The subsequent application of the variational circuit VQC 2 transforms this into ρ 2 pθ, ϕq " U 2 pϕqρ 1 pθqU : 2 pϕq, " after which a measurement of the energy is performed.
The von Neumann entropy is then calculated as S " ´ i p i log p i , with p i « n i {n s , where n i represents the counts for state i in the computational basis and n s represents the number of shots.By varying the parameter set ξ " pθ, ϕq in an optimal way, one obtains a variational approximation to the Gibbs state ρGibbs " e ´β Ĥ .The fact that the circuit is broken up into two parts is intended for noisy intermediate-scale quantum (NISQ) devices, as the depth of each variational circuit can be kept relatively shallow.This approach, however, naturally leads to a truncation in the number of states in Eq.( 5), as both the number of shots, n s , and the memory necessary to store the histogram describing the counts scales with the dimension of the Hilbert space of the underlying problem.We pursue this discussion in further detail in Sect.III B 4. In [11], the algorithm was demonstrated to reach convergence for the transverse field Heisenberg-Model in d " 1, 2. Applying qVQT to gauge theories will introduce the need for a gauge-invariant algorithm implementation.We show below for the case of Z 1`1 2 how this can be achieved.

B. Z2
LGT in 1 `1d This section briefly outlines our two approaches to Z 2 lattice gauge theory in the Hamiltonian formulation.In the first approach, introduced in [13] and referred to as the gauge-redundant formulation, local gauge invariance must be enforced on the states to define the physical Hilbert space.This can be achieved, for example, by introducing a penalty term as in [6].In the second approach, used in [14] and referred to as the resourceefficient formulation, the matter states are decoupled, LGT with staggered fermions on a lattice with linear extent L. The links are all in the eigenstate of X with eigenvalue `1, the even sites are empty (fermions), and the odd sites are filled (anti-fermions).We note here that the links at the boundary are set to unity.
Depiction of a single-layer of a variational quantum circuit in the gauge-redundant formulation.The layer is composed of one sublayer of single-qubit rotation gates and a sublayer containing the decomposed three-qubit gate which is formed from the gauge-matter term in the gauge-redundant formulation.All relevant gates contain variational parameters ∆i. and the constraint imposed by Gauss' law is built into the Hamiltonian.

Gauge-Redundant Formulation
In the so-called group basis, Z 1`1

2
LGT with a single flavour of staggered fermion is formally given by the Hamiltonian H " H E `HGM `HM `Hµ with In the above equations, the fermionic creation and annihilation operators have been replaced by Pauli matrices living at the sites of the lattice via the usual Jordan-Wigner transformation.Here we have adopted the familiar convention of denoting Pauli operators acting on the link Hilbert space of the gauge bosons by {X, Y, Z}, and those acting on the fermionic matter space by {σ x , σ y , σ z }.The occurrence of both gauge-bosonic and fermionic degrees of freedom implies a gauge redundant Hilbert space.With our choice of basis, a gauge transformation at site n is defined by [16] Θ g pnq " X n X n´1 e iπ[ψ : n ψn`1 2 pp´1q n ´1q] , We illustrate a gauge-invariant state |ψ⟩ of the lattice system in Fig. 2, where the fermionic states are in the computational basis, and the gauge bosons are in the Xeigenbasis, X |˘⟩ " ˘|˘⟩.This subspace of the total Hilbert space is needed to compute physical observables.

Resource-Efficient Formulation
It turns out that one can eliminate the fermionic matter fields while simultaneously removing the Hilbert space redundancy of the theory.This is accomplished in a twostep process whereby the fermions are first mapped to hard-core bosons by a unitary transformation and are then decoupled from the gauge bosons by a second unitary transformation.This method is, in fact, a very powerful formalism which can also be applied in arbitrary dimensions to continuous gauge groups such as SU pN q and U pN q [17].For Z 2 LGT in 1 `1 dimensions, the Hamiltonian in this formulation is as follows One notices that the Hamiltonian obtained by performing the above-mentioned transformations only vaguely resembles the original 1 `1-dimensional Kogut-Susskind Hamiltonian.For example, the gauge-matter interaction in eq. ( 12) couples sets of three adjacent links.Similarly, for the term originating from the staggered mass in Eq. ( 13), one obtains a coupling between two adjacent links.With this resource-efficient formulation, one can simulate larger systems without worrying about the elements of the variational circuit preserving gauge invariance.

III. IMPLEMENTATION AND RESULTS
The Z 2 lattice gauge theory has a chiral symmetry which is spontaneously broken at low T .At large temperatures as well as at large values of the chemical potential, this symmetry is restored.To benchmark our approximation to the Gibbs state, we calculate the chiral condensate, ⟨ ψψ ⟩, which is sensitive to restoring chiral symmetry.As this is a local, time-independent observable, one would like a second, time-dependent observable that measures correlations in the system.We chose to study the unequal-time density-density correlator, which provides information about charge dynamics in the system.

A. Gauge redundant formulation
In this formulation of the theory, gauge invariance has to be considered explicitly.Apart from the gaugeinvariant initial state (c.f.Fig. 2), this not only applies to the variational circuits VQC 1{2 but also to the intermediate measurement.A schematic depiction of the circuit illustrating the above-mentioned components is shown in Fig. 1.For VQC 1{2 , a gauge-invariant gate set can be constructed by following the general considerations laid out in [18].This consists of simply taking the Pauli strings in our Hamiltonian Eq.( 6) as elementary building blocks.Due to the simplicity of the model, the gate set will hence decompose into parameter-dependent single qubit R x and R z gates, together with a 3-qubit gate " e iσ x Zσ x `iσ y Zσ y stemming from the hopping term in Eq. (6).Following the procedure presented in [19], the latter can be easily decomposed into single-qubit and entangling gates.In Fig. 3, one observes the resulting circuit layer consisting of a sublayer of single-qubit gates along with a sublayer of the decomposed multiqubit gate for a three-qubit system where the ∆ i represents generic variational parameters.Several such layers then yield the VQC i , with parameter sets θ and ϕ for i " 1, 2, respectively (Eq.( 5)).With the composition of the variational quantum circuits at hand, we still need to define the figures of merit quantifying the quality of the optimization step depicted in the classical block of Fig. 1.In addition to the relative difference in free energy, ∆F " F ρvar {F ρ Gibbs ´1, we take the fidelity [20] as a distance measure between two quantum states represented by ρ 1 and ρ 2 , respectively.In Fig. 4, one sees the results of the classically simulated variational optimization process for our gauge-redundant theory for N " 4 matter sites corresponding to n q " 7 qubits using open boundary conditions with a two-layer ansatz in the VQC i .Simulating at aϵ " am " 0.5, aµ " 0 and n s " 10 4 for varying T {m, chiral symmetry is explicitly broken in our model.On the left, we plot the performance measures defined above, where each point corresponds to the minimum in free energy over up to 10 2 samples, resulting in optimal parameter sets pθ ˚, ϕ ˚q.The quality of the approximation by the variational ansatz decreases slightly in the crossover region until, with increasing T {m, the entropic term in Eq.( 2) becomes dominant, yielding a high-fidelity solution.On the right, the chiral condensate ⟨ ψψ⟩ evaluated using ρ 2 pθ ˚, ϕ ˚q is shown as a function of T {m.The displayed errors were obtained from the variation of the observable over solutions corresponding to the 20th percentile of the free energy distribution for each set of samples.

B. Resource efficient formulation
Turning now to the Hamiltonian Eqs.(11) with matter degrees of freedom eliminated and gauge invariance builtin, the choice of variational ansatz circuits for the VQC i is in principle unrestricted, leaving us with the freedom to use, e.g.hardware efficient circuits such as EfficientSU2 from qiskit's standard circuit library [21].As for the set of physical gates, by inspecting the Hamiltonian one sees that the single and multi-qubit terms can be employed as elementary ansätze.As an exaample, the decomposed three-qubit term is depicted in Fig. 5.An alternation of a single-qubit sublayer containing R x and R y gates, followed by a staggered multi-qubit sublayer, will hence serve as our variational ansatz for both VQC 1{2 .

Observables at T, µ ą 0
As was done in the gauge-redundant formulation, we quantitatively investigate the quality of our variational approximation to the Gibbs state and calculate the chiral condensate.This has been done as a function of temperature T {m and varying chemical potential µ{m at ϵ{m " 0.5 with n s " 10 4 .Owing to the resource efficiency of the formulation, n q " 7 classically simulated qubits now correspond to twice the system size (N " 8) when compared to the gauge-redundant formulation.While simulations for varying T {m across the cross-over region, Fig. 6, show good agreement between exact and simulated results, entering the transition region for fixed temperature and increasing µ{m, we see that our algorithm has difficulty to match the exact result.One notes that as we lower the temperature, the transition becomes more and more discontinuous.In Fig. 8  variation of ⟨ ψψ⟩ over the top 20th percentile of variational solutions.We expect an increase in performance can be achieved here by changing from the gradient-free optimization algorithm (COBYLA) used in this study to gradient-based methods.This is an approach which we have not pursued further.

Thermal Unequal-Time Correlators
In addition to the study of quantum systems at finite fermion density, quantum simulation in the NISQ era and beyond will allow the investigation of real-time observables.These include thermal unequal-time correlation functions.In general, an unequal-time correlator of two operators takes the form where H is the Hamiltonian of the system and Aptq is the operator in the Heisenberg picture.In this study, the expectation value in ( 16) will be taken with respect to our variational density matrix, ρ 2 pθ ˚, ϕ ˚q.
The question of how one measures such a correlation function naturally arises when studying lattice gauge theories.Fortunately, a general procedure has already been devised to measure the correlation function of two arbitrary unitary operators A and B on a universal quantum computer [22,23].This procedure is commonly referred to as Ramsey interferometry and only involves adding a single ancillary qubit entangled with our system's original n q qubits.The time evolution can be performed with a Trotterization of the time evolution operator whose circuit depth increases linearly with the number of steps.A variational approach also exists to obtain the time evolution of an arbitrary quantum state with a fixed circuit depth [15].In this work, a simple Trotterization has been used.
In Fig. 9 we show the unequal-time thermal densitydensity correlation function ⟨O 0 p0qO i ptq⟩ at T {m " 1 for N " 8, where O i " ε i Z i´1 b Z i .The different data sets represent different spatial separations i " 1, . . ., 5 of the density operator on the lattice plotted as a function of time with a remarkable agreement between simulated and exact results for the time range displayed.In these classical simulations of our variational quantum circuit, both VQC 1 and VQC 2 contain two layers of alternating single and multi-qubit sublayers.For simplicity's sake, we have assumed infinite statistics in the computation of ρ 2 pθ ˚, ϕ ˚q, hence excluding shot noise.One is ultimately interested in obtaining transport coefficients from realtime dynamics.As the latter involves time integrals of unequal-time correlators, it remains to be seen if such simulations yield sufficient accuracy.

qq meson screening
It is useful to check that our algorithm correctly reproduces physical observables, which are more easily accessible by standard Eulcidean methods.As an example, we consider the thermal expectation value of a gaugeinvariant qq meson operator, i.e. a quark anti-quark pair connected by a flux tube [14,24], We have written the above operator in the gaugeredundant version for simplicity but have used the corresponding resource-efficient expression in our calculations.As demonstrated recently with the help of DMRG [25], at T " 0, this equal-time two-point function shows an exponential decay for any ϵ ą 0. The exponents of a spatial decay are screening masses, which in this case give the energy of the flux tube and its excitations, whose linear increase with separation indicates confinement.We expect this behaviour to persist for small temperatures until the onset of string breaking, when the flux gets increasingly screened and its energy no longer grows with separation.
Using the same variational setup as described in Sect.III B 2, Fig. 9 (right) shows |⟨M i ⟩| at T {m " 0.5 and T {m " 1.25.Whereas the low-temperature correlator exhibits marked similarities to the zero-temperature decay (c.f.[25]), which is accurately captured by the variational Gibbs state, at higher temperatures, the algorithm fails to reproduce the exact correlator accurately, possibly due to the low complexity of the ansatz circuits in our simulations (two layers of variational gates for the VQC i ).While this precludes a detailed picture of the long distance behavior with increasing temperature, our variational calculation does exhibit two qualitative features known from other techniques: i) the screening masses generally grow with temperature as expected (faster decay); ii) the difference between high and low temperature is small at short distances, but pronounced at large distances, because of the flux screening.

Resources and Entropy Estimation
We now turn to the question of resources in runtime and classical memory.Considering the limited system sizes currently studied in the NISQ era, estimating the von Neumann entropy S " ´ i p i log p i from the intermediate measurements via p i « n i {n s is a viable approach.Applying the algorithm with the information flow depicted in Fig. 1 for realistic system sizes will necessarily involve a certain level of approximation.This stems from the fact that the number of p i 's grows with the dimension of the Hilbert space, 2 nq .One possible way proposed in [11], is the partitioning of the full n q qubit system into n q {n ss independent subsystems of size n ss in the entropic part of the algorithm (VQC 1 ) that prepares the latent distribution ρ 1 pθq.In the version of the algorithm which employs a classical estimation of the entropy [7] this corresponds to a particular choice of the latent distribution known as factorized latent space models.It remains to be seen if such an approximation (Left) The unequal-time thermal correlation function ⟨O0p0qOiptq⟩ for varying distance i " 1, . . ., 5 (shifted in magnitude for displaying purposes) at T {m " 1 for N " 8, where Oi " εiZi´1 b Zi, i.e. the relevant term of the number operator.The solid lines are the exact solutions while the crosses represent the values calculated with the Trotterized timeevolution operator (aδt " 0.2) for the variatonal Gibbs state ρvar produced with infinite statistics.(Right) The absolute value of the meson string operator as a function of the string length (lattice units).The results for two temperatures are plotted using both the variational Gibbs state and the exact density matrix.One notices that at larger temperatures the current variational ansätze give states which deviate from the exact result.
is a valid approach for gauge theories where thermalization and subsystem entanglement are non-trivially connected [26].
There are, however, alternative methods that can be used to estimate the entropy efficiently.For example, it turns out that an evaluation of the Taylor series approximation to the entropy can be useful [27].This method requires the addition of n q ancillary bits to perform a qubitefficient, repeated SWAP-test evaluation of the trace of various powers of the latent density matrix, trpρ n 1 pθqq.The results of this approach for the simple gauge theory used in this study are displayed in Fig. 10.For temperatures up to T {m " 1, only a few orders are needed to give reasonable values of the fidelity.

IV. CONCLUSION
We have explored the use of variational quantum algorithms in producing thermal states.These approaches were applied to Z 2 lattice gauge theory in 1 `1 dimen-sions using two independent formulations, and several observables of interest were computed.We find this variant of the VQT a suitable algorithm to approximate the thermal states of quantum systems in the NISQ era as estimating the entropy of larger systems will necessitate the introduction of truncations in the estimation.A particular aspect we have left for future work is the performance of intermediate measurements for mixed-state creation on actual quantum devices compared to the usual approach of entangling ancillary systems followed by a terminal measurement.A further possible extension of our work is the formulation of the VQT to study non-Abelian gauge theories at finite temperatures.As the latter will generally come with a gauge-redundant Hilbert space, even after reformulating in a resource-efficient manner, this amounts to respecting gauge invariance or equivalent local constraints throughout the variational quantum circuit. .The fidelity of the Gibbs state produced using the Taylor-series approximated von Neumann entropy as a function of T {m at aµ " 0 for various truncation orders K on a system with N " 8 matter sites.One sees that for low temperatures this approximation works well.

Figure 2 .
Figure 2. Depiction of a gauge-invariant state for Z 1`1 2

Figure 4 .Figure 5 .
Figure 4. (Left)The fidelity of the approximation to the Gibbs state produced by VQT on the classical simulator and as a function of T {m at µ " 0 on a system with N " 4 matter sites, corresponding to nq " 7 qubits in the gauge-redundant formulation.(Right) The chiral condensate ⟨ ψψ⟩ calculated on the same data sets where for clarity the exact result has also been included.

Figure 6 .Figure 7 .
Figure 6.(Left) The fidelity of the approximation to the Gibbs state produced by VQT on the classical simulator and as a function of T {m for µ " 0 on a system with N " 8 matter sites, corresponding to nq " 7 qubits in the resource-efficient formulation.(Right) The chiral condensate ⟨ ψψ⟩ calculated on the same data sets.

Figure 8 .
Figure 8. (Left) The fidelity of the approximation to the Gibbs state produced by VQT on the classical simulator and ∆f as a function of µ{m for T {m " 0.4 on a system with N " 8 matter sites, corresponding to nq " 7 qubits in the resource-efficient formulation.(Right) The chiral condensate ⟨ ψψ⟩ calculated on the same data sets.

8 Figure 10
Figure 10.The fidelity of the Gibbs state produced using the Taylor-series approximated von Neumann entropy as a function of T {m at aµ " 0 for various truncation orders K on a system with N " 8 matter sites.One sees that for low temperatures this approximation works well.