Real-time Spin Systems from Lattice Field Theory

We construct a lattice field theory method for computing the real-time dynamics of spin systems in a thermal bath. This is done by building on previous work of Takano with Schwinger-Keldysh and functional differentiation techniques. We derive a Schwinger-Keldysh path integral for generic spin Hamiltonians, then demonstrate the method on a simple system. Our path integral has a sign problem, which generally requires exponential run time in the system size, but requires only linear storage. The latter may place this method at an advantage over exact diagonalization, which is exponential in both. Our path integral is amenable to contour deformations, a technique for reducing sign problems.


I. INTRODUCTION
In this paper we develop a lattice field theory method for computing the real-time dynamics of spin systems in a thermal bath.By "spin-system", we mean a collection of SU (2) generators, s(x), assigned to sites on a lattice, interacting through a Hamiltonian H(s).This lattice could correspond to an actual lattice in space, but may be also be an arbitrary graph.No assumptions are made about the Hamiltonian except that it is a function of a finite number of spins.
Our purpose is twofold.First, we wish to provide an alternative to exact diagonalization for reproducing the realtime dynamics of digital quantum computers.Today, the only way to check a digital quantum computer's simulation of real-time physics for a general hamiltonian is to diagonalize it or perform a Trotterized task of equal complexity [1][2][3][4].The cost of such computations is exponential in the number of qubits in both storage and run-time.The method we develop is exponential in runtime, but is linear in storage.Second, we wish to provide an alternative to tensor networks for computing the real-time properties of spin systems, especially in more than one spatial dimension, where tensor networks are generally less efficient [5,6].
Our main result is a path integral discretization for spin systems which provably reproduces the exact real-time dynamics for any Hamiltonian of a finite number of spins.The path integral is discretized in both space and time, and the desired dynamics are recovered in the "time continuum limit".The field variables of this path integral are points on the two-sphere which are in one-to-one correspondance with points on the Bloch sphere.What is shown is that correlation functions of points on the Bloch sphere in the lattice theory converge exactly to spin correlation functions in the quantum theory.What is perhaps unconventional about our formulation is the appearance of a "Schwinger-Keldysh action" [7,8], which encodes both real-time and thermal information.
The most important property of our discretization is that has the correct continuum limit.The same cannot be said of all discretizations on the market.Demonstrated analytically in [9] then corroborated numerically in [10], even textbook path integrals for spin systems, such as those found in [11,12], do not in general have time continuum limits.In spite of being explicitly mentioned [12], this fact appears little appreciated today.
Our discretization builds upon a foundation layed by Takano in the 1980's [13], where he developed a Monte Carlo method for computing thermal properties of spin systems.The basic building block of Takano's path integral are "spin coherent states" 1 .Developed in [14][15][16][17], spin coherent states are to spin operators what bosonic coherent states are to bosonic operators: they are an overcomplete set of states providing a resolution of the identity that can be used to build a path integral.The resolution of the identity involves an integral over the Bloch sphere, and this is ultimately why the path integral has field variables on the sphere.
What we add to Takano's work is real-time dynamics and simplification.In his original work the possibility of real-time dynamics was not considered.Furthermore, the original formulation relies on a certain "product formula" for inserting observables into the thermal trace [18].Implementing the product formula quickly becomes prohibitively difficult as both both the hamiltonian and observables increase in complexity.We skirt the product formula altogether by using functional differentiation, a common technique in lattice field theory for observable insertions.This allows for the analysis of more general Hamiltonians and observables than would have otherwise been practical.In the following sections we derive the path integral, demonstrate that it works, then offer a discussion.

II. DERIVATION
Consider a spatial lattice of spins s(x) interacting through a hamiltonian H and satisfying the canonical commutation relations [s i (x), s j (y)] = iϵ ijk s k (x)δ xy .In this section we develop a path integral representation of correlation functions where x, x ′ are lattice sites, t, t ′ are real times, and β = T −1 is the inverse temperature of the system.This path integal representation is obtained in two steps.First, a trotterized, Hilbert-space object is formed which reproduces the correlation functions in the time continuum limit.Second, this Hilbert-space object is approximated by a path integral with controlled errors.While our presentation is for two-point functions it will be clear the procedure generalizes to n-point functions.Throughout we assume a lattice of finite extent but arbitrary spatial dimension.Consider a time extent t max larger than t & t ′ and the following object: where where the sums run over lattice sites, 1 is the identity operator, N is an integer number of "timeslices", ∆τ = β/N , ∆t = t max /N , and j + , j − , j E are classical sources.For brevity, we collect all sources into a symbol j and define Z(j) ≡ Z(j + , j − , j E ).The short-time propagators P + , P − , P E build up a discretized "Schwinger-Keldysh contour" in the complex time plane, shown in Fig. (1).This contour defines a "Schwinger-Keldysh action" [7,8], which we will elaborate on later.We presently show that derivates of Z with respect to sources produce the real-time correlation functions Eq. (1) when N → ∞ at fixed β, t max .We borrow language from lattice field theory and call this the "time-continuum limit".
FIG. 1. Schematic of the discretized Schwinger-Keldysh contour.The path followed is indicated by the arrows: beginning at zero the contour sweeps right to positive tmax, then back to zero, then down an amount β.The vertical splitting between the real-time segments is only drawn so that the forward and backward contours do not overlap; an iϵ prescription is neither needed nor employed.
First, since the Hilbert space is finite dimensional, the Hamiltonian is bounded and N can always be taken large enough that where the errors are small as desired.Therefore Z(0) converges linearly to tr e −βH in the time continuum limit.Next, consider two spacetime points (x, t) & (x ′ , t ′ ) with t < t ′ < t max , and define integers ( t, t′ , tmax ) = (t, t ′ , t max )/∆t.Then where P a ≡ P a (0) for a = +, −, E. Since P t−1 + = e iHt + O(∆t), and similarly for the other products of short time propagators, one obtains Taking a time-continuum limit while holding the physical time separation between derivatives fixed, a real-time correlator is obtained.While this construction produced a time-anti-ordered correlator, any time ordering can be obtained.One finds: We will now relate Z(j) to a lattice path integral using Takano's framework [13].The key step is to express the short-time propagators in terms of spin coherent states, which is accomplished using a "diagonalization theorem" [14,15,19].We first set notation: to every spatial lattice site x we associate a sphere whose points we denote Ω x = sinθ x cosϕ x , sinθ x sinϕ x , cosθ x .A spin coherent state is then defined as where |s, s⟩ is the highest weight state in the spin-s representation of SU (2).Then the diagonalization theorem states that any operator on spin space may be written as an integral over spin coherent states: where the function f O depends on the choice of operator and The theorem allows to write x and similarly for P − , P E2 .Here h is related to the Hamiltonian through H = d ¯Ω h(Ω) |Ω⟩⟨Ω|.We now introduce the operators3 which are equal to P + (j + t ), P − (j − t ), P E (j E t ) up to errors of O(∆t 2 , ∆τ 2 ).Therefore, is equal to Z(j) up to errors of O(∆t, ∆τ ).Computing the trace, one finds: where: x x x We call S SK (Ω, j) the Schwinger-Keldysh action.To clarify the notation, the Ω + t appearing within Z(j) arise from the product of P + operators, while the Ω − t and Ω E t arise from the products of P − and P E .We collect these three into a composite variable Ω which is a "path ordered" concatenation of Ω + , Ω − , Ω E , in that order.Referring to the Schwinger-Keldysh contour in Fig. (1), Ω + t , Ω − t , Ω E t are assigned to the forward, backward, and Euclidean segments, respectively.Since Z(j) and Z(j) differ at O(∆t, ∆τ ), so do their derivatives.One therefore has the following relations, up to errors of O(∆t, ∆τ ), which vanish in the time continuum limit: Here ⟨O⟩ SK ≡ Z(0) −1 DΩ e S SK (Ω,0) O(Ω) and Ω a t x i is the i th component of Ω a on lattice site x and timeslice t.Eq. ( 14) and Eq. ( 15) are our main results.
We conclude with some details about the Schwinger-Keldysh action.The inner product of coherent states appearing in Eq. ( 14) is the analog of the "Berry phase" in our formulation [12].Each overlap has the explicit expression The only dependence of the path integral on the system in question is through the structure of the lattice and the function h(Ω).As long as the hamiltonian does not involve products of spins on the same site, h is obtained from the Hamiltonian through the simple relation: where on the right hand side we have distinguished the operator ŝ from the spin representation s.If the Hamiltonian does contain products of spins on the same site, Eq. ( 17) no longer holds, and it becomes necessary to use the product formula of [13,18] to derive h.We emphasize however that the only complication is in the relation between H and h.Correlation functions retain the form Eq. ( 15).

III. DEMONSTRATION
We now engage in a numerical demonstration.We consider the Hamiltonian for two-qubits in one spatial dimension with periodic boundary conditions.We will obtain the real-time correlation functions ⟨s i (x, t)s i ′ (x ′ , t ′ )⟩ by computing the Schwinger-Keldysh path integral Eq. ( 13) and taking a time-continuumlimit.The h function corresponding to this hamiltonian is Usually, lattice path integrals can only be estimated stochastically using Monte Carlo methods.However, by considering a small system, we are able to elide the Monte Carlo and compute the lattice path integral exactly.Without the burden of stochastic noise we demonstrate that the time continuum limit of the Schwinger-Keldysh path integral Eq. ( 13) both exists and converges to the correct result.Since the method works for multiple qubits, albeit two, we see no reason it should fail for any finite number of them.
To circumvent the Monte Carlo, we compute matrix representations of the P + , P − , P E appearing in Eq. (12).For this two-qubit theory, these are 4 × 4 matrices which, at fixed ∆t and ∆τ , can be computed numerically by evaluating the integrals over spheres in Eq. (11).To form correlation functions, we also numerically compute the matrices In the following example we take β = 3.0j −1 and t max = 10j −1 .Due to spacetime and internal symmetries, there are only four independent two-point correlation functions in this theory: ⟨s i (x, t)s i (x, 0)⟩, ⟨s i (x, t)s i (x + 1, 0)⟩ for i = x, y.We show results only for i = x correlation functions; the qualitative conclusions drawn from i = y are the same.
Our numerical results are found in Fig. (2).In the left panel we plot the exact correlation functions ⟨s x (x, t)s x (x ′ , 0)⟩ in dotted lines, and the unordered correlation functions (s + 1) 2 ⟨Ω + t x i Ω − tmax x ′ i ′ ⟩ SK in transparent lines.The lattice theory is taken to have N = 5000 timeslices.The visual agreement between the dashed and transparent curves is Error is defined as the difference between the lattice and exact theories, and each point is obtained from a lattice calculation with N timesteps.The lines are linear fits to the finest three lattice spacings.evidence the lattice method works.This evidence is corroborated in the right panel of Fig. (2), which illustrates time continuum limits.The vertical axis is the difference between the lattice and exact results at time separation t = 5.0j −1 .The points are the results of lattice calculations at different number of timeslices N , and the lines are linear fits to the finest three lattice spacings.Linear convergence to zero as N −1 → 0 demonstrates both the existence and correctness of a continuum limit of the lattice theory.Clearly the fits in the right panel of Fig.

Fit window
(2) converge very close to zero.However, our fits to not converge to exactly zero.This is simply because at any non-zero lattice spacing there are always some O(∆t 2 , ∆τ 2 ) errors contaminating the fit.We have therefore performed fits at extremely fine lattice spacing.As tabulated in Table (I), we find a strictly shrinking error as we approach the continuum.We find this to be evidence beyond reasonable doubt that the method works.While we have presented detailed results only for the ⟨Ω + Ω − ⟩ SK correlator, which is able to access twice the time separation as the ⟨Ω + Ω + ⟩ SK and ⟨Ω − Ω − ⟩ SK correlators, we have computed all three and find they all converge to the correct continuum, limit.We have also computed the ⟨s y s y ⟩ correlation functions and obtain correct continuum limits in all cases.

IV. DISCUSSION
In this paper we developed a lattice field theory method for computing the real-time dynamics of spin systems.This method extends previous work by Takano with Schwinger-Keldysh and functional differentiation techniques.Our main results are Eq.( 14) and Eq. ( 15), which respectively give the Schwinger-Keldysh action for generic hamiltonians and formulae for spin correlation functions.We then demonstrated the method on a simple two-qubit system.
Our method can be easily extended in several ways.First, convergence to the continuum can be accelerated with higher-order approximations to the short-time propagators.The cost of this acceleration is the need to write products like H 2 , H 3 , ... in terms of coherent states which requires the product formula.Such expressions are cumbersome though not impossible to write.Second, three-point and higher-order spin correlation function are simple to obtain by taking more functional derivatives.In fact, n-point functions of arbitrary observables O are simple to obtain too.From the beginning one simply couples sources to O rather than to spins.The same Schwinger-Keldysh path integral results, but rather than Ω correlators one requires f O correlators (f O is defined through Eq. ( 8)).Combining this fact with the Schwinger-Keldysh path integral's ability to compute any time ordering, "out of time ordered correlators" of arbitrary observables can be computed with our method.This may prove to be useful in studies of quantum chaos [20,21].Third, though we have chosen a particular approach to the continuum (i.e.we have taken the number of timeslices to be idential on all three legs of the Schwinger-Keldysh contour), infinity other approaches exist and may be more efficient.Indeed, in highly asymmetric cases such as high-temperate & long real-time, or low-temperature & short real-time, it may be profitable to choose a different number of timeslices on each of the legs.Finally, we emphasize that no "+iϵ" prescription ever appears in our method.It is not needed.This is important: it is unwise to add unnecessary parameters to lattice simulations requiring extrapolation.There are already enough to do.
A Monte Carlo evaluation of the path integral will be unavoidable for larger systems.While the storage cost for Monte Carlo evaluations of lattice path integrals is linear in the number of space time points, our path integral has a sign problem.This generally requires exponential runtime to resolve [22].However, our path integral has a holomorphic path integrand for any Hamiltonian.This means it is anemable to "path integral contour deformations", a general method for taming sign problems where the integration manifold is deformed into complex field space to a surface with reduced phase oscillations [23].That our path integral is over many spheres makes it particularly straightforward to apply contour deformations, and existing machine-learned manifolds developed for other systems may prove useful [24].While the real-time dynamics of spins is a challenging problem, it is encouraging that other real-time problems have been solved with contour deformations [25,26].

FIG. 2 .
FIG. 2.Left: Real-time correlation functions.The transparent curves are unordered correlators obtained from the disctetized theory at fixed N = 5000, while the dashed curves are the exact result obtained by matrix exponentiation.Right: Time continuum limits at fixed t = 5.0j −1 .Error is defined as the difference between the lattice and exact theories, and each point is obtained from a lattice calculation with N timesteps.The lines are linear fits to the finest three lattice spacings.

TABLE I .
Table of continuum extrapolation errors.Printed within each row is the error in a linear continnum extrapolation of three lattice calculations with the number of timeslices indicated by the fit window.The column labels indicate the observable whose error is printed.