Real-time lattice gauge theory actions: unitarity, convergence, and path integral contour deformations

The Wilson action for Euclidean lattice gauge theory defines a positive-definite transfer matrix that corresponds to a unitary lattice gauge theory time-evolution operator if analytically continued to real time. Hoshina, Fujii, and Kikukawa (HFK) recently pointed out that applying the Wilson action discretization to continuum real-time gauge theory does not lead to this, or any other, unitary theory and proposed an alternate real-time lattice gauge theory action that does result in a unitary real-time transfer matrix. The character expansion defining the HFK action is divergent, and in this work we apply a path integral contour deformation to obtain a convergent representation for U(1) HFK path integrals suitable for numerical Monte Carlo calculations. We also introduce a class of real-time lattice gauge theory actions based on analytic continuation of the Euclidean heat-kernel action. Similar divergent sums are involved in defining these actions, but for one action in this class this divergence takes a particularly simple form, allowing construction of a path integral contour deformation that provides absolutely convergent representations for U(1) and SU(N) real-time lattice gauge theory path integrals. We perform proof-of-principle Monte Carlo calculations of real-time U(1) and SU(3) lattice gauge theory and verify that exact results for unitary time evolution of static quark-antiquark pairs in (1 + 1)D are reproduced.


I. INTRODUCTION
Lattice quantum field theory methods have long been used to calculate equilibrium properties of strongly coupled quantum systems. Discretization of quantum field theories on a spacetime lattice provides a UV regulator that allows expectation values of observables to be computed after renormalization and extrapolation to the continuum limit. Lattice gauge theory (LGT) in particular allows observables to be computed in strongly coupled gauge theories, for example including observables in quantum chromodynamics (QCD) that are relevant for understanding the dynamics of the strong force in the Standard Model.
LGT calculations typically rely on Monte Carlo methods to evaluate Euclidean path integrals by stochastically sampling gauge field configurations from a probability distribution proportional to e −S E , where the Euclidean action S E is assumed to be real. This approach is suitable for determinations of many equilibrium gauge theory properties, which can be related to expectation values in Euclidean spacetime.
Non-equilibrium properties of gauge theories are more difficult to obtain from Euclidean correlation functions; for example, transport coefficients of the quark-gluon plasma and neutron stars, properties of early universe phase transitions, and inclusive hadron scattering crosssections involving resonance production have been inaccessible to ab initio Euclidean LGT approaches despite significant theoretical and phenomenological interest. In principle, out-of-equilibrium gauge theory properties defined by correlation functions of timelike separated operators can instead be accessed using the Schwinger-Keldysh formalism [1,2], which represents out-of-equilibrium observables with path integrals involving regions of both Minkowski and Euclidean signature. In practice, Monte Carlo methods cannot directly be used for efficient nonperturbative calculations of these lattice-regularized realtime or Schwinger-Keldysh LGT path integrals, since the path integral weights involve a pure phase factor e iS M , where S M is the action for the Minkowski region, and the path integral weights thus cannot be interpreted as a probability distribution to be used for importance sampling.
In some cases, path integral contour deformations have been used to tame this sign problem associated with realtime path integrals; for example, real-time calculations have been performed using path integral deformations in (0 + 1)D quantum-mechanical theories [3][4][5][6][7]. Suitably chosen path integral contour deformations in complexified field space can exactly preserve the total path integral while reducing phase fluctuations of the integrand and introducing fluctuations in the magnitude of e iS M that are amenable to importance sampling. Path integral contour deformations have also been applied to sign problems for finite-density systems and other observables in a variety of quantum field theories; see Ref. [8] for a recent review. This promising approach to mitigating sign problems has recently been extended to generic observables with signal-to-noise problems in SU (N ) gauge theories [9], but the construction of contour deformations for real-time LGT path integrals has not yet been addressed.
It was pointed out by Hoshina, Fujii, and Kikukawa (HFK) in Ref. [10] that another obstacle facing calculations of real-time LGT path integrals is the determination of a suitable discretized action S M . In particular, Ref. [10] argues that applying the same discretization as the Wilson gauge action [11] to the real-time continuum theory results in a real-time LGT that is non-unitary.
Although it is formally possible to construct a unitary time evolution operator for LGT through analytic continuation of the eigenvalues of the imaginary-time transfer matrix [12], this formal construction cannot be practically implemented without first solving for the LGT spectrum. Instead, HFK propose an alternative real-time LGT transfer matrix obtained by analytically continuing the character expansion of the kinetic term in the Wilson action (i.e. all terms involving timelike plaquettes). The resulting real-time transfer matrix is unitary and formally allows the definition of a discretized real-time action, called the HFK action below, and of associated real-time LGT observables.
The HFK action, however, is defined by an infinite series and we demonstrate below that it does not converge for some or all values of the gauge field in U (1) and SU (N ) gauge theory, making it impossible to calculate the weights necessary for an importance sampling approach. This divergent representation is a generic feature of real-time LGT actions that are defined by local kinetic and potential energy terms and that give rise to unitary transfer matrices, as discussed below and detailed in Appendix C. In order to remedy this situation, this work introduces path integral contour deformations that explicitly depend on summation indices appearing in the definition of the action. These path integral contour deformations can change the convergence properties of the action for fixed gauge field values, though singularities still arise for gauge field configurations corresponding to endpoints of the path integral contour that must be kept fixed during deformation. Rotating the prefactor of the kinetic term in the action starting from −1 and approaching i (analogous to a Wick rotation of the continuum theory) regularizes these singularities and makes the path integral absolutely convergent everywhere outside of the limit. By deforming the path integral contour as a function of the prefactor and exactly cancelling contour segments that are related by periodicity, the remaining path integral is made absolutely convergent even after analytically taking the limit to i. A simple contour deformation that renders the U (1) HFK path integral convergent in this sense is introduced below. Unfortunately, it is challenging to construct an analogous contour deformation that would lead to an absolutely convergent representation of the SU (N ) HFK path integral.
This challenge motivates the introduction of another class of actions based on analytic continuation of the heat-kernel action of Menotti and Onofri [13]. The key feature of these actions is a kinetic term that is exactly unitary at all values of the lattice spacing, much like the HFK action; in comparison to the HFK action, however, a simpler summation is involved in their definition. The heat-kernel kinetic term can be combined with any real potential without violating unitarity. Possible choices of the potential include either a term with the same heat kernel form applied to the spacelike plaquettes or the Wilson potential term (i.e. the terms in the Wilson action involving spacelike plaquettes). In the following,  I. Unitarity in the continuum limit, convergence of the definition of the action, and existence of convergent path integral representations using contour deformations for the four actions considered in this work -the Wick-rotated realtime Wilson action, the Hoshina-Fujii-Kikukawa (HFK) realtime action [10], the real-time heat-kernel action (HK) derived from the Euclidean action given in Ref. [13], and the real-time modified heat-kernel (HK) action introduced in the main text. A truncated heat-kernel action which is convergent and non-unitary is also studied in Appendix D, including details on the continuum limit of this action. Based on the contour deformations introduced in this work, the HK action is the only unitary choice for which a convergent Monte Carlo calculation can be performed in real-time SU (N ) LGT, but it is possible that more sophisticated contour deformations could be discovered to provide convergent deformations for the cases marked with .
these two options are respectively termed the real-time heat-kernel (HK) action and the modified real-time heatkernel (HK) action. The summations defining the HK and HK actions are also divergent, but a simple path integral contour is shown to provide an absolutely convergent representation of the HK action for U (1) and SU (N ) gauge theory in Minkowski spacetime with any dimension. This is achieved using a prescription similar to the one applied for the HFK action: parameterizing a "Wick rotation" of the kinetic term prefactor, exactly cancelling pieces of the deformed contour related by symmetry, and then analytically taking the limit to prefactor i. Thus real-time path integrals using the HK action can be evaluated after contour deformation using Monte Carlo techniques, and real-time LGT observables can be numerically computed. The comparison of the HK and HK actions to the HFK action and the real-time Wilson action is summarized in Table I and is analyzed in detail below.
Although the real-time Wilson action is not suitable for quantum LGT, the derived equations of motion do give well-defined deterministic evolution of the fields involved and result in a well-formed classical theory. In these approaches, the lattice spacing is a free parameter that is independent of the gauge field coupling, rather than a dynamical quantity whose value in physical units must be determined by tuning the gauge coupling, as described for example in Ref. [37]. The non-existence of a unitary continuum limit for the real-time Wilson action in quantum LGT is thus irrelevant for calculations of solutions to the classical equations of motion associated with the real-time Wilson action.
In the quantum setting, the (non-unitary) real-time Wilson actions for U (1) and SU (2) LGT have previously been used in complex Langevin calculations [56,57]. Complex Langevin methods are not guaranteed to reproduce exact results for real-time LGT, and Ref. [56] finds that complex Langevin results only reproduce analytically calculable results for simple real-time LGT observables for values of the Langevin evolution time that are not too large. Methods based on reweighting and gauge fixing are found to increase the region in which complex Langevin reproduces exact results in Ref. [57]. The exact results for the one-plaquette model calculated using the real-time Wilson LGT action in this reference agree with the exact results for unit area Wilson loops in (1 + 1)D Minkowski spacetime presented in Sec. III E below. The non-unitarity of time evolution in the oneplaquette model with the real-time Wilson action implied by these results is not mentioned in Ref. [57]. The lack of a well-defined continuum limit for the SU (2) oneplaquette model with the real-time Wilson action is discussed in the reference; however, Ref. [57] assumes that the corresponding continuum limit of (3 + 1)D LGT with the real-time Wilson action exists and can be used to calculate physical observables in real-time gauge theory. It is demonstrated below that the real-time Wilson action is not unitary in arbitrary spacetime dimensions, even in the small gauge coupling limit, and that the continuum limit of real-time LGT with the Wilson action in (3 + 1)D either does not exist or is non-unitary, depending on the choice of gauge group SU (N ). In either case, the real-time Wilson action does not provide a suitable starting point for calculating physical observables in realtime gauge theory.
Finally, it is possible to simulate real-time gauge theory dynamics using the Hamiltonian formalism rather than the path integral, as is often considered in the context of tensor network and quantum computing approaches; see Ref. [58] for a recent review. Though unitarity of time-evolution in the Hamiltonian formalism is also a key condition, the focus of the present work is on the construction of actions suitable for classical simulation of discrete real-time path integrals, and the discussed actions and contour deformations are not immediately relevant to Hamiltonian quantum simulation. However, the analysis of unitarity explored here may have relevance for similar analyses for quantum computing approaches. For example, Ref. [59] shows how real-time path integrals in lattice field theory can be described using quantum circuits with Trotterization errors described as discretization effects, but asserts that a time-evolution operator for quantum simulation of SU (N ) gauge theory obtained using the real-time Wilson action is unitary. 1 It would be interesting to consider analogous quantum circuit descriptions of real-time transfer matrices for the unitary real-time LGT actions considered in this work.
The remainder of this work is structured as follows. The (non)-unitarity of real-time transfer matrices defining discretized path integrals of compact and non-compact variables, including the non-unitarity of the real-time Wilson action, is given in Sec. II. The unitary real-time LGT actions shown in Table I are introduced in detail and discussed in Sec. III. In Sec. IV, path integral contour deformation techniques are used to construct absolutely convergent representations of path integrals involving the HFK action for U (1) LGT, the HK action for U (1) and SU (N ) LGT, and the associated Schwinger-Keldysh action; proof-of-principle Monte Carlo calculations in (1 + 1)D are also discussed. Our conclusions and outlook are summarized in Sec. V.

II. REAL-TIME EVOLUTION OF COMPACT AND NON-COMPACT VARIABLES
In a quantum field theory (QFT) defined on a Minkowski spacetime background, the continuum action S M provides a useful starting point for understanding and perturbatively calculating expectation values of quantum operators and corresponding physical observables. The continuum action is manifestly Lorentz and translation invariant, ensuring that the physics encoded in the action satisfies the Poincaré symmetry observed in nature. A Hamiltonian and related quantum operators can be defined using the construction of a Hilbert space on a co-dimension-one submanifold of spacetime, which superficially breaks Poincaré invariance [60]. The path integral formalism [61], however, provides a way to relate expectation values of quantum operators to integrals over fluctuations of classical fields weighted by the manifestly Poincaré-invariant function e iS M . Path integrals in imaginary time can be used to calculate expectation values for systems in thermal equilibrium when suitable temporal boundary conditions are used for bosons and fermions and the inverse temperature β is set by the length of the imaginary-time direction.
1 The real-time transfer matrix elements computed in Ref. [59] can be represented as integrals with pure-phase integrands, as shown in Eq. (H10) of that work. The real-time transfer matrix can also be written as an integral over unitary operators, as given in their Eq. (18). However, these features do not imply that the eigenvalues of the real-time transfer matrix are pure phases and therefore do not establish unitarity.
Formally, the imaginary-time path integral can be related to a real-time path integral by Wick rotation x 0 ↔ ix 0 of the time coordinate [62]. The action S M to be used in real (Minkowski) time consists of the kinetic minus potential energy, whereas the action S E in imaginary (Euclidean) time consists of the kinetic plus potential energy, with the corresponding path integral weights related by Wick rotation as e iS M ↔ e −S E . Although these path integral relations are straightforward in continuum QFT, subtleties arise in the connections between real-and imaginary-time path integrals with lattice-regularized actions involving compact variables as discussed below.

A. The SHO and quantum rotator
The continuum path integral is ill-defined without a prescription for regularizing UV divergences and extrapolating to remove the regulator. Discretizing spacetime provides such a regulator with the desirable properties of being non-perturbative, gauge-invariant, and amenable to numerical simulation. This discretization is well-understood for gauge fields in imaginary-time lattice field theory, but as discussed in Ref. [10] and below subtleties arise when considering real-time path integrals or path integrals involving real-time components, such as the path integrals required in the Schwinger-Keldysh approach to computing out-of-equilibrium observables [1,2]. These difficulties can be simply demonstrated even in (0 + 1)D quantum mechanical systems. Though these systems do not possess a Lorentz symmetry, the challenges in these (0 + 1)D path integrals are investigated as simple analogues to the challenges that arise in (3 + 1)D LGT.
We first consider the simple harmonic oscillator (SHO), the (0+1)D theory of a single noncompact variable x ∈ R constrained by a quadratic potential. The continuum actions for the SHO in real and imaginary time can be written respectively as where x(t) and x(τ ) are the position histories in real and imaginary time respectively, and we work in units for which the mass is set to 1. Transition amplitudes between states |x and |x after time evolution by L T are given in real time by x |e −iĤL T |x and in imaginary time by x |e −ĤL T |x . These transition amplitudes can be used to extract the energy spectrum and other physical properties. A discretized path integral can be used to write the transition amplitudes using time steps of size a as gives the prefactor associated with time evolution on the real-and imaginary-time contours, respectively. Factoring out 1/s in the exponent allows one to collect the exponentials into a temporally-discretized version of the continuum action in Eq. (1), where Dx ≡ dx a . . . dx L T −a and the discretized action is given by with the kinetic and potential energy terms K and V given by Noting that s 2 = −1 for real time and s 2 = +1 for imaginary time, we arrive at the usual conclusion that the discretized action for real or imaginary time can be related by simply replacing the sign in front of the potential term and using the appropriate prefactor of 1/s in the path integral weights in Eq. (3). In either real or imaginary time, the discretized path integral Dx e 1 s S(x;s) gives an approximation to x |e sĤL T |x which is accurate to O(a 2 ) under the assumption that the integral exists as a → 0 and can be expanded in powers of a. The existence of such a limit is critical to extrapolating physical quantities of interest to the continuum. For the SHO, the a → 0 limit is welldefined and can be studied using the transfer matrix. The matrix elements of the real-and imaginary-time transfer matrices respectively are given by one factor of the path integrand, where t = na and τ = na in real and imaginary time. The corresponding SHO Hamiltonian is given byĤ =K +V , where x |V |x = ω 2 x 2 2 δ(x − x) andK =p 2 /2 is defined in term of the momentum operatorp satisfying [x,p] = i.
Noting that x |e sap 2 /2 |x ∝ e (x −x) 2 /(2as) , the real-and imaginary-time transfer matrices can be expressed aŝ By the Lie-Trotter product formula [63], products of the real-and imaginary-time transfer matrices approximate products of the corresponding time evolution operators in real and imaginary time, e sĤa , with errors that vanish as a → 0.
The continuum real-and imaginary-time transfer matrices e −iaĤ and e −aĤ are unitary and positive definite respectively for HermitianĤ. These properties allow a spectral decomposition in either case and are desirable to maintain in the discretized theory. Unitarity of the SHO real-time transfer matrix (up to a constant normalization factor suppressed in Eq. (7)) can be demonstrated by direct computation, and can be shown to be positive definite by the fact that the Gaussian integral kernel itself is positive definite [64], (9) Fundamentally, these desired properties of the transfer matrix emerge for the SHO because the kinetic portion of the discretized action is respectively a unitary integral kernel or a positive-definite integral kernel when the prefactor is i or −1. Specifically, in Eq. (8), it is exactly unitarity of the kinetic integral kernel that produces the term δ(x − x ). The potential factors multiplied by this delta function are inverses which cancel, regardless of the specifics of the potential. Similarly, in Eq. (9) only the positive-definiteness of the Gaussian kinetic term is required to show overall positive definiteness, given that the potential appeared in the same way on each side of the transfer matrix (as long as the potential is real and bounded from below so that f (x)e −aV (x) ∈ L 2 (R)). For other theories of noncompact variables, the kinetic term in the action can also generically be discretized as a Gaussian integral kernel that satisfies these properties, and these properties therefore extend to lattice field theories of noncompact variables.
We next consider the planar quantum rotator in (0 + 1)D in both real and imaginary time. This system can be physically interpreted as the SHO on a compact domain, which can be chosen for example to be the circular domain x ∈ [0, 2π], with x = 0 identified with x = 2π. This choice normalizes the length of the compact domain such that x can be considered as the angular variable of the quantum rotator. Since the microscopic description of the theory is identical to the SHO, the continuum action given for the SHO in Eq. (1) also describes the physics of the quantum rotator. We are also free to set ω = 0 in the potential because the compact domain ensures convergence of path integrals in the free theory, which simplifies the analytical manipulations below.
Although the continuum actions for the ω = 0 SHO and the quantum rotator are the same, the definitions of the discretized action and path integral for the quantum rotator require care. In particular, the derivatives used in the kinetic operator should be compatible with the identification of x = 0 and x = 2π. A typical choice is to write the path integral using a cosine for the discrete kinetic term, where as above s ∈ {−i, −1} gives the appropriate prefactor for real or imaginary time, respectively. The corresponding discretized action for the quantum rotator is The Taylor expansion of 1−cos( (11) is equivalent to the free theory SHO action (ω = 0) for small fluctuations of the position in lattice units. The two actions should therefore be perturbatively equivalent in the continuum limit.
The differences in behavior at non-zero lattice spacing become apparent when considering the transfer matrix description. The transfer matrices in real and imaginary time associated with the path integral in Eq. (10) are defined respectively bŷ The lack of any potential terms reflects our choice of working with the free theory. As argued above, the unitarity and positive definiteness of the transfer matrix in real and imaginary time depends only on the behavior of this kinetic integral kernel, as any (real) potential terms will cancel from the relations in Eq. (8) and (9) if the kinetic term satisfies the desired properties. For this choice of discretization, the imaginary-time transfer matrix does satisfy positivity, and can be expanded in terms of Fourier eigenfunctions using the Jacobi-Anger expansion [65] where I k is the modified Bessel function of the first kind with rank k. The eigenvalues e −1/a I k (1/a) all vanish as a → 0, but physical observables are determined by ratios I k (1/a)/I 0 (1/a) which converge to 1 in the continuum limit, allowing renormalization and extraction of quantities of interest in the continuum. On the other hand, the non-Gaussian nature of the kinetic integral kernel results in non-unitarity of the real-time transfer matrix, which can be seen by direct calculation, The breakdown of unitarity at non-zero a is an undesirable feature, but could be considered acceptable if the ratios of transfer matrix eigenvalues converged to unitnorm values in the continuum limit. Instead, many ratios of eigenvalues simply do not have a continuum limit, as can be seen by performing a similar Jacobi-Anger expansion in real time, As a → 0, the ratio I k (−i/a)/I l (−i/a) → 1 for k ≡ l mod 2 but the limit does not exist for k ≡ l + 1 mod 2, demonstrating that observables that depend on these ratios of eigenvalues do not have either a well-defined continuum limit or a spectral representation consistent with unitary time evolution. If a potential is included, these free theory eigenfunctions can generically be expected to mix, and the non-existence of a unitary continuum limit for certain eigenstates of e −iaK can be expected to spoil the existence of continuum limits for generic observables. This simple exploration of the SHO and quantum rotator highlights a concern that must be addressed if attempting to work with discretized real-time path integrals in a position-space representation. For compact variables, a position-space kinetic term in the action that satisfies periodicity may not be compatible with the replacement K+V → K−V in moving from the imaginarytime action to the real-time action, as this simple replacement may result in non-unitarity of the transfer matrix and prevent extrapolating to continuum physics.

B. Lattice gauge theory in real and imaginary time
The Standard Model of particle physics involves the gauge groups U (1) and SU (N ), and thus we focus on these two groups in explorations of suitable real-time actions for lattice gauge theory. As an Abelian group, the continuum limit of U (1) gauge theory can be accessed by extrapolating to zero lattice spacing using either the compact gauge group U (1) or the non-compact gauge group R; see for example Ref. [66]. Using the compact gauge group is susceptible to the subtleties discussed above and is the focus of our U (1) studies. In order to describe continuum Euclidean SU (N ) gauge theory, standard formulations of LGT [11,67] use SU (N ) variables in the fundamental group representation, which consists of N × N unitary matrices with unit determinant. These formulations are therefore also susceptible to challenges associated with path integrals involving compact variables when moving to real time.
A lattice gauge theory for gauge group G in D spacetime dimensions is defined in terms of a set of gauge fields U x,µ ∈ G, where x = (x 0 , x 1 , . . . x D−1 ) are the spacetime lattice points, and µ = 0, . . . , D − 1 labels the lattice axes with the (real or imaginary) time direction specified by µ = 0. We assume a lattice with an extent of L/a sites in each spatial direction and L T /a sites in the temporal direction, where a is the lattice spacing in physical units. When applying the Schwinger-Keldysh formalism for out-of-equilibrium observables, the lattice extent in the temporal direction is divided into regions of Euclidean, forward Minkowski, and reverse Minkowski time evolution; discussion on extending real-time actions to this setting is deferred to Sec. IV D and in other sections the spacetime signature is assumed to be uniform throughout the lattice. Each component U x,µ of the gauge field is associated with an edge connecting neighboring sites x and x +μ. For all actions under consideration an exact gauge symmetry holds: transforming all gauge field components by U x,µ → Ω x U x,µ Ω † x+aμ , for any field Ω x ∈ G, does not modify the value of the action.
Equilibrium properties of gauge theories can be determined from the continuum limits of expectation values of "observables" in LGT, defined as generic functions of the gauge field, O(U ), with Euclidean path integral representations where DU ≡ x,µ dU x,µ is the product of the Haar measure dU x,µ for each gauge field degree of freedom, S E is the Euclidean action, and Z E = DU e −S E (U ) is the partition function. We restrict to considering Euclidean actions that can be expressed as a sum of potential and kinetic energy functions involving the lattice gauge fields on individual timeslices and pairs of adjacent timeslices respectively, where U τ = {U x,µ | x 0 = τ } and N T = L T /a is the number of lattice sites in the (imaginary) time direction.
Restricting to actions of this form allows the construction of transfer matrices and explicit analyses of unitarity. Constructing unitary actions that violate this form, e.g. due to Symanzik improvement [68] or inclusion of matter fields, is left as the subject of future work. The Hilbert space for pure gauge theory on a fixed timeslice can be represented as a product of Hilbert spaces for group-valued quantum rotators [67]. Gauge field operatorsÛ x,k can be defined for k = 1, . . . , D − 1, i.e. they are associated with spacelike links, and are analogous to position operatorsx na for the quantum rotator discussed in Sec. II A. Assuming temporal gauge, the Hilbert space for pure gauge theory is spanned by polynomial functions of these operators. States |U x,k are defined by the eigenvalue relationÛ x,k |U x,k = U x,k |U x,k and are normalized to satisfy U x,k |U x,k = δ(U x,k , U x,k ). This Hilbert space can equivalently be described using a basis of L 2 -normalizable complex-valued functions f (U ) as detailed in Ref. [12]. Hilbert space states |f can be associated with this function basis by (⊗ x,k U x,k |) |f = f (U ), and the actions of Hilbert space operators in the function basis can be represented using integral kernels.
The imaginary-time transfer matrixT E can be concretely defined by an integral kernelT E (U τ +a , U τ ) = U τ +a |T E |U τ , where |U τ and |U τ +a are arbitrary tensor product basis states given by ⊗ x,k |U x,k for particular gauge field configurations U (τ, x),k and U (τ +a, x),k . The integral form of the action of this operator on a functionbasis state |f is then The integral kernelT E (U τ +a , U τ ) is analogous to the coordinate space matrix elementsT E (x τ +a , x τ ) for the quantum rotator. For a LGT action of the form in Eq. (17), the imaginary-time transfer matrix is defined in a general gauge by the integral kernel [69] where U ∂(τ +a,τ ) = {U x,0 | x 0 = τ }. The imaginary-time transfer matrix describes discretized imaginary-time evolution in LGT corresponding to the generic action in Eq. (17). For example, Euclidean correlation functions involving a pair of temporally separated operators A(U τ ) and B(U 0 ) have a transfer-matrix representation It is possible to formally (although not practically) construct a Hermitian HamiltonianĤ and unitary real-time evolution operator e −iĤt directly from the imaginary-time transfer matrix. Assuming that the transfer matrix for a particular choice of Euclidean action is positive definite, then it is possible to construct the Hamiltonian operatorĤ defined bŷ without encountering singularities of the logarithm [12]. Formally, a perfect real-time transfer matrix can then be constructed,Û By constructionÛ is unitary, with eigenvalues e −iEna related to the eigenvalues e −Ena ofT E . The existence of a positive-definite transfer matrix therefore guarantees the existence of a unitary time-evolution operator with the same energy spectrum. The positivity of the transfer matrix associated with the Wilson action [11] for Euclidean LGT was established early on in the study of lattice QFT through proofs of reflection positivity [70,71], and the transfer matrix was then explicitly constructed [72] and explicitly demonstrated to be positive [12]. These results and their generalizations to other actions crucially allow the energy spectra of Euclidean gauge theories obtained from the continuum limits of LGT results to be identified with the energy spectra of the corresponding Minkowski continuum gauge theories relevant for experiments involving real-time dynamics. This approach to determining the energy spectra of continuum gauge theories corresponds to first taking a → 0 and then subsequently taking τ → it in Fig. 1 and is expected to be valid for any Euclidean LGT with a positive-definiteT E . However, the calculation of correlation functions with timelike separated operators in Minkowski spacetime, relevant for example for inclusive scattering cross-sections and transport coefficients, is an ill-posed and practically challenging problem when using analytic continuation of numerical results for Euclidean correlation functions. For these and other applications, it may be advantageous to consider the opposite order of limits shown in Fig. 1, in which a real-time transfer matrix is constructed for Minkowski LGT and physical results are obtained by subsequently taking the continuum limits of real-time observables obtained in Minkowski LGT. The operatorÛ defined in Eq. (22) is not a suitable starting point for practically computing Minkowski LGT observables because it cannot be constructed without explicitly diagonalizingT E and working in the energy eigenbasis. A seemingly promising alternative approach is to replace the sum of kinetic and potential terms K + V with a difference K −V to move from the LGT Euclidean action to the LGT Minkowski action, with the goal of recovering the physics encoded inÛ in the continuum limit. The Minkowski action obtained by such a replacement is Minkowski expectation values described by this action are defined by where Z M = DU e iS M (U ) . A real-time transfer matrix T M can be defined for this action in analogy to Eq. (19) This real-time transfer matrix can be used to equivalently write matrix elements of products of operators separated in Minkowski time, otherwise given by real-time LGT path integrals involving products of operators. For example, Minkowski correlation functions involving a pair of temporally separated operators A(U t ) and B(U 0 ) are given analogously to Eq. (20) by The energy spectrum for Minkowski LGT with action S M can therefore be obtained from the eigenvalues ofT M . The real-and imaginary-time transfer matrices can be decomposed into products of potential-energy evolution operators that are diagonal in the coordinate basis and kinetic-energy evolution operatorsR E andR M , The real-and imaginary-time kinetic-energy evolution operatorsR M andR E are defined bŷ In temporal gauge the integrals over U ∂(τ +a,τ ) appearing in Eq. (28) are trivial, and the kinetic-energy evolution operators can be simplified tô The real-and imaginary-time kinetic-energy evolution operators satisfy a relationR M (U τ +a , U τ ) = R E (U τ +a , U τ ) −i with superficial similarities to Eq. (22); however, this relation between coordinate-basis matrix elements ofR M andR E does not imply that the eigenvalues ofR M andR E satisfy analogous relations.
In general the real-time transfer matrixT M in Eq. (27) is distinct from the unitary time-evolution operatorÛ defined by Eq. (22). The fact thatT M =Û is not by itself problematic, as the validity of results obtained using real-time LGT calculations only requires that a continuum limit exists in which the energy spectra associated withT M andÛ agree. However, ifT M is non-unitary for all limits of the action parameters then it is not possible to achieve a continuum limit in whichT M coincides witĥ U and the real-and imaginary-time energy spectra will contain unphysical differences. This undesirable scenario corresponds to the non-commutation of limits shown in Fig. 1 and can arise in practice for theories with compact variables such as the quantum rotator discussed in Sec. II A. For non-compact scalar field theory, this scenario does not arise and the unitarity of time-evolution associated with the standard finite difference discretization of the continuum Minkowski action has been demonstrated in 0+1D calculations in Refs. [4][5][6][7]. A proof of the unitarity of the real-time transfer matrixT M for scalar field theory proceeds analogously to the SHO case in Sec. II A and is detailed in Appendix A.
Below, the unitarity and continuum limits ofT M will be studied for specific LGT actions. Proofs of (non-)unitarity ofT M are facilitated by the observation that U |T MT † M |U , which must be proportional to a delta function δ(U, U ) forT M to give rise to unitary physics, satisfies Similarly, Assuming that V (U ) is real here and below, it follows

C. The real-time Wilson action
For the gauge groups U (N ) or SU (N ), the (Euclidean) Wilson action is defined by [11] S E,W (U ) = 2 g 2 x µ<ν Re Tr (1 − P x,µν ) , For the case of G = SU (N ), the bare coupling g is normalized so that the "naive continuum limit" obtained by defining U x,µ = e iaAµ(x) and taking the a → 0 limit of Eq. (32) agrees with the continuum Yang-Mills action for A µ (x) ∈ g with the same gauge coupling g. For the case of G = U (1), the rescaled coupling e = g/ √ 2 has the correct normalization for the naive continuum limit of Eq. (32) to match the continuum U (1) gauge action with gauge coupling e. The eigenvalues of the "plaquette" P x,µν can be represented as x,µν = 0 mod 2π. The Wilson action can be represented in terms of these eigenvalues as The quantum rotator discussed in Sec. II A is in fact equivalent to U (1) LGT using the Wilson action in (1 + 1)D, and the arguments in that section explicitly demonstrate thatT M is non-unitary in this case. As discussed in Ref [10] and reviewed here, character expansion methods can be used to demonstrate non-unitarity of the Wilson gauge action more generally for G = U (1) or G = SU (N ) in arbitrary spacetime dimensions. A similar check for (non-)unitarity can be applied to any LGT action that only depends locally on the P x,µν in a way that can be decomposed into a kinetic energy piece that is a function of "timelike plaquettes" P x,0k on each timeslice and a potential energy piece that is a function of "spacelike plaquettes" P x,ij on each timeslice. Since timelike plaquettes are given in temporal gauge by U x,k U † x+a0,k , the corresponding kinetic-energy evolution operatorR M (U, U ) only depends on the products U x,k U † x,k . The kinetic-energy evolution operator associated with any plaquette LGT action can therefore be represented using a character expansion of the form where r labels the representations of the gauge group G, d r denotes the dimension of representation r, and χ r (U ) is the character of the group element U ∈ G in representation r; see for example Refs. [69,73] for further discussion in the context of LGT. The coefficients of the character expansion c M r (g 2 ) are functions of only the gauge coupling g. Using the character orthogonality properties and the characters χ r (U ) ≡ x,k χ r x,k (U x,k ) can be seen to be eigenfunctions ofR M , The eigenvalue associated to each eigenfunction χ r (U ) is therefore the product of the character coefficients c M r x,k (g 2 ) associated with the representation r x,k specified by χ r (U ) for each gauge link. This implies that the kinetic-energy evolution operatorR M will be unitary if and only if all c M r (g 2 ) are unit-norm complex numbers. The character expansion coefficients can be obtained using character orthogonality, where N L = (D − 1)(L/a) D−1 is the number of spatial links U x,k at fixed t. Unitarity can then in principle be checked for particular actions.
The kinetic and potential energy terms K W and V W for the Wilson action are given by where k, l = 1, . . . , D − 1 denote spatial indices. The Minkowski Wilson action obtained by combining K W and V W with a relative minus sign and multiplying by a to remove the factors of 1/a above is given by The kinetic-energy evolution operatorR M,W associated with the Minkowski Wilson action is then given in temporal gauge bŷ Applying the character expansion in Eq.
where the character expansion coefficients for the Wilson action are given by The kinetic-energy evolution operatorR M,W is exactly unitary if and only if |c M,W r (g 2 )| = 1 for all r. However, a change in the overall normalization of the path integral uniformly rescales the magnitudes of all character expansion coefficients while leaving all operator expectation values invariant. Therefore if there is a particular choice of overall path integral normalization that leads to a unitary real-time transfer matrix then the corresponding action is compatible with unitary real-time evolution. 2 To demonstrate the non-unitarity of the Wilson gauge action, we therefore explicitly analyze ratios of the Wilson action character expansion coefficients for the groups U (1) and SU (N ) with N ∈ {2, . . . , 9} below. These ratios are independent of the overall path integral normalization. In all cases, the ratios of Wilson action character coefficients are found to have non-unit magnitudes almost everywhere in g 2 and in particular in the would-be continuum limit g 2 → 0, which is sufficient to establish that R M,W is non-unitary and cannot be made unitary by adjusting the path integral normalization. As discussed below Eq. (30), it follows thatT M,W is non-unitary and further cannot be made unitary by adjusting the path integral normalization.
In the simplest case of G = U (1), group elements U ∈ U (1) can be represented as U = e iφ , the Haar measure is simply dφ 2π , and there are only one-dimensional representations with characters χ r (e iφ ) = e irφ . Eq. (43) therefore gives 2 We thank Henry Lamm for bringing this point to our attention.  (2) includes an infinite number of singularities that accumulate as g 2 → 0 and is replaced by a gray background for g 2 < 0.1. The Wilson action for both SU (2) and SU (3) can be seen to have a non-unitary coefficient ratio almost everywhere in g 2 . In the SU (2) case, the limit g 2 → 0 is ill-defined due to the repeated singularities, while for SU (3) there is an apparently well-defined limit (explained by a stationary phase expansion in Appendix B) which is non-unitary. The HFK action for both SU (2) and SU (3) gauge theory has a unitary ratio for all choices of g 2 .
The modified Bessel functions I r (−i/e 2 ) are oscillating functions of e 2 with vanishing magnitude and increasingly rapid oscillations in the e 2 → 0 limit associated with the continuum limit of Euclidean U (1) LGT. The form of these functions immediately gives that |c (e 2 )| = 1 almost everywhere in e 2 and for r > 0 -these ratios do have unit norm for particular choices of e 2 and r as seen for the fundamental representation in Fig. 2, but not for all r at any e 2 ). It is also possible to directly confirm non-unitarity in the continuum limit e 2 → 0, which is not proportional to the identity integral kernel as e 2 → 0.
Analogous explicit results can be obtained for SU (2) gauge theory, where the irreps are labeled by half-integers j or equivalently integers r = 2j and have dimension d r = r + 1. Denoting the eigenvalues of group elements U by e iφ and e −iφ , the Weyl character formula gives χ r (U ) = sin((r + 1)φ)/ sin(φ). For functions of the eigenvalues, the Haar measure corresponds to dφ π sin 2 (φ), and the character expansion coefficients are therefore given by Explicit results for the SU (N ) character expansion coefficients for the imaginary-time Wilson action are known [73] and can be used to derive analogous results for the general SU (N ) real-time Wilson action. For the imaginary-time Wilson action, the kinetic-energy evolution operatorR E,W has a character expansion given bŷ The coefficients of this expansion are given analogously to Eq. (38) by The character identity χ(U ) = χ(U † ) * along with invariance of the trace and Haar measure under the transfor- It is proven in Ref. [12] that c E,W r (g 2 ) > 0 by expanding the exponential in Eq. (47) and observing that c E,W r (g 2 ) is equal to a positive number (for g ∈ R) times a counting factor related to the number of times the irrep r appears at a given order of the expansion. Comparing the definition of c M,W r in Eq. (43) and the coefficients c E,W r , one finds The explicit forms of c E,W r (g 2 ) for generic SU (N ) gauge groups presented in Ref. [73] can thus be used with Eq. (49) to obtain results for c M,W r (g 2 ). We numerically compute results for c M,W these numerical results are shown in Fig. 3 and indicate that, although c (g 2 )| = 1 for all g 2 and for the limit g 2 → 0, and therefore the real-time SU and non-unitarity for general N . In the appendix the g 2 → 0 limit is also analyzed using the stationary phase approximation.

III. UNITARY REAL-TIME LGT ACTIONS
Ref. [10] directly constructs a unitary real-time evolution matrix, here denotedT M,HF K , whose spectrum in the naive continuum limit is related to the spectrum of the usual Euclidean Wilson transfer matrix. Unitarity is guaranteed by defining the kinetic energy evolution operator,R M,HF K , in terms of an explicitly unitary spectrum. The spectrum is determined by the character expansion with coefficients denoted c M,HF K r (g 2 ) that are constructed to satisfy |c M,HF K r (g 2 )| = 1. Unitarity for any choice of potential, including the Wilson potential V W , then follows from Eq. (30).
An action S HF K (U ) can be formally defined using the character expansion forR M,HF K [10], as reviewed below in Sec. III A. The resulting series, however, is not absolutely convergent in all gauge field configurations as discussed below and cannot be numerically evaluated using either systematically improvable truncations or Monte Carlo sampling techniques. For U (1) gauge theory, a simple path integral contour deformation can be found that provides an absolutely convergent representation of path integrals involving the HFK action (see Sec. IV), but for SU (N ) gauge theory it is challenging to find an analogous contour deformation that gives convergence.
In this work, an alternative real-time LGT action based on analytic continuation of the Euclidean heat-kernel LGT action [13] is therefore studied; the Euclidean construction is reviewed in Sec. III B, and the Minkowski version is introduced in Sec. III C. This Minkowski heatkernel, or HK, action is also formally defined by a divergent series, and it is difficult to find path integral contour deformations that result in absolutely convergent representations of this action. However, a real-time modified heat-kernel, or HK, action is introduced in Sec. III D for which path integral contour deformations can be used to construct absolutely convergent representations of the U (1) and SU (N ) unitary real-time LGT actions. The real-time HK action is thus in principle suitable for evaluation of unitary real-time LGT dynamics for these gauge groups. Analytic results for unitary and non-unitary actions computable in (1 + 1)D are discussed in Sec. III E.

A. The real-time HFK action
We begin by reviewing in this section the construction of the HFK action given in Ref. [10]. Positivity of c E,W r (g 2 ) permits the definition of a set of character coefficients that satisfy |c M,HF K r (g 2 )| = 1 by construction. Explicitly,R M,HF K is then defined from the character expansion in temporal gauge bŷ (51) A gauge-invariant action whose kinetic energy term is related toR M,HF K by Eq. (29) can then be defined as The real-time transfer matrix associated with the HFK action iŝ and by unitarity ofR M,HF K , the entire transfer matrix is unitary.
Positivity of the Euclidean coefficients c E,W r (g 2 ) guarantees the existence of an operatorK W satisfyinĝ whereV W (U, U ) ≡ V W (U )δ(U, U ) and the imaginarytime Wilson transfer matrix is given similarly bŷ The relationship betweenT M,HF K andT E,W in Eqs. (54)-(55) is the expected analytic continuation relating discretized real-and imaginary-time evolution operators. Proving the existence of the continuum limit for LGT in either real or imaginary time is outside the scope of this work, but from the forms ofT M,HF K andT E,W in Eqs. (54)-(55) and the Lie-Trotter product formula one can expect thatT E,W andT M,HF K will satisfy the commutative diagram in Fig. 1 in the continuous time limit. The HFK action therefore provides a theoretically suitable action for real-time LGT.
A practical complication associated with the HFK action is that the character expansion in the definition of the action does not provide an absolutely convergent function of the gauge field U across the entire group domain. In particular, the constraint |c M,HF K r (g 2 )| = 1 and the fact that each character is normalized by dU |χ r (U )| 2 = 1 and χ r (U ) cannot vanish for all U ∈ G implies that the r-th term in the sum does not vanish as r → ∞ using any enumeration of the representations of G, at least for a set of non-zero measure in G. 3 Since |c M r (g 2 )| = 1 is also precisely the condition required for unitarity of a real-time LGT transfer matrix, it is clear that this divergence is a generic feature of the character expansions for unitary real-time LGT actions. A simple proof that this character expansion diverges for any 3 Stronger statements can be proven for particular choices of the gauge group. For G = U (1) the characters χr(U ) = e irφ satisfy |χr(U )| = 1 and therefore the sum in Eq. (52) diverges for all U . For G = SU (2), there are combinations of r and U satisfying χr(U ) = 0; however, the Weyl character formula can be used to explicitly show that limr→∞ χr(U ) is oscillatory and not equal to zero for any U and therefore that the sum in Eq. (52) diverges for all U . The explicit SU (3) character formula in Ref. [74] can be used to analogously prove that limp,q→∞ χ (p,q) (U ) is oscillatory and not equal to zero for any U where the irreps are enumerated using p, q ∈ Z.
real-time LGT action with kinetic and potential energy densities that only depend locally on timelike and spacelike plaquettes, respectively, is presented in Appendix C. Path integral contour deformations can be used to improve this convergence as discussed in Sec. IV, where a simple contour deformation is obtained for which the HFK action for U (1) LGT is represented by an absolutely convergent character expansion. Obtaining an analogous contour deformation providing an absolutely convergent representation of the SU (N ) HFK action suitable for numerical calculations is challenging. For this reason, alternate real-time LGT actions are introduced below for which the kinetic-energy evolution operator takes a simpler form for both G = U (1) and G = SU (N ), and a contour deformation that leads to an absolutely convergent path integral representation can be constructed in Sec. IV.

B. The imaginary-time heat-kernel action
The central role of the kinetic-energy evolution operator in establishing unitarity of the real-time transfer matrix suggests that it may be easier to construct unitary real-time LGT actions in the eigenbasis of the LGT kinetic energy operator. A kinetic energy operator that generalizes the Laplacian operator for non-compact coordinates to compact variables including the quantum rotator as well as SU (N ) gauge fields was used by Kogut and Susskind to construct a lattice gauge theory Hamiltonian [67] which includes the same potential energy term as the Euclidean Wilson action but a different kinetic energy term that is expected to be equivalent to the Wilson action kinetic term in the continuous-time limit. Defining oper-atorsL A x,k by the commutation relation where the t A are Hermitian generators of g normalized by Tr(t A t B ) = 1 2 δ AB , and the operator∆ appearing in Eq. (56) is defined bŷ This can be recognized as the quadratic Casimir operator (up to a sign) and for gauge groups U (1) and SU (N ) acts on functions of the gauge field by∆ where ∆ x,k is the Laplace-Beltrami differential operator as shown in Ref. [13]; see also Refs. [69,75]. Denoting the character expansion for f (U x,k ) by , the action of the Laplace-Beltrami operator on f (U x,k ) is given by where C (2) r is the eigenvalue of the quadratic Casimir operator for representation r of G.
Ref. [13] constructs eigenfunctions of the kineticenergy time-evolution operator e −aK by solving the diffusion or "heat-kernel" equation obtained by analytic continuation of the Schrödinger Omitting the prefactor in Eq. (56) for simplicity, the heat-kernel solution K E (U, τ ) is defined as the solution to the differential equation with boundary condition K E (U, 0) = δ(U, 1). The heatkernel solution provides an integral kernel for∆, where the state |U x,k ⊗ 1 is defined by assigning U x,k to link ( x, x +k) and the identity to all other links and the first equality follows from commutativity of the Laplace-Beltrami operator with group multiplication [13]. An integral representation for the (temporal-gauge) Kogut-Susskind kinetic-energy time-evolution operator is therefore given by The factor U x,k U † x,k can be identified with P x,0k in temporal gauge in order to construct a gauge-invariant expression. In Euclidean spacetime it is convenient to identify the kinetic and potential energy terms as identical functions of timelike and spacelike plaquettes respectively in order to obtain a LGT with D-dimensional rather than (D − 1)-dimensional hypercubic symmetry. Such an isotropic Euclidean heat-kernel action is defined by [13] which is gauge-invariant by the gauge-invariance of the Laplace-Beltrami operator (as well as more explicit arguments below). The kinetic-energy piece of the transfer matrix associated with the heat-kernel action is given by construction asR E,HK = e −aK . Its character expansion in temporal gauge can be obtained from Eq. (59) and Eq. (62) aŝ from which the character expansion coefficients can be identified as Positivity of the quadratic Casimir eigenvalues C r gives c E,HK r (g 2 ) > 0, from which the positivity ofR E,HK and T E,HK and formal existence of a unitary time-evolution operatorÛ HK = [T E,HK ] i follows.
The explicit construction of the Euclidean heat-kernel action requires an explicit form for K(U, τ ) for particular G. For the non-compact group G = R the heat-kernel equation reduces to the usual diffusion equation which has the well-known solution with the normalization constant N = 1/ √ 2πτ fixed by the boundary condition K E,R (x, 0) = δ(x) at τ = 0 and the evolution equation in Eq. (60) for later τ . In the heat-kernel solutions derived for various groups below, normalizing constants are not distinguished for conciseness; distinct normalizing constants apply to each case with the appropriate normalization clear from context. For G = U (1), the heat-kernel equation takes an analogous form and has a solution given by a sum of Gaussian terms of the form in Eq. (67) over coordinates x = φ + 2πn for all possible integers n. Noting that the appropriate coupling constant normalization is given from Eq. (63) by τ = e 2 , the form of the U (1) heat-kernel required for Euclidean LGT calculations is given by with the normalization N = 1/ √ 2πe 2 fixed by K E,U (1) (U, 0) = δ(U, 1). This heat kernel weight corresponds to the Villain action [76]. For G = SU (N ), additional factors arise in coordinate descriptions of the heat-kernel equation from the non-trivial metric of the Riemannian manifold associated with G. The construction of the heat-kernel for this case is detailed in Ref. [13], and the solution is given in terms of the eigenvalues e iφ1 , . . . , e iφ N of U by where {n} describes a sum of infinite sums ∞ n A =−∞ and N is fixed by K E,SU (N ) (U, 0) = δ(U, 1) as in the U (1) case. The integers n A , A = 1, . . . , N are subject to a constraint A n A = 0 analogous to the eigenvalue phase constraint A φ A = 0 that ensures that det(U ) = 1. Note that φ A +2πn A is treated as a non-compact variable in the heat-kernel action in the sense that A φ A = 0 and A n A = 0 are enforced rather than A φ A = 0 mod 2π [13]. Eq. (70) is invariant to permutation of eigenvalues e iφ1 , . . . , e φ N ensuring that the kernel does not depend on the (unphysical) ordering of eigenvalues. Gauge transformations act by matrix conjugation on P x,µν and therefore do not affect the unordered set of eigenvalues e iφ1 , . . . , e iφ N , leaving Eq. (70) invariant.

C. The real-time heat-kernel action
The analog of the heat-kernel equation that describes the real-time evolution of U ∈ G in the absence of a potential is the free Schrödinger equation on G, This Schrödinger equation is related to Eq. (60) through analytic continuation using the identification τ = it. The solution to the Schrödinger equation in Eq. (71) can therefore be immediately obtained through analytic continuation of the solution to the Euclidean heat-kernel equation, Analytic continuation of Eq. (62) gives The Minkowski heat-kernel action defined by analytic continuation of Eq. (63), therefore has a kinetic-energy time-evolution operator with integral kernel which therefore satisfiesR M,HK = e −iaK . The unitarity ofR M,HK follows immediately from the Hermiticity of the Laplace-Beltrami operator andK.
Unitarity can also be directly verified through analytic continuation of the Euclidean heat-kernel character expansion in Eq. (64), which gives The coefficients c M,HK r (g 2 ) of the character expansion of R M,HK analogous to Eq. (43) are given by It follows that |c M,HK r (g 2 )| = 1, which implies that R M,HK is unitary. Unitarity ofT M follows for any choice of real potential from Eq. (30). In particular, this implies that the Minkowski heat-kernel action defined in Eq. (74) leads to a unitarity LGT time-evolution operatorT M,HK .
Comparing Eq. (65) and Eq. (77), the Minkowski and Euclidean heat-kernel actions are further seen to satisfy It follows from this that the real-and imaginary-time kinetic-energy evolution operators associated with the heat-kernel action satisfyR M,HK = e −iaK =R i E,HK . The corresponding real-and imaginary-time transfer matrices can therefore be expected to satisfy the commutative diagram shown in Fig. 1. Since the Euclidean heatkernel and Wilson actions are expected to be equivalent in the continuum limit, the Minkowski heat-kernel action is a unitary real-time LGT action that should be equivalent to the HFK action in the continuum limit.
For the U (1) gauge group, the explicit form of the Minkowski heat-kernel solution is given by where the normalizing constant N is related to the Euclidean normalizing constant by the analytic continuation e 2 → ie 2 . For the SU (N ) gauge group, the Minkowski heat-kernel solution is where J ({φ}, {n}) is given in Eq. (70) and {n} again denotes a set of integers n A with A = 1, . . . N subject to the constraint A n A = 0. The SU (N ) normalizing constant is similarly related to the Euclidean case by the analytic continuation g 2 → ig 2 . The infinite sums in Eqs. (79)-(80) are divergent since the summand does not vanish in the n → ±∞ limit. This non-convergence of oscillating sums defining the path integral weights mirrors the problems faced in the HFK action. Motivated by a saddle point expansion, a possible approach to this problem is explored in Appendix D in which sums over n A are truncated to n A = 0. This truncation breaks unitarity at non-zero lattice spacing, but the dominance of the n A = 0 terms in the saddle point approximation might suggest that this unitarity breaking is removed in the continuum limit. However, it is shown in the appendix that unitarity is not in fact recovered in the analytically solvable case of (1 + 1)D U (1) LGT, indicating that this truncated version of the heat-kernel approach is not a useful starting point for real-time LGT and underscoring the importance of preserving unitarity in real-time LGT actions.

D. The modified real-time heat-kernel action
The convergence issues of the HFK and Minkowski heat-kernel actions and the lack of a scaling limit for the truncated heat-kernel action motivate the definition of a real-time modified heat-kernel action that includes the heat-kernel kinetic term and the Wilson potential term, where for the U (1) case J ({φ x,0k }, {n x,k }) is replaced by unity, and P † x,ij = P −1 x,ij has been used to express Eq. (81) in a form suitable for discussions of contour deformations in Sec. IV. Since the kinetic term is the usual heatkernel one,R M,HK =R M,HK = e −iaK is unitary and regardless of the different choice of potential in S HK (U ) and S HK (U ) both are therefore unitary real-time LGT actions. The real-time transfer matrix associated with S HK (U ) is given bŷ and is manifestly unitary. This is the expected form for a discretized real-time evolution operator and by the Lie-Trotter product formulaT t/a M,HK should converge to e −it(K+V W ) in the continuous-time limit. Since T τ /a W is expected to converge to e −τ (K+V W ) in the continuoustime limit,T E,W andT M,HK are expected to satisfy the commutative diagram shown in Fig. 1. It can also be seen at the level of the action that the real-time HK action is equivalent to the (truncated) Minkowski heatkernel action in the g 2 → 0 limit where in the stationary phase approximation the Wilson potential approaches The real-time HK action therefore provides another unitary real-time LGT action with the same naive continuum limit as the Minkowski heatkernel and HFK actions.
It is noteworthy that all of the Minkowski actions above include at least a sign difference between kinetic and potential terms arising from Eq. (23) and are therefore isotropic in (D − 1) spatial dimensions but not D dimensions. There is no subgroup of the Lorentz group that provides a symmetry of any of the Minkowski LGT actions described above, and the real-time HK action shares the same symmetries as the Minkowski heat-kernel and HFK actions. Although the real-time HK action also shares the same downside as the HFK and Minkowski heat-kernel actions -a definition involving a divergent sum -it is demonstrated below in Sec. IV that path integral contour deformations can be used to construct an alternative representation of S M,HK in which the sum defining the kinetic energy is absolutely convergent. This representation of the real-time HK action provides a welldefined starting point for numerical investigations of realtime LGT with unitary time-evolution at non-zero lattice spacing.

E. Exact results in (1+1)D
In (1 + 1)D with open boundary conditions (OBCs), both the Wilson and heat-kernel Euclidean actions do not include potential terms and SU (N ) gauge theory is analytically solvable [13,77,78]. These solutions can be straightforwardly analytically continued in order to compare results for observables constructed using Minkowski actions leading to non-unitary and unitary time evolution.
A simple, non-trivial observable in (1 + 1)D is the Wilson loop, where x,µ∈∂A U x,µ denotes an ordered product along the boundary of the two-dimensional rectangular region A with spatial extent L and temporal extent denoted τ in Euclidean spacetime and t in Minkowski spacetime. Wilson loops can be interpreted as propagators for static quark-antiquark pairs separated by a distance L. In Euclidean spacetime, for a lattice gauge theory with positive-definite transfer matrix, Wilson loops satisfy the spectral representation where E n > 0 is the energy of the n-th energy eigenstate with appropriate quantum numbers for describing the static quark-antiquark system and Z n is the overlap factor onto the n-th energy eigenstate. An ideal spectral representation for the corresponding Minkowski theory is obtained by analytic continuation with τ = it, A suitable real-time LGT action should give rise to a spectral representation of the form Eq. (85) with energies and overlap factors that may differ at non-zero lattice spacing but should agree in the continuum limit for states with energies much below the lattice cutoff. With a local action S E = x L E (P x ) that depends only on the plaquettes P x ≡ P x,01 , Wilson loop expectation values in (1+1)D with OBCs take a simple factorized form [77,78] This can be simplified using the character expansion of e −L E (P ) , where r = f denotes the antifundamental representation and has dimension d f = N . The condition c E r (g 2 ) > 0 required for positivity ofT E guarantees that this can be expressed as which is the expected spectral representation for a single state with energy E 0 = σL. An analogous relation holds in Minkowski spacetime, The condition c M r (g 2 ) = [c E r (g 2 )] i is sufficient for Eq.
while the real-time Wilson action expectation value becomes This does not correspond to a real-time spectral representation of the form Eq. (85) for any value of g 2 . In particular lim g 2 →0 W A

SU (2)
M,W does not exist, as shown in Fig. 4. This demonstrates that even though the g 2 → 0 limit of the real-time transfer matrix exists in the SU (2) case, well-defined continuum limits for real-time observables do not necessarily exist. Similar results can be obtained for SU (N ) gauge groups in (1 + 1)D, and again the lack of a continuum limit where the real-time transfer matrix becomes unitary is associated with the lack of a well-defined continuum limit for demonstrating that the HFK action leads to unitary results that correspond to the τ = it analytic continuation of the corresponding Euclidean Wilson action LGT results. In higher dimensions where LGT potentials are non-trivial, this exact t = iτ correspondence will not hold in LGT but should emerge in the continuum limit.
The Euclidean heat-kernel action leads to a different form for the (1 + 1)D LGT string tension [13] σ SU (N ) It can be explicitly seen that demonstrating that W A M,HK has the expected spectral representation associated with unitary time evolution in LGT. This demonstrates that lattice artifacts leading to differences between the energies appearing in the Minkowski and Euclidean spectral representations are also absent for the heat-kernel action in (1 + 1)D, a feature which does not persist in higher dimensional LGT where potential operators are present. The same result can be derived directly by inserting the character expansion coefficients from Eq. (77) into Eq. (89) and using C (2) In the U (1) case C .
The real-time HK action S M,HK defined in Eq. (81) uses the Wilson potential and heat-kernel kinetic terms. Since there is no potential term present for LGT in (1+1)D, S M,HK coincides with the Minkowski heat-kernel action S M,HK and leads to Wilson loop results identical to Eq. (95). In higher dimensions S M,HK will not coincide with the S M,HK exactly, but S M,HK still leads to a unitary time-evolution operator and therefore Wilson loop expectation values with a unitary spectral representation. "Lattice artifacts" take the form of modifications to the energy spectrum and overlap factors that vanish as g 2 → 0 and are not expected to spoil the scaling properties of the continuum limit (unlike the non-unitary truncated heat-kernel action investigated in Appendix D).

IV. REAL-TIME LGT PATH INTEGRAL CONTOUR DEFORMATION
Real-time LGT path integrals in more that two dimensions cannot be calculated analytically with current techniques. In order to perform Monte Carlo calculations of real-time LGT path integrals, convergent representations of path integrands involving the real-time actions under study must be constructed. Path integral contour deformation techniques previously used to tame sign problems are used here to construct such con-

A. Sign problems and path integral contour deformations in real time
For a compact Lie group the absolute value of the Minkowski path integral weight and measure, provides a well-defined probability measure for performing Monte Carlo sampling, unlike in the case of noncompact scalar fields in real time discussed in Refs. [3][4][5][6][7]. However, the "reweighting" approach defined by sampling with respect to DU generically leads to an exponential sign-to-noise (StN) problem: the variance of Monte Carlo estimates of the average phase factor e iS M required to determine the full partition function grows exponentially as the size of the system is increased. Standard arguments for sign problems in Euclidean spacetime compare the "phase-quenched" partition function, defined by the integral of the absolute value of the path integrand, to the full partition function; if ignoring phase fluctuations leads to a free energy f Q distinct from the free energy f in the full theory, then the variance of estimates of the partition function will be given for large volumes by e −2f Q V β − e −2f V β and positivity of the variance requires f Q ≤ f . The signal-to-noise (StN) ratio associated with computing the ratio of phase-quenched to full partition function will therefore scale as e −(f −f Q )V β and will decrease exponentially with increasing Euclidean spacetime volume [79][80][81][82][83]. Analogous arguments can be extended to Schwinger-Keldysh correlation functions [7].
For the case of Schwinger-Keldysh LGT in particular, the phasequenched theories describing the Minkowski segments are defined by trivial path integral weights and can be given a thermodynamic interpretation as a Euclidean action that is exactly zero, corresponding to the strong coupling limit. The phase-quenched Schwinger-Keldysh partition functions will thus scale with volume as Numerical calculations using the unitary Minkowski actions studied above -the HFK, real-time heat-kernel (HK), and modified real-time heat-kernel (HK) actions -face the additional challenge that in all cases e iS M (U ) is formally defined by an infinite sum that does not converge for fixed U and its average over the distribution DU cannot be calculated using Monte Carlo methods even in principle. In the (1 + 1)D examples discussed above, where exact results are available, it is clear that performing the gauge field integral before the infinite sums provides convergent results that match the desired spectral representations obtained by analytic continuation of Euclidean results. One could imagine ordering the summation outside integration and explicitly performing Monte Carlo integration over e iS M (U ) for each combination of terms in the HFK character expansion or integers in the heat-kernel sum below a specific cutoff, but the need to perform O(V L T ) Monte Carlo calculations and the systematic uncertainties in the combined result arising from truncating the infinite sums make this approach undesirable. If the combined sum-integral defining the path integral over e iS M for a unitary action could be rendered absolutely convergent, it would instead be possible to perform the sums and integrals together in a joint Monte Carlo calculation. This is achieved below using path integral contour deformation methods.
The foundation of path integral contour deformations is the complex analysis result that for holomorphic integrands Oe iS M , the integration contour of the path integral can be deformed in order to affect the StN properties of the integral without modifying the total integral value. Previously, path integral contour deformations have been used to improve the sign and associated StN problems affecting real-time (0 + 1)D quantum mechanics models [3][4][5][6][7], as well as imaginary-time theories of scalars and fermions , U (1) gauge theory [109][110][111][112][113][114], dimensionally reduced (single-or few-variable) non-Abelian gauge theory [104,[115][116][117][118][119][120], and recently large Wilson loops in (1 + 1)D SU (N ) gauge theory [9]. By modifying the integrand magnitude and phase, contour deformations also have the potential to improve the convergence problems highlighted above.
A lattice gauge theory path integral can be interpreted as an iterated integral over a set of compact, groupvalued variables. For the U (1) gauge group, these integrals can be written in terms of one angular variable φ ∈ [0, 2π] per U (1) gauge link. For SU (N ) gauge groups, deformations can be defined in terms of an angular parameterization of each SU (N ) gauge link, given by a set Ω = {φ a , θ b } of azimuthal angles φ a ∈ [0, 2π] and zenith angles θ b ∈ [0, π/2]. To be a valid contour deformation, endpoints must be handled properly for both φ and θ angles: for periodic φ angles (appearing in both the U (1) and SU (N ) parameterizations) any deformation that keeps the endpoints identified will be valid, while for non-periodic θ angles the endpoints must be held fixed. Further details on angular parameters and deformations of SU (N ) variables can be found in Ref. [9]. The Haar measure appearing in the path integral can be related to the natural measure on the relevant angular coordinates by DU = a dφ a b dθ b H(Ω) = DΩ H(Ω), where H(Ω) can be straightforwardly computed for particular angular parameterizations [9,121].
After rewriting the path integral in terms of the chosen coordinates, i.e. replacing the measure as above and replacing instances of U with U(Ω), contour deformation can be directly applied to the compact integration paths of each real-valued angular variable, potentially conditioned on other angular variables. For a valid deformation Ω → Ω(Ω) describing integration on a new manifold M, a generic integral can be deformed as where Cauchy's theorem gives the equality between the first and second line, the Jacobian J(U ) = J(Ω)H( Ω)/H(Ω) accounts for the change in Haar measure arising from the deformation, and the deformed gauge fieldŨ ≡ U ( Ω(Ω)) is a member of the complexified group. 4 To be valid, the deformation map must also be continuously connected to the identity map Ω Id (Ω) ≡ Ω. By deforming and then writing integration on the deformed manifold in terms of coordinates Ω in the same domain as the original integration, the net effect of contour deformation is to replace the integrand f (U ) with J(U )f ( U ) without modifying the integral value. Applying Eq. (98)  A similar approach can be applied when the path integral is defined on a Schwinger-Keldysh contour consisting of both Euclidean and Minkowski spacetime regions. Sign problems in such path integrals arise from the weights associated with the Minkowski region, suggesting that useful deformations can therefore be largely restricted to this region, up to boundary effects. Path in- 4 For example, U ∈ SL(N, C) for U ∈ SU (N ). tegral contour deformations can therefore be used to improve the sign and StN problems associated with calculations of e iS M (U ) in real-time and Schwinger-Keldysh LGT in complete analogy to previous applications to quantum mechanical models.

B. Convergent U(1) HFK path integrals
By changing the integrand magnitude, contour deformations can also render a divergent sum-integral absolutely convergent in certain cases, giving a prescription for evaluation. In particular, if integrating first then summing gives a well-defined result, it may be possible to apply contour deformations to each integral within the sum such that the joint sum-integral becomes absolutely convergent. Doing so preserves the value defined by the integrate-then-sum order of operations while enabling reordering or joint integration/summation. However, because the endpoints of the integration contour must remain fixed to define a valid contour deformation, there will always be a neighborhood of the endpoint that has similar convergence properties to the original integration contour. In order to define an absolutely convergent prescription for jointly performing sum-integrals below, an additional regularization of the infinite sum appearing in the HFK kinetic term analogous to a continuum Wick rotation is introduced below. In particular, the phases of U (1) gauge fields and plaquettes posses a 2π shift symmetry that is used to enforce cancellations between integral contributions from segments near both endpoints before taking the limit in which the Wick rotation regularization is removed. The resulting sum-integral is absolutely convergent for all points on the deformed contour after this cancellation is enforced and the limit is subsequently taken. This can be thought of as a rigorous coordinate-based approach to dealing with the identification of points related by the 2π shift symmetry, which is sufficient to avoid these endpoint singularities.
A simple contour deformation can be constructed that in this sense provides an absolutely convergent representation of the HFK action for G = U (1). Real-time LGT path integrals involving the HFK action can be represented as where Z M,HF K = DU {r} e iS M,HF K (U,r) and the indices {r} ≡ {r x,k : x ∈ sites, k ∈ {1, . . . , D − 1}} label the representations appearing in the character expansion for the HFK kinetic term involving temporal plaquettes P x,0k = e iφ x,0k . These indices can be considered auxiliary variables in the path integral, with a local action S M,HF K (U, r) defining the weights over U and r simultaneously. For G = U (1) this action is given by which is a convergent function of U x,µ and r x,k . As discussed in Sec. III A, however, the path integral in Eq. (101) involves sums over r x,k and when these are evaluated inside the integral (i.e. holding the gauge configuration U fixed) the large-r behavior is divergent. This prevents the sum-integrals defining O M,HF K from being evaluated using Monte Carlo methods because the absolute value of the weights do not give a well-defined probability measure. On the other hand, exchanging the order of summation and integration in Eq. In the G = U (1) action above, integration over the gauge variables produces strong cancellation in the term e ir x,k φ x,0k when any of the r x,k are taken large, suppressing these terms in the sum. A similar suppression occurs for G = SU (N ), though it takes a more complex form due to the non-Abelian characters χ r appearing in the integrand. Although this doesn't enable joint Monte Carlo sampling of U and r (the ordering of summation/integration would be lost), it does provide a starting point for contour deformations. In order to define an absolutely convergent representation of real-time LGT path integrals involving the HFK action, we first regularize the divergent sums by replacing the factor of [c E,W r (g 2 )] i appearing in the HFK action with [c E,W r (g 2 )] e iθ where θ ∈ [0, π/2] parameterizes a transformation between the Euclidean Wilson action with θ = 0 and the Minkowski HFK action with θ = π/2. For generic θ, the kinetic term appearing in the U (1) HFK action becomes For θ ∈ [0, π/2) the first factor I r x,k (1/e 2 ) cos(θ) vanishes faster than exponentially at large |r x,k | [65], leading to absolute convergence of path integrals of the form Eq. (103) involving this transformed kinetic term. It is only exactly at the Minkowski limit corresponding to θ = π/2 that this term is equal to unity and absolute convergence is lost. Before taking the Minkowski limit, we perform an r x,kdependent deformation of the φ x,0k integration contour defined by the map which can be interpreted as a distinct contour deformation for each term in the sum in Eq. (103). "Constant vertical deformations" analogous to this one have previously been used to reduce sign problems in several applications, see for example Refs. [102,103], typically including one or more free parameters that define the imaginary shift and can be optimized to determine an integration contour with minimal phase fluctuations. 5 Constant vertical deformations like Eq. (105) are only valid for periodic compact variables: the contours parallel to the imaginary axis required to complete a closed path involving the original and deformed integration contours differ by a shift equal to the domain of periodicity and have opposite orientation, causing integrals along these contours to cancel [9]. For θ ∈ [0, π/2), the absolute convergence discussed above allows this cancellation to be performed. After enforcing the cancellation of the perpendicular contours to leave only the integral along the shifted contour corresponding to Eq. (105), the Minkowski limit can be taken. For all points on the shifted contour, exponential damping factors appear in the transformed kinetic energy term from e i φ x,0k r x,k = e iφ x,0k r x,k e −|r x,k | . As detailed below, these exponential damping factors are sufficient to render path integrals involving the HFK action absolutely convergent for all points on the shifted contour, even in the Minkowski limit. In order to apply the transformation defined in Eq. (105) for timelike plaquette phases, it is necessary to first perform a change of variables from the original set of path integral variables {U x,µ } to a set of variables that includes the timelike plaquette phases. Although it is not possible to perform a one-to-one change of variables from links to plaquettes, it is possible to transform to the set of variables {U x,0 , P x,0k } and an additional fixed gauge field U (L T , x),k acting as a boundary condition in the real-time direction. For the OBC case used in Sec. III E as well as Secs. IV E-IV F below, U (L T , x),k = 1. Periodic boundary conditions in real-time can be viewed as infinite-temperature Schwinger-Keldysh contours and include singular observables not suitable for numerical simulation. 6 For Schwinger-Keldysh contours with nontrivial Euclidean extent, U (L T , x),k can be taken to be the first link on the Euclidean side of the boundary with each Minkowski region as detailed in Sec. IV D. For U (1) gauge theory, it is convenient to represent the boundary field as U (L T , x),k = e ib x,k , the gauge field as U x,µ = e iax,µ , and the plaquette as P x,µν = e iφx,µν as above. The relation can be used to solve for the spacelike components of the gauge field phase as The spacelike plaquette phase can then be expressed as a function φ x,ij (φ x,0k , A x,0 , b x,k ) by inserting Eq. (107) into Eq. (106) and the full spacelike plaquette further ob- Path integrands of the form in Eq. (103) can therefore be transformed from base coordinates where the measure is proportional to x,µ dA x,µ into coordinates with measure proportional to x,k dφ x,0k dA x,0 through a linear transformation with unit Jacobian. In this basis, the contour deformation defined by Eq. (109) can be simply applied to φ x,0k with A x,0 held fixed. The spacelike plaquette phase for the deformed contour is therefore given (with θ = π/2 taken after enforcing cancellation of perpendicular contours) by sign(r x+t0,0i ) + sign(r x+t0+î,0j ) − sign(r x+t0+ĵ,0i ) − sign(r x+t0,0i ) .
The HFK action evaluated for the transformed gauge field can therefore be expressed using P x,ij ≡ e iφx,ij ( φ x,0k ,Ax,0,b x,k ) and P −1 The r-dependent kinetic term [I r x,k (1/e 2 )] i e ir x,k φ x,0k has unit magnitude. The factors of P x,ij in the potential term only depend on {sign(r x,k )} and cannot have asymptotically diverging magnitude in the large-r limit. Any observable that is a function of gauge fields similarly can only depend on {sign(r x,k )}. The summand in Eq. (103) The absolute value weight | e iS M,HF K ( U (U,r),r) Z M,HF K | now defines a valid probability measure over the combined space (U, r) and Monte Carlo sampling techniques can be applied.
For G = SU (2), the sum over representations can be converted to a sum over integers r A x,k with A = 1, 2 satisfying A r A x,k = 0 of the form Eq. (110). The undeformed action in this case is given by where φ 1 x,0k = −φ 2 x,0k are the phases of the eigenvalues e iφ A x,0k of P x,0k and r x,k ≡ r 1 x,k . These weights cannot be brought to an absolutely convergent form through either a constant vertical deformation or through the more general affine transformations explored in Sec. IV C below. The exploration of more sophisticated contour deformation to transform the HFK action to an absolutely convergent representation for G = SU (N ) is left to future work.

C. Convergent SU(N) and U(1) HK path integrals
We now turn to the construction of a path integral contour deformation that provides absolutely convergent representations of path integrals using the real-time HK LGT action S M,HK (U ) that are suitable for numerical calculations using Monte Carlo techniques. Deformations for the case of both a U (1) and SU (N ) gauge group are explored.
Since the divergent sums appear in the heat-kernel kinetic term, which is expressed in terms of the eigenvalues of the plaquette, it is convenient to begin with a change of path integral variables to a set that includes P x,0k . For real-time systems, one is generally interested in fixed or open boundary conditions in time or Schwinger-Keldysh contours as discussed in Sec. IV D rather than periodic boundary conditions in real time. As discussed above in Sec. IV B, a boundary field configuration U ( x,L T ),k is available for each of these cases that be used to perform a change of variables from spatial links to timelike plaquettes. For both U (1) and SU (N ) gauge theory, it is possible to invert the relation and solve for the spatial gauge field in terms of the plaquettes, boundary field configuration, and temporal links as It follows that once U (L T , x),k is specified the path integral variables U x,µ are in one-to-one correspondence with the set of variables {P x,0k , U x,0 }. Further, the Jacobian for performing the change of variables to {P x,0k , U x,0 } is unity by the group multiplication invariance of the Haar measure. The plaquettes P x,0k can be diagonalized and written in terms of eigenvalues e iφ A x,0k and corresponding unitary eigenvector matrices V x,0k satisfying where φ A x,0k ∈ [−π, π] and, for G = SU (N ), A = 1, . . . , N with A φ A x,0k = 0 mod 2π. Together these provide a new set of path integral variables where V x,0k and U x,0 are parameterized according to some coordinate description of G. These variables are overcomplete in comparison to {P x,0k , U x,0 } since V x,0k can be multiplied by any unitary diagonal matrix without changing P x,0k , but this redundancy only affects the normalization of the path integral measure since the unitary phases can be "integrated in" without changing the integral value. In general an integral over an SU (N ) plaquette can thus be rewritten using these variables as [122, where dV x,0k is the Haar measure for U (N ). Using the variables in Eq. (115) as independent path integral variables with the measure appearing on the right-hand-side Eq. (116) permits the use of deformations of the φ A x,0k integration contours that leave V x,0k and U x,0 unaffected.
In analyzing the convergence of HK path integrals, it will be useful to consider the limit as the prefactor in the kinetic energy term of the HK action is rotated from −1 (Euclidean) to i (Minkowski), and for intervening points we use the prefactor −e −iθ , with θ ∈ [0, π/2], and represent a generic heat-kernel kinetic term for link (x, k) as For the gauge group U (1), this term is already in a suitable form for deforming as below; however, for the gauge group SU (N ) further rewriting of this term is helpful. The determinant constraint for SU (N ) implies that φ N x,0k and n N x,0k are not independent variables that are integrated over, and instead can be related to the other variables by φ N x,0k = − x,k . In terms of the N − 1 independent phases and integers, G x,k is a correlated Gaussian, where (120) New variables can be defined for this eigenbasis by This permits changing variables from φ A x,0k to the newly introduced ψ A x,0k with a Jacobian given by | det(S)| = (N −1)!. The heat-kernel term G x,k can thus be expressed as an uncorrelated Gaussian, To change variables in the path integral, the particular linear combinations involved in the definitions of m A x,k and ψ A x,k affect the domain of summation/integration. Unlike the set of n A x,k , for m A x,k it is necessary to write an iterated summation m 1 x,k m 2 x,k . . . with bounds and increments determined by the linear combinations involved in Eq. (121). A convenient choice of order is to sum over m 1 x,k outermost, then successively sum over m N −1 x,k , . . . , m 2 x,k in decreasing order. In this order, each m A x,k can be related to variables in outer sums by the (123) The domain of the m 1 x,k sum is thus Z, the domain of the m N −1 x,k sum is m 1 x,k + (N − 1)Z, and the domain of the remaining m A x,k sums are m A+1 x,k + AZ, indicating that the increment of the sum over m A x,k is always A though the offsets are generally functions of variables in outer sums. A similar iterated integral is required to properly cover the original domain when integrating over the new ψ A x,0k . These variables satisfy relations analogous to Eq. (123), and from these relations we find that each ψ A x,k should be integrated over a domain of width 2πA, with offsets determined by outer integration variables. Based on the form of the diagonalized heat-kernel term in Eq. (122) and the domains of the sums above, the pair . The shift in m A x,k can be absorbed into the infinite summation, and thus we are free to use this symmetry to shift the domain of integration of each ψ A x,k to the simple range [−Aπ, Aπ]. In terms of these new variables, Fig. 5 shows the m A x,kdependent deformation of ψ A x,0k that is used to construct a convergent representation of HK path integrals. 7 The contour for each value of m A x,k and a fixed angle θ consists of the union of three segments. The middle segment is The sign function is defined by sign(x) = x/|x| with the midpoint convention sign(0) = 0. The two outer segments of the contour are defined to be the linear intervals [−Aπ, α(m A x,k )] and [β(m A x,k ), Aπ] in the complex plane. To understand the convergence properties of this deformed path integral, it is useful to split each integration over ψ A x,k into an integral over a coordinate parameterization of the middle segment and a second integral over a coordinate parameterization of the outer segments. The first integral can be defined using a coordinate transformation ψ (0),A x,0k (ψ A x,0k , m A x,k ) relating the base coordinate interval [−Aπ, Aπ] to the middle segment of the deformed contour, The second integral can be defined using a similar coordinate transformation ψ (1),A x,0k (ψ A x,0k , m A x,k ) relating the base coordinate interval [−Aπ, Aπ] to the union of the two outer deformed contour segments, The transformed coordinates ψ are indexed by c A x,0k ∈ {0, 1}, abbreviated c in the ψ label, which in the full deformed path integral must be summed over. The eigenvalue phases of each deformed temporal plaquette can then be obtained by inverting Eq. (121), with the value of φ (c),N x,0k and n N x,k obtained by solving the constraints φ (c),N x,0k x,0k and n N x,k = − x,k . These deformed eigenvalue phases can be further used to define deformed timelike plaquettes as and from these the deformed gauge fields using Eq. (113), This contour deformation therefore defines a map U x,µ → U (c) x,µ (U, n) that depends non-locally on the gauge field as well as the integers n A x,µ ensuring 2π shift invariance of the heat-kernel action. Transformations of functions of U x,µ can be obtained immediately from the definition of U (c) x,µ (U, n), for example transformed spacelike plaquettes are given by This transformation does not leave the Haar measure invariant. The coordinate piece of the Jacobian is given by Including the explicit change of the Haar measure in the coordinates Eq. (116), the full Jacobian J (c) for each deformed contour segment is given by For G = U (1) the explicit change of the Haar measure is trivial and h (c) x,k should be replaced by unity in Eq. (133). The HK action given in Eq. (81) is a holomorphic function of φ A x,0k , with the independent variables V x,0k and U x,0 left undeformed. The spacelike plaquettes appearing in the action are functions of these variables, as defined in Eq. (131), and are also holomorphic in φ A x,0k .
These eigenvalue phases φ A x,0k are themselves holomorphic functions of the linear combinations ψ A x,0k introduced in Eq. (121), and results for expectation values of observables that are similarly holomorphic functions of ψ A x,0k (including in particular Wilson loops and other polynomial functions of U x,µ ) are therefore unchanged by deformations of the ψ A x,0k integration contours, as long as endpoints of each contour are held fixed and the unit determinant constraint for G = SU (N ) is not violated. Expectation values using the deformation defined by Eq. (126)-(127) of the HK path integral can therefore be expressed as where Z M,HK = DU {n} e iS M,HK (U,n) and {c} denotes a collection of sums 1 c A x,k =0 . Using the general prefactor −e −iθ , the modified heat-kernel action applied to the transformed variables takes the form where J is the SU (N ) heat-kernel factor defined in Eq. (70) applied to the transformed coordinates φ A x,0k , This factor can be naturally analytically continued using the complex sin function. The poles at 1 2 φ A − φ B + 2π(n A − n B ) = πk for k = 0 appear to cause problems with holomorphy, but these are exactly the cases when φ A = φ B mod 2π and these poles are cancelled by zeros of the deformed Haar measure appearing as the numerator of Eq. (133).
To analyze absolute convergence, we consider the phase-quenched path integral given for general rotation angle θ by The summation {c} is over a finite set and the integral DU is over a compact domain, so showing boundedness of the integrand and absolute convergence of the sum {n} is sufficient for convergence of the absolute sum-integral. The term in Eq. (137) involving the spatial plaquettes P (c) x,ij and the term DU |J (c) (U, n)| determining the deformed Haar measure are functions of only the variables {sign(m A x,k )} and therefore do not diverge asymptotically. To understand the large-n behavior we focus on the remaining terms. The heat-kernel factor J on the second line grows polynomially with n A x,k (the singularities in the denominator are cancelled by the Haar measure, as discussed above). The kinetic-energy term on the third line contains non-trivial dependence on n A x,k , and it is helpful to expand the kinetic-energy term for a general decomposition into real and imaginary components of the deformed Gaussian eigenbasis variables ψ For θ < π/2, the first term provides a Gaussian cut- for large m A x,k , which equivalently corresponds to large-n behavior. This provides absolute convergence for Wick rotated path integrals with θ < π/2 for all points on all contour segments.
In the desired θ → π/2 limit, only the last term survives. If lim m→±∞ sign(z A x,k ) = sign(m A x,k ), that is if the sign of the imaginary part of φ (c),A x,0k matches the sign of m A x,k asymptotically, then this surviving term provides an exponential cutoff for large m A x,k sufficient to overcome the polynomial growth in J and provide convergence. 8 This is true of all points on the middle contour segments in Fig. 5 identified by c A x,k = 0. However, this does not hold for the endpoints of the outer segments identified by c A x,k = 1.
The deformation employed here is therefore convergent everywhere except a set of measure zero for the action and deformation with θ = π/2, but this is not sufficient to guarantee absolute convergence. In fact, these singularities do cause the phase-quenched partition function Z pq θ,HK to diverge for θ = π/2. Fundamentally, the fixed endpoints of the deformation contour prevent convergence at these points for any choice of contour deformation. This problem can be resolved, and the path integral rendered absolutely convergent, by noting that in the full path integral these outer segments can be identified and cancelled using 2π shift symmetry. As discussed above, path integrals using the HK action are invariant under the transformation ( . This symmetry relates the outer contour endpoints β(m A x,k ) = α(m A x,k + A) + 2πA for integrals along the blue contours shown in Fig. 5 for ψ A x,0k integrals with m A x,k and m A x,k + A. The other endpoints of these contours at ±Aπ are also identical because of this shift symmetry, and therefore the blue contour segments in Fig. 5 describe oppositely oriented contour integrals of holomorphic functions between identified endpoints and therefore cancel. Strictly, this identification and cancellation should only be considered valid for an absolutely convergent integral. For all θ < π/2, the term cos θ(2πm A x,k ) 2 gives absolute convergence for the deformation under study, and for all θ < π/2 the outer contour segments in the neighborhood of the endpoints can be identified and exactly cancelled. This defines a limiting procedure for defining HK path integrals in which integrals are defined for a generic Wick rotation angle θ and the cancellation of the blue contours in Fig. 5 is performed before taking the θ → π/2 Minkowski limit.
Finally, any observable O(U ) without explicit ndependence leads to an O( U (c) (U, n)) whose only n-dependence can be described as dependence on sign(m A x,k ) and is therefore finite in the n A x,k → ∞ limit. This demonstrates that the exponential convergence provided by the kinetic term is not spoiled and that HK path integrals including observables are absolutely convergent.
It is noteworthy that in addition to providing an exponential convergence factor for the integer sums in the are amenable to importance sampling, and in particular are the same weights that appear in the Euclidean heat-kernel action. Phase fluctuations arise from terms with n A x,k = 0; however, the contributions from these terms are exponentially suppressed for small g 2 . This contour deformation can therefore be expected to significantly improve the sign problem associated with the real-time heat-kernel kinetic term for sufficiently weak couplings in addition to providing absolute convergence.

D. Convergent LGT Schwinger-Keldysh path integrals
For applications involving the Schwinger-Keldysh formalism, an action S SK must be constructed for the time integration contour The real-time transfer matrices for x 0 ∈ T + M and x 0 ∈ T − M correspond toT M,HK andT † M,HK respectively and are therefore unitary. The imaginary-time transfer matrix for x ∈ T E corresponds toT E,W and is therefore positive. With the contour ordering shown in Fig. 5, the Schwinger-Keldysh partition function associated with S SK,HK is equal to the Euclidean Wilson partition function (even before taking the continuum limit), It is noteworthy that either the Wilson or heat-kernel kinetic or potential terms could be used for x 0 ∈ T E , and in all cases the imaginary-time transfer matrix for x 0 ∈ T E would be positive. 9 However, the choice of Wilson kinetic term for x 0 ∈ T E , as well as the choice of Wilson potential term for x 0 ∈ T ± M , will be seen below to be important for establishing convergence of Schwinger-Keldysh path integrals using path integral contour deformations analogous to those in Sec. IV C.
In order to apply a generalization of this contour deformation to Schwinger-Keldysh path integrals, auxiliary variables n A x,k corresponding to the indices for all of the sums in the heat-kernel kinetic terms appearing Eq. (139) can be introduced and an n-dependent action can be defined as truncated heat-kernel potential is used in place of the Wilson potential, holomorphy of the potential as a function of φ A x,0k is no longer manifest because path integrals employing the truncated heat-kernel action do not possess a φ A x,0k → φ A x,0k +2π shift symmetry and are therefore sensitive to branch cuts in arg(e iφ A x,0k ) that could cause further complications. 9 It is also possible to construct a more general action with the same symmetries by introducing independent gauge couplings g 1 for the Wilson action for x 0 ∈ T E and g 2 and g 3 for the kinetic and potential terms, respectively, of the HK action for x 0 ∈ T ± M . Choosing a different trajectory in this three-dimensional coupling space besides the choice g 1 = g 2 = g 3 in Eq. (139) would modify the LGT spectrum in the Euclidean and Minkowski regions (which in general differ by lattice artifacts) and the approach to the continuum limit defined by lim g 1 ,g 2 ,g 3 →0 . The investigation of lattice artifacts and the existence of a continuum limit for the trajectory g 1 = g 2 = g 3 is deferred to future work.
Nonequilibrium LGT observables with Schwinger-Keldysh contour path integral representations can be defined using the HK action as where Z SK,HK = {n} DU e iS SK,HK (U,n) . A straightforward generalization of the contour deformation introduced in Sec. IV C can be used to provide an absolutely convergent representation of O SK,HK as described below.
For Schwinger-Keldysh applications the path integrals involving gauge fields on T E are convergent before applying contour deformations and do not have sign problems arising from the action. The gauge field variables on T E can therefore be held fixed, and only the gauge field variables on T ± M need to be actively deformed (induced transformations to plaquettes on T E near the boundaries with T ± M are discussed below). The fixed gauge field variables on T E can then be used as boundary field configurations to perform a change of variables from {U x,µ } to {P x,0k , U x,0 } for x 0 ∈ T ± M . The original gauge field variables {U x,µ } will continue to be used a path integral variables for x 0 ∈ T E . For x 0 ∈ T + M , this change of variables is given by the function U x,µ (P x,0k , U x,0 ) defined as where η = −i(β − a) denotes the imaginary part of x 0 for x 0 ∈ T ± M . An analogous change of variables valid for Further introducing the timelike plaquette eigenvalues e iφ A x,0k and eigenvalues V x,0k defined by Eq. (114) as well as the variables ψ A x,0k with A = 1, . . . , N − 1 defined by Eq. (121) for which SU (N ) heat kernel takes an uncorrelated Gaussian form, the variables {ψ A x,0k , V x,0k , U x,0 } can be used as a set of independent path integral variables on T ± M . The redundancy introduced by treating V x,0k as independent variables is irrelevant by Eq. (116). Deformed contours for ψ A x,0k can then be introduced that are equal to the contour shown in Fig. 5 for T + M and equal to its reflection about the horizontal axis for T − M . The central segments of these contours correspond to the co- where the ± signs in e ±iπ/4 correspond to x 0 ∈ T ± M . Defining φ A x,0k in terms of ψ A x,0k and m A x,0k using Eq. (128), this transformation turns the n A x,k dependent oscillatory factors in Eq. (141) into exponential damping factors that render the sums in Eq. (142) exponentially convergent as in Sec. IV C. Integrals over the additional contour segments required to complete the contours in x,k after regularizing the sums with transformations of kinetic prefactor exactly equivalent to the "Wick rotation" procedure described in Sec. IV C for x 0 ∈ T + M and with θ → −θ for x 0 ∈ T − M . Formally, Schwinger-Keldysh path integrals can then be defined by performing this cancellation with a generic Wick rotation angle and subsequently taking the Minkowski limit.
A noteworthy difference between Schwinger-Keldysh contour deformations as compared with the purely realtime case discussed in Sec. IV C is that the contour deformation associated with Eq. (145) induces a nontrivial transformation of the timelike plaquette P x,0k with x 0 = −i(β − a) ∈ T E . If a heat-kernel kinetic term were used for the Euclidean segments, this induced transformation could introduce a divergence to the Euclidean heat-kernel sum for this boundary plaquette; however, the choice of the Wilson kinetic term for x 0 ∈ T E ensures that the only infinite sums appearing are associated with x 0 ∈ T ± M and therefore that the transformations of fields on T E cannot spoil convergence.
An absolutely convergent representation for Schwinger-Keldysh path integrals with the HK action is therefore given by where the full Jacobian J(U, n) is given in analogy to Eq. (133) by where the contributions from each link are given by and LGT actions in (1 + 1)D Euclidean spacetime (left) and the HFK action in Minkowski spacetime (right) analogous to Fig. 7. The sign problem for the HFK action using the deformed contour studied in this work is much more severe than the corresponding sign problem for the HK action and contour deformation, and identical sized statistical ensembles were used for the HFK action with a volume of LT L = 4 and the HK action with a volume of LT L = 32.
action, the only two actions for which the contour deformations above result in well-defined Monte Carlo sampling schemes. In (1 + 1)D the remaining gauge links after gauge fixing are in one-to-one correspondence with the plaquettes P x ≡ P x,01 . The factorization of Wilson loop path integrals into products of one-dimensional integrals for each P x can be exploited to obtain a simple Monte Carlo simulation strategy in which plaquette variables P x and integer-valued auxiliary variables are drawn independently from the probability distribution Results for the expectation values of Minkowski Wilson loops as a function of area tL are shown in Fig. 7 for an ensemble of N cfg = 50, 000 configurations for a lattice with total area L T L = 32a 2 and gauge coupling e = 0.6. In (1 + 1)D with OBCs these results are only sensitive to the total area L T L and the total Wilson loop area tL, independent of the shape, due to factorization of the path integral. Results for Euclidean Wilson loop expectation values obtained using the heat-kernel action for an ensemble with the same number of configurations, lattice area, and gauge coupling are also shown for comparison in Fig. 7. Good agreement between numerical Monte Carlo and exact analytic results is obtained in both Euclidean and Minkowski cases. The inclusion of integer-valued auxiliary variables is not found to significantly increase autocorrelation times or introduce other undesirable numerical features in comparison to calculations with the Wilson action.
In this simple theory, similar statistical precision is achieved in the Euclidean and Minkowski cases. This can be understood by noting that e and therefore that the probability distribution is identical to p M,HK (P x , n x ) for the case of G = U (1) from which W A ( P (P, n)) can be explicitly constructed as  Fig. 7. This level of precision for Minkowski observables is remarkable given the severity of the sign problem before applying the contour deformation in Eq.  is identical to the exact string tension σ U (1) E,HK = e 2 /2 for the heat-kernel action with e = 0.6. Results for Monte Carlo calculations using this gauge coupling and the same ensemble size, N cfg = 50, 000, and lattice area as the heat-kernel calculations are shown in Fig. 8. Comparison to the Euclidean heat-kernel results in Fig. 7 shows that the two Euclidean actions achieve similar precision as well as good agreement with the corresponding (identical) exact results. Significantly lower precision is achieved by the Minkowski HFK action in comparison to HK action calculations above, which can be attributed to the presence of a reweighting factor e iRe[S U (1) M,HF K ( U ,r)] with severe phase fluctuations whose StN is observed to decrease exponentially with increasing lattice area L T L. The results shown in Fig. 8 use a lattice with area L T L = 4a 2 . A statistical ensemble many orders of magnitude larger would be required to extend these results to a lattice with area comparable to the L T L = 32a 2 areas explored in the HK case. Alternatively, a more sophisticated contour deformation than the simple shift in Eq. (109) might be able to significantly improve the sign/StN problem associated with e iRe[S U (1) M,HF K ( U ,r)] while maintaining the convergence properties. Studying the best choice of action and deformation for practical simulations is deferred to future work.

F. Monte Carlo simulations of two-dimensional SU(3) gauge theory
Exactly as in the U (1) case, an exploratory study of the feasibility of numerical Monte Carlo calculations of real-time non-Abelian gauge theory observables can be performed in (1 + 1)D with OBCs and validated against exact results from Sec. III E. The HK action is the only unitary SU (3) LGT action studied in this work for which a contour deformation is obtained that results in a welldefined distribution for Monte Carlo sampling. Factorization of path integrals into products of one-dimensional integrals is again exploited to obtain a Monte Carlo algorithm in which P x ∈ SU (3) is sampled along with integervalued auxiliary variables n A x satisfying the constraint A n A x = 0 from the probability distribution defined in Eq. (150). Expectation values are then approximated using Monte Carlo ensembles averages defined by Eq. (151). Wilson loops factorize as and the corresponding transformed observables are given by where P x (P, n) is defined by Eq. (129). Fig. 9 shows results from Monte Carlo calculations using the SU (3) HK action with an ensemble of N cfg = 50, 000 gauge field configurations, using a lattice with area L T L = 32a 2 , and gauge coupling g = 0.6. Results for the SU (3) Euclidean heat-kernel action for an ensemble with identical parameters are also shown in appearing in the denominator of Eq. (151) can be calculated with sub-percent-level precision for the lattice area and ensemble size used here. The precision obtained for Minkowski and Euclidean SU (3) Wilson loop results is similar, and the familiar exponential StN degradation with loop area is seen in both cases (note that this is distinct from, and much milder than, the expected sign problem extensive in total lattice volume associated with undeformed real-time simulation). In higher dimensions, these phase fluctuations may become more severe and additional phase fluctuations will arise from potential terms in the action. Further studies are needed to explore the feasibility of Monte Carlo calculations of realtime SU (N ) LGT in higher dimensions, and more sophisticated contour deformations or other approaches to improving the sign problem may be required to achieve precise results for four-dimensional real-time LGT observables.

Naive Wick rotation of the Euclidean Wilson
LGT action leads to a Minkowski LGT with non-unitary time evolution, as pointed out in Ref. [10]. It is demonstrated in this work that there is no way to recover a unitary time-evolution operator even in the continuum limit using the real-time Wilson LGT action. Further, it is demonstrated that for exactly solvable examples in (1+1)D the weak-coupling limit associated with the continuum limit of Euclidean LGT does not exist using the real-time Wilson LGT action. Real-time LGT calculations must therefore use an alternative Minkowski action.
One alternative action leading to unitary timeevolution for real-time LGT is provided by HFK in Ref. [10]. The character expansion defining this action is divergent and requires regularization. Here, a path integral contour deformation that renders the character expansion absolutely convergent is provided for the gauge group U (1), making it possible to use Monte Carlo methods for this action. Exploratory Monte Carlo calculations of U (1) Wilson loops in (1 + 1)D Minkowski spacetime were performed using the HFK action and found to be consistent with exact results obtained through analytic continuation of results for the Euclidean Wilson action. Finding a similar contour deformation is difficult for the SU (N ) HFK action because the character expansion involves sums over irreducible representations of SU (N ), and understanding the convergence properties of these sums is non-trivial. However, no fundamental obstacle was encountered that would prevent the construction of such a deformation, and it may be an interesting subject of future work to find a convergent representation of the SU (N ) HFK path integral because of its closeness to the commonly employed Wilson gauge action.
Another class of actions leading to unitary timeevolution in Minkowski LGT is obtained in this work through analytic continuation of the Euclidean heatkernel action. Unitarity of the real-time transfer matrix depends only on using the heat-kernel or HFK kinetic term in the action and does not depend on the form of the potential term (besides the assumption that the potential term is real). In particular, this permits the definition of a modified heat-kernel, or HK, action using the potential term from the Wilson action. The real-time HK action is defined by a divergent series required to ensure periodicity of angular variables, but by exploiting the Gaussian form of the heat kernel we obtain path integral contour deformations leading to absolutely convergent representations of HK path integrals for both U (1) and SU (N ) real-time LGT. Exploratory Monte Carlo calculations of U (1) and SU (3) Wilson loops in (1 + 1)D are performed and found to agree with exact results obtained through analytic continuation of results for the Euclidean heatkernel action.
The contour deformations of the HFK and HK path integrals are crafted to provide a complete cancellation of the sign problem associated with the kinetic energy term in the weak coupling limit, beyond giving a convergent representation of these path integrals amenable to Monte Carlo sampling. In the (1 + 1)D proof-of-principle calculations, this is apparent as an exponential reduction in the sign problem associated with the reweighting factors e iS M,HF K and e iS M,HK , respectively. For the HFK action, the sign problem is not completely eliminated and the deformed theory still suffers from a signal-to-noise problem that grows exponentially with the lattice size, making numerical results less precise than in analogous Euclidean LGT calculations. On the other hand, the HK deformation completely removes the sign/StN problem for the U (1) case and exponentially improves it for the SU (3) case, in both cases resulting in comparable precision for Minkowski and Euclidean Wilson loop results. These results should however be considered specific to (1 + 1)D; gauge theories in (1 + 1)D do not include a potential term and therefore a careful analysis of the kinetic energy term is sufficient to directly find good contour deformations. It can be expected that sign/StN problems that are approximately or entirely solved in (1 + 1)D will be significantly worse in (3+1)D. Exploration of the contour deformations needed to practically tame sign/StN problems in higher dimensions is left to future work. It is possible that more sophisticated contour deformations must be constructed or that numerical optimization of a variational ansatz for the contour deformation will be necessary to make progress.
This work focused exclusively on pure gauge theory, but the inclusion of fermionic matter fields is not expected to spoil unitarity [10]. Convergence of the HK action is not affected by the presence of fermionic matter since the fermion determinant does not explicitly depend on the heat-kernel sum index and therefore will only depend on sign(m A x,k ) after applying our contour deformation. Care should be taken with real-time fermion doublers and the continuum limits of lattice gauge theories including fermions, and a study of these subtleties is deferred to future work.
The absolutely convergent representation of HK path integrals defined here provides a suitable starting point in principle for Monte Carlo calculations of SU (3) LGT in (3 + 1)D Minkowski spacetime and calculations of Schwinger-Keldysh path integrals in QCD. If phase fluctuations arising from sign problems in (3 + 1)D can be tamed for some real-time LGT observables, significant strides could be made towards predicting real-time gauge theory observables relevant for phenomenology. has a corresponding time-evolution operator T M (ϕ, ϕ ) = e −iaV (ϕ)/2R M (ϕ, ϕ )e −iaV (ϕ )/2 , where the kinetic term is included in Defining the measure Dϕ = x dϕ x and δ(ϕ − ϕ ) = x δ(ϕ x − ϕ x ), it can be explicitly verified thatR M is unitary after multiplying by a normalization factor of The unitarity ofT M follows immediately from Eq. (A5).  (6) includes an infinite number of singularities that accumulate as g 2 → 0 and is replaced by a gray background for g 2 < 0.1. Numerical precision limits the approach to g 2 → 0 for other choices of N , and the results are replaced by a gray background for g 2 < 0.02. The Wilson action can be seen to have a non-unitary coefficient ratio almost everywhere in g 2 for all N . In the SU (6) case, the limit g 2 → 0 is ill-defined due to the repeated singularities, while for other SU (N ) there is an apparently well-defined limit which is non-unitary. The limit values are explained by a stationary phase expansion in the text. The HFK action has a unitary ratio for all choices of g 2 and N .
to give a universal constant. The Haar measure over eigenvalues must also be included, giving in total c M,W,SU (N ) ) .

(B6)
The constant of proportionality is independent of the representation r, and will cancel from ratios.
Finding the stationary points of the action requires addressing the unit determinant constraint for the N eigenvalues. One approach to including this constraint is through the use of a Lagrange multiplier term, α(2πm − A φ A ) with arbitrary m ∈ Z. The stationary points of the action with respect to the N − 1 free directions are then given by the solutions to This is solved by the group elements for which the imaginary component of every eigenvalue is identical; for example, this includes the identity element, which is also the global minimum of the Wilson action with Euclidean signature. For specific choices of N , other solutions are also possible. The Haar measure in square brackets in Eq. (B6) suppresses eigenvalue degeneracy, so the stationary points with lowest eigenvalue degeneracy will dominate the integral. Table II gives these lowest-degeneracy stationary points for N ∈ {2, . . . , 9}. For a generic N , we can label the set of k dominant N Stationary points (lowest degeneracy) χ f 2 (1, 1) / (−1, −1) 2 / −2 3 (1, −1, −1) −1 4 (e iθ , e iθ , e i(π−θ) , e i(π−θ) ) 4i sin(θ) 5 (1, 1, 1, −1, −1) 1 6 (1, 1, 1, 1, −1, −1) / (1, 1, −1, −1, −1, −1) 2 / −2 7 (1, 1, 1, −1, −1, −1, −1) −1 8 (e iθ , e iθ , e iθ , e iθ , e i(π−θ) , e i(π−θ) , e i(π−θ) , e i(π−θ) ) 8i sin(θ) 9 (1, 1, 1, 1, 1, −1, −1, −1, −1) 1 TABLE II. The group elements that are stationary points of the Wilson action for the group SU (N ) with N = 2, . . . , 9; only the elements with lowest eigenvalue degeneracy are listed, as these dominate the character coefficient integrals due to the suppression of degenerate eigenvalues by the Haar measure. Group elements are reported in terms of their spectra (i.e. sets of eigenvalues) because the Wilson action is invariant under transformations of the eigenvectors alone. Stationary points always include the center elements of the group, but additional elements may also solve the constraints. The center elements have maximum eigenvalue degeneracy, so when other solutions exist they will dominate the integrals. The fundamental representation character is given for all cases with a single dominant stationary phase.
stationary points P 1 , . . . , P k ∈ SU (N ). Expanding the evaluation of character expansion coefficients in the limit g 2 → 0 gives c M,W,SU (N ) (B8) where the eigenvalues of P are denoted by e iφ A . The Haar measure factor in brackets in Eq. (B8) vanishes exactly at the stationary points, and to be precise one should expand to higher order in the stationary phase approximation and then determine the limiting values of character expansion coefficient ratios using these higherorder results. However, these higher-order terms are determined by the structure of the Haar measure and the action about the stationary points, which is independent of the representation r, and coefficients of these terms will cancel in ratios of character expansion coefficients. Ratios of the approximate expressions in Eq. (B8) will therefore give correct g 2 → 0 results for ratios of character coefficients under the assumption that the stationary phase expansion has a non-trivial radius of convergence about g 2 = 0.
For some choices of N , it is apparent from Fig. 3 and Fig. 10 that the g 2 → 0 limits of character expansion coefficient ratios do not exist. These values align with the cases of N ≡ 2 (mod 4) where multiple stationary points contribute to the value of the integral as shown in Table II, and one expects the divergences to correspond to cancellations between these contributions to c M,W,SU (N ) 0 (g 2 ) for particular values of g 2 . For other N , the character expansion coefficients ratios shown in Fig. 3 and Fig. 10 appear to approach a limiting value as g 2 → 0, suggesting that the stationary phase expansion can be reliably used to understand the g 2 → 0 limit in these cases. For N ≡ 0 (mod 4) a class of group elements parameterized by an angle θ ∈ [−π/2, π/2] contributes to the value of the integral, as detailed in Table II these group elements are weighted by the fundamental character, which is proportional to sin(θ). Integrating over θ ∈ [−π/2, π/2] gives cancelling contributions from this set, and indeed the ratios can be seen to vanish as g 2 → 0 for these cases in Fig. 10.
In the remaining case of odd N , character expansion coefficient ratios are dominated by a single stationary point, and labelling the corresponding group element P 1 , the would-be continuum limit of these ratios is given by The fundamental character χ f (P 1 ) * in these cases corresponds to ±1, and the ratio therefore approaches ±1/N with the sign depending on the particular choice of N . This behavior is observed for SU (3) in Fig. 3 as well as SU (5), SU (7), and SU (9) where P x,µν is the plaquette obtained from U x,µ . In temporal gauge P x,0k = U x,k U † x+a0,k and therefore the kinetic-energy evolution operatorR M defined bŷ x,k r d r c r (g 2 )χ r (P x,0k ) .
(C3) If S M (U ) is a unitary real-time LGT action, that is ifR M and therefore the corresponding real-time transfer matrix T M are unitary, then it follows from the character decomposition in Eq. (34) that |c r (g 2 )| = 1. The fact that each character is normalized by dU |χ r (U )| 2 = 1 implies that χ r (P x,0k ) is non-vanishing for a set of non-zero measure in path integrals involving e iS M (U ) , and it therefore follows from |c r (g 2 )| = 1 that the r-th term in the sum in Eq. (C3) using any enumeration of the representations of G does not vanish as r → ∞. It immediately follows that the character expansion for e iS M (U ) in Eq. (C3) diverges for a set of gauge fields with non-zero measure in real-time LGT path integrals.
Appendix D: The truncated real-time heat kernel action In the classical limit, the real-time LGT path integral is dominated by the stationary phase solutions. Though the real-time heat kernel path integral introduced in Sec. III C involves weights given by non-convergent sums over integers n A , this classical limit is dominated by n A = 0. This section explores the path integral defined by a truncation of these sums to n A = 0, and demonstrates that ultimately it results in non-unitarity in the continuum limit in the analytically tractable case of U (1) LGT in (1 + 1)D.
Explicitly, the truncated real-time heat-kernel action is defined by for G = SU (N ) with J ({φ x,0k }, {0}) replaced by 1 for G = U (1). In the case of G = U (1), the normalizing constant is explicitly N = 1/ √ 2πie 2 . The truncated realtime heat-kernel action corresponds to the n = 0 term in the real-time heat-kernel action defined by Eqs. (74), (79), (80). The properties of the truncated real-time heat kernel action can be analyzed most simply for G = U (1), and we specialize to this case below. The kinetic-energy evolution operator associated with this action for G = U (1) is given bŷ which establishes that unitarity is recovered for a fixed lattice volume in the e 2 → 0 limit. However, the scaling limit of interest to continuum QFT is defined by taking e 2 → 0 with some dimensionful observable held fixed. It is discussed in Ref. [13] and Sec. III E that, in (1 + 1)D where LGT is analytically solvable, correlation lengths diverge ∝ 1/e in the e 2 → 0 limit. The appropriate scaling limit in (1 + 1)D is Lt . (D7) The factor in brackets approaches unity exponentially faster than in the Minkowski case, and the truncated Euclidean heat-kernel action leads to O(e −1/e 2 ) lattice artifacts. Further, the limit of W A M,n=0,U (1) at fixed Lτ σ U (1) HK exists and equals the Euclidean heatkernel result e −Lτ σ U (1) HK in the same limit. This example demonstrates that "lattice artifacts" associated with non-unitary Minkowski actions can be difficult to remove, even if unitarity is restored in the would-be continuum limit. It is unlikely that the truncated heat-kernel action provides a useful starting point for studying the continuum limits of higher-dimensional lattice gauge theories.