Imaginary time plays a prominent role in multiple branches of physics, including cosmology, statistical mechanics, and quantum field theory. The seemingly simple replacement of real time, t, with its imaginary counterpart, τ = -it, leads to fundamental connections between quantum theory and statistical mechanics [1]. Such connections enable the efficient simulation of many quantum systems using quantum Monte Carlo techniques [2][3][4]. However, for many physically interesting models, these methods suffer from the prohibitive "negative sign problem" [5,6], which requires an exponential amount of computational resources to obtain reasonable accuracy for quantum many-body systems. Many outstanding problems in condensed matter, such as those involving high-temperature superconductors or topologically ordered phases, require an understanding of complex interacting models which are unsolved with present techniques.
One class of Monte Carlo methods that can avoid the sign problem are so-called variational Monte Carlo (VMC) methods [7][8][9][10][11]. In VMC, one assumes a sufficiently general trial state that depends on adjustable parameters. These parameters are then chosen to minimize the energy with respect to the given Hamiltonian. Finding an effective trial state such as Jastrow [11], matrix product states [12][13][14], or neural network states [15][16][17][18], can result in efficient simulation of interacting quantum systems. The key to the success of these techniques is a well-chosen ansatz that reflects the properties of the target phase and the existence of a viable optimization scheme [19,20].
The recent advent of quantum computers and simulators has motivated the development of new variational approaches [21]. Such variational quantum algorithms involve applying a sequence of unitary operators, parametrized by several variables onto a easy-to-prepare initial state. The variables are chosen to optimize a given cost function involving the resulting wave function. For example, in the quantum approximate optimization algorithm (QAOA) [22][23][24][25][26] the cost function is a classical Hamiltonian encoding a combinatorial optimization problem, and the variational wave function is prepared by alternating between evolving with the Hamiltonian and a transverse field. The evolution times constitute variational parameters that are optimized to minimize the Hamiltonian cost function. This variational approach has been generalized for preparing both strongly correlated and highly entangled states on near-term quantum devices [27][28][29].
Motivated by the success of such variational approaches, in this paper we propose a variational ansatz for ground states of quantum many-body systems which involves sequentially evolving with different Hamiltonians in imaginary time. In contrast to real-time evolution with local Hamiltonians, which is limited by Lieb-Robinson bounds on the growth of correlation functions, imaginary-time evolution does not have this constraint and can exhibit remarkable efficiency in traversing Hilbert space.
As proof of concept, we demonstrate the efficiency of our ansatz in representing the ground state of the transverse-field Ising model at criticality. Whereas standard projector methods require imaginary-time scaling with system size L to reach the critical ground state, we show numerically that our ansatz requires time scaling logarithmically with L. Furthermore, we analyze how entanglement grows after each imaginary-time operation in our ansatz, and we find an exponential growth that is a unique feature of imaginary-time dynamics. We conclude by demonstrating that the ansatz continues to perform well in the presence of integrability-breaking perturbations, and we mention generalizations of our ansatz to other models, including those with sign problems. We envision the main purpose of this ansatz to be an efficient trial wave function for quantum many-body physics on (classical) computers; however, it is possible that one can also implement such imaginary-time evolution natively on a quantum computer [30,31].
Many Hamiltonians are naturally a linear combination of two components H = H A + gH B , where H A,B are individually tractable to analyze. Examples include transverse-field Ising, Hubbard, and the J 1 -J 2 model. Motivated by the QAOA procedure, we consider the following variational imaginary-time ansatz (VITA) for the ground state of H:
where P is the number of pairs of variational parameters α = (α 1 , . . . , α P ), β = (β 1 , . . . , β P ), |ψ 0 is an initial state, and N is a normalization factor. We further define the total imaginary time τ = 1 2 P p=1 (α p + β p ). A circuit representation is shown in Fig. 1. The bang-bang, or square-pulse, style of the ansatz is optimal for quantum control in real-time QAOA as per Pontryagin's principle [32].
While VITA is applicable to any Hamiltonian, in specific cases there is explicit physical motivation for considering such an ansatz. For example, for the fermionic Hubbard model, the P = 1 ansatz reduces to Otsuka's generalization of the Gutzwiller variational wave function [7,9,10], which seeks to balance single occupancy per site with itinerancy. The P 3 case has been considered in Ref. [33] for the twodimensional Hubbard model, but a systematic analysis of how its performance scales with system size and P was not carried out. A related variational approach for the Hubbard model has also been considered in Ref. [34].
The standard projector method for attaining the ground state of H involves evaluating e -τ H |ψ 0 for τ 1/ where is the many-body spectral gap. This can be decomposed via Trotterization into a sequence of the VITA form, with parameters α p = β p = τ /2P for large P. This guarantees that VITA can exactly represent the ground state in the P → ∞ limit. However, the projector method is especially expensive for critical systems where ∼ 1/L, and hence τ scales polynomially with L. One can consider Eq. ( 1) as a nonuniform Trotterization with large (and variable) time steps. We will show that remarkably high fidelities can be attained even with τ that is exponentially smaller compared to the aforementioned estimate from the standard projector method.
We first present some general considerations of why such an ansatz may be efficient. It is useful to first compare with
1. Circuit representation of the trial state |ψ P (α, β) for P = 1. Each box denotes imaginary-time evolution with the enclosed Hamiltonian. the real-time analog, which are QAOA-type circuits involving alternating real-time evolution between two Hamiltonians. The Lieb-Robinson bound dictates that real-time evolution with local Hamiltonians can generate correlations only within a light cone, and thus there are lower bounds on the time it takes to prepare highly correlated states starting from unentangled product states. For example, in one dimension, the total time to prepare the Greenberger-Horne-Zeilinger (GHZ; "cat") state 1
scales at least linearly with system size L [28,35]. In contrast, by evolving with the GHZ parent Hamiltonian H ZZ = -L i=1 Z i Z i+1 in imaginary time (Z is the Pauli-Z matrix), the GHZ state can be prepared with imaginary-time scaling as log L [36].
We test VITA on the transverse-field Ising model (TFIM)
with periodic boundary conditions on a system with N spins. Z, X are the Pauli matrices, and h is the transverse-field strength. Our ansatz in this case starts from the paramagnetic ground state of H X , |+ , and alternates between
The Ising chain can be mapped to free fermions via the Jordan-Wigner transformation [37] which allows for the efficient evaluation of Eq. ( 1). We can thus optimize our ansatz for very large system sizes and several pulses. This allows us to properly characterize how efficient the VITA ansatz is without introducing sampling error. For a fixed P, optimization involves finding the minima of the energy cost function E P (α, β) = ψ P (α, β)|H |ψ P (α, β) .
We first focus on approximating the critical ground state of H at h = 1. Figure 2(a) shows the relative error in energy rel = |(E P (α, β) -E exact )/E exact | where E exact is the exact ground state energy at the critical point, for various P. Evidently, increasing P dramatically improves the accuracy in the energy, even for large system sizes. Optimized parameters for various P are provided in the Supplemental Material [38]. Since the exact ground state for the TFIM is known, we also compare the fidelity, f , of the optimized trial state with the target state. The error in fidelity, 1
is shown in Fig. 2(b) for various P, L. The efficiency is quite remarkable; for example, for L = 64, P = 2 is already sufficient to approximate the critical state to within around 10 -4 in relative energy and 10 -2 in fidelity. Recall that the error in fidelity provides an upper bound for the error in any observable [38].
While there is no Lieb-Robinson bound limiting the rate for generating long-range correlations in our ansatz, entanglement considerations provide lower bounds on the iterations P required to prepare the critical state. Ignoring normalization, imaginary-time evolution with a local Hamiltonian can be represented by a (nonunitary) quantum circuit; each iteration of our ansatz corresponds to three layers shown in Fig. 1. After P pulses, the bipartite entanglement entropy between the left and right halves of the system, hereafter abbreviated EE, can be attained by bisecting the circuit through P bonds. Hence, EE after P iterations is at most P log D, where D is the Schmidt rank (number of singular values) upon decomposing a single two-qubit imaginary-time operator.
In order to generate the EE of the critical state, which scales as S ∝ log L, we need the number of pulses scaling at least as P ∝ log L. We observe that Fig. 2(b) is consistent with this scaling form. For example, for a target fidelity error of 10 -10 , each additional P in the ansatz can represent a system approximately twice as large.
The entanglement dynamics in imaginary-time evolution can be considerably different from its real-time counterpart. For real-time evolution, EE across a bipartition can increase only by acting with an operator supported on both sides of the partition. If the circuit in Fig. 1 were unitary, any increase in EE from one layer to the next would be bounded by a constant depending on the two-qubit unitary but not on the state being acted on (see the "small incremental entangling theorem" [39]). In contrast, even imaginary-time operators acting on one side of the bipartition can generate entanglement across the cut. As a very simple example involving two spins, the action of e -βZ 1 on η + |11 + η -|00 can increase EE as long as |η + | > |η -| > 0. This illustrates that the more entangled the initial state, the more imaginary-time evolution can change the entanglement. The change is not simply bounded by a state-independent constant. This allows in principle an exponential growth of EE, as long as the total EE after P steps lies below P log D.
Our ansatz exhibits such dynamics. We take the P = 5 ansatz and analyze the EE of the states at intermediate steps, p, of our protocol using the technique of [40,41]. We find that the EE increases exponentially with imaginary time [Fig. 3(a)]. Moreover, for every intermediate state, we plot the mutual information (S A + S B -S AB ) between two spins A, B as a function of their separation [Fig. 3(b)]. The power-law decay for every step is in stark contrast to any local realtime evolution and illustrates the ability of imaginary-time evolution to generate long-range correlations [42]. We find that under imaginary-time evolution, entanglement starts to grow immediately following a local quench, in contrast to real-time evolution [43] where it takes a time proportional to ℓ before growing, ℓ being the distance between the location of the local quench and the entanglement cut.
To compare directly with the projector method, we also investigate the total imaginary time τ ≡ 1 2 P p=1 (α p + β p ) required to achieve a target fidelity. Motivated by the P ∝ log L scaling for achieving a target fidelity, we propose a scaling form of 1f = G[τ (log L) -ν ] for some exponent ν. In Fig. 4, we perform a scaling collapse for L ∈ [4, 6, . . . , 1024], and P ∈ [1, . . . , 7] and find the optimal exponent ν = 2.3 ± 0.1. This logarithmic scaling is an exponential speedup compared to the linear scaling of the projector method, 1f = F (τ L -1 ) [44].
While the TFIM model admits a dual representation as free fermions, for a general model sampling methods are crucial for estimating the energy cost function. As a proof of concept, we also use Monte Carlo sampling for stochastically optimizing VITA.
The quantum-classical correspondence maps quantum observables to dual classical observables of a classical anisotropic Ising model on an L × (2P + 1) lattice. We denote the classical spin configurations by {s} and the spatial and imaginary time by (i, p), respectively 0 0.5 1 1.5 2 h 10 -6 10 -5 10 -4 10 -3 10 -2 rel (a) L = 64 P = 1 P = 2 2.5 3 3.5 4 h (b) N = 10 × 10 FIG. 5. Relative error in energy, rel , for VMC using our ansatz on the TFIM: (a) one-dimensional model with L = 64 spins. Solid lines denote the results from the free fermion approach. (b) Twodimensional model on a 10 × 10 square lattice. Energies are compared with those from zero-temperature stochastic series expansion [47].
The expectation value of a quantum observable O is
where Õ are dual classical observables, and p α,β (s) is the Boltzmann weight corresponding to Ising model with couplings J x (p) = α p , J τ (p) = 1 2 ln coth β p between nearest neighbors in space and imaginary time, respectively. The couplings vary with imaginary time, p, but are uniform in space. In this way, the energy of the trial wave function can be sampled efficiently with Monte Carlo.
We use this scheme with P = 1, 2 to target the ground states for various values of the transverse field h in both the (integrable) one-dimensional and the (nonintegrable) twodimensional TFIM. We use stochastic natural gradient descent (stochastic reconfiguration) to optimize the parameters [20]. We find rapid convergence for P = 1, while higher P becomes more difficult, especially for with a noisy objective function.
The relative error in energy achieved is shown in Fig. 5, with the free fermion results for comparison. For P = 1 the VMC achieves the same accuracy as the free fermion method. However, for P = 2 away from the critical point h = 1, the VMC performance is limited by sampling error [46].
We have introduced a variational technique that is motivated by both projector methods and recently developed quantum algorithms. It provides substantial shortcuts to the usual Trotterization of imaginary-time evolution, at the expense of making the procedure variational. Using TFIM as a first testbed, we have demonstrated that this ansatz is viable for sampling methods and highly efficient. In particular, the number of variational parameters required to represent the TFIM critical state scales as ∼log L, in contrast to other variational methods such as density matrix renormalization group, which in this critical case requires bond dimension scaling with L and thus the number of parameters scaling with L 2 . One reason for this efficiency is the fact that imaginarytime evolution is not subject to many bounds for real-time evolution; despite being generated by local Hamiltonians, our ansatz exhibits an exponential growth of entanglement entropy and rapid generation of long-range correlations, features unique to imaginary-time evolution.
Our variational approach is potentially useful in the many situations where imaginary-time Trotterization involves a prohibitively large number of steps. For example, in models with a sign problem, the computational cost scales exponentially with space and imaginary time O(τ L d ). Our ansatz provides a variational shortcut that significantly reduces τ [from τ ∼ L to τ ∼ (log L) 2.3 in the critical one-dimensional TFIM] which could enable the study of larger systems even with a sign problem. For few variational parameters (P = 1, 2), the optimization may be feasible, and we leave these investigations to future work.