Real-time dynamics in 2+1d compact QED using complex periodic Gaussian states

We introduce a class of variational states to study ground state properties and real-time dynamics in (2+1)-dimensional compact QED. These are based on complex Gaussian states which are made periodic in order to account for the compact nature of the $U(1)$ gauge field. Since the evaluation of expectation values involves infinite sums, we present an approximation scheme for the whole variational manifold. We calculate the ground state energy density for lattice sizes up to $20 \times 20$ and extrapolate to the thermodynamic limit for the whole coupling region. Additionally, we study the string tension both by fitting the potential between two static charges and by fitting the exponential decay of spatial Wilson loops. As the ansatz does not require a truncation in the local Hilbert spaces, we analyze truncation effects which are present in other approaches. The variational states are benchmarked against exact solutions known for the one plaquette case and exact diagonalization results for a $\mathbb{Z}_3$ lattice gauge theory. Using the time-dependent variational principle, we study real-time dynamics after various global quenches, e.g. the time evolution of a strongly confined electric field between two charges after a quench to the weak-coupling regime. Up to the points where finite size effects start to play a role, we observe equilibrating behavior.


I. INTRODUCTION
Gauge theories are of paramount importance in fundamental physics. Its most prominent example, the standard model of particle physics, describes electromagnetic, weak and strong interactions. In some regimes, interactions can be treated in terms of perturbative expansions. However, since the coupling in quantum field theories is typically scale-dependent, there are regimes (e.g. lowenergy QCD) where non-perturbative methods are required [1,2].
Lattice gauge theory is a gauge invariant lattice regularization of gauge theories, in which either spacetime [3] or space [4] is discretized. This has allowed to uncover many interesting features of non-perturbative quantum field theories, in particular using Monte-Carlo simulations [5]. Nevertheless, certain aspects are difficult to study within this framework, e.g. fermionic theories with finite chemical potentials may suffer from the sign problem [6] and time dynamics are difficult to access as Monte-Carlo simulations require a formulation in Euclidean spacetime.
One class of approaches to these problems is based on a Hamiltonian formulation of lattice gauge theories, first proposed by Kogut and Susskind [4]. Other formulations in the Hamiltonian picture include the quantum link model [7][8][9][10] or the prepotential approach [11]. It has been shown that these Hamiltonians or truncations [12] thereof can be mapped to Hamiltonians of quantum devices (e.g. ultracold atoms, trapped ions or superconducting qubits) in order to study such theories by quantum simulation [13][14][15][16]. Another option is to study the Hamiltonian by designing appropriate variational ansatz states which are both efficiently tractable and capture the most relevant features of the theory.
Both ideas have been successfully applied to one-dimensional theories. The implementation of quantum simulators has been demonstrated using trapped ions [17] and ultracold atoms [18][19][20][21]. On the numerical side, there has been a lot of success in applying matrix product state (MPS) methods to (1+1)-dimensional Abelian and non-Abelian lattice gauge theories [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36], enabling the study of finite chemical potential scenarios and out-ofequilibrium dynamics which would not have been accessible in Monte-Carlo simulations of Euclidean lattice gauge theory. Also some generalizations of Gaussian states have proven to be suitable for these theories [37]. The situation becomes more challenging in higher spatial dimensions, in particular due to appearance of magnetic interactions, leading to four-body plaquette terms on the lattice. There have been ideas on how to overcome this problem in quantum simulators (either by employing a digital [38][39][40][41] or an analog simulation scheme [42]) but so far they are out of experimental reach. On the numerical side, tensor network methods in 2+1d have been applied to pure gauge theories [43] and for studying U(1) ground states in quantum link models [44,45]. It has also been shown that fermionic Gaussian PEPS can be gauged [46] and serve as numerical ansatz states for lattice gauge theories, admitting a sign-problem free Monte-Carlo contraction scheme [47].
In this work, we study (2+1)-dimensional compact quantum electrodynamics (compact QED). It is a good starting point for the study of higher dimensional lattice gauge theories since it shares some features with (3+1)-dimensional Quantum Chromodynamics (QCD), e.g. that it is in a confined phase for all values of the coupling constant [48]. To access physics which is difficult to simulate with Monte-Carlo simulation of Euclidean lattice gauge theories, we not only study ground state properties but also non-equilibrium physics, namely real-time dynamics after a quantum quench.
Since exact diagonalization methods become infeasible in higher dimensions for reasonable system sizes, in particular due to the infinite local Hilbert space of the gauge field, it seems unavoidable to use variational techniques (in 1+1d the infinite dimension can be avoided either by integrating out the gauge field nonlocally [29,49,50] or by using the natural restriction of gauge symmetry which makes the dimensions finite [51]).
We choose to work with complex periodic Gaussian states, a generalization of periodic Gaussian states, first proposed in [52] to prove confinement in the weakcoupling limit of 2+1d compact QED, thus establishing the existence of one confining phase for all couplings also in the Hamiltonian picture (after it had been proven in the action formalism [48]). As expectation values with respect to periodic Gaussian states cannot be evaluated analytically, the authors of reference [52] used Feynman diagram techniques to evaluate all relevant quantities in the weak-coupling regime. In contrast to that approach, we develop a numerical approximation scheme to evaluate these states for the whole coupling region. By extending the variational manifold to complex periodic Gaussian states we are also able to account for real-time dynamics. One appealing feature of these states is that they do not require any truncation in Hilbert space which allows us to study truncation effects which are required in other approaches and give estimates in which coupling regimes they are justified.
The manuscript is structured as follows: In Sec. II, we introduce the model and the variational ansatz including a scheme for its numerical evaluation. In the first part of Sec. III, we study ground state energy density and string tension over the whole coupling region. In the second part, we investigate truncation effects by comparing the variational ground state energy with exact diagonalization results where the local Hilbert space is truncated in the electric basis. In Sec. IV, we study real-time dynamics after a quantum quench using the time-dependent variational principle. In Sec. V, we conclude.

A. (2+1)-dimensional compact QED
We define the theory of (2+1)-dimensional compact QED on a square lattice of extent L × L with periodic boundary conditions. The gauge fields reside on the links; U x,i denotes the gauge field operator on the link emanating from site x in direction e i . The Hamiltonian in lattice units takes the following form, originally proposed by Kogut and Susskind [4]: with g 2 being the coupling constant and U p ≡ U x,1 U x+e1,2 U † x+e2,1 U † x,2 where x is the bottom left corner of plaquette p. U x,i is in the fundamental representation of U (1), it can also be written in terms of an angle θ x,i , U x,i = e iθx,i with −π < θ x,i ≤ π. The restriction of the gauge field to this compact interval is the reason why the model is called compact QED and why it exhibits interesting features such as confinement in contrast to the non-compact theory [53]. E x,i is the electric field operator fulfilling the following commutation relations: Since we work in the temporal gauge, there is a residual spatial gauge symmetry defined by the Gauss law operators G x . All physical states need to be eigenstates of them: where the eigenvalue Q x gives the static charge configuration at x. These local constraints put quite severe restrictions on the choice of variational states. Following [52], we thus change to variables where gauge invariance is already incorporated (at least up to a global constraint). This can be achieved by splitting the electric field E x,i into its transversal part E T x,i , which is dynamical, and a longitudinal part E L x,i which is fixed by the static charge configuration. Since the transversal part of the electric field can be expressed by a plaquette field L p (the lattice analogue of a solenoidal vector field), the remaining dynamical degrees of freedom {L p , U p = e iθp } reside on plaquettes, having the same Hilbert space structure and fulfilling the same commutation relations as the link variables: The operator U p creates an electric flux excitation around plaquette p. However, to construct all possible gauge-invariant flux configurations two global noncontractible flux loops around the torus (one for each spatial direction) are required, their operators are denoted as {θ 1 , L 1 } and {θ 2 , L 2 } specifying the topological sector of the flux configuration. L 1 and L 2 commute with the Hamiltonian and we will restrict ourselves to the topological sector with L 1 = L 2 = 0 which corresponds to no electric flux loops winding around the torus. For more details see [54] or Appendix A. Writing the Hamiltonian in terms of these new variables, reads where E C is an energy offset given by the lattice Coulomb energy and p accounts for the transversal part of the electric field caused by the static charges only, i.e p = 0 in case of no static charges. Even in this formulation there is one remaining global constraint left which is intuitively clear since raising the electric flux around all plaquettes should return the same state due to the periodic boundary conditions. Thus, For details on this formulation, we refer the reader to [52,54]. A rigorous derivation of eq. (5) from eq. (1) and an explicit formula for the calculation of p and E C can be found in Appendix A.

B. The variational ansatz
We formulate our variational ansatz states in terms of the θ p -variables defined above such that it only needs to fulfill the global constraint (6). Starting from periodic Gaussian states introduced in [52], we extend the variational wavefunction to have an imaginary part in order to account for real-time dynamics. The ansatz is based on a complex Gaussian state: with x p ∈ R and p = (p 1 , p 2 ), p 1 , p 2 ∈ [0, .., L − 1]. The linear part in the exponent, i.e. p , is fixed by the static charge configuration (see section II A and Appendix A) and is defined by the variational parameters γ R k and γ I k . In the following, we will use the shorthand notation pk ≡ 2π p1k1+p2k2 L . Since the disorder introduced by static charges is incorporated in p , the quadratic part A is assumed to be translationally invariant. The factor of 1/π is chosen for later convenience.
Written in terms of Fourier components x k = 1 L p e ipk x p , the quadratic part in the exponential be- Thus, to guarantee convergence of Ψ CG we need to require γ R k > 0 ∀k. Since |x k | 2 = |x −k | 2 , the variational parameters γ R/I k and γ R/I −k are redundant. We define the equivalence relation and With the quotient set K ≡ [0, .., L − 1] 2 \ (0, 0) /∼ k we can define a set of independent variational parameters, . Choosing a set of independent parameters will be important later on for applying the time dependent variational principle (see section IV A).
To construct a suitable ansatz state for compact U (1) gauge fields (θ p ∈ [−π, π]) we sum over complex Gaussian states, thus ensuring periodicity: The delta function needs to be included in order to satisfy condition (6) for physical states. To shorten notation, we will denote the product over infinite sums p +∞ Np=−∞ by {Np} and the product over integrals p π −π dθ p by π −π Dθ. The Gaussian nature of the wavefunction is exploited when evaluating expectation values of observables O by combining the integral over 2π with one of the two infinite sums to an integration over the real axis The integral f O ({N p }) can be carried out analytically and the remaining infinite sum needs to be evaluated numerically.
Exemplary, we show this procedure for the norm of the variational state, Ψ CP G |Ψ CP G . The computation of observables follows analagously; details on their exact form can be found in Appendix C. After carrying out the integrals, the remaining function The γ k parameters determine how fast contributions to the sum in eq. (11) decrease exponentially with increasing |N k | 2 .
We group the configurations N p of this sum in different orders such that within one order the configurations only change up to permutations. Since all relevant configurations will contain mostly zeros, we will denote orders by its non-zero elements, e.g. {N } 1 is the set of all permutations of the configuration N defined by N p=0 = 1 and N p =0 = 0, i.e. {N } 1 ≡ S N . If the parameters γ k are large enough, the sum can be approximated by orders having small Euclidean norm, The higher number of permutations in orders with larger norm cannot compensate for the exponential suppression (this would not be the case if the γ k were arbitrarily small). Using this scheme, the constraint δ p N p is useful since it excludes many orders, e.g. {N } 1 or {N } −1 . The order with the lowest non-zero norm is therefore {N } 1,−1 . In fact, the sum in eq. (11) can be expanded in orders containing only pairs of 1, −1: {N k=0 =0} denotes the sum over the set of all N p configurations with N k=0 = 0, i.e. fulfilling the global constraint. For sufficiently large γ k higher orders of the type {N } 2,−2 or {N } −2,1,1 are exponentially suppressed as well as orders with a large number of 1, −1 pairs. Thus, the above expansion can be truncated after the first few terms. Each of the remaining orders is evaluated numerically. The fact that configurations only change up to permutations within one order can be used to highly parallelize the computation. On an 8 × 8 lattice we are able to compute the first three orders exactly. This procedure is sufficient for most configurations of variational parameters with γ k 1. However, in the intermediate regime γ k ≈ 1 more orders are required to obtain good convergence. In these cases, higher orders are computed using uniform sampling. Since for all our purposes the different γ k parameters were of the same order of magnitude and the N p configurations only change up to a permutation within an order, a uniform probability distribution is a suitable ansatz for the exponential in eq. (13). This is only the case for sampling within one order; it would fail if one tried to sample the whole sum. This combined approach of exact evaluation and uniform sampling has the advantage that it introduces almost no error for most of the variational manifold (up to truncated orders which are exponentially suppressed) and even for regions where uniform sampling is required the error is still suppressed since it only occurs in higher orders. For a detailed error analysis due to truncating orders and uniform sampling see Appendix B.
When the γ k become small, the above approximation fails. In that case, one can exploit the fact that Ψ CP G |Ψ CP G can be written as a multidimensional Riemann theta function [55] which is defined as where z ∈ C g , Ω ∈ C g×g , such that Ω = Ω T and Im(Ω) is strictly positive definite. To bring Ψ CP G |Ψ CP G into this form one can rewrite the delta function as the limit of a Gaussian and exchange the limit with the infinite sum due to uniform convergence. One can now exploit invariance of the Riemann theta function under modular transformations, in particular the following relation holds (for details see [55]): If we insert this relation and take the limit, we obtain: with γ −1 0 = 0. The exponential weight depends now on γ −1 k which allows in principle to approximate the sum with only a very limited number of orders for sufficiently small γ k . However, the sum is not well defined since all constant configurations N p = c(1, 1, ..., 1) have weight one for c ∈ Z. Fortunately, since all f inv,O ({N k =0 }) are independent of N k=0 (as a result of the global constraint on physical states), all these configurations can be factored out such that they cancel when calculating expectation values. This can be formulated rigorously by defining an equivalence relation for N p configurations: When calculating expectation values of observables only a sum over representatives of this equivalence relation is required: If we choose the representative to be the one closest in norm to the N p = 0 configuration, we can expand the sum again in orders having mostly 0's. In this case we have no constraint so that all orders must be taken into account. For more details see Appendix B.
A nice way to check the validity of both numerical approximation schemes presented above is to see whether they agree in the parameter region γ k ≈ 1. This check has been carried out throughout this work since it also indicates that the whole variational manifold can be accessed which is required in order to study the whole coupling region.
To illustrate that both approximation schemes complement each other, we give the variational energy of Ψ CP G with respect to the Kogut-Susskind Hamiltonian given in eq. (5), written both in the infinite sum representation for high and for low γ k : The brackets denote the infinite sums: If we set γ I k = 0 ∀k the expressions for high γ R k , i.e. with the sums {N k=0 =0} , agree with the results given in [52] up to redefinitions. It is important to emphasize that the convergence of infinite sums is determined by For real-time evolutions, e.g. a quantum quench, γ I k 2 will typically become large and so will γ k , irrespective of the real part γ R k . This allows to truncate the expansion in eq. (14) already after the first term such that everything can be evaluated without resorting to sampling. This property makes the ansatz well suited for real-time evolution compared with other methods where sampling at all times often makes it difficult to reach long times.

III. STATIC PROPERTIES
In this section, we study the variational ground state of 2+1d compact QED over the whole coupling region. To minimize the energy we applied a gradient descent algorithm (the formula for the gradient can be found in Appendix C). We used different initial seeds to prevent the possibility of getting stuck in local minima. To make sure that our variational state can approximate the ground state, we compare it first to known exact results. One should note that exact diagonalization methods cannot be applied to the full theory since the local Hilbert space is infinite. However, for the case of a single plaquette exact analytical solutions are known, namely the Mathieu functions.

A. Benchmark for one plaquette
For benchmarking our variational ansatz, we will restrict ourselves to the sector without static charges. The Hamiltonian given in the formulation of the previous chapter, written in the basis of θ, reads: The corresponding Schroedinger equation for ξ(θ) can be written as a Mathieu equation: andξ(z) ≡ ξ(θ/2).ξ is therefore not 2π-periodic but π-periodic. The π-periodic solutions are usually separated into even ce 2r (z, q) (r ≥ 0) and odd se 2r (z, q) (r ≥ 1) solutions. The lowest energy, i.e. the lowest characteristic value a, corresponds to the solution ce 0 (z, q). In Fig. 1, this exact ground state energy is plotted against the minimized variational energy. They agree very well over the whole coupling region, even in the regime where the difference is maximal (g 2 ∼ 0.7) the relative error is still around 0.5%.

B. Ground state properties
In this section, we study the properties of the varational ground state for an extended lattice and investigate its finite size effects. We start by studying the ground state energy density e 0 (L) for lattice sizes up to 8 × 8 plaquettes without static charges. We see that for couplings g 2 1.0 this size is already enough to get a linear scaling with 1 L 2 . The thermodynamic limit e 0 (L = ∞) is then extracted with the following fit For large couplings the thermodynamic limit can be reached with even smaller lattice sizes. The region which limits the evaluation of our variational state to 8 × 8 is around g 2 ∼ 1.1 since the variational parameters are of order one (γ R k ∼ 1 , γ I k = 0) and thus both approximation schemes agree (see Appendix B). Hence, for couplings below this transition region we can simulate larger lattices, namely 14 × 14 for g 2 = 0.8, 0.9 and 20 × 20 for 0.1 ≤ g 2 ≤ 0.7. For such lattice sizes, the finite size effects become again small enough to extrapolate to the thermodynamic limit. The result for the ground state energy density in the thermodynamic limit over the whole coupling region is shown in Fig. 2. To illustrate, we show the extrapolation to the thermodynamic limit for g 2 = 0.5 and g 2 = 2.0 in Fig. 3.
In the next step, we study the string tension over the whole coupling region. We can measure it in two ways: First, we place static charges and analyze the scaling of the ground state energy depending on the distance between static charges. We will fit the potential with the following function: where σ is the string tension and V Coul is the lattice Coulomb potential in two dimensions which becomes a logarithmic potential in the continuum limit. The values for V (d) are computed as the difference between the ground state energy with static charges separated by a distance d and the ground state energy without static charges. Exemplary, we show the fit of the potential for g 2 = 2.0 in Fig. 4.
In the second approach we use the scaling of spatial Wilson loops to extract the string tension. This works at zero temperature since on the Euclidean lattice spatial and temporal Wilson loops are related by O(4) symmetry. At finite temperature this symmetry is broken due to a compactified temporal dimension [56]. The formula to calculate Wilson loops of arbitrary size with complex periodic Gaussian states in both the low and high γ k approximation can be found in Appendix C. On 8 × 8 lattices, we consider all rectangular loops R 1 × R 2 with R 1 , R 2 ≤ 4 (four is the maximal physical length due to the periodic boundary conditions). Furthermore, we require |R 1 − R 2 | ≤ 1 to avoid additional finite size effects coming from an asymmetry in the edges. For weak couplings, where larger lattices are accessible, we extend the allowed maximal edge length to 7 and 10 (for 14 × 14, resp. 20 × 20). We fit the Wilson loop scaling according to the following formula:  The first term corresponds to area law scaling with string tension σ and the second term to perimeter law scaling.
To illustrate the procedure, we show the fit for the ground state at g 2 = 0.5 in Fig. 5. We also tried to extract the string tension via Creutz ratios [57] but the results were less reliable than the Wilson loop fits. The result for both approaches is shown in Fig. 6. For large values of the coupling constant, the fit for the static potential works well and agrees with the strong-coupling prediction g 2 2 . Since a large coupling implies a signifi- Since both methods complement each other we chose to make the string tension data for the static potential transparent for couplings g 2 ≤ 1.5 and the ones extracted by Wilson loops scaling for g 2 > 1.5. The remaining full data points in Fig. 6 are the most reliable estimates for the string tension. For small couplings an exponential decay of the string tension is expected according to the formula [58]: If we fit this formula to the string tension data of the Wilson loop fits between 0.5 ≤ g 2 ≤ 0.9 we obtain c = 23.53 and ν 0 = 0.318 which is close to the theoretical prediction (ν 0,theo = 0.321) [59].

C. Truncation effects
Since our wave function does not require a truncation, we can study truncation effects of other methods. Here, we will focus on a truncation in the electric basis. To see these effects we will study the variance of the electric field operator. For simplicity, we will look at this effect without static charges, since they only introduce -shifts (−1/2 < < 1/2) in the electric field. Since the expectation value of the electric field vanishes in the absence of static charges, we can write the variance in terms of the electric energy The variance is plotted in the inset of Fig. 8 for the ground state which was computed in the last section. To quantitatively show the difference, we compare our variational state to an exact diagonalization calculation of a Z 3 lattice gauge theory. To reduce the required Hilbert space dimension, we formulate it in terms of plaquette variables, in the same style as we did for the U (1) theory. The Hilbert space is truncated in the eigenbasis of L p to three states (corresponding to the eigenvalues m = 0, 1, −1). To make this a consistent theory we define the gauge field operators cyclically: This is equivalent to a Z 3 lattice gauge theory formulated in link variables: x is the vertex at the bottom left corner of plaquette p and Q x,i the cyclic raising operator of the electric field on link (x, i), such that (see [60] for details) The maximal lattice size we can achieve in our ED calculation for a reasonable amount of time is 3×3 plaquettes. We calculate the ground state energy density for this lattice size with ED and our variational ansatz. The result is shown in Fig. 8. The two approaches exhibit good agreement in the strong coupling regime. For intermediate couplings differences becomes more pronounced leading to qualitatively different results in the weak-coupling limit g → 0 .
Since the electric Hamiltonian becomes bounded in the truncated theory, it does not contribute in the weak coupling limit. In the U (1) theory, however, the electric Hamiltonian is unbounded and the growth in electric energy leads to a finite result for the ground state energy in the continuum limit.

IV. REAL-TIME DYNAMICS
In this section, we study out-of-equilibrium dynamics by applying the following quench protocol: We prepare the ground state for the compact QED Hamiltonian at some coupling g 2 , quench to a Hamiltonian with a different coupling constant g 2 quench and observe the subsequent time evolution. The observables we track during the evolution are Wilson loops and the electric field (their expectation values in terms of the variational parameters can be found in Appendix C). In addition we check whether the energy is conserved throughout the whole time evolution.

A. Time-dependent variational principle
To study dynamical phenomena, we employ the timedependent variational principle. The equations of motion are projected onto the tangent plane of our variational manifold. For every variational parameter γ R/I k we define a corresponding tangent vector Ψ where P Ψ ensures orthogonality to |Ψ CP G : If we restrict the momenta k of the variational parameters to the set K defined in eq. (9), all tangent vectors become linearly independent. This allows to invert the Gram matrix G k k ≡ Ψ R k Ψ R k with k, k ∈ K. Since our variational manifold is Khler, we can express the time evolution of the variational parameters γ R/I k (k ∈ K) in the following way [61]:

B. Benchmark of variational ansatz
Since we are dealing with a variational ansatz, one should try to test it against exact results. For a comparison, we use the exact diagonalization results of the Z 3 theory. Since the truncation in the electric basis led to significant differences in the ground state energy already for intermediate coupling and time-dynamics increase the variance in the electric field, we can only expect reasonable agreement for a quench within the strong coupling region. We choose to quench the Hamiltonian from g 2 = 2.5 to g 2 = 4.0. The result is shown in Fig. 9. Even though truncation effects might still play a minor role in that quench, the comparison shows that the variational state can approximate amplitude and frequency of the oscillation.

C. Quench dynamics
We start with quenches in the weak-coupling regime where finite-size effects are most pronounced. We are interested in the maximal time up to which we can extract physics in the thermodynamic limit before boundary effects due to our finite lattice start to play a role. To compute that point in time, we perform the same quench on different lattice sizes and check where they start to deviate from each other. In order to easily compare observables for different lattice sizes, we restrict ourselves to the  10. Variational time evolution after a quench from g 2 = 0.8 to g 2 = 0.5 for lattice sizes of 8 × 8, 10 × 10 and 12 × 12. The inset shows the relative error in energy E with respect to the initial energy E0 after the quench. sector without static charges. We will focus on tracking the 1 × 1 Wilson loop during time evolution. We probed two different quenches, one from g 2 = 0.8 to g 2 = 0.5 for an 8 × 8, 10 × 10 and 12 × 12 lattice (shown in Fig. 10) and another one from g 2 = 0.6 to g 2 = 0.3 for lattice sizes of 16 × 16, 18 × 18 and 20 × 20 (shown in Fig. 11). The time evolution on the 8 × 8 lattice agrees with the 12×12 lattice up to t max,8 ∼ 3.8, the 10×10 lattice up to t max,10 ∼ 4.8. The energy is conserved for all lattice sizes up to a relative error of the order 10 −3 . During the time spans where we can reliably extract the time evolution, the Wilson loops indicate equilibrating behavior. This statement is supported by the second quench, where the smaller coupling constants allow us to reach larger lattices. The 16 × 16 and 18 × 18 lattice agree with the 20 × 20 lattice up to t max,16 ∼ 8.5 and t max,18 ∼ 9.5. The energy is conserved up to a relative error of 10 −6 . We can only make a statement about the equilibration of Wilson  [62,63].
In the next step, we look at a quench from weak to strong coupling (g 2 = 0.5 to g 2 = 4.0) for an 8 × 8 lattice without static charges. We track the time evolution of quadratic Wilson loops with edge sizes ranging from one to four. The result is shown in Fig. 12. Although all Wilson loops equilibrate at zero on short time scales (between t eq,4 ∼ 0.2 for the 4×4 Wilson loop and t eq,1 ∼ 0.5 for the 1 × 1 Wilson loop), we carried out the same evolution on a 7 × 7 lattice and found the same behavior.
The coupling constant at g 2 = 4.0 is large enough to approximate the spectrum by the strong-coupling limit g 2 → ∞, where the eigenstates |n become diagonal in the electric basis (this can be seen e.g. in the spectrum of the Z 3 theory which is available due to exact diagonalization). In this limit the thermal expectation value of Wilson loops vanishes trivially: For this special quench, we can thus verify that the Wilson loops equilibrate at their thermal expectation value. The next quench we will study is from strong to weak coupling. We quench on an 8 × 8 lattice from g 2 = 4.0 to g 2 = 0.5 with static charges horizontally separated by four links. Besides the 1 × 1 Wilson loop at the origin, we observe how the electric field of the ground state at g 2 = 4.0, a strongly confined fluxtube, evolves after the quench, in particular the electric field E 1 (x 1 = 2, x 2 = 4) (one of the links inside the fluxtube, see Fig. 13). It starts close to one, the strong-coupling value of the electric field, and decreases rapidly to E C 1 (2, 4) = 0.322, the value of the Coulomb electric field on that link (shown in the red dashed line). The Wilson loop seems to equilibrate on longer time scales.
The energy is conserved up to a relative error of 10 −2 . The larger error compared to previous quenches can be explained by the fact that around t ∼ 0.25 the approximation method of the infinite sums appearing in the evaluation of expectation values changes from the low γ k to the high γ k approximation (see section II B). In that transition region higher orders need to be calculated using uniform sampling (see Appendix B) which introduces additional errors. However, the relative error is still small and observables have no visible jump in this region, indicating that the two approximation schemes work. After the transition region the energy is well conserved due to the fact that the variational parameters γ I k increase, making the approximation of the infinite sums involved in the calculation of expectation values very easy (see section II B).
The spreading of the electric field from inside the flux tube between the two charges towards the Coulomb configuration of the electric field is illustrated in Fig. 14. An interesting question is whether the state becomes deconfined at long times. We cannot use the scaling of spatial Wilson, this only serves as an indicator for confinement in the ground state [56]. Since in our formulation the value of the longitudinal (Coulomb) part of the electric field is fixed and only the transversal part is dynamical (see Appendix A), we can measure precisely how much an electric field configuration differs from the Coulomb configuration. At t = 2.0, in the last of the three pictures in Fig. 14, the difference to the Coulomb configuration is of order 10 −12 for the whole lattice, with no remnant of an electric flux tube between the two charges. This is a strong indication that the state becomes deconfined, corresponding possibly to a thermal state with a temperature above the confinement-deconfinement transiton [64,65].

V. CONCLUSION
We introduce a new class of variational states, complex periodic Gaussian states, to study ground state properties and real-time dynamics in a (2+1)-dimensional U (1) lattice gauge theory. The evaluation of expectation values can only partially be done analytically, an infinite sum remains to be computed numerically. We present a scheme to approximate them for all variational parameters on an 8 × 8 lattice and for the weak-coupling regime up to 20 × 20. This allows us to study the variational ground state of these states over the whole coupling region and extract the thermodynamic limit. We benchmark our ansatz against the exact ground state for the one-plaquette case. We also compute the string tension  13. Variational time evolution on an 8 × 8 lattice after a quench from g 2 = 4.0 to g 2 = 0.5 lattice with a positive charge placed at (x1 = 2, x2 = 4) and a negative charge at (x1 = 6, x2 = 4). The 1 × 1 Wilson loop at the origin and the electric field on a link between the two charges (E1(2, 4)) is measured. The red dashed line represents the Coulomb value of the electric field. The inset shows the relative error in energy E with respect to the initial energy E0 after the quench.
using two different methods: First, by fitting the static potential between two charges with a 2d Coulomb potential and a linear potential. Secondly, we fit the exponential decay of Wilson loops with an area and a perimeter law. The two approaches are complementary since in the strong-coupling regime Wilson loops become difficult to fit due to the tiny value of large Wilson loops while the static potential approach works well as energy differences become larger. In the weak-coupling regime, however, the string tension becomes too small to extract the linear part of the potential on the given lattice sizes while Wilson loops decay only modestly allowing reliable fits. We are able to observe the expected exponential decay of the string tension in the weak-coupling regime.
Since our variational states do not need a truncation in the local Hilbert space, we compare our U (1) ground state data with exact diagonalization results for a Z 3 theory to study trunctation effects in the electric basis. The results agree for strong couplings and start to differ significantly for intermediate couplings. While the ground state energy of the truncated theory goes to zero in the continuum limit g 2 → 0 (since the electric energy is bounded), the variational ground state energy tends towards a finite value due to the variance of the electric field growing unboundedly.
In the last section, using the time-dependent variational principle, we probe out-of-equilibrium dynamics after a quench of the coupling constant. As a benchmark, we compare the variational time evolution after a quench within the strong-coupling regime with exact diagonalization results of the Z 3 theory. We then start by studying quenches within the weak-coupling regime where we ex-FIG. 14. Variational time evolution of the electric field on an 8 × 8 lattice after a quench from g 2 = 4.0 to g 2 = 0.5 with a positive charge placed at (x1 = 2, x2 = 4) (blue dot) and a negative charge at (x1 = 6, x2 = 4) (red dot). The color of the charges is only for graphical illustration (not related to the colorbar). The expectation value of the electric field is shown at t = 0.0, t = 0.2 and t = 2.0. At t = 0.0, the state is in the variational ground state for g 2 = 4.0 where the electric flux is confined between the two charges. After the quench, the electric field starts to spread over the lattice (t = 0.2) and equilibrates at the Coulomb value for this charge configuration (t = 2.0).
pect finite size effects to be significant. We compare the time evolution of a Wilson loop after the same quench for different lattice sizes to estimate at which time scales smaller lattices deviate from the thermodynamic limit. The times we can reach are large enough to indicate equilibration of Wilson loops.
In the next step, we perform a quench from weak (g 2 = 0.5) to strong coupling (g 2 = 4.0) and track the time evolution of differently sized Wilson loops. They all equilibrate at zero which is the thermal expectation value in the strong-coupling limit (g 2 → ∞). Since the spectrum at g 2 = 4.0 is close to the strong-coupling limit, this indicates that the Wilson loops equilibrate at their thermal expectation values. We also study a quench from strong to weak coupling in the sector of two static charges. We observe that the electric flux, which is perfectly confined for the strong-coupling ground state, spreads over the whole lattice and equilibrates at the Coulomb value for the electric field to very high accuracy, leaving no trace of confinement.
In all considered quenches, we see equilibrating behavior of observables up to the times where boundary effects start to play a role. It would be interesting to compare the equilibrated expectation values to thermal expectation values which can be computed by Monte-Carlo simulations [62,63]. Another interesting application for Monte-Carlo methods would be in the numerical evaluation of the variational ansatz by approximating the infinite sums. This could potentially allow to reach larger system sizes. The accuracy of the simulations would need to be high in order to carry out the evolution over reasonable time scales while ensuring energy conservation. Another natural extension of this work is the treatment of (3+1)-dimensional compact QED. By generalizing an idea in [52] to complex Gaussian states, a variational ansatz can be designed for 3+1 dimensions. However, due to additional local constraints appearing in 3+1 dimensions (compared to one global constraint in 2+1 dimensions) a new numerical approximation scheme would be required. Another interesting idea is to include dynamical matter. In order to couple the gauge degrees of freedom to matter, it is essential to find a transformation of such a theory to variables which are well-suited for complex periodic Gaussian states. Work in this direction is in progress. Extending the ansatz to non-Abelian gauge theories is more difficult since they do not allow a translationally invariant formulation in terms of gauge-invariant plaquette variables. However, other gauge-invariant variables could be used to construct similar ansatz states [66].
In this section, we want to give a short review on the separation of gauge fields into (almost) gauge invariant plaquette variables and a static part corresponding to the longitudinal Coulomb field. For simplicity and since this is the charge configuration used throughout the paper, we will focus on a situation with two static charges placed vertically at x 2 = d 2 separated horizontally by a distance d. Other charge configurations follow analogously. We want to split the electric flux line between the two charges into a transversal component, generated by the lattice curl of a field on the plaquettes and into a longitudinal component, generated by the lattice gradient of a scalar field φ on the vertices. All other electric flux configurations can be created on top of it by exciting an electric flux loop around a plaquette or around the whole lattice.
The longitudinal part is by definition of the form where ∇ (+) is the lattice forward derivative. Using Gauss law, with ∇ (−) the lattice backward derivative, we arrive at a lattice version of Poisson's equation: The solution for φ is x 2 ) and x 1 , x 2 ranging from 0 to L−1. The same applies to y and k. There is no k = 0 contribution since the total charge on a periodic lattice needs to be zero because of gauge invariance. E L i (x) then follows straightforwardly from (A1).
We write the transversal part as the curl of an -field on the plaquettes, where the plaquette corresponding to x is the one having x as its the bottom left corner. We take the curl of the above expression and use a lattice analog of the vector identity ∇ × ∇ × A = ∇(∇ · A) − ∆A, here in two dimensions, to obtain a Poisson equation for the -field: where we sum over repeated indices. This equation is solved by In our case, ij ∇ (+) i E j (x) is 1 on the plaquettes above the electric string connecting the charges and −1 below the string. However, since the sum over this expression will always be zero, we cannot generate a constant electric field with the -field which is required sinceẼ 1 (k = 0) = d L . It is important to also consider the Polyakov loop winding horizontally around the lattice. We choose it to wind around the lattice at x 2 = d 2 and the electric field along it to be poly,1 = d L . We define an additional -field const on the plaquettes, on top of : It is defined in such a way that giving us the k = 0 component of the electric field. We can now rewrite the electric field operator as: as the Coulomb electric field, we obtain the electric Hamiltonian given in eq. (5): . Besides orthogonality between longitudinal and transversal component of the electric field, we also used Plancherel's theorem which ensures orthogonality between the constant part and the other two since their k = 0 component is zero. In this section, we review the numerical evaluation of complex periodic Gaussian states in more detail. We saw in section II B that the region with γ k = γ R k + γ I k 2 γ R k −1 ≈ 1 is the most difficult to evaluate. Since the variational ground state (for which γ I k = 0) varies from high γ R k for low couplings to low γ R k for large couplings, there is a transition region at g 2 ∼ 1.1 where γ k approaches one. We therefore want to study the approximations to all infinite sums involved in the computation of the variational energy in eq. (20) on an 8 × 8 lattice without static charges for the ground state at g 2 = 1.1 which is the highest coupling where the high γ k approximation is used and g 2 = 1.2 which is the lowest coupling for which the low γ k approximation is used. For all other couplings the contributions to infinite sums decay faster with higher orders compared to one of the two examples discussed below. The variational parameters γ R k for these states (rounded to three digits) are shown in Table I and II. The values are displayed not only for the independent γ R k (k ∈ K), but are split between the dependent parameters γ R k and γ R −k , to illustrate that they lie in the transition region γ k ≈ 1 between the two approximation methods.
Using γ I k = 0 and p = 0, the expressions we need to compute for the variational ground state at g 2 = 1.1 simplify significantly: for the electric energy, The orders are organized in groups. P denotes orders which contain a growing number of 1's. Since Np and −Np are evaluated together, P also represents orders with a growing number of −1's. M (1)P contains orders whose non-zero elements are a single −1 and a growing number of 1's. The first order in M (1)P contains a pair of (1,-1) as non-zero elements. M (2)P is structured in the same way as M (1)P but with two −1's. Analogously for the other groups. The Np = 0 configuration (denoted as 0) has vanishing contribution to J el but a non-zero contribution to Jmag.
for the magnetic energy and the normalization We include orders with N p configurations of up to 8 pairs of {1, −1} and the rest zeros. The first three are computed exactly and the remaining five by uniform sampling. Additionally, we compute exactly the orders {N } 2,−1,−1 and {N } −2,1,1 to show they have negligible contributions. Orders like these, whose N p configura-tions differ only by a minus sign can be evaluated together by evaluating for every permutation not only the contribution of N p but also of −N p . Therefore, from now on orders which are not closed under reflection will also include all their permutations multiplied by minus one. This will be heavily used in the low γ k approximation.
The exact evaluation of orders is based on an algorithm which generates all permutations of a multiset in O(1) time [67], i.e. the time to generate a new permutation is independent of the permutation size. It is much smaller than the time needed to do computations with a permutation which allows to highly parallelize the process and reach higher orders. The evaluation of an observable with respect to a set of permutations {N } with uniform sampling is based on the approximation: where S is a set of s randomly drawn N p configurations from {N } and p the number of permutations within this order. For all orders which are computed with uniform sampling we use s = 10 8 in the high γ k approximation and s = 10 7 in the low γ k approximation. The contributions to I el and I mag for the high γ k approximation are displayed in Fig. 15. We do not show this analysis for the normalization since its contributions decay faster than the ones for I el and I mag . The errors due to uniform sampling are too small to be shown in the plot, the biggest error occurs in the order with four pairs of {1, −1} which has a contribution of 347.54 (15) to I el and of 622.70 (24) to I mag .
For the variational ground state at g 2 = 1.2, the infinite sums in eq. (20) reduce to (B5) for the computation of the electric energy, with 1 2 p k = 1 2L e −ipk for the computation of the magnetic energy and for the normalization. Since we do not have a global constraint in the low γ k approximation, more orders contribute to the infinite sums. The contributions to J el and J mag of different orders are given in Fig. 16. The errors are again too small to be displayed, the biggest one occurs in the order {N } −1,1,1,1,1,1 with contributions of 15.22 (1) to J el and of 44.07(4) to J mag .
Both approximation schemes decay reasonably well with higher orders and the truncation of even higher orders can be justified. Moreover, the errors introduced due to uniform sampling are small, in particular since the lowest orders were still calculated exactly. The algorithm we applied during computations to decide with which approximation method an expectation value should be evaluated was to select higher orders and compute them by uniform sampling with a low sample size of s = 10 5 . This allowed us to choose the scheme which had a better decay with higher orders.

Appendix C: Observables
In this section, we provide formulas for important quantities which are too lengthy to fit into the main body of the manuscript. This includes formulas for expectation values of observables, namely Wilson loops and electric field, and a formula for the gradient of the energy with respect to the variational parameters which is essential to minimize the variational energy and carry out the timedependent varational principle. For the latter, we present additionally a formula for the Gram matrix. For every infinite sum appearing in the expressions, we provide both the high and low γ k approximation.
We start with the expectation value of a Wilson loop along a contour C, where p ∈ C denotes all plaquettes within this contour: 2L p∈C e −ipk and γ k = γ R k + γ I k 2 γ R k −1 . Another observable which is used in the manuscript is the electric field. We present for simplicity the expectation value of a horizontal link emanating from vertex x, E x,1 . The plaquette above the link is denoted asp. The expectation value for a vertical link follows analogously.
The next quantity we present is the gradient of the variational energy with respect to the independent parameters γ R k and γ I k (k ∈ K). We split the energy into an electric and a magnetic part to make the expressions less cumbersome. We start with the derivatives of the electric energy with respect to γ R k and γ I k (k ∈ K): We denote by m k the number of elements in the equivalence class k ∈ K which is two if k = −k and one if k = −k.
The expression for |N k | 2 for both high and low γ k approximation can be found in eq. (21). The gradient of the magnetic energy with respect to γ R k and γ I k takes the form: