Towards a variational Jordan-Lee-Preskill quantum algorithm

Rapid developments of quantum information technology show promising opportunities for simulating quantum field theory in near-term quantum devices. In this work, we formulate the theory of (time-dependent) variational quantum simulation of the 1+1 dimensional $\lambda \phi^4$ quantum field theory including encoding, state preparation, and time evolution, with several numerical simulation results. These algorithms could be understood as near-term variational quantum circuit (quantum neural network) analogs of the Jordan-Lee-Preskill algorithm, the basic algorithm for simulating quantum field theory using universal quantum devices. Besides, we highlight the advantages of encoding with harmonic oscillator basis based on the LSZ reduction formula and several computational efficiency such as when implementing a bosonic version of the unitary coupled cluster ansatz to prepare initial states. We also discuss how to circumvent the"spectral crowding"problem in the quantum field theory simulation and appraise our algorithm by both state and subspace fidelities.

Rapid developments of quantum information technology show promising opportunities for simulating quantum field theory in near-term quantum devices. In this work, we formulate the theory of (time-dependent) variational quantum simulation of the 1+1 dimensional λφ 4 quantum field theory including encoding, state preparation, and time evolution, with several numerical simulation results. These algorithms could be understood as near-term variational quantum circuit (quantum neural network) analogs of the Jordan-Lee-Preskill algorithm, the basic algorithm for simulating quantum field theory using universal quantum devices. Besides, we highlight the advantages of encoding with harmonic oscillator basis based on the LSZ reduction formula and several computational efficiency such as when implementing a bosonic version of the unitary coupled cluster ansatz to prepare initial states. We also discuss how to circumvent the "spectral crowding" problem in the quantum field theory simulation and appraise our algorithm by both state and subspace fidelities.

I. INTRODUCTION
Quantum information science is currently an important direction of modern scientific research across several subjects, including quantum physics, computer science, information technology, and quantum engineering. The rapid development of quantum technology brings us evidence that quantum computers in the near-term are able to perform some specifically scientific computations using dozens of qubits, but errors appearing in the noisy quantum circuits might set certain limits of the computational scale [1,2]. At the current stage, it makes sense to assume a reasonable quantum device exists and study potential scientific applications of such a device. This forms one of the main topics in the modern research of quantum information science.
Among numerous quantum applications, physicists, in particular, might care about how quantum devices could enlarge the range of computational capacity on certain problems in fundamental physics. In modern physics, quantum field theory is a general language or paradigm for describing almost all phenomena existing in the world, * junyuliu@uchicago.edu † lizm@mail.sustech.edu.cn ‡ hanz98@uchicago.edu § jinzhao.sun@physics.ox.ac.uk from sub-atomic particle physics, string theory and gravity, to condensed-matter and cold-atomic physics. If we could imagine the existence of powerful quantum computers, it will be natural, important, and interesting to consider if quantum computation could address open problems appearing in the study of quantum field theories, where many of them are at strong coupling and strong correlation. In fact, simulating quantum field theories in quantum devices is one of the earliest motivations of quantum computation [3], and becomes an important new research direction recently in the physics community, see the references [4][5][6][7] as examples. When simulating quantum field theories, or more generally, solving some well-defined computational tasks using quantum computation, theorists will either assume a universal, fault-tolerant quantum computer, or a noisy, near-term quantum circuit without enough quantum error correction. Both of them are wise choices and important scientific directions. Using fault-tolerant quantum computing is helpful for theoretical, conceptual problems or development of quantum devices usually appearing in the long-term, while studying near-term, early quantum computation will allow us to use existing machines and do experiments. In this paper, we will focus on the second direction, by exploring how far quantum simulation could go using near-term devices, with the help of specific problems in quantum field theories. It is helpful to see the usage and limitations of the currently existing, or future possible quantum hardware to simulate quantum field theories, and benchmark our quantum devices using interesting problems in fundamental physics [8]. Eventually, we believe that a universal, fault tolerant quantum device will come true, and we believe that our work might be helpful to speed up the process.
Here, we are specifically looking at the Jordan-Lee-Preskill scattering problem [6,7] in the 1+1 dimensional λφ 4 quantum field theory. The research about scattering problems has a long history in physics, from the scattering experiment of alpha particles by Rutherford to the modern discovery of the Higgs boson. Performing scattering experiments and determining scattering matrices are important themes in particle physics and quantum field theories. In Refs. [6,7], Jordan, Lee, and Preskill designed a full algorithm running in a universal quantum computer to perform particle scattering in quantum field theories, containing initial state preparation, time evolution, and measurement, where the proof of polynomial complexity is presented. In this work, we will construct closely-related algorithms that are more suitable for nearterm quantum computers.
We will be most interested in the circumstance where we have a machine to perform variational quantum simulation and hybrid quantum-classical calculations (see, for instance, Refs. [9][10][11][12][13][14][15][16][17][18]). In those algorithms, we will imagine that quantum gates or states are parametrized by a few parameters, and we iteratively perform measurements from quantum states and construct variational algorithms to optimize those parameters. We believe that those algorithms realized in the laboratory might be able to perform useful computations and could tell us something unknown about fundamental physics. In this work, we will systematically evaluate the possibility of variational quantum simulation in the context of λφ 4 quantum field theory. The paper is organized as following: • Basis choice. We will make a detailed comparisons between the field basis and the harmonic oscillator basis, momentum space, and coordinate space in Section II. All those choices have pros and cons. The field basis will cause the field correlations to be easy to measure and make the Lagrangian density local in coordinate space, but it will not be directly connected to the Feynmann rules and scattering calculations in the momentum space. Moreover, finding the eigenstates (for instance, the vacuum and low-lying one-particle states) might be not easy. It requires non-trivial digital quantum algorithms (with truncation error) for encoding. On the other hand, the harmonic oscillator basis is easy to formulate, track, and identify the energy levels of states, but may not be easy to identify field profiles. The Hamiltonian is non-local but still sparse in the momentum space. In this paper, we will work exclusively with harmonic oscillator basis (HO basis for brevity) as opposed to the field basis used in the original Jordan-Lee-Preskill algorithm. Besides the aforementioned advantages, free theory eigenstates are defined naturally under the HO basis through LSZ formalism as asymptotically far away in-states. We summarize the results in Proposition 1 and explain more details in Section.II F and in Appendix.
• Initial state preparation. In order to prepare the interaction wave packets, Jordan-Lee-Preskill algorithm uses adiabatic state preparation to turn on the coupling from free theory constructed under field basis. In the variational setup, alternative strategies could be directly used and solve the initial scattering directly. In Section III, we show how to prepare the state by variational algorithms. We introduce a bosonic version of the unitary coupled cluster (UCC) ansatz, which can be efficiently implemented on an NISQ device to test our simulation algorithms, and show the optimization using the imaginary time evolution [19]. Our initial state preparation strategies admit a simple interpretation according to the LSZ formalism: the UCC ansätze act on the free far-past particle eigenstate to some wave-packet in the interaction region at t i . The use of quantum imaginary time evolution further evolves this unknown wave-packet to low excited states at t f . We then theoretically and numerically investigate the spectral crowding phenomena in quantum field theories in both weakly-coupled and strongly-coupled theories.
• Real-time evolution and scattering. Besides digital quantum simulation algorithms (see Refs. [20][21][22][23][24]), the real-time evolution algorithms could also be tracked by variational methods, see Refs. [12,13,25,26]. During the real time evolution, variational errors might be hard to control especially for the non-perturbative regime and violent scattering processes. Nonetheless, we can track the simulation error during the time evolution, and we can adaptively construct the quantum circuit to achieve the desired accuracy within a polynomial circuit depth [16]. A theoretical framework for the dynamics simulation and comments on the challenges of the scattering process is provided in Section III F.
• Simulation fidelity. Simulation errors in the variational setting will not only be limited to the digital simulation error (like the Trotter error) but also the variational error from the ansatz, measurement error, and noise in the devices. In this work, we observe that in the initial scattering state preparation, as long as the total particle number and type are not changed significantly, the scattering experiment could still be performed, even starting with imprecise wave packets. Thus, the task of scattering state preparation could tolerate more noise. It could be qualified by particle subspace fidelity and suitable for NISQ devices. In Section IV, we show the numerical simulation for the ground state and excited states preparation using variational quantum algorithms, and compare it with adiabatic evolution by fidelities. We also provide a resource analysis of our method. Finally, in Section V, we highlight a list of future directions.
We summarize error and efficiency analyses in the following Proposition. Proposition 1. Let us consider a discrete λφ 4 theory on a spatial lattice Ω = aZ N with total length L, lattice spacing a and number of site N = L a . We also define its momentum space dual lattice Γ = 2π L Z N . Let us put n q qubits at each site (n q = log 2 (1 + 2φ max /δ φ ) for field basis and n q = log n cut for HO basis with more explanations on the notation in Section II). Then the following facts hold: 1. Let be an acceptable truncation error of truncated wavefunctions of free theory eigenstates under basis field. Then n q = O(log log 1 ). While eigenstates under HO basis are simulated by single computational basis elements without truncating any wavefunction. On the other hand, both simulated field operators φ(x), π(x) (Eq. (12)) and ladder operators a † k , a k (Eq. (13)) do not satisfy the canonical commutation relation perfectly. Let denote the corresponding error, then n q = O(log 1 )) (see Section II F & Appendix).

The full theory Hamiltonian H is sparse under
HO basis with O(N 3 ) nonzero entries for each of its row/column and can be complied by O(N 3 n cut ) Pauli operators (Section II F). The UCC ansatzT 2 (Eq. (22)) proposed in Section III B can be be compiled by O(N n 2 cut ) Pauli operators acting on 2n q qubits. To judge our algorithm, an n-particle subspace fidelity can be measured in O(N n ) times (see Section III E).

II. FORMALISM AND STATE ENCODING
At the starting point, we show how to encode our Hamiltonian from quantum field theory to a quantum device. In Section. II A, we first give a review of the λφ 4 theory in 1+1 dimension, and we point out the use of harmonic oscillator basis in calculating scattering amplitudes in the scalar field theory by Lehmann-Symanzik-Zimmermann (LSZ) reduction formula. Then, we provide a detailed comparison on various versions of bases, including the field basis and the harmonic oscillator basis in the coordinate space and momentum space, respectively. Some similar discussions can be found in the quantum chemistry context [27].
A. λφ 4 theory and the LSZ reduction formula In this theory, we have a scalar quantum field φ with the Hamiltonian Moreover, we discretize it in the lattice, The theory is defined on the spatial lattice Ω with dual lattice Γ (see Proposition 1). The field momentum π(x) is defined as the Fourier conjugate of the field φ(x) with the following commutation relation, The discretized version of the derivative is given by where m 0 is the (bare) mass term in the free theory, and the λ0 4! φ(x) 4 term represents the coupling. When λ 0 = 0, we call it the free theory. In the case of the free theory, we could diagonalize the Hamiltonian by the following mode decomposition in the continuum, The discrete version can be defined similarly, where the canonical algebra of φ and π leads to the commutation relation as [a k , a † l ] = Lδ k,l for lattice and [a k , a † l ] = 2πδ(k − l) for continuum, respectively. The energy dispersion is given by In such a basis, the Hamiltonian is diagonalized as where An important physical observable we could measure, is the (Wightman) two-point function as G(x − y) = Ω|φ(x)φ(y)|Ω , where |Ω is the ground state of the theory. In the free theory case where λ 0 = 0, one can compute the two point function explicitly The two-point function of the scalar defines the scalar mass of the theory. In the continuum limit, when we turn on the interaction λ 0 , in weakly-coupled regime one could compute the correction to the mass through Feynman diagrams. The theory will experience a second-order phase transition at strong coupling, where the universal behavior belongs to the 2D Ising universality class and the use of perturbation theory is difficult. In Appendix, we address the relation between lattice models and their field theory description, emphasizing the importance of simulating quantum field theories. From a non-perturbative perspective, computing two-point functions will tell us the information about masses of particles through Källén-Lehmann spectral representation [28,29]. We also review necessary backgrounds through LSZ reduction formula in Appendix that is particularly suited to the use of HO basis in this paper.
B. The field basis in the coordinate space One of the simplest considerations is the field basis. For a scalar quantum field theory discretized in a lattice Ω, we could define the state decomposition Here, N is the total number of sites. The state decomposition for an arbitrary state ψ gives the above wavefunction ψ (φ 1 , . . . , φ N ), where we abuse the notation φ i chosen as an arbitrary number in R to denote an eigenvalue of the local field operator φ(x i ). The corresponding eigenstates |φ i form a basis of the local Hilbert space. This definition is similar to the coordinate basis in quantum mechanics. Now, since we are using a quantum computer, we need to truncate the local Hilbert space. Moreover, we want an increment δ φ in discretization of the spectra of each φ(x i ) such that we do not need to choose variables in a continuous interval. The states that we are interested in is truncated and discretize as and thus the number of qubits we need to encode at each site x i is n q = log 2 (1 + 2φ max /δ φ ) . There are bounds on φ max , δ φ and n q from the scattering energy E derived in Refs. [6,7] which are useful to prove the polynomial complexity of the Jordan-Lee-Preskill algorithm. However, the original bound n q = O(log 1 ) of qubit number with respect to truncation error is not tight, we compute rigorously by several properties of Hermit-Gauss functions (eigenfuncations of quantum harmonic oscillator) and show in Appendix that this bound can be refined as n q = O(log log 1 ).

C. The HO basis in the coordinate space
There is another important basis, the harmonic oscillator basis to define a digital representation of states in lattice quantum field theories. We firstly consider the following transformation, with π(x) being transfomred similarly. Here m x is a free parameter we could choose. Then the canonical algebra of φ and π leads to a x , a † y = a −1 δ x,y . Now, the creation operator a † x and its conjugate could define the number states |n x at the site x. On the lattice Ω of N sites, let us say that we are mostly interested in the maximal energy level n cut , so we cut the Hilbert space and define Now we introduce the dual field basis in the momentum space. Remember that we define the dual lattice Γ based on the spatial lattice Ω. Thus, one can directly write the Hamiltonian in terms of the momentum coordinate. To be more specific, consider the free theory mode expansion Eq. (1) with Fourier transformation of φ(x) and the dual field momentum π k being defined similarly. Then we can discretize the interaction piece of the Hamiltonian in the momentum space by We can make a truncation on the field range in the momentum space and discretize a state like Eq. (7).

E. The HO basis in the momentum space
Similarly, we could consider the HO basis in the momentum space. Under this basis, a discretized state decomposition is written as |ψ cut = n k i ≤ncut ψ(n k1 , . . . , n k N ) |n k1 , . . . , n k N . (11) The number of qubits needed at each momentum mode k i is thus n q = log n cut . The number states |n ki now are generated by the creation operator from Eq. (1). The above state has a very clear physical meaning: the basis directly show the scalar particle numbers in different momenta. This also provides a good initial guess for the excited states in the interacting theory. In Section IV, we discuss the particle excitations in the momentum space in more details.

F. A comparison
Besides a brief comparison mentioned in the Introduction, we discuss more details on discretization of scalar field under different bases. For each position site x i , a truncated field operator φ(x i ) is given by Pauli Zmatrices: with π(x i ) being defined as the discrete Fourier transform of φ(x i ). The corresponding free theory vacuum is then the discrete Gaussian prepared by Kitaev-Webb Algorithm. It is shown in [30][31][32] that they can be efficiently constructed in a quantum circuit. On the other hand, the creation/annihilation operators on momentum space are defined as where k i ∈ Γ specifies a momentum mode with |s being the computational basis of n q = log n cut qubits at each mode. The corresponding free theory vacuum in this case is simply |s = 0 k1 ⊗ · · · ⊗ |s = 0 k N and hence can be prepared at a constant circuit depth. This is also true for other excited initial state used in our numerical simulation (e.g., Eq. (59) & (58)). Easily constructible initial states is the first advantage when working with HO basis. Besides, these states are taken without truncating any wavefunction and hence there is no need to consider the truncation error. Even though, we should mention that the aforementioned discrete commutation relations false for both truncated energy level of HO basis and truncated field strength of field basis. The corresponding error only decays exponentially with n q . More details can be seen in Appendix.
Furthermore, as the computational basis encodes particle numbers of momentum modes, the free theory Hamiltonian Eq. (3) is automatically diagonalized. To check the implementation efficiency of interaction Hamiltonian Eq. (10), we first examine the sparsity. By Eq. (9), each row/column of the matrix representation of φ ki contains at most 2 nonzero entries. Even H int is non-local, any of its terms φ k1 φ k2 φ k3 φ −k1−k2−k3 is a four-fold tensor product and hence contains at most 2 4 nonzero entries in each row/column. Comparing with the 2 nqN -dimensional total Hilbert space, each term is sparse. Because H int has N 3 terms when summing over the dual lattice Γ with momentum conservation, nonzero matrix elements scales cubically with the number N of momentum mode. This makes H applicable under most existing quantum algorithm, especially the imaginary time evolution employed in this work. [33] To count the number of Pauli operators needed to compile this Harmonization, we first consider how to expand creation/annihilation operators by the Pauli basis {I, σ x , σ y , σ z } ⊗nq at each mode. To calculate the needed number, we transform the Pauli basis into matrix unit basis {E ij }: consider the n q = 1 (onequbit) case: Each matrix unit E ij can be written as a linear combination of 2 Pauli operators. Hence the corresponding transformation matrix M has 2 nonzero entries in each of its columns. Since the creation operator a † on one-qubit space is colinear with E 12 , it decomposes into 2 Pauli operators. For two-qubit space, the transformation matrix simply equals M ⊗2 with 2 2 nonzero terms. As a † is now a linear combination of E 12 , E 23 , E 34 and one can check that E 12 , E 34 are transformed from the same subcollection of Pauli operators, a † decomposes into 2 · 2 2 pieces in total. By induction, the decomposition of a † on n q -qubit space has n q ·2 nq pieces. With respect to the energy cut-off, we need n cut +log n cut Pauli operators which scales linearly with n cut . By the same method, the Hermitian operator a + a † can be built by 1 2 (n cut + log n cut Pauli operators with the factor 1 2 coming from cancellation of anti-Hermitian terms when we sum a and a † together. As a simple example to verify this point, let n q = 1, then a + a † is colinear with a single Pauli operator σ x . On the other hand, by Eq. (10) & (9), H int is built by at most 8N 3 (n cut + log n cut ) Pauli operators. We will apply this method to verify gates efficiency of the UCC variational ansatz in Section III and Section IV.
One can check that the above analysis automatically holds for general d+1 dimensional theory where N stands for the number of momentum nodes in d dimension. Except the efficiency of preparing initial states and implementing the Hamiltonian, the HO basis is also useful to keep track of the simulation results in real-time, since one could quickly identify the basis overlap and find the particle number and their momenta. Indeed, the harmonic basis specifies the momentum sectors of the asymptotically far past in-states, where the fields satisfy the on-shell condition. For the interacting theory, when the interaction is turned on, one could specify the momentum sectors again by the adiabatic state preparation and we could use this method to define the wave packets in the given momentum sectors from adiabetically preparing a † F (t), as shown in Appendix. Thus, in this paper, we will mainly work on the HO basis in the variational setup.

A. The variational ansatz
Variational quantum simulation is a useful technique especially for the near-term quantum computer. The variational algorithm starts by preparing the quantum state by a quantum circuit as Here U s are some unitary operators that could be realized in the quantum device, for instance, U (θ ) = e −iθ X , with the variational parameters θ = (θ ). X s are some Hermitian operators, for instance, elements in the Pauli group, and |ψ 0 could be some simple initial states that could be easily prepared. The target state will, in principle, be approximated by some optimal choices of θ, say θ * , which could be found using the variational principles. For example, a typical problem in quantum simulation is to find the ground state, then we could minimize the energy with respect to the variational parameters H θ ≡ ψ(θ)| H |ψ(θ) . The general strategy for the ground state searching is by updating the parameters as where θ µ (t) represents the optimization dynamics with step t, and the learning rate is given by η µ (t). Here, we use A(θ(t)) to represent the metric matrix at the parameter θ(t). The metric matrix in the gradient descent algorithm is simply the identity matrix. In the following section, we will show its explicit form during the optimization. One can ask if there exists a regime where there is a convergence guarantee and, if so, the rate of convergence for these variational paramterization. One can study this question from over-paramterization using quantum neural tangent kernel (QNTK) [34]. Further taking H θ(t) ≡ z(θ(t)), Eq.(15) implies: Assuming A is identity matrix and η to be paramterization-independent, the resultant is precisely the QNTK defined in [35]. Note that we can interpret A −1 µν as the learning rate tensor as part of definition of NTK in classical neural networks [31]. In particular, for the circuits that form at least approximate 2-design that satisfies certain concentration conditions (See in [34]), the average convergence is of the form where ≡ z(θ) − E 0 , E 0 is the ground state energy and with L being the total count of variational parameters and H be full Hamiltonian. The dimension of the Hilbert space in our case is n N cut . The exponential convergence rate is guarantee on average in the over-parametrization regime where L ≈ dim(H) 2 / tr(H 2 ). When A fails to be an identity matrix such as in the case of quantum imaginary time evolution used in the following, no precise analytical convergence guarantee is known. However, it seems to be reasonable to extrapolate the hypothesis that such methods, due to its more physical/geometric nature, would have convergence rates lowered-bounded by the naive gradient descent methods. The dependence of the square of size of Hilbert space would imply the above analysis only is suitable to small size qubit system (See in Appendix for more details).
The next question is how to choose the initial state |ψ 0 and U s? The precise strategy of choosing |ψ 0 , U and the optimization scheme will specify the variational quantum algorithm we use. There are many variational algorithms (see Ref. [14,15] for a recent review). In this work we will discuss the following bosonic unitary coupled cluster (UCC) ansatz and imaginary time evolution where we practically find the best in our physical system. Different from the quantum computational chemistry literature, where the UCC ansatz consists of the fermionic excitations in the active space, our algorithm expresses the UCC ansatz directly with the bosonic mode.

B. Bosonic UCC ansatz
As is mentioned before, the variational algorithm may not be very sensitive to the locality of the Hamiltonian. With known implementation efficiency, we will focus on the HO basis in the momentum space. Prior work has extensively investigated the coupled cluster methods to solve the electronic energy spectra and vibrational structure in the chemistry and materials science, and the quantum version, unitary coupled cluster ansatz, has been suggested and further experimentally demonstrated to solve the chemistry problems on a quantum computer [36,37]. Other prior works [38,39] discussed the usage of bosonic UCC in studying vibronic properties of molecules.
The general form of unitary coupled cluster is given by whereT is the sum of symmetry preserved excitation Hermitian operators truncated at finite excitations aŝ The key ingredient of UCC ansatz is to search for the true ground state of the interacting fermionic theory by considering the particle-conserving excitations above a reference state.
In our quantum field theory setup, a bosonic version of the UCC ansatz [38,39] could be natural to capture types and particle numbers for scalar particles. In the momentum space, the effective action preserves the momentum reflection symmetry (k → −k =k). Therefore, we may express the lth excitation Hermitian operators of the bosonic UCC ansatz in the momentum space aŝ . (20) Here, k 1 , ..., k are distinct momentum modes taken from the lattice Γ. Since the λφ 4 field could lift the excitation up to 4 level, we impose the energy constraint |s i − t i | ≤ 4 forT at each momentum mode which makes each term ofT spares like the λφ 4 Hamiltonian. To count Pauli operators, we first note that the concerned local Hilbert space is defined by n q qubits (recall that n q = log n cut ) and then apply the same from Section II F toT . With energy constraint, it can be built by O(n cut ) Pauli operators and each of which acts nontrivially on at most n q qubits. The total ansatz is thus compiled by O(N n cut ) Pauli operators with the same order of number of parameters. When = 4, expanding H int by Eq. (10), (9) & (13), we can set parameters ofT 4 being the expansion coefficients. We can even vary these parameters as θ(s), s ∈ [0, 1] such that H int =T 4 (θ(1)) and hence is tantamount to adiabatic turn-on of the interaction. To simulate the adiabatic evolution, we have to divide the time interval [0, 1] into M pieces with M large enough (depending on the energy gap). We then apply Trotter formula to approximate the time evolution using O(N 4 n 4 cut ) Pauli operators for each product term. In NISQ devices however, we wish to further reduce the computational cost. Thus we focus on variational ansätze and restrict to use single excitation operatorT 1 which can be constructed by at most 4N 1 2 (n cut +log n cut ) Pauli operators. The second term is obtained like expanding a + a † in Section II F. We also employ a modified double excitation operatorsT 2 aŝ which considers the pairing correlations of the momentum k andk. It makes k the only momentum variable when taking summation. With the requirement to be Hermitian and the energy constraint, this ansatz can be built by at most 32N 1 2 (n cut + log n cut Pauli operators such that each of which acts nontrivially on at most 2n q qubits and hence reduces a large number of parameters comparing Eq. (20). We may even discard the second term in Eq. (22) to further reduce the gate count in the variational quantum circuits. For small n cut , the number of used Pauli operators would be even fewer (see Fig. 9 in Section IV). As an inevitable consequence, these variational ansätze cannot replace adiabatic evolution in searching ground state. We will remedy this problem by employing the quantum imaginary time evolution in the next section.

C. Variational state preparation
We now discuss how to use variational quantum algorithms for finding the ground state and the low-lying excited states. We first briefly review the variational quantum simulation algorithm of imaginary time evolution [19]. The normalized imaginary time evolution at imaginary time τ is given by The population of the energy eigenstate |e j will decay exponentially with the energy E j , and the ground state can be obtained in the long time limit ψ (0) = lim τ →∞ |ψ(τ ) . While the nonunitary imaginary time evolution cannot be directly implemented on a quantum computer, one could still simulate imaginary time evolution on a quantum computer by using the hybrid quantum-classical algorithm. Instead of simulating the imaginary time evolution directly, we assume that the time-evolved state can be approximated by a parametrized trial state |ψ(θ(τ )) , with variational parameters θ(τ ) = (θ µ (τ )). As mentioned in [19], by minimizing the distance between the ideal evolution and the evolution of the parametrized trial state, the evolution of the target state |ψ(τ ) under the Schrödinger equation can be mapped to the trial state manifold as the evolution of parameters θ.
Using McLachlan's variational principle, we have and the evolution of the parameters under the imaginary time evolution could be determined by with the matrix elements of A and C given by Here, |ψ = ψ|ψ is the norm of the quantum state, we denote ∂ i ≡ ∂/∂θ i , and we assume the parameters are real. By tracking the evolution of the variational parameters, we can effectively simulate imaginary time evolution. This actually serves as an optimizer to update the parameters in Eq. (15). It is worth mentioning for the pure state imaginary time evolution, this approach is equivalently to the quantum natural gradient descent method, and the matrix A is indeed the Fisher matrix [40].
The quantum imaginary time evolution minimizes the energy loss function: where its total variation [19]: where the fact that the ij C i A −1 ij C j is nonnegative follows that the fact that A is non-negative definite. The nonnegativity of ε follows from the variational theorem where E 0 is the smallest eigenvalue.
Moreover, having found the ground state ψ (0) , we can construct a modified Hamiltonian , where α is the regularization term that lifts the ground state energy, and is sufficiently large comparing to the energy scale of the system.
The ground state of the modified Hamiltonian H (1) becomes ψ (1) , the first excited state ψ (1) of the original Hamiltonian H. As ψ (0) is an excited state of the modified Hamiltonian, we can evolve the system under H (1) in the imaginary time to suppress the spectral weight of ψ (0) and obtain the first excited state ψ (1) . This process can be repeated to obtain the higher-order excited states. More specifically, the (n + 1)th excited state is the ground state of effective Hamiltonian Instead of preparing the Hamiltonian directly, we can simulate the imaginary time evolution under H (n+1) by tracking the evolution of the parameters, which are now modified as while the matrix A keeps the same as in Equation (25). These addition terms in C i can be evaluated using the low-depth swap test circuit. Other variational excited state preparation techniques can be found in a recent review paper [14].
We wish to remark that the circuit ansatz for the imaginary time evolution does not have to be fixed. Instead, the circuit ansatz could be adaptively determined by tracking the distance of the ideal evolved state and the variational state. In the extreme case, we could construct the circuit by approximating the normalized state at every single time step, which reduces to the quantum imaginary time evolution, firstly proposed in Ref. [41]. Suppose the Hamiltonian has the decomposition H = L l=1ĥ l , where the Hamiltonian contains L local terms and eachĥ l acts on at most k neighboring qubits. Using the first-order Trotterization, the evolved state after applying nonunitary operator e −δτĥ l within imaginary time δτ by where c is the normalization factor andÂ is a Hermitian operator that acts on a domain of D qubits around the support ofĥ l . The unitary operator e −iδτÂ can be determined by minimizing the approximation error in Eq. (30), which is similar to the derivation in Eq. (23). For a nearest-neighbor local Hamiltonian on a d-dimensional cubic lattice, the domain size D is bounded by O(C d ), where C is the correlation length. More details about the algorithm complexity can be found in Ref. [41].
This circuit construction strategy can be regarded as a special case in the variational imaginary time evolution given by Eq. (23) and Eq. (25). If we fix the old circuit ansatz θ(τ ) constructed before imaginary time τ , and determine the new added unitary operator θ(δτ ) that approximates the effect of e −iHδτ , this is exactly the same as Eq. (30). However, to further reduce the circuit depth, we can jointly optimize the parameters θ(τ )⊕θ(δτ ), making it more compatible for the near-term quantum devices.

D. Spectral crowding
Before we start to apply variational algorithms, we will make a short investigation on the spectrum of the λφ 4 quantum field theory. In the momentum space, HO basis, one might have a large number of degeneracies in the energy eigenstates (similar problems appear in other bases as well), bringing potential problems for quantum simulation. We will borrow the terminology "spectral crowding" that has been used in the ion trap systems [42] referring to this situation.
For excited states, degeneracy might happen even in the free theory in our construction. For instance, say that in the free theory, it might be the case where Here, we have states represented in the HO basis in the momentum space, with particle numbers and momenta Spectral crowding for λφ 4 theory: Blue/left: free theory; Red/right: interacting theory. We use m0 = 0.369 and λ0 = 0 (free) or λ0 = 0.481 (interacting), with maximally three excitations ncut = 3, system size N = 4, and the lattice spacing a = 1, for the HO basis in the momentum space. n i , p i , orn j ,p j , and their energies are precisely identical. A typical example is that considering the continuum limit, we might have where n ∈ Z >0 . In those cases, their states are degenerate. Another typical case the role of parity which anticommutes with the momentum: since we are not able to distinguish the left-moving and right-moving states only by their energies. Figure 1 provides an example for spectral crowding, where we fix m 0 = 0.369 and λ 0 = 0 (free) or λ 0 = 0.481 (interacting), with maximally three excitations n cut = 3, system size N = 4, and the lattice spacing a = 1. The choice of parameters is aiming on avoiding the situation in the Eq. (32). Spectral crowding might bring us difficulties on identifying states in the output, and defining different directions of momenta for particles, especially when states are excited. Instead of looking at the general structure of density of states, we start with the maximally oneparticle states in this simple system. In the free theory, we have the ground state with the energy 2.662. Moreover, the single-particle states have the energies: We know that this degeneracy is made by the boundary condition of the momentum p ∼ 2π/L − p, which is the parity Z 2 [43]. In general, for an n-particle state, since we could freely choose the direction of momentum, the spectral crowding will be enhanced at least O(2 n ). Now, we consider to turn on the interaction. In the adiabatic process where we slowly turn on the φ 4 Hamiltonian as Eq. (21). Since the interacting Hamiltonian is invariant under the parity transformation, we could use the adiabatic process to define the direction of the momentum. In Figure 2, we show an example of the adiabatic evolution numerically, with the number of adiabatic steps T = 100 (which means that we are dividing the interval s ∈ [0, 1] to 101 steps). We find all single-particle eigenstates could agree with the corresponding energy eigenstates with high fidelities (we only show p = 2π L (1, 3) example in the plot, but all four adiabatic state preparations are also verified). Note that this operation specifies the direction of momenta in the interacting theory. This is an advantage of our basis, where we could specify the direction of momenta in this way.
The above algorithm could also be made variationally. Recall that in the variational process, we are starting from a wave packet state |ψ initial , and we slowly turn on the interaction λ from the free theory λ = 0. Thus, during this process, the Hamiltonian is time-dependent. Instead of considering Lie-Trotter-Suzuki decomposition in a digital quantum computer, one could perform the above calculation in a quantum computer with a variational form. We will use the variational approach of time evolution introduced in Refs. [13,19]. Similar to the imaginary time proposal, we will use the McLachlan's variational principle and take care of the time-dependent global phase. Now consider the situation where we adiabatically turn on the coupling of the Hamiltonian. We restrict our solution inside the variational form similar from before, Note that we are starting from the corresponding momentum eigenstates of the free particle. The differential equation of θ based on the McLachlan's variational principle is given by where M i,j = ReA i,j + ∂ i ψ(θ(s))|ψ(θ(s)) ∂ j ψ(θ(s))|ψ(θ(s)) , V i = ImC i + i∂ i ψ(θ(s))|ψ(θ(s)) ψ(θ(s))|H(s)|ψ(θ(s)) , and A and C are similarly defined in Eq. (25). One could show that the solution of θs are always real in our variational form, and the variational answer is consistent with the actual answer up to a time-dependent global phase. More detailed discussions on the simulation error during the dynamics can be found in Section III F. Similar to the above variational adiabatic state preparation algorithm, the imaginary time evolution could also start from the corresponding free theory states. Practically, we find in our example, the imaginary time evolution algorithm performs better (this is intuitively because we are looking for low-lying states with low energies).
Moreover, these methods can be integrated together. We can turn on the interaction, similarly as in Eq. (21), but with much fewer steps. We could then use the variational algorithms to find the ground state of the intermediate Hamiltonian H(s), using which as the initial state for the next time step until s reaches 1. Compared to finding the ground state of H I , this method may avoid the local minimum.
Finally, we comment on other methods to snake around the spectral crowding problem. A useful trick to find the state with both fixed momentum and energy is through measurement in quantum devices. We could consider measuring the momentum operator during the variational process, making sure that it keeps the sign when the interaction is turning on. However, the momentum operator only has its meaning in the free theory, so we only expect the above algorithm to be useful in the sense of weakly-coupled theory. Another useful trick for keeping the momentum is similar to the idea of the tangent space method in the language of matrix product state (see a review [44]). Usually, we expect that our momentum−p eigenstate could have the following form Here T x is the translation operator with the vector x. If the state |Φ is already a momentum-p eigenstate, we have Moreover, if the state |Φ is a linear superposition of the momentum-p state and the momentum-(−p) state where Note that the c 2 term is suppressed because it sums over a pure numerical phase. Thus, for the state we obtained from the variational quantum simulation, we could make a linear superposition weighted by e ipx to obtain a momentum eigenstate with a fixed momentum direction, at least in the case of the single-particle scattering experiment. However, the above method seems to be mostly useful when we know how to construct the translation operator. It is manifest in the coordinate space, but not easy in the momentum space.
E. State fidelity, one-particle subspace fidelity and generalizations Here we discuss some concepts about fidelities that are useful for the variational, scattering-state preparation setting. Say that we originally have a wave packet centered around a given momentum, and it is a one-particle state in the free theory. Now we could turn on the interaction slowly. Ideally, as we discussed before, a oneparticle state will still remain a one-particle state in the interacting theory. In fact, if we consider momentum eigenstates of a single particle, |p , we could define the one-particle subspace by V one-particle,free = span p (|p ).
Now, if we are adiabatically turning on each state |p towards the coupling λ 0 , the space will become V one-particle,λ0 = span p (adiabatic evolution λ0 • |p ).
In fact, if the adiabatic evolution is slow enough, the above expression will define the one-particle space in the interacting theory. This makes our number eigenstate definition more precise. Counting the one-particle eigenstates in free theory on the lattice Ω, they span an N -dimensional subspace. The dimension of a n-particle space is equal to the number of compositions (n 1 , ..., n N ) of n. It is mathematically analogous to calculate the dimension of nth symmetry tensor power of R N and the answer is n+N −1 n . It should be noted that since we also truncate particle numbers by n cut at each momentum mode, an n-particle space with n > n cut cannot be fully realized and thus has a lower dimension. Even though, low energy subspaces can be fully defined with dimension being bounded polynomially in N . Now, say that we are doing the state preparation using the variational algorithm (which is not the ideal adiabatic process). Due to the limitation induced by the variational ansatz, we will have some systematic errors (or other errors). However, the resulting state, although suffering from the noise, might still have a large overlap with the one-particle subspace V one-particle,λ0 . In fact, we could define the state fidelity which is an overlap between the accurate state from an ideal adiabatic simulation without any error, and the state obtained from the variational algorithm. We could also define the one-particle subspace fidelity F one-particle,adiabatic = | ψ variational |Λ|ψ variational | . (45) Here, Λ is the projector of the space V one-particle,λ0 . By definition, F one-particle,adiabatic should be no less than F state,adiabatic (e.g., Fig. 4 (b), (c) or (e), (f )). In principle, we wish our fidelities to be always high enough, which means that we are performing high-quality state preparations. Ideally, we wish the state fidelity to be high. But in practice, when we do not really care about the form of the wave packet, and we only care about if the state is still approximately a one-particle state, we could only use the one-particle subspace fidelity. As a summary, the level of requirements we want about fidelities is closely related to the actual physical motivation we have in the simulation experiment. Let us end this subsection by making a final comment on the fidelities. The definition of fidelities is indeed related to the physical task we want when doing the experiment. The definition of one-particle subspace fidelity corresponds to the choice when we wish to maintain the one-particle subspace during the initial state preparation. When we have other requirements, we could demand other versions of fidelities be high. For instance, we could define the momentum fidelity by measuring the momentum center of the wave packet. We could also define the wave packet profile fidelity by measuring the wave packet form. Stronger definitions on fidelities would require higher quality when we are doing the variational state preparation.
In the most general setting, we could define the projector to the subspace Λ as where |q i is the basis and D Λ is the dimension of the subspace. Thus, for a variational state |ψ we have Here c i and c e are the expansion coefficients, and |e is the perpendicular component of the subspace projector Λ. Say that the state is normalized, we have At the same time, the ideal state is given by the expan- thus the state fidelity is given by This equation illustrates the relation between the state fidelity and the subspace fidelity. In the limit where εs are small, those two fidelities are almost equal [45]. We verify numerically in Section IV a high one-particle subspace fidelity (Fig. 4) as well as state fidelity (Fig. 5, 6 7 & 8) for a lower number of clean qubits as a successful benchmark for our variational algorithm about the adiabatic state preparation before particle scattering. State fidelities can be easily measured by Hadamard test with the help of adiabatic quantum computing. For n-particle subspace fidelities, as we analyzed above, the projector Λ defined in Eq. (46) can be built by at most n+N −1 n eigenstates evolving adiabatically from the free theory. Thus measuring O(N n ) times like the case of state fidelity, we can compute the subspace fidelity. Especially for n = 1, one particle subspace fidelity can be obtained by N -time measurements.

F. Particle scattering
Similar to the variational state preparation, we could also make the variational version of particle scattering. Now, when we are considering the variational time evolution, the only difference comparing to the adiabatic state preparation, is that now the Hamiltonian is static, and the bare coupling λ 0 is fixed. The evolution of the variational parameter is given by where M and V can be similarly expressed as M i,j = ReA i,j + ∂ i ψ(θ(t))|ψ(θ(t)) ∂ j ψ(θ(t))||ψ(θ(t)) , V i = ImC i + i∂ i ψ(θ(t))|ψ(θ(t)) ψ(θ(t))|H|ψ(θ(t)) .
Here the Hamiltonian H does not depend on the time t. Now we show how to approximate the ideal evolved state during the dynamics up to given error ε. Suppose at time t the ideal state Φ(t) is approximated by Φ(t) ≈ |Ψ(θ(t)) . Within time step δt, we evolved the state |Ψ(θ(t)) by updating the parameters θ(t) to θ(t + δt), which introduces an approximation error at t + δt as δε = |Ψ(θ(t + δt)) − e −iHδt |Ψ(θ(t) .
Minimizing the error will give the similar results as that from the McLachlan's variational principle in Eq. (51).
In the extreme case, if we choose the unitary operators in the ansatz U from all the Hamiltonian terms (ĥ l ), for instance the Trotterization, this error could be reduced to zero. This indicates that for single step, we can guarantee the simulation error at each time t up to certain threshold. We can also track the accumulated error during the whole scattering process. Starting from an initial state |Ψ 0 , the accumulated error until time t + δt can be bounded by where we have used triangle inequality and the distance invariance under the unitary transformation in the second line. The single step error can be bounded by where the first-order order error is with the matrix A and C defined in Eq. (25). The total error during the simulation can be bounded by where max ∆ is the maximum error during the evolution. In practice, we could add the operators from the Hamiltonian term (ĥ l ) into the circuits to decrease the error to a certain threshold ε 0 by setting ∆ cut = ε 0 /T . Therefore, by tracking the simulation error at each step, we can ensure the simulation accuracy. If at time t, the error ∆(t) is measured to be above the threshold, i.e., ∆(t) > ∆ cut , we repeat to add new operators from the Hamiltonian term (ĥ l ) until ∆ ≤ ∆ cut The efficiency of the adaptive strategy is guaranteed by the following theorem.
1. The first-order error ∆ strictly decreases at each iteration until 0; 2. In each circuit construction process (if ∆(t) > ∆ cut ), each Pauli term,ĥ l , in the Hamiltonian is only needed to appear once; 3. We can achieve an error ∆ ≤ ∆ cut in at most L iteration for any ∆ cut ≥ 0. Here, L is the number of terms in the Hamiltonian The key idea of the proof is that in the circuit construction subroutine, there always exists an operation h k ∈ (ĥ l ), by appending which to the old circuit, the distance strictly decreases if ∆ = 0. Theorem 1 indicates that circuit construction process will terminate in a finite number of steps during the total time evolution. In an extreme case, we can optimize the parameters directly to make ∆ ≤ ∆ cut , such that no additional gates are required to be added. This reduces to the conventional variational algorithms in Eq. (51).
We also remark that ∆ is a measurable quantity. The additional measurement cost for the adaptive circuit construction comes from H 2 , which could be measured efficiently by using the compatibility of the Pauli operators in the Hamiltonian. For instance, ifĥ l andĥ k qubit-wise commute with each other, we can simultaneously measure them within one Pauli basis, which can significantly reduce the measurement cost.
We conclude this section by making the following comments about the variational realization of the particle scattering algorithm.
• The bosonic ansatz Eq. (20) in the momentum space allows the creation of new particles during the scattering process, and they could capture the particle excitations along time evolution. • Since we are considering the scattering process of the wave packet states, some challenges might appear because of the limitation of the variational ansatz: we cannot cover the full space during the variational simulation. Furthermore, besides the error appearing in the near-term quantum devices, we might have some other errors in the variational process due to the level crossing phenomena among different excited states. For a given theory, lots of tests need to be done to obtain some numerically satisfying results, and we leave those opportunities to future research.
• During the scattering process, we might wish to read off some explicit results for the S matrix elements. Thus, the result should be sensitive to the error, from the adiabatic state preparation to the scattering process. In this situation, we don't want uncontrolled errors from the quantum noise or some systematic errors from the assumption of the variational ansatz. However, if we only want some collective, statistical properties of the output states, for instance, some macroscopic quantities or random averages that could contain some intrinsic noises (for instance, the jets), we might have fewer constraints on the fidelity of the variational algorithm.
• Other hybrid-classical quantum simulation methods, such as hybrid tensor networks [46], could be leveraged to simulate this scattering process with fewer quantum resources. Moreover, perturbative quantum simulation methods that do not rely on the circuit ansatz could be applied to this task [47].

IV. NUMERICAL RESULTS
In this section, we demonstrate how the techniques can be used to find the ground state and excited states of the interacting lattice field. We also discuss the spectrum of the lattice field with different bare mass and coupling strength.
Similar to the analysis before, we consider a lattice Ω with total length L = 4 and lattice spacing a = 1, and its dual lattice Γ has 4 sites in the momentum space. We use the HO basis in the momentum space with n cut = 4 energy levels. To benchmark the performance of the variational algorithms, we consider finding the ground state using the bosonic UCC ansatz with an increasing coupling of the interacting field. We prepare the initial state as the ground state, |0 ⊗8 , in the free theory. Here, we use the compact mapping for the creation and annihilation operators in Eq. (13). We truncate the highest energy level to be 3 in the double excitation operator in Eq. (22) to reduce the number of parameters and the quantum circuit depth. In order to find the variational state, we use the imaginary time evolution to evolve and identify the low-lying spectra in the interacting theory. The regularization term in the excited state search is fixed as α = 8. In the numerical simulation, the error of the results is compared with exact diagonalization.
The relative error of the ground state energy and the state associated with different coupling strengths of the interacting field, λ 0 , is shown in Figure 3. We characterize the relative energy error by (E variational − E ideal )/E ideal , where E ideal is the correpsonding eigenenergies calculated by the exact diagonalization. The state fidelity is defined by the overlap between the variational state and the ideal excited state ψ ideal |ψ variational . We can see that the simulation error increases when the interaction strength λ 0 increases, ranging from λ 0 = 0.5 ∼ 4!, but even for a large interaction strength λ 0 = 4!, we can achieve a high simulation accuracy below 10 −3 both for the energy error and the state fidelity, which indicates a strong representation capability of the quantum circuit ansatz. In the following, we will choose two coupling strength λ 0 = 1 and λ 0 = 10 to test the performance of the variational algorithms in several regimes [48].
Moreover, we extend the discussions to the excited state preparation, which could be more complicated due to the spectral crowding and degeneracy of the lattice field. We first consider the static mass m 0 = 1.27, such that the energy of single-particle excitation is lower than that of multi-particle excitations in both the free field and interacting field. In this case, we prepare the initial state for the excited states, searching in the corresponding single-particle excited-state space of the free Hamiltonian λ 0 = 0.
We show the relative error of the energy and the fidelity of the ground state and low-lying excited states towards the iteration, see Figure 4, (a-c) for λ 0 = 1, and (d-f) for λ 0 = 10, respectively. As is noticed before, the single-particle excitation has a two-fold degeneracy for excitation at momentum p = 2π/L(1, 3) due to the boundary condition of the momentum. For the degenerate states, we compare the state fidelity in the subspace of the degenerate states. From the simulation results, we can find that the eigenstates obtained from the varia- tional algorithms can be found with a high state fidelity, verifying the effectiveness of the variational algorithms. Figure 4, (c,f) shows the state fidelity in the oneparticle subspace. Here, the one-particle subspace is obtained by adiabatically evolving the one-particle state in the free theory. In the adiabatic evolution, we set the time step dt = 0.01 and total time T = 50 to ensure the state fidelity error below 10 −4 . The results indicate a high state overlap in this one-particle space in the presence of interaction λ 0 = 1 and λ 0 = 10, consistent with our analysis.
We then discuss the simulation in the spectral crowding regime with a relatively small static mass m 0 , where the many-particle state occupied at zero momentum p = 0 will have lower energies compared to the single-particle state. In this regime, we should note that the excitations in the interacting field will not be local and may not have a well-defined particle number as that in the free theory. Therefore, searching for the excited state could be difficult in general, even when we could have access to adiabatic evolution. In what follows, we will discuss the low-lying excited states for two static mass m 0 = 0.5 and m 0 = 0.37, and show the search for eigenstates using the variational algorithms.
(58) Compared to the case of m 0 = 1.27, the two-particle states have lower energies than the single-particle state, and the state |0, 0, 1, 0 is indeed a highly excited state (higher than three-particle states). The excited states of the interacting Hamiltonian with λ 0 = 1 follows the same order as that of the free Hamiltonian. We compare the simulation results for λ 0 = 1 using the variational methods and adiabatic evolution in Figure 5. Figure 5, (a,b) and (c,d) shows the results using variational methods and adiabatic evolution, respectively. Overall, the low-lying eigenstates obtained from the variational algorithms can be found with a high state fidelity.
Then, we use the variational quantum algorithms to  search for the low-lying excited states. We show the simulation results with variational methods and adiabatic evolution for λ 0 = 10 in Figure 6. Figure 6 shows that we could still use the excited state of the free Hamiltonian as the initial guess and obtain the target state with relatively high fidelity.
In the interacting theory, the single-particle excitation may not be well defined, and the eigenstate of the free theory may not be adequate for finding the ground state in the interacting theory, especially for the large interacting field λ 0 = 10. Note that the problem for the choice of the initial state exists in the adiabatic evolution.
(60) We first show the convergence towards iteration in Figure 7 for λ 0 = 1 for both the variational state preparation and adiabatic state preparation.
Similar to the case of m 0 = 0.5, the two-particle excitation for m 0 = 0.37 has a small energy in the free theory, but it has a much higher energy compared to all the single-particle excitations and even higher than the three-particle excitation in the strong coupling regime. For instance, for a large interaction strength λ 0 = 10, the energy has the following relation E (|3, 0, 0, 0 ) < E (|0, 0, 1, 0 ) < E (|2, 0, 0, 0 ) , (61) which indicates the energy single-particle excitation is between the multi-particle excitation in the interacting field. However, in the interacting field, we can find that the single-particle state is actually much close to the excited states in terms of state fidelity. Therefore, we similarly choose the initial states as the single-particle excited states. We show the convergence towards iteration in Figure 8 for λ 0 = 10. As shown in Figure 8 (b,d), we can find that the state fidelity of the second excited state and the third excited state are relatively lower than the others. This is because these two states are actually evolved from the two degenerate states due to the boundary condition in the momentum space. However, the state fidelity in the subspace spanned by the degenerate states is numerically tested to be over 99%.
Finally, we show a detailed resource analysis of our method. In the simulation, we can reduce the number of  Pauli operators by restrict the higher order transitions.
Here, we consider to fix the energy cutoff as n cut = 4 and we restrict the higher order transition to be less than 3 in T 1 and the modifiedT 2 operator. Then the energy constraint |s i − t i | ≤ 4 is trivially hold and thus we need 6 Pauli operators to construct a single term fromT 1 (they are σ x σ x , σ y σ y , σ x σ z , σ x I, σ z σ x , Iσ x ). A linear combinations of tensor products of all pairs of these Pauli operators to form a single term from Eq. (22). In practice, we may only use the second term in Eq. (22), which saves half resources. Therefore, the number of Pauli rotation operators required inT 2 is upper bounded by 18N . Figure. 9 shows the number of Pauli rotation operators in the N -qubit variational circuitT 1 +T 2 .

V. OUTLOOKS
In this paper, we discuss constructions of a variational version of the Jordan-Lee-Preskill algorithm in the nearterm quantum computer. We justify the validity of the algorithm by several numerical simulations. We believe that our hybrid quantum-classical algorithm will eventually benefit possible solutions to open problems in fundamental physics, and benchmark tasks of near-term quantum devices. Here, we summarize some potential research directions along our path.

A. Relation to the physical observables
In the previous discussions, we demonstrate the numerical simulation with fixed lattice spacing and lattice sites. To obtain the expectation value of physical observables in the real scalar field, say φ , we can first measure the expectation value with a series of increasing number of sites N , φ (0) a,N , and extrapolate to the infinite volume limit N → ∞, φ (0) a,N →∞ . Then, we further extrapolate these results to the continuum limit a→0,N →∞ . Finally, we renormalize the expectation value as φ (R) = Z φ (0) with the renormalization constant Z. We leave the discussions to future work.

B. Simulating quantum field theories in the NISQ era
Our work opens up a new direction of formulating several quantum field theory tasks in the setup of variational quantum simulation. In the era of noisy intermediatescale quantum (NISQ), we expect that hybrid, variational quantum simulation algorithms might be one of the most accessible ways regarding near-term quantum hardware.
There is a landscape of quantum field theories with a full basket of open problems, who are looking for the potential computational capacity of quantum devices. Our work about strongly-coupled λφ 4 theory is one of the simplest examples, whose non-perturbative nature is not fully understood by quantum field theorists. One could consider generalizing the scattering paradigm and its relevant techniques towards other quantum field theories. Specifically, lattice gauge theories in the four dimensions are particularly important for particle physicists, since it is related to quantum chromodynamics (QCD) and the Standard Model in particle physics. We refer to Refs. [22,23,25,26,[49][50][51][52] for recent theoretical and experiment advances in the quantum simulation of lattice gauge field theories. One could look for other strongly-coupled quantum field theory problems, for instance, phase transitions in the finite-temperature quantum field theories that are closely related to nuclear physics [4,53].

C. Identifying possible quantum advantages
Variational algorithms running on near-term devices might have further advantages for fundamental studies in quantum information science. Specifically, since we could design hybrid quantum-classical algorithms, it is easy for us to diagnose which classical or quantum steps have advantages practically. Although in this work, we do not focus on this comparison, we expect that similar comparisons could be performed in future studies. In the future, people might work out practically, which steps in the whole algorithms might have the quantum advantage, and if so, how much advantage they will have. Those studies might be helpful to construct the most useful quantum algorithms using practical experiences, and use those experiences to benchmark near-term devices. Quantum simulation of quantum field theories is a field that is still young, but we expect that finally, more techniques and hardcore developments could be formulated (see some similar analysis in computational quantum chemistry [54]).
VI. ACKNOWLEDGEMENTS so that we can express the two-point function with time-ordered a free part and a spectral function ρ(M 2 ): where now we integrate over d + 1 energy-momentum and the spectral function defined: The spectrum function leads to isolated peak of Driac Delta functions for the single-particle state where M 2 is less than 4m 2 or m 2 bound at the rest frame |p 0 ⟩: The field-strength renormalization is given by: where the last equality follows the fact that the spectrum function is Lorentz invariant. One can prove that Z 3 ≤ 1 and equals one if and only if it is in the free theory. As a result, to even write down the two-point function, we require the field to be renormalized: In the latter part, we always assume ϕ(x) is properly renormalized. In simulating scattering amplitude of the ϕ 4 theory in quantum circuits, one should, therefore, always first measure the two-point function in term of the self-energy to obtain the proper field-strength renormalization factor [? ]. This is an important caveat to measuring the scattering amplitude in quantum circuits using the Harmonic oscillator basis. In scattering theory, the use of harmonic oscillator basis is to prepare in-and-out states (particularly in-states). It is conventional to assume that the in-states start far past as a wave packet and then evolve in time as free fields to interaction, as a consequence of the locality of interaction. Therefore, the use of harmonic oscillator in quantum circuits is to encode the far past in-states, where we now use the Heisenberg picture to denote the fields ϕ H in (t, x) (ϕ H out (t, x)). Then we can express the creation and annihilation operators at some arbitrary time instance t 0 : (I.12) Note that at t 0 → −∞, the above equation simply implies the mode expansion of the scalar fields in the free theory where we can express the scattering amplitudes: ⟨q 1 , · · · , q m ; out |p 1 , p 2 , · · · , p n , in⟩ = ⟨Ω| a q1 (t f ) · · · a qm (t f )a † p1 (t i ) · · · a † pn (t i ) |Ω⟩ . (I.13) The variational ansätze constructed in this paper can simulate initial states |i⟩ specified momentum sector as well as a wavepacket: The specific form of the wave-packets should not concern us here; however, it is an independent issue to analyze the encoding of F (k) in the quantum circuits if we were to prepare a wave-packets in the Harmonic oscillator basis. We cam analogously define the wave-packet creation operator: Notice that the wave-packets must preserve the inner products; thus, in a quantum circuits with harmonic oscillator basis, one can write: For some unitary matrix V also defined on the harmonic oscillator basis. Note that: where the LHS can be directly applied to Eqn.(I.13) where ⟨Ω| a † (∞) = 0 and similarly a(−∞) |Ω⟩ = 0. The RHS can be evaluated from the wave-packets Eqn.(I.15) and Eqn.(I.12): where we use the convention dtdx = d 2 x and the integration by parts in: Note that if we are describing a free (on-shell) theory, (∂ 2 + m 2 0 )ϕ = 0 implies there is no quantization, i.e., no particles being created or annihilated in the process. In the interaction, the on-shell condition is replaced by some external source j(x) such that: (I.20) Analogously, let assume that the final state is described by some wave-packets: We would can derive the LSZ formula by putting everything together: Therefore, we can express the more general version of the scattering amplitudes from Eq.(I.13): ⟨q 1 , · · · , q m ; out |p 1 , p 2 , · · · , p n , in⟩ = (i) m+n d d+1 y 1 · · · d d+1 y m d d+1 x 1 · · · d d+1 x n e iq1y1 · · · e iqmym e −ip1x1 e −ipnxn (∂ 2 y1 + m 2 ) · · · (∂ 2 ym + m 2 )(∂ 2 x1 + m 2 ) · · · (∂ 2 xn + m 2 ) ⟨Ω| T ϕ(y 1 ) · · · ϕ(y m )ϕ(x 1 ) · · · ϕ(x n ) |Ω⟩ = (i) m+n (q 2 1 + m 2 ) · · · (q 2 m + m 2 )(p 2 1 + m 2 ) · · · (p 2 n + m 2 ) ⟨Ω| Tφ(q 1 ) · · ·φ(q m )φ(p 1 ) · · ·φ(p n ) |Ω⟩ , (I. 23) where at the last step we assume the wave-packets to be the planar wave functions localized at δ d (q 1 ) · · · δ d (q m )δ D (p 1 ) · · · δ D (p n ).
The m + n-point function can be calculated from the partition function pertubatively from Feynman Diagrams. Combing with that, we can derive rigorously the Feynman rules. In particular, the LSZ reduction formula implies that it, to obtain the scattering amplitude, is sufficient to apply using Eq. (I.2): (I.24) Simulating the above evolution from free states would ensure to prepare the proper in states for the scattering amplitudes, where the factor (⟨Ω|0⟩) −1 normalizes the vacuum bubbles.
The UCC ansätze used in this work can be seen to create a wave-packet in interaction region with some fixed t. Recall Eq.(I. 16 Note thatT 4 (θ(t)) in the adiabatic picture can be used to simulate the full Hamiltonian (Eq. (21) in the main text); hence, in which case V is the identity operator. Therefore, Eq.(I.24) provides a simple interpretation of the use of variational UCC ansätze: to create a wave-packet at some intermediate time t, for which we can denote as initial time t i of the optimization process using variational quantum imaginary time evolution [? ] or other optimization method such as [? ]. The optimization would produce a final time t f (or τ f ≡ it f ), where we further evolve the wave-packet into lowest energy eigenstates where one can utilize physic prior to recognizing these states and compute corresponding scattering amplitudes.

II. TRUNCATION ERROR OF FIELD BASIS
We provide a rigorous analysis on the truncation error of free theory eigenstates under field basis and compare with the case of harmonic oscillator basis. Let |ψ⟩ be any state expanded under field basis as with truncated approximation Then the truncation error ϵ can be defined as For simplicity, we assume N = 1 which corresponds to the case of a free quantum harmonic oscillator. The eigenstates are given by Hermit-Gauss functions: where H n (x) are physicists' Hermit polynomials defined as The derivative is nonnegative in [0, 1 2 ] while negative in ( 1 2 , +∞). Accordingly, f is nonnegative in (0, 1]. We only need to make sure that f will not turn into negative later. Note that by L'Hôpital's rule, f (x) → 0 when x → +∞. If f (x) could be negative, it must admit a global minimum at x 0 ∈ ( 1 2 , +∞) and at which f "(x 0 ) = 0. However, this is impossible and hence the inequality holds on R.
Now the problem is analyzing how would the remaining term containing H 2j−1 (x) decay with respect to both x and n. As Eq. II.1 converts H n (x) 2 from the integrand into H 2j−1 (x), it allows us to perform Cramer inequality [? ]: |ψ 2j−1 (x)| ≤ π − 1 4 for all x ∈ R. With this fact we have Explicitly expansion shows that Therefore, For fixed energy level n, an exponentially small error can be obtained when x > O( √ n). In the field context, this x is though as ϕ max which is given as 2 nq (multiplying with a chosen increment δ ϕ ) on a quantum computer with n q qubits. Thus we can conclude that n q = O(log log 1 ϵ ) of the error for approximating low-lying excited states of free theory and the number of qubits. As a comparison, the original truncation error bound of Jordan-Lee-Preskill algorithm is obtained by Chebyshev's inequality. As they make no presupposition on the concerned wavefunction, Chebyshev's inequality only yields a polynomial decay on the error.
The scaling O( √ n) coincides with what indicates by WKB approximation: the classically forbidden region is given by the square root of energy. If one takes x ≈ O( √ n), the above inequality (II.2) should be further refined and numerical results from [? ] indicate that the error only decays polynomially for such cases.
When we change to harmonic oscillator basis, there is no such a truncation error on the number eigenstates because, as we explained in the main text, each of them is simply a single computational basis element |s (k) ⟩ and there is no need to talk about truncating Hermit-Gauss functions in this setting. This is another advantage to work directly with harmonic oscillator basis. However, the error still stems from elsewhere. For instance, the creation/annihilation operators do not satisfy the commutation relation [a † , a] = iI rigorously (this also happens for the field operator ϕ(x) and its conjugate momentum π(x)). Expanding the Lie bracket explicitly by matrices, one can check that this error ϵ ′ decays exponentially with the dimension of the local Hilbert space, and hence we have n q = log 1 ϵ ′ where n q is the number of qubits per site.

III. THE RATE OF CONVERGENCE OF VARIATIONAL ALGORITHMS
It is of interest to determine the rate of convergence of the variational quantum imaginary evolution [? ] the variational UCC bosonic ansätze, as it determines the number of measurements one needs to perform for the parameter updating until reaching the optimal parameters θ * . In particular, having an exponential rate of convergence would give an optimal scaling of both the subspace and state fidelity seen in Section IV.E. Analyzing the rate of convergence for heuristic optimization algorithms is currently an active research front [? ? ? ] and thus beyond the scope of this paper. We here give a only heuristic remark on convergence rate using quantum neural tangent kernel (QNTK).
The general strategy for the ground state searching is by updating the parameters as where θ µ (t) represents the optimization dynamics with step t, the learning rate is given by η µ (t), and ⟨H⟩ θ ≡ ⟨ψ(θ)| H |ψ(θ)⟩.
Here, we use A(θ(t)) to represent the metric matrix at the parameter θ(t). The metric matrix in the gradient descent algorithm is simply the identity matrix. In the following section, we will show its explicit form during the optimization. One can ask if there exists a regime where there is a convergence guarantee and, if so, the rate of convergence for these variational paramterization.
One can study this question from over-paramterization using quantum neural tangent kernel (QNTK) [? ]. Further taking ⟨H⟩ θ(t) ≡ z(θ(t)), Eq.(III.1) implies: Assuming A is the identity matrix and η to be parameterization-independent, the resultant is precisely the QNTK defined in [? ]. Note that we can interpret A −1 µν as the learning rate tensor as part of the definition of NTK in classical neural networks [? ]. In particular, for the circuits that form at least approximate 2-design that satisfies certain concentration conditions (See in [? ]), the average convergence is of the form ϵ(θ(t)) ≈ e −γt ϵ(θ(0)), where ϵ ≡ z(θ) − E 0 , E 0 is the ground state energy and With L being the total count of variational parameters and H be full Hamiltonian. The dimension of the Hilbert space in our case is n N cut . The exponential convergence rate is guarantee on average in the over-parametrization regime where L ≈ dim(H) 2 / tr(H 2 ). When A fails to be an identity matrix such as in the case of quantum imaginary time evolution used in the following, no precise analytical convergence guarantee is known. However, it seems to be reasonable to extrapolate the hypothesis that such methods, due to its more physical, geometric nature, would have convergence rates lowered-bounded by the naive gradient descent methods.
Clarifying the gradient descent dynamics would enable us to determine the query complexity of the quantum circuits during the optimization phase; that is the number of measurements the variational quantum circuits needed to perform or simply the number of parameter updates. Let us consider the k-particle subspace fidelity: f (k) (τ ) = ⟨ψ(θ(τ ))| Λ (k) |ψ(θ(τ ))⟩ , (III. 5) where the N -dimensional projection operator: where D Λ (k) is bounded by k+N −1 k . As an intuitive example, we consider the projection to Λ (0) = |Ω⟩ ⟨Ω| and f (0) (τ ). Suppose we initially at f (0) (τ i ) and we are aiming at fidelity f (0) (τ f ) = 1 − δ, an important question to ask is what is the query complexity (τ f − τ i )/∆τ ? We wish this value to be bounded by Poly(N log n cut ). Recall from quantum field theories we know that the overlap between physical vacuum state and the free vacuum state ⟨Ω| |0⟩ is non-vanishing even in the continuum limit. Therefore, having initializing the UCC ansätze |ψ(θ(τ i ))⟩ = c 0 |0⟩ + |⊥⟩ for some constant c 0 would ensure that the variational algorithms to converge to the physical vacuum state |Ω⟩. Note that the physical vacuum state is unique so we expect that we can approximate |Ω⟩ with fidelity f (0) (τ f ) = 1 − δ if ϵ(θ(τ f )) ≤ δ where we can further write from Eq. (III.3): for some error δ and initial error ϵ 0 . The log of inverse error looks promising but unfortunately as we computed 1/γ scales as the dimension of Hilbert space so not super ideal. If it 1/γ were to be bounded by Poly(N log n cut ), it then follows that the query complexity is polynomially bounded.
To obtain the scaling of γ we first assume that the variational ansätze generated byT 1 andT 2 (Eq.(25) & (26) in Main text) form an at least approximate 2-design and satisfies the concentration condition [? ] in the large n cut and N limit. Then For eacĥ T ℓ , we have O(N ℓ n ℓ cut ) variational parameters, corresponding to each Pauli operator. As a result, the rate of convergence can be computed: (III.8) Thus the over-parameterization regime occurs at: log k + log η + (3 + ℓ) log N − (2N − ℓ − 1) log n cut = 0, (III.9) where the trace norm of the Hamiltonian is bounded by its number count of Pauli operators O(N 3 n cut ) and k some constant. Note that this regime fails to hold at large N ; so the QNTK could only provide some analytic guarantee of exponential convergence rate at some small momentum modes where N scales linearly with the fixed ℓ. In this regime, we set the learning rate log kη ≈ ℓ log ncut log ℓ . As a result, in general, it is difficult from QNTK alone to analyze the rate of convergence in the ansätze for a fixed and lower ℓ as the constant kη has to be linearly increasing with respect to n cut . For NISQ applications, where the number of qubits is within tens of qubits, however, it might be possible to achieve an exponential rate of convergence. For instance, if we set the constant k = 10000 and the learning rate η = 0.01 with N = 4, then we have roughly n cut = 32. ThusT 4 can achieve an exponential rate of convergence, provided the aforementioned assumptions, up to the 20-qubit system.
It should be noted that we expect that NISQ era to have roughly between a few tens to few hundred qubits, clearly above the regime where the rudimentary QNTK theory can provide a theoretical guarantee on logarithms or polynomial query complexity. To this end, to obtain any polynomial query complexity in the variational setting is, to our best knowledge, difficult and largely unexplored; as a result, we leave it to as an open remark. The above analysis could also be understood as a motivation for the significance on researching the direction of QNTK or similar methods, where it could enable us to obtain a rigorous algorithmic prediction instead of relying on heuristic arguments.

IV. ϕ 4 THEORY AS AN EFFECTIVE FIELD THEORY
In this appendix, we simply review how the effective description of the classical Ising model in the long distance could be given by field theories, for the convenience of the quantum information science community. The observation where field theory could serve as an effective description of many-body physics has an important significance in theoretical physics, sharpening the importance of simulating quantum field theories.
The ϕ 4 field is one of the simplest interacting field theories related to the emergent behavior in the physical systems. In this section, we will review that the d-dimensional classical Ising model can be described by the ϕ 4 field at the long distance, following the discussions in [? ] rather closely.
The classical magnetism can be described by the Hamiltonian where J ij is a translationally invariant correlation matrix describing the interaction of the spins. The classical partition function can be written as where we defineJ ≡ −βJ andh i ≡ βh i . By introducing a continuous auxiliary field ψ i , the partition function is given by where the Gaussian integral has been evalued as The equation above is under the saddle point approximation, which we preserves the second order of the functional Note that the interaction between spins is decoupled by the continuous auxiliary field. In this case, we could work with the multi-dimensional integral Dψ instead of the summation over discretized spins where we have used S i = ±1, and the normalization factor has absorbed the trivial prefactor. Notice that ln cosh(x) = 1 2 x 2 − 1 12 So the Ising model could be described by the ϕ 4 theory at the long distance.
To make the fourth-order term more clear, we change the integration variable ψ to the field variable via and we have Z = N Dϕe − ij ϕiJij ϕj + i ϕihi+ i ln cosh(2 jJ ij ϕj ) .
(IV. 8) We note that the amplitude of the field ϕ indicates a measure of magnetism. If the coupling is transnational invariant, which indicates a spatially smooth field, we can capture this feature by the Fourier transformation: which can be determined from the interaction coupling. We can thus see its correspondence to (??), the ϕ 4 field in the momentum space.
In the real space, the effective action admits the form of ϕ 4 field as When the coupled external field is zero, i.e., h = 0, it becomes the prototypical ϕ 4 field discussed in the main text. Suppose c 2 > 0, for c 1 > 0, the action has a global minimum at ϕ = 0, and thus the ground state is in the the paramagnetic phase. While for c 1 < 0, the action has two degenerate minima. Thus the phase is spontaneously broken with a ferromagnetic ordering. If the interaction strength J ij < 0, the constant c 1 is positive at sufficiently low temperature, and it becomes zero at the transition temperature. Therefore, the ϕ 4 field could capture the critical behavior of the paramagnetic phase to the ferromagnetic phase.