# 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, explicitly designed for quantum simulation of quantum field theory. We develop hybrid quantum-classical algorithms for crucial ingredients in particle scattering experiments, including encoding, state preparation, and time evolution, with several numerical simulations to demonstrate our algorithms in the 1+1 dimensional λϕ^4 quantum field theory. These algorithms could be understood as near-term analogs of the Jordan-Lee-Preskill algorithm, the basic algorithm for simulating quantum field theory using universal quantum devices. Our contribution also includes a bosonic version of the Unitary Coupled Cluster ansatz with physical interpretation in quantum field theory, a discussion about the subspace fidelity, a comparison among different bases in the 1+1 dimensional λϕ^4 theory, and the "spectral crowding" in the quantum field theory simulation.

## Authors

• 11 publications
• 1 publication
• 3 publications
09/16/2019

### Near-term quantum algorithms for linear systems of equations

Solving linear systems of equations is an essential component in science...
09/22/2020

### Control dynamics using quantum memory

We propose a new quantum numerical scheme to control the dynamics of a q...
03/16/2021

### Variational Quantum Algorithms for Euclidean Discrepancy and Covariate-Balancing

Algorithmic discrepancy theory seeks efficient algorithms to find those ...
10/03/2018

### Simulating the weak death of the neutron in a femtoscale universe with near-Exascale computing

The fundamental particle theory called Quantum Chromodynamics (QCD) dict...
06/17/2019

### Time-dependent Hamiltonian simulation with L^1-norm scaling

The difficulty of simulating quantum dynamics depends on the norm of the...
12/10/2020

### Variational Quantum Algorithms for Trace Distance and Fidelity Estimation

Estimating the difference between quantum data is crucial in quantum com...
02/08/2021

### Long-time simulations with high fidelity on quantum hardware

Moderate-size quantum computers are now publicly accessible over the clo...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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

preskill2018quantum; arute2019quantum. 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, 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 feynman1982simulating, and becomes an important new research direction recently in the physics community, see the references Preskill:2018fag; cyber; Jordan:2011ne; Jordan:2011ci 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 Milsted:2020jmf. 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 Jordan:2011ne; Jordan:2011ci in the 1+1 dimensional 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. Jordan:2011ne; Jordan:2011ci, 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 near-term 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. Peruzzo2014; farhi2014quantum; McClean2016; li2017efficient; yuan2019theory; sugurureview21; cerezo_variational_2021; zhang2020low; Endo20Variational; XU2021). 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 quantum field theory. We summarize our contribution as follows.

• Basis choice. We will make comments on the comparisons between the field basis and the harmonic oscillator basis, momentum space, and coordinate space. All those choices have pros and cons. Coordinate space setup will make the Lagrangian density local, but it is not directly connected to the Feynmann rules and scattering calculations in the momentum space. The momentum space setup is not local but still sparse, leaving a door open for variational (and also digital) simulations. 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 field basis will cause the field correlations to be easy to measure, but finding the eigenstates (for instance, the vacuum and low-lying one-particle states) might be not easy (which requires computing correlation functions at hand and non-trivial digital quantum algorithms for encoding).

• Initial state preparation. In order to prepare the initial wave packets, the original Jordan-Lee-Preskill algorithm uses adiabatic state preparation to turn on the coupling. In the variational setup, alternative strategies could be directly used and solve the initial scattering directly. In this work, we will discuss possible variational algorithms in detail, especially on the treatment of the momentum sectors. Specifically, we perform numerical simulations using the imaginary time evolution mcardle2019variational and a bosonic version of the unitary coupled cluster (UCC) ansatz with physical interpretation in quantum field theory C9SC01313J; O_Malley_2016; Shen_2017 to test our simulation algorithms. We also theoretically and numerically investigate the spectral crowding phenomena in the quantum field theories in both weakly-coupled and strongly-coupled theories.

• Real-time evolution and scattering. Besides digital quantum simulation algorithms (see Refs. tro; low2016hamiltonian; shaw2020quantum; chakraborty2020digital; bender2018digital), the real-time evolution algorithms could also be tracked by variational methods, see Refs. li2017efficient; yuan2019theory; kokail2019self; Paulson_2021. In this work, we present the formulation of the corresponding algorithms, and discuss how to guarantee a certain simulation accuracy. We also comment on the particle excitations, circuit ansatz, and the challenges.

• Simulation fidelity. Simulation errors in the variational simulation will not only be limited to the digital simulation error (like the Trotter error) but also the variational error from the variational 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 (where we could define particle subspace fidelity to quantify it), thus suitable for noisy intermediate-scale quantum (NISQ) devices. We numerically verify the fidelity in the one-particle subspace in the interacting theory. The real-time evolution task might require higher precision to make sure consistent scattering results, and we analyze the simulation error during the scattering process.

The paper is organized as the following. In Section II, we discuss the theoretical formulation of the theory and the Jordan-Lee-Preskill algorithm, commenting on its variational extension. In Section III, we discuss different bases for state encoding and their comparisons. In Section IV, we discuss the variational initial state preparation in the Jordan-Lee-Preskill algorithm. In Section V, we discuss variational quantum simulation of the scattering dynamics, and analyze the simulation error during time evolution. We comment on the possibility of realizing the particle scattering process. In Section VI, we show the numerical simulation for the ground state and excited states preparation using variational quantum algorithms, and compare it with adiabatic evolution. We discuss the simulation results for different strength of the interacting field. Finally, in Section VII, we highlight a list of future directions.

## Ii λϕ4 theory and the Jordan-Lee-Preskill algorithm

### ii.1 Basics

In this section, we give a review of the theory in 1+1 dimensions. In this theory, we have a scalar quantum field with the Hamiltonian

 H=∫dx(12π2+12(∂xϕ)2+12m20ϕ2+λ04!ϕ4),

and the Lagrangian

 L=∫dx(12(∂tϕ)2−12(∂xϕ)2−12m20ϕ2−λ04!ϕ4).

Moreover, we discretize it in the lattice,

 H=∑x∈Ωa[12π2+12(∇aϕ)2+12m20ϕ2+λ04!ϕ4],

with the Lagrangian

 L=∑x∈Ωa[12(∂tϕ)2−12(∇aϕ)2−12m20ϕ2−λ04!ϕ4].

The theory is defined on a lattice with total length and lattice spacing as , and , and we could also define its dual lattice We use to denote the sites in the lattice. The field momentum is defined as the Fourier conjugate of the field ,

 lattice: [ϕ(x),π(y)]=ia−1δx,y, continuum: [ϕ(x),π(y)]=iδ(x−y).

The discretized version of the derivative is given by where is the (bare) mass term in the free theory, and the term represents the coupling. When , 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,

 ϕ(x)=∫dp2π√12ωp(ap+a†−p)eipx, π(x)=−i∫dp2π√ωp2(ap−a†−p)eipx, (1)

and in the discrete system we have,

 ϕ(x)=∑p∈Γ1Leipx√12ω(p)(ap+a†−p), π(x)=−i∑p∈Γ1Leipx√ω(p)2(ap−a†−p), (2)

where the canonical algebra of and leads to the commutation relation as for lattice and for continuum, respectively. The energy dispersion is given by

 ω(p)=√m20+4a2sin2(ap2)a→0−−→ωp≡√m20+p2. (3)

In such a basis, the Hamiltonian is diagonalized as

 H0=∑p∈Γ1Lω(p)a†pap+E0a→0−−→∫dp2πωpa†pap+E0, (4)

where

 E0=∑p∈Γ12ω(p)a→0−−→∫dp2π12ωp×2πδ(0). (5)

An important physical observable we could measure, is the (Wightman) two-point function as where is the ground state of the theory. In the free theory case where , one can compute the two point function explicitly

 G0(x−y)=∑p∈Γ1L12ω(p)eip(x−y)a→0−−→∫dp2π12ωpeip(x−y). (6)

The two-point function of the scalar defines the scalar mass of the theory. In the continuum limit, when we turn on the interaction

, in the weakly-coupled regime one could compute the correction to the mass through the Feynman diagrams. If the coupling is strong enough, studying the theory becomes challenging with only analytic perturbative techniques. The theory will experience a second-order phase transition at the strong coupling, where the universal behavior belongs to the 2D Ising universality class. In general, computing the two-point function will tell us the information about the mass of the particle in the lattice discretization. In the Appendix, we address the relation between lattice models and their field theory description, emphasizing the importance of simulating quantum field theories.

### ii.2 Jordan-Lee-Preskill algorithm and extension

In this section, we discuss the Jordan-Lee-Preskill algorithm Jordan:2011ne; Jordan:2011ci. The algorithm requires the Hamiltonian of the field theory, say the theory in the lattice, and the input state data specified by some wave packets.

###### Algorithm 1 (Jordan-Lee-Preskill Jordan:2011ne; Jordan:2011ci).

The algorithm consists the following smaller steps.

• State encoding. We need to formulate a qubit system in a device with a proper number of qubits. In the Jordan-Lee-Preskill algorithm, we use the field basis in the coordinate space. The number of required qubits is bounded by the scattering energy we need.

• Preparing initial states in the free theory. This could be done by the Kitaev-Webb algorithm in the quantum computer kitaev2008wavefunction. In the scattering experiment, we need some single-particle wave packet states with given momenta as initial states, which could be prepared from the creation operator acting on the Gaussian vacuum.

• Adiabatic state preparation towards the interacting theory. This treatment is to convert the wave packets of single particles in the free theory towards reasonable single-particle wave packets in the interacting theory. To achieve this goal, one could use the adiabatic state preparation to slowly turn on the interaction, which is closely related to the mass gap in the theory.

• Time evolution of the Hamiltonian. The evolution in a universal quantum computer is using the product formula (Lie-Trotter-Suzuki formula). Note that the same trick is also used in the adiabatic state preparation procedure.

• Measurement. One could measure quantum observables in the system using quantum algorithms (see discussions in Ref. Jordan:2011ne; Jordan:2011ci).

In this work, we discuss how to perform the above analysis in the variational platform. Before a more detailed discussion, we will make a brief summary on the variational realization of the Jordan-Lee-Preskill algorithm.

• State encoding. There are many choices of bases for bosonic lattice quantum field theory: coordinate space or momentum space, field basis or harmonic oscillator basis. The choice of bases is depending on the purpose of our simulation and the hardware we are going to use. One of the main advantages of the coordinate space is the locality in the action, which will be ruined when we move to the momentum space. However, the Hamiltonian is still sparse in the momentum space, and the variational algorithms may not be sensitive to the locality at all, since the computation is completely done by introducing variational ansatz. Thus, all four combinations of basis choices could be considered, which might provide different levels of convenience in the ground state, wave packet state construction, and measurement, where we will discuss in more detail in later sections, especially Section III (see other formulations of quantum simulation of quantum field theories in the matrix models, especially in the bosonic form Liu:2020eoa; Gharibyan:2020bab; Buser:2020cvn; Rinaldi:2021jbg; Kreshchuk:2020dla; Kreshchuk:2020aiq; Kreshchuk:2020kcz).

• Preparing initial states for the scattering problem. Different choices of bases will lead to different methods for state preparation. In the field basis, one could construct the ground state using the Kitaev-Webb algorithm, which unfortunately is not suitable for the near-term. One could also use variational forms to find the ground state as a Gaussian distribution (see a related work

Macridin:2018oli, and an earlier discussion in Klco:2018zqz). Moreover, the ground state and the single-particle states in the free theory are defined naturally in the harmonic oscillator basis in the momentum space. Even if we switch to the harmonic oscillator basis in the coordinate space, the true vacuum state will be slightly different. We will address the above issue in the later sections, especially in Section IV.

• Time evolution of the scattering states. The real time evolution of the algorithm could be done by the variational algorithm in li2017efficient; yuan2019theory. 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 zhang2020low. We provide a theoretical framework for the dynamics simulation and make comments on the challenges of the scattering process in Section V.

## Iii State encoding

At the starting point, we wish to encode our Hamiltonian from quantum field theory to a quantum device. Thus, we have to choose some computational bases. This section is devoted to having a detailed review and comparison on different versions of bases (see some simular discussions in the quantum chemistry context sawaya2020resource).

### iii.1 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

 |ψ⟩=∫∞−∞dϕ1⋯∫∞−∞dϕNψ(ϕ1,…,ϕN)|ϕ1,…,ϕN⟩. (7)

Here, is the total number of sites. The state decomposition for an arbitrary state gives the above wavefunction , where could choose from an arbitrary number in in principle. The basis is constructed such that a local field operator

hits on the basis will return the real number as its eigenvalue. 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. The states that we are interested in is truncated as

 |ψcut⟩= ∫ϕmax−ϕmaxdϕ1⋯∫ϕmax−ϕmaxdϕNψ(ϕ1,…,ϕN)|ϕ1,…ϕN⟩. (8)

Moreover, we want an increment in discretization such that we don’t need to choose variables in a continuous interval. Thus, the total number of qubits we need to encode is Assuming the cutoff error such that the state fidelity satisfies , it is shown in Refs. Jordan:2011ne; Jordan:2011ci that there are bounds on from the scattering energy , and we have

 ϕmax=O⎛⎝√NEam20ϵcut⎞⎠,    δϕ=O(√ϵcutaNE), nb=O(log(NEm0ϵcut)), (9)

where we demand

The bound is derived from the Chebyshev inequality relating the probability outside the range of the field variable cutoff to the field fluctuations, which are bounded by energy based on the expression of the Hamiltonian. The bound is useful to prove the polynomial complexity of the Jordan-Lee-Preskill algorithm, but one should note that the bound itself may not be tight in the practical calculations.

### iii.2 The harmonic oscillator 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. Since we are dealing with a bosonic field theory, the local degrees of freedom are given by harmonic oscillators.

We firstly consider the following transformation,

 ϕ(x)=1(2mx)1/2(ax+a†x), π(x)=−i(mx2)1/2(ax−a†x). (10)

Here is a free parameter we could choose. Then the canonical algebra of and leads to the commutation relation

Now, the creation operator and its conjugate could define the number states at the site . Thus, in the harmonic oscillator basis, we could define the state decomposition

 |ψ⟩=∞∑n1=0…∞∑nN=0ψ(n1,n2,…,nN)|n1,n2,…,nN⟩. (11)

Here, the number states are given by creation operators hitting on the vacua for the location . Let’s say that we are mostly interested in the maximal energy level , so we could instead cut the Hilbert space and define

 |ψcut⟩=ncut∑n1=0…ncut∑nN=0ψ(n1,n2,…,nN)|n1,n2,…,nN⟩. (12)

### iii.3 The field basis in the momentum space

Now we introduce a basis that is dual to the field basis: the field basis in the momentum space. Remember that we define the dual lattice based on the real lattice . Thus, one could directly write the Hamiltonian in terms of the momentum coordinate. To be more specific, consider the free theory solution

 ϕ(t,x) =∫dk(2π)√2ωk(ak+a†−k)e−iωkt+ikx ≡∫dk2πϕk(t)eikx, (13)

with the lattice version

 ϕ(t,x) =1L∑k∈Γ1√2ωk(ak+a†−k)e−iω(k)t+ikx ≡1L∑k∈Γϕk(t)eikx. (14)

We could write the system in terms of the modes . The free theory part is given by

 L0=∫dk(2π)[12(∂tϕk)(∂tϕ−k)−12ω2kϕkϕ−k], H0=∫dk(2π)[12πkπ−k+12ω2kϕkϕ−k], (15)

and in the lattice we have,

 L0=1L∑k∈Γ[12(∂tϕk)(∂tϕ−k)−12ω(k)2ϕkϕ−k], H0=1L∑k∈Γ[12πkπ−k+12ω(k)2ϕkϕ−k], (16)

where we similarly define . Now the commutation relation is

 discrete: [ϕk(t),πp(t)]=iLδk,−p, continuum: [ϕk(t),πp(t)]=2πiδ(k+p),

in the continuum. We could also write the interaction piece of the Hamiltonian in the momentum space by

 discrete: λ04!∑x∈Ωϕ4= λ04!1L3∑k1,k2,k3∈Γϕk1ϕk2ϕk3ϕ−k1−k2−k3, continuum: λ04!∫dxϕ4= λ04!∫dk1dk2dk3(2π)3ϕk1ϕk2ϕk3ϕ−k1−k2−k3. (17)

Thus, we could encode the system in the following basis

 |ψ⟩=∫∏pi∈Γdϕpiψ(ϕp1,ϕp2,…)∣∣ϕp1,ϕp2,…⟩, (18)

and we could similarly make a truncation on the field range in the momentum space.

### iii.4 The harmonic oscillator basis in the momentum space

Similarly, we could consider the harmonic oscillator basis in the momentum space. In this basis, the state decomposition is

 |ψ⟩=∑pi∈Γ,npi∈[0,ncut]ψ(np1,np2,…)∣∣np1,np2,…⟩. (19)

The number states now is generated by the creation operator specified by

 ϕk=1√2ωk(ak+a†−k), (20)

and the dual field momentum. 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 VI, we discuss the particle excitations in the momentum space in more details.

### iii.5 A comparison

Comparisons on the various aspects of the encoding schemes need a detailed analysis.

The momentum space is widely used in particle physics. For instance, it is useful to write the Feynman rules and Feynman diagrams in the momentum space. In the momentum space picture, it is also natural to interpret particles as momentum eigenstates in the theory, and the S-matrix should have a clear definition in the momentum space. However, unfortunately, in the momentum space, the Hamiltonians are usually not local, even if we are studying the free theory. So in order to make the Hamiltonian in a local way, we wish to simulate the Hamiltonian in the coordinate space when we are computing the time evolution dynamics.

The field basis is useful for estimating field correlation functions. For instance, if we wish to read the field distribution for a given state, it is straightforward to use the field basis. However, if we wish to see the state in the space of energies and particle numbers, we might consider using the harmonic oscillator, although those quantities are defined in the sense of the free theory. S-matrix might also have a better interpretation under the harmonic oscillator basis. The paper

Klco:2018zqz makes a detailed discussion and numerical comparison between those two bases.

Here, we will give a summary of the number of terms in the Hamiltonian for each basis.

• The field basis in the coordinate space. If we only count the number of terms evaluated purely in the field and field momentum operators, the total number of terms is .

• The harmonic oscillator basis in the coordinate space. As we have described before, we will replace the field and field momentum operators in the position space by introducing local harmonic oscillators as the computational basis. Since this happens completely in the position space, the number of total terms in the Hamiltonian is , where each local term is made by local creation and annihilation operators and their nearest neighbour products.

• The field basis in the momentum space. As we have discussed before, the Hamiltonian in the momentum space is written in the non-local form. For the free theory Hamiltonian, the number of terms will scale as when we count the field and field momentum operators. When we add interactions, the number of terms will scale as .

• The harmonic oscillator basis in the momentum space. Similarly, we get terms in the free theory, and terms in the interacting theory.

Moreover, the problem of choosing the basis is closely related to the initial state construction before scattering. Going back to the Jordan-Lee-Preskill situation, in the field basis, the ground state is Gaussian for the free theory, so we could construct the Gaussian state with the help of the Kitaev-Webb algorithm as mentioned before. Then, adiabatic state preparation is used to simulate the single-particle wave packet in the interacting theory. Note that the wave packets centered at a given momentum are in fact defined by the above way.

Back to the harmonic oscillator basis, it will be very convenient to work in the momentum space, since the momentum sectors are naturally defined in the free theory by the harmonic oscillator. The 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. 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.

In summary, one of the main advantages of the momentum space, especially for the harmonic oscillator basis, is to quickly identify the momenta and the particle number of output states during the simulation. Thus, in this paper, we will mainly work on the harmonic oscillator basis in the variational setup.

## Iv Variational state preparation

### iv.1 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

 |ψ(θ)⟩=(∏Lℓ=1Uℓ(θℓ))|ψ0⟩. (21)

Here s are some unitary operators that could be realized in the quantum device, for instance, with the variational parameters . s are some Hermitian operators, for instance, elements in the Pauli group, and 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 .

The general strategy for the ground state searching is by updating the parameters as

 θℓ(t+1)=θℓ(t)−λℓ(t)A−1(θ(t))∂∂θℓ⟨H⟩θ(t), (22)

where represents the optimization dynamics with step , and the learning rate is given by . Here, we use to represent the metric matrix at the parameter

. 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.

The next question is how to choose the initial state and s? The precise strategy of choosing , and the optimization scheme will specify the variational quantum algorithm we use. There are many variational algorithms (see Ref. sugurureview21; cerezo_variational_2021 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.

### iv.2 Bosonic UCC ansatz

As is mentioned before, the variational algorithm may not be very sensitive to the locality of the Hamiltonian. Thus we will focus on the harmonic oscillator basis in the momentum space.

For a single harmonic oscillator, the space of the truncated eigenstates with energy levels can be encoded by a compact mapping using qubits. The creation operator and annihilation operator can be represented as

 (23)

The ladder operators can be decomposed into Pauli operators with terms. Thus, we can encode the Hamiltonian with a quantum computer 111The Hamiltonian with truncated energy levels is represented in a low-energy subspace.. 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 O_Malley_2016; Shen_2017. Other prior works C9SC01313J; ollitrault2020hardware discussed the usage of bosonic UCC in studying vibronic properties of molecules.

The general form of unitary coupled cluster is given by

 |ψ(θ)⟩=exp(^T−^T†)|ψ0⟩, (24)

where is the sum of symmetry preserved excitation operators truncated at finite excitations as 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 C9SC01313J; ollitrault2020hardware could be natural to capture types and particle numbers for scalar particles. In the coordinate space, the single and double excitation operators preserving the excitation mode can be expressed as

 ^T1= L∑l=1ncut∑s,t=0θs(l),t(l)|s(l)⟩⟨t(l)|, ^T2= L∑l1,l2n% cut∑s1,s2,t1,t2=0θs(l1)1,t(l1)1,s(l2)2,t(l2)2|s(l1)1s(l2)2⟩⟨t(l1)1t(l2)2|, (25)

where , in the superscript denote the site in the coordinate space, i.e. the mode of the harmonic oscillator basis, and denote the energy level of the th site, and is the variational parameter in terms of certain energy level and site. Here, and represent the creation and annihilation of the particle excitation carrying energy and in the th harmonic oscillator mode, respectively. The field could lift the excitation up to levels at each mode so we could constrain the energy level as . Moreover, since we are interested in the low-energy eigenstates, the excitations to high energy levels can be truncated to and by considering the dominant low-energy excitations. We note that this ansatz is especially useful for lattice model with strong static correlations. For relatively small , the discretized theory in the coordinate space have strong locality. Therefore, we expect this ansatz would be suitable for this problem.

In the momentum space, the effective action preserves the momentum reflection symmetry (). Therefore, we may express the single and double excitation operators of the bosonic UCC ansatz in the momentum space as

 ^T1= ∑kncut∑s,t=0θs(k)t(k)|s(k)⟩⟨t(k)|, ^T2= ∑k1,k2ncut∑s1,s2,t1,t2=0θs(k1)1,t(k1)1,s(k2)2,t(k2)2× (|s(k1)1s(k2)2⟩⟨t(k1)1t(k2)2|+|s(~k1)1s(~k2)2⟩⟨t(~k1)1t(~k2)2|). (26)

Here, we use the same notation as in Eq. (25) and () is the (reversed) momentum. We can similarly have the constraint for the field in the momentum space.

This type of ansatz could be useful for the quantum field theory since the four-point coupling could induce the

scattering process, which represents two excitations coming in and two excitations moving out in the perturbation theory and the Feynman diagram (note that this argument is relatively heuristic, since the expansion from the exponential will induce more complicated terms).

Other heuristic ansatz that takes account of the symmetries of the field could be designed. For instance, one may construct the ansatz in the momentum space as

 ^T1= (27) ^T2= ∑kncut∑s1,s2,t1,t2=0θs(k)1,t(k)1,s(~k)2,t(~k)2× (|s(k)1s(~k)2⟩⟨t(k)1t(~k)2|+|s(~k)1s(k)2⟩⟨t(~k)1tk)2|), (28)

which considers the pairing correlations of the momentum and . Compared to ansatz in Eq. (26), this ansatz reduces a large number of parameters in the ansatz, and could be efficient for the implementation and optimization in practice. We may even discard the second term in Eq. (28) to further reduce the gate count in the variational quantum circuits. Especially, we find this ansatz works well numerically.

Now, say that we can decompose the operator into Pauli operators as with elements of the Pauli group and the variational parameters . Then the unitary coupled cluster operator can be decomposed in terms of the imaginary part of as We can use the first-order Trotter formula, where each term can be efficiently implemented on a quantum computer. It is worth mentioning that since the ansatz in Eq. (25) and Eq. (26) is local in terms of the mode, and therefore the Pauli terms could be easily calculated.

### iv.3 Variational quantum algorithms

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 mcardle2019variational. The normalized imaginary time evolution at imaginary time is given by The population of the energy eigenstate will decay exponentially with the energy , and the ground state can be obtained in the long time limit 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 mcardle2019variational, 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

 δ∥(dτ+H−Eτ)|ψ(θ(τ))⟩∥=0, (29)

and the evolution of the parameters under the imaginary time evolution could be determined by

 ∑jAi,j˙θj=−Ci, (30)

with the matrix elements of and given by

 Ai,j =Re(∂i⟨ψ(θ(τ))|∂j|ψ(θ(τ))⟩), Ci =Re(∂i⟨ψ(θ(τ))|H|ψ(θ(τ))⟩). (31)

Here, is the norm of the quantum state, we denote , 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. (22). 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 is indeed the Fisher matrix stokes2019quantum.

Moreover, having found the ground state , 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 becomes , the first excited state of the original Hamiltonian . As is an excited state of the modified Hamiltonian, we can evolve the system under in the imaginary time to suppress the spectral weight of and obtain the first excited state . This process can be repeated to obtain the higher-order excited states. More specifically, the ()th excited state is the ground state of effective Hamiltonian

 H(n+1)=H+αn∑j=0|ψ(j)⟩⟨ψ(j)|. (32)

Instead of preparing the Hamiltonian directly, we can simulate the imaginary time evolution under by tracking the evolution of the parameters, which are now modified as

 Ci= Re(∂i⟨ψ(θ(τ))|H|ψ(θ(τ))⟩+ αn∑j=0∂i⟨ψ(θ(τ))|ψ(j)⟩⟨ψ(j)|ψ(θ(τ))⟩), (33)

while the matrix keeps the same as in Equation (31). These addition terms in can be evaluated using the low-depth swap test circuit. Other variational excited state preparation techniques can be found in a recent review paper sugurureview21.

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. motta2020determining. Suppose the Hamiltonian has the decomposition where the Hamiltonian contains local terms and each acts on at most neighboring qubits. Using the first-order Trotterization, the evolved state after applying nonunitary operator within imaginary time by

 |Ψ(τ+δτ)⟩=c−1/2e−δτ^hl|Ψ(τ)⟩≈e−iδτ^A|Ψ(τ)⟩, (34)

where is the normalization factor and is a Hermitian operator that acts on a domain of qubits around the support of . The unitary operator can be determined by minimizing the approximation error in Eq. (34), which is similar to the derivation in Eq. (29). For a nearest-neighbor local Hamiltonian on a -dimensional cubic lattice, the domain size is bounded by , where is the correlation length. More details about the algorithm complexity can be found in Ref. motta2020determining.

This circuit construction strategy can be regarded as a special case in the variational imaginary time evolution given by Eq. (29) and Eq. (31). If we fix the old circuit ansatz constructed before imaginary time , and determine the new added unitary operator that approximates the effect of , this is exactly the same as Eq. (34). However, to further reduce the circuit depth, we can jointly optimize the parameters making it more compatible for the near-term quantum devices.

### iv.4 Spectral crowding

Before we start to apply variational algorithms, we will make a short investigation on the spectrum of the quantum field theory. In the momentum space, harmonic oscillator 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 landsman2019two 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

 ∑iniω(pi)=∑j¯njω(¯pj). (35)

Here, we have states represented in the harmonic oscillator basis in the momentum space, with particle numbers and momenta , or , and their energies are precisely identical. A typical example is that considering the continuum limit, we might have

 n|m0|=√m20+p2, (36)

where . In those cases, their states are degenerate. Another typical case the role of parity which anti-commutes with the momentum:

 ω(p)=ω(−p), (37)

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 and (free) or (interacting), with maximally three excitations , system size , and the lattice spacing . The choice of parameters is aiming on avoiding the situation in the Eq. (36).

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 one-particle 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:

 p=2πL(0,1,2,3):E=2.754,3.027,3.170,3.027. (38)

We know that this degeneracy is made by the boundary condition of the momentum , which is the parity  222However, for this set of parameter choices, the energies of two-particle and three-particle zero-momentum states are lower than the single-particle excited states. In general, for an -particle state, since we could freely choose the direction of momentum, the spectral crowding will be enhanced at least .

Now, we consider to turn on the interaction. In the adiabatic process where we slowly turn on the Hamiltonian,

 H(s)=H0+sλ0HI, (39)

and , 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 (which means that we are dividing the interval to 101 steps). We find all single-particle eigenstates could agree with the corresponding energy eigenstates with high fidelities (we only show 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 , and we slowly turn on the interaction from the free theory . 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. yuan2019theory; mcardle2019variational. 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,

 |ψ(θ)⟩=(∏Lℓ=1Uℓ(θℓ))∣∣ψfree% theory states⟩. (40)

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

 ∑jMi,jdθjds=Vi, (41)

where

 Mi,j=ReAi,j+∂i⟨ψ(θ(s))|ψ(θ(s))⟩∂j⟨ψ(θ(s))|ψ(θ(s))⟩, Vi=ImCi+i∂i⟨ψ(θ(s))|ψ(θ(s))⟩⟨ψ(θ(s))|H(s)|ψ(θ(s))⟩,

and and are similarly defined in Eq. (31). 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 V.

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. (39), but with much fewer steps. We could then use the variational algorithms to find the ground state of the intermediate Hamiltonian , using which as the initial state for the next time step until reaches . Compared to finding the ground state of , 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

 P=a∑x∈Ωπ∇aϕ, (42)

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 vanderstraeten2019tangent). Usually, we expect that our momentum eigenstate could have the following form

 |Φp⟩=∑x∈ΩeipxTx|Φ⟩. (43)

Here

is the translation operator with the vector

. If the state is already a momentum- eigenstate,

 |Φ⟩=∑y∈ΩeipyTy|Ψ⟩, (44)

we have

 ∑x∈ΩeipxTx∑y∈ΩeipyTy|Ψ⟩=∑x,y∈Ωeip(x+y)Tx+y|Ψ⟩ (45)

Moreover, if the state is a linear superposition of the momentum- state and the momentum- state where

 ∑x∈ΩeipxTx⎛⎝c1∑y∈ΩeipyTy|Ψ⟩+c2∑y∈Ωe−ipyTy|Ψ⟩⎞⎠ =#×c1∑z∈ΩeipzTz|Ψ⟩+c2∑x,y∈Ωeip(x−y)Tx+y|Ψ⟩ =#×c1∑z∈ΩeipzTz|Ψ⟩+c2∑ueipu∑vTv|Ψ⟩ (46)

Note that the 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 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.

### iv.5 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 one-particle state will still remain a one-particle state in the interacting theory. In fact, if we consider momentum eigenstates of a single particle, , we could define the one-particle subspace by

 Vone-particle,free=spanp(|p⟩). (47)

Now, if we are adiabatically turning on each state towards the coupling , the space will become

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 momentum eigenstate definition more precise.

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 . 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

Here, is the projector of the space . By definition, if is high enough, should also be high enough.

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. A successful benchmark for our variational algorithm about the adiabatic state preparation before particle scattering will require a detailed check about fidelities in a lower number of clean qubits before scaling the system up.

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

 (51)

where is the basis and is the dimension of the subspace. Thus, for a variational state we have

 |ψ⟩=DΛ∑i=1ci|qi⟩</