Quantum Simulation of Bound State Scattering

The last few years have seen rapid development of applications of quantum computation to quantum field theory. The first algorithms for quantum simulation of scattering have been proposed in the context of scalar and fermionic theories, requiring thousands of logical qubits. These algorithms are not suitable to simulate scattering of incoming bound states, as the initial-state preparation relies typically on adiabatically transforming wavepackets of the free theory into wavepackets of the interacting theory. In this paper we present a strategy to excite wavepackets of the interacting theory directly from the vacuum of the interacting theory, allowing the preparation of states of composite particles. This is the first step towards digital quantum simulation of scattering of bound states. The approach is based on the Haag-Ruelle scattering theory, which provides a way to construct creation and annihilation operators of a theory in a full, nonperturbative framework. We provide a quantum algorithm requiring a number of ancillary qubits that is logarithmic in the size of the wavepackets, and with a success probability vanishing at most like a polynomial in the lattice parameters and the energy of the wavepacket. The gate complexity for a single iteration of the circuit is equivalent to that of a time evolution for a fixed time. Furthermore, we propose a complete protocol for scattering simulation using this algorithm. We study its efficiency and find improvements with respect to previous algorithms in the literature.


I. INTRODUCTION
The successes of last few years towards the implementation of quantum algorithms on real platforms are creating growing expectation regarding about the opportunities that quantum computation may open in different fields.One prominent area that has been particularly fruitful in providing examples of the potential advantage offered by quantum computation is high-energy physics [1], especially in what concerns data analysis [2][3][4][5][6][7][8][9][10][11] and simulations of lattice quantum field theories .Such theories have been used for many decades as a tool to numerically investigate several aspects of quantum field theory through the Euclidean path integral approach, which is a powerful tool in its range of applicability, but does not cover the whole class of phenomena that require the study of real-time evolution.For such problems the Hamiltonian formulation in Minkowski spacetime is the natural framework, but it is hardly approachable with classical computation.For this class of phenomena -and arguably in all the theory of fundamental interactions -scattering events are of special interest since they are essentially the only means we have to access those regimes of physics where quantum field theory is necessary.This work concerns state preparation for digital quantum simulation of scattering.
The prospect of large-scale, fault-tolerant quantum computers, however far they may be in the future, has already started to change our approach to lattice field theory through the seminal papers of Jordan, Lee, and Preskill (JLP) [12,13].To exploit the potential advantage offered by quantum computers, it is better to formulate the theory on a space lattice with continuous time rather than on a spacetime lattice, as traditionally done.In Refs.[12,13], an efficient quantum algorithm for simulation of scattering in the scalar theory ϕ 4 , requiring thousands of logical qubits, is provided and analyzed.The present work aims to contribute to this long-term perspective by making the first step towards digital quantum simulation of scattering events with incoming composite particles.Preparation of bound states is an essential task to ultimately perform simulations of important, real-life collider events, such as protonproton scattering at the Large Hadron Collider (see Ref. [50] for recent developments in the context of tensor networks and Refs.[51][52][53][54][55][56] for recent developments in the context of analog quantum simulation).
The same problem was treated in Refs.[57][58][59], with use of different, albeit equivalent, formulations such as the wavelet basis or the multiparticle decomposition of the Hilbert space.A common feature of these studies is that they all rely on excited states of the associated free theory to prepare wavepackets of the interacting theory.An immediate consequence is that these approaches can be used only for states of the interacting theory that can be obtained by smooth interpolation (typically by an adiabatic transformation) from states of the free theory.This excludes the case of bound states.
Having in mind the framework in Refs.[12,13], we provide in the present work a general strategy to prepare singleparticle wavepackets of elementary or composite particles, with lower and upper mass gaps, on a quantum computer.We assume that the preparation of the vacuum state of the interacting theory is available and that we have access to an interpolating operator between the vacuum and the particle we want to create.Once the incoming wavepackets are created, the time evolution and measurements steps of the quantum simulation algorithm proceed as in Refs.[12,13].
The key idea is to use the Haag-Ruelle scattering theory, which is an alternative and complementary approach to the Lehmann-Symanzik-Zimmerman (LSZ) theory.The Haag-Ruelle formalism is of great conceptual importance in the context of axiomatic quantum field theory since it provides a link between the LSZ framework and the Wightman axioms (see, for instance, Ref. [60] for more details).From an operative point of view, its success has been rather limited, on one hand because it is outperformed by the LSZ formalism in the context of perturbation theory and on the other hand because real-time evolution, essential to any scattering theory, is unmanageable with traditional lattice techniques.Quantum computers may offer a new boost to the Haag-Ruelle formalism in terms of operativeness.
The method we propose requires a number of ancillary qubits that is logarithmic in the size of the wavepacket.We provide a quantum circuit with a gate complexity that is equivalent to that of a time evolution for a fixed time, and a certain probability of success.We argue that this probability does not vanish faster than a polynomial in the lattice parameters in the continuum limit.For definiteness, we work here with a single scalar field, because it is an illustrative case and is directly comparable with what is available in the literature.However, the idea holds, mutatis mutandis, with other theories as well.
State preparation is typically the most difficult step in digital quantum simulation of scattering, both technically and in terms of complexity.This work provides an innovative approach to the topic and sets the route to the preparation of composite particles.Here we consider the problem of preparing the vacuum only briefly.It is an interesting problem on its own and has already been addressed in other papers [61][62][63][64][65][66].For example, one may use the free vacuum preparation and the adiabatic transformation in Ref. [13], with the simplification that no backward time evolution is needed to contain premature wavepacket propagation.We consider this possibility and study its efficiency, with the purpose of providing an estimation of the total scaling of state preparation using our approach.We also compare our approach with the approach of Jordan, Lee, and Preskill, and find that our protocol has a comparable or better scaling than theirs in some cases.If one is interested in only scattering amplitudes, a simpler approach not involving quantum simulation would be the one used in Ref. [67].However that approach works only for fixed final states and for processes with a small total number n of ingoing and outgoing particles, as the complexity scales exponentially with n.
We mentioned previously that we rely on the existence of lower and upper mass gaps, but in general we can also have bound states immersed in the continuum of multiparticle states under the condition that they are protected by some symmetry.In this case our strategy is still suitable with some extra caveats.Fermionic theories do not present extra problems apart from the ones related to encoding of anticommuting degrees of freedom on a qubit system.Gauge theories, as usual, require special attention, not only because of the the well-known issues related to quantum simulation of these theories, but also because of the formulation of the Haag-Ruelle theory in the presence of massless particles such as photons.Nevertheless, we assume that with proper care these problems can be solved.For ultrarelativistic particles, the gap closes, and correspondingly the difficulty of creating a particle increases with its momentum.This approach is not suited for the creation of massless particles, for which there is no mass gap.
The rest of this paper is divided into three main sections and an appendix.In Sec.II we very briefly introduce the axiomatic approach to quantum field theory, where the Haag-Ruelle scattering theory is developed.We list the Wightman axioms and introduce basic concepts of the theory.We discuss first a theory with a single elementary particle, and then proceed to a theory in one space dimension with both an elementary particle and a composite one.We end Sec.II with some remarks on the applicability of our work.In Sec.III we describe how the Haag-Ruelle scattering theory can be used for state preparation in digital quantum simulation.We provide a quantum algorithm for particle creation from the interacting vacuum, and we analyze its complexity and probability of success.In Sec.IV we consider a full protocol for initial-state preparation based on adiabatically preparing the interacting vacuum and then creating wavepackets from it.In the Appendix we discuss the truncation of the field operator in a generic site of the lattice.

II. THE HAAG-RUELLE SCATTERING THEORY
The Haag-Ruelle scattering theory is best developed in the framework of axiomatic quantum field theory.The axiomatic approach provides a rigorous framework to construct a quantum field theory.It is based on a set of axioms formulated by Wightman and incorporates at the same time the principles of quantum mechanics and special relativity.
An important feature of this approach is that free and interacting theories are treated on equal terms.In particular, interacting theories are not seen as extensions of free theories obtained by adding interaction terms to a quadratic Hamiltonian.The difference between free theories and interacting theories is mainly of a pragmatical nature, because we can solve and construct explicitly free theories, but for most of the interacting theories the same is not true.For this reason it is quite uncommon to use this approach operatively.In the next two sections we review the Haag-Ruelle formalism, first in the context of a theory with only an elementary particle, and then in the context of a theory with bound states.

A. Scalar Field Without Bound States
In this section we give a brief review of the axioms following essentially chapter 9 in Ref. [60].We consider d = D −1 space dimensions, and we consider here a theory containing a single scalar field, whose dynamics give rise to a single kind of particle of mass m.Examples of such theories include the free scalar field theory and the ϕ 4 theory at weak coupling.The latter has been shown to satisfy the axioms for d = 1 and d = 2.In the physical case d = 3, it is believed to be trivial, i.e. equivalent to a free theory, while for d ≥ 4, triviality has been proven [68][69][70][71][72][73] (see also the discussion at the beginning of Sec.2.1 in Ref. [13]).The free theory is known to satisfy the axioms.The framework can be adapted straightforwardly to the case of a theory with fermionic fields.All that is required is to change commutation relations into anticommutation relations.
It is convenient to divide the whole set of axioms into a few families according to their content.The first one concerns the space of states and the spectral properties of the theory: 2. Axiom Ib The infinitesimal generators P µ of the translation subgroup T (x) = U (1 1, x) of the Poincaré group have a spectrum p µ restricted to the forward light cone, p 0 ≥ 0, p 2 ≥ 0.

Axiom Ic
There is a unique state |Ω⟩, the vacuum, with the isolated eigenvalue p µ = 0 of P µ .

Axiom Id
The theory has a mass gap: the squared-mass operator has an isolated eigenvalue m 2 > 0, and the spectrum of P 2 is empty between 0 and m 2 .The subspace of H corresponding to the eigenvalue m 2 carries an irreducible spin-0 representation of the homogeneous Lorentz group.These are the single-particle states of the theory.The remaining spectrum of P 2 is continuous, and it begins at (2m) 2 .
Clearly the specific form of the generators P µ depends on the theory at hand and the generator of translations in time is the Hamiltonian, P 0 = H.One-particle states can be labeled according to their momentum and can be written as wavepackets in momentum space with some wave function ψ.
The next family of axioms concerns the operator content of the Hilbert space and establishes what kind of fields appear in the theory: 5. Axiom IIa An operator-valued (tempered) distribution φ(x) exists such that for any Schwartz test function f (x) the smeared field is an unbounded operator defined on a dense subset D ⊂ H.Moreover, ϕ f D ⊂ D, allowing the definition of arbitrary (finite) products of smeared fields.
We recall that a Schwartz function f (x) is infinitely differentiable (i.e. is C ∞ ) and, together with all its derivatives, falls faster than any power of x as x goes to infinity.
6. Axiom IIb Under the unitary representation of the Poincaré group U (Λ, x) introduced in Axiom Ia, the smeared fields transform as 7. Axiom IIc Let f 1 and f 2 be Schwartz functions of compact support: thus, if f 1 vanishes outside a compact spacetime region v 1 and f 2 vanishes outside the compact region v 2 , and if 8. Axiom IId The set of states obtained by applying arbitrary polynomials in the smeared fields ϕ f (with all possible Schwartz functions f ) to the vacuum state |Ω⟩ is dense in the Hilbert space H.
In Axiom IIa we introduce an important difference with respect to canonical quantization, namely smeared operators.
In a free theory we can identify the operator-valued distribution φ(x) with the familiar field operator.This is not a well-defined operator because the state φ(x) |Ω⟩ has an infinite norm.To avoid this problem, it is necessary to introduce smearing and treat φ(x) as a distribution.In this way the state ϕ f |Ω⟩ has a finite norm for any Schwartz function f and ϕ f is a well-defined (unbounded) operator.Then, by Axiom IIb, we can define There are two more axioms of great importance to develop a satisfactory scattering theory: 9. Axiom IIIa For some one-particle state with discrete eigenvalue m 2 of the squared-mass operator, the smeared field ϕ f (x) has a nonvanishing matrix element from this single-particle state to the vacuum, 10. Axiom IIIb (asymptotic completeness) The Hilbert space H in (or H out ) corresponding to multiparticle states of far-separated, freely moving stable particles in the far past (or far future) are unitarily equivalent, and may be identified with the full Hilbert space H of the system.
It should be noted that Axiom IIIb plays a crucial role in the derivation of the LSZ reduction formula, but here it is somewhat superfluous.The joint set of eigenvalues of P µ , labeled by p µ , is composed of three disconnected subsets (see Figure 1).There is the vacuum subset, containing only the origin p µ = 0, the one-particle mass hyperboloid, containing all the p µ points such that p 2 = m 2 , and the multiparticle continuum with all the points such that p 2 ≥ 4m 2 .In the two-particle subspace, for instance, the squared-mass operator gives With this in mind we define an operator ϕ 1 (x) exactly as in Eq. ( 7) but with a smearing function f 1 (x) chosen as the Fourier transform of a function f1 (p) with support in the region sandwiching the one-particle mass hyperboloid.This guarantees that the state ϕ 1 (x) |Ω⟩ is a one-particle (and oneparticle only) state, by Axiom IIIa.
Next consider a positive-energy solution of the Klein-Gordon equation, FIG. 1: Structure of the spectrum of P µ with one space dimension.The blue line represents the one-particle mass hyperboloid, the red region is the multiparticle continuum, and the green region is the region defined by and we define the operator where the double derivative is defined as in The state ϕ 1,g (τ ) |Ω⟩ can be shown to be independent of τ , This is no longer true in general if we consider multiple applications of such operators at the same time τ , |Ω⟩, but in this case we can rely on the following theorem (see chapter 9 in Ref. [60] for more details).
Theorem 1 (Haag Asymptotic).The time-dependent state vector converges strongly in the limit τ → −∞ to the n-particle in-state with momentum wave functions The convergence rate of the limit |Ψ, τ ⟩ → |Ψ⟩ in is |τ | −d/2 in the general case, which includes the case where some or all of the wavepackets g 1 , . . ., g n overlap with each other and come from the same direction.In this case, the convergence is guaranteed by spreading of the wavepackets to the point where they are so broad that they cease to interact.In practice we consider only wavepackets coming from different directions, and in this case the fast decrease of the Schwartz functions ensures that the convergence rate is faster than any inverse power of |τ |.
The state |Ψ⟩ in is a Heisenberg state, which implies that it is not to be visualized in general as a state made of n (spatially) well-separated wavepackets.Its form depends on the reference time at which the Heisenberg state and the Schrödinger state coincide.As the reference time, we can choose a moment well before the collision between the wavepackets occurs, in which case we indeed have well-separated wavepackets, or a moment during the collision or later, in which case we can expect to have a complicated state more or less spread in space.The strong convergence of the theorem is to be taken at the same reference time, whether in the far past or not, for both |Ψ, τ ⟩ and |Ψ⟩ in .This reference time should not be confused with the parameter τ appearing in the theorem.But since we will work in the Schrödinger picture this subtlety is not relevant to us.
Let us see more explicitly what we have stated so far, starting from the state where is obtained from (12) after some simple manipulations.For our purposes it is convenient to move to the Schrödinger picture by writing x = (t, x) and plugging into (18).Furthermore, shifting the integration variable in (18) by t → t + τ , we get Let us have a closer look at ψ(t + τ, x; τ ).We can rewrite it as If, without loss of generality, we assume f 1 (t, x) peaked at around t = 0, so is f ′ 1 (p; t).Roughly speaking, the effect of f ′ 1 (p; t) is to spread ψ(t + τ, x; τ ) without changing its position.Therefore, by (23) we see that ψ(t + τ, x; τ ) is essentially a solution of the Klein-Gordon equation moving through space with time τ .When τ < 0, a † ψ (τ ) creates a wavepacket at some point along the past trajectory of ψ(t + τ, x; τ ).Then, the time evolution operator in (21) moves it to where ψ(t + τ, x; τ ) lies at τ = 0.This accounts for (14).
Similarly, in the case of two incoming particles we can write We can choose g 1 (τ, y 1 ) and g 2 (τ, y 2 ) such that their wavepackets ψ 1 and ψ 2 are always well separated from each other for τ ≤ 0, and are on a collision course for some τ > 0.Then, from (25), the action of the operators ϕ 1,g1 (τ )ϕ 1,g2 (τ ) on the vacuum is clear: as τ → −∞, the two wave functions ψ 1 (t 1 +τ, x 1 ; τ ) and ψ 2 (t 2 +τ, x 2 ; τ ) are sent far from each other where the two creation operators a † ψ1 (τ ) and a † ψ2 (τ ) act undisturbed by each other.Then, the time evolution operator e iHτ evolves the system forward making the two wavepackets approach and interact with each other.Let us call r the distance between the regions where ψ 1 (t 1 + τ, x 1 ; τ ) and ψ 2 (t 2 + τ, x 2 ; τ ) are concentrated at τ = 0. Provided that interactions between particles are short-ranged in the theory, as is the case for gapped theories, we can ignore the interactions occurring between the two wavepackets during the evolution from τ = −∞ to τ = 0, with error vanishing faster than any inverse power of r.Thus, at the cost of slightly increasing r, we can take τ = 0 in (25) with excellent precision.

B. Scalar Field with Bound States
The theory with a single scalar field with interactions λ( φ6 − φ4 ), in a single space dimension, is a simple example of a quantum field theory displaying bound states.This theory is known to satisfy the Wightman axioms [73], and at weak coupling, it has a single composite particle below the two-particle threshold 2m.Its mass m b can be computed perturbatively [74][75][76] as The bare mass m 0 (λ) can be set such that m = m(m 0 (λ), λ) is fixed.We notice at this point that Axiom Id should be modified in an obvious way to accommodate the composite particle with mass m b .This is a well-controlled model that is perfect to see in a simple and explicit way how to treat bound states in the Haag-Ruelle formalism.The last ingredient we need to know is that the field : φ2 : (x) interpolates between the vacuum and one-particle states of the composite particle [76].The Wick ordering : • : is necessary in the continuum to avoid divergences of φ2 (x), but on the lattice it amounts to a constant, finite shift.This is completely irrelevant to our purposes as can be easily seen.
Let Φ(x) denote either φ(x) or : φ2 : (x).Then, if we shift it by a constant A and smear it with the function ψ in (19), we obtain The second term is proportional to f1 (0), which is zero by construction.
The joint spectrum of P 0 = H and P 1 is depicted in figure 2. Starting from this picture, we find that the Haag-Ruelle theory discussed in the previous section holds equally well to obtain wavepackets of the elementary particle with mass m, or wavepackets of the composite particle with mass m b .If we want to have an elementary particle, we use the interpolating field φ(x), a smearing function f 1 (x), with Fourier transform f1 (p) sandwiching the hyperboloid of mass m, and a solution g(x) of the Klein-Gordon equation to build the operator ϕ 1,g (τ ) of equation ( 12) as before.If, on the other hand, we want to obtain a composite particle, we use the interpolating field : φ2 : (x), a smearing function f b (x), with Fourier transform fb (p) sandwiching the hyperboloid of mass m b , and a solution h(x) of the Klein-Gordon equation to obtain the operator where Then everything proceeds as before.A useful remark is that the field φ(x) does not interpolate between the vacuum and one-particle states of the composite kind, and vice versa, : φ2 : (x) does not interpolate between the vacuum and one-particle states of the elementary kind.In the formulae, if |α 1 ⟩ is an elementary one-particle state and |α b ⟩ is a composite one-particle state, we have The first equality holds because the one-particle sector of the elementary kind is spanned by vectors of the form ϕ 1 (x) |Ω⟩.Then the matrix element ⟨α 1 | : φ2 : (x) |Ω⟩ is an integral containing ⟨Ω| φ φ φ |Ω⟩, which is zero because the Hamiltonian contains only even powers of the field φ.A similar argument leads to ⟨α b | φ(x) |Ω⟩ = 0.As a consequence of the equalities in Eq. ( 30), f1 and fb are allowed to have supports intersecting both hyperboloids of mass m and m b .
We conclude this section with some remarks on the validity of the theory just described.In general, we can have bound states whose mass falls above the two-particle threshold of lighter particles, on the condition that they are protected by some symmetry (internal quantum numbers).In this case the symmetry selects a sector of the Hilbert space, and when we restrict the mass-squared operator to this sector, the mass of the composite particle appears as a discrete point in the spectrum of P 2 again.
The Haag asymptotic theorem critically depends on two assumptions: 1.There exists a lower and an upper mass gap for the particle we want to create in such a way that it is possible to sandwich the mass hyperboloid corresponding to such a particle.
2. We have access to an operator interpolating between the vacuum and one-particle states of the particle we want to create.
Clearly, studying these two conditions strongly depends on the theory under consideration and can be very difficult, but we can use standard techniques of lattice quantum field theory to study these properties model by model, or the Bethe-Salpeter equation as done in Refs.[74][75][76].Smearing a field operator in time can be avoided if we have at our disposal an operator that does not couple the vacuum to multiparticle states (as happens for free theories).For a bound state whose mass is immersed in the continuum of other particles, we also need to ensure that the interpolating operator couples only to the sector where the bound state lives, in a way similar to the reasoning leading to (30).
As an example, consider a theory where we have a composite particle with spin 0 and mass m (a pion), and another composite particle with spin 1/2 and mass M > 2m (a proton).These two particles are made of the same underlying elementary fields, but the difference in spin should make it easy to build interpolating fields from the vacuum to each of the two particles without crossing, even though the heavier particle has mass in the continuum of the lighter one.

III. STATE PREPARATION USING THE HAAG-RUELLE THEORY
In the following we assume that all of this holds as an approximation on the lattice.The Haag-Ruelle scattering theory was developed for spin systems in Ref. [77] and for Euclidean lattice field theory in Ref. [78].Detailed studies of the effects of latticization on applications for quantum simulation will be the subject of future work.
From now on we take τ = 0, ψ(t, x) = ψ(t, x; 0), as given in (19), and We want to prepare the state a † ψ1 a † ψ2 |Ω⟩, with ψ 1 and ψ 2 well separated.Because of this, we can focus on one wavepacket at a time, so in the following we will see how to implement the operator a † ψ on a quantum computer.In (19), we choose g peaked at around p with support of size δ p and f1 (p) peaked at around ( Ē, p) with support of size δ E in the p 0 direction and δ p in the other directions, with Ē = E(p).Also, we can take δ E = O(m).First we truncate the integration over t and the summation over x around the spacetime region where ψ is significantly different from zero, which, since ψ is a Schwartz function, introduces an error vanishing faster than any power as the hypervolume of the region is increased.We label the space points in this region by x 1 , . . ., x S and we approximate the integral with a sum over time points t 1 , . . ., t N with spacing δ t .By the uncertainty principle, the linear size of ψ over space is proportional to 1/δ p , and t N − t 1 is proportional to 1/δ E = O(1/m).Thus, as an approximation of a † ψ , we have anc. / V Φ(t1) Φ(t2) . . .FIG. 3: Description of the circuit implementing ā † ψ .The slash at the beginning of each line means that the line represents a register of qubits: "anc." is the ancillary register of N a qubits, x i is the register of k qubits dedicated to the site x i , and Γ is the set of qubits dedicated to the rest of the lattice.
We work in the field basis [12,13,79], where the operator φ(x) is diagonal.If k qubits are dedicated to the lattice site x, then φ(x) is implemented by a linear combination of Z Pauli matrices, The effects of truncating and discretizing the spectrum of φ(x) were studied in Refs.[12,13,79], where it is also shown that it is enough to take k quite small in a broad range of situations.The operator (33) is not unitary but can be implemented with linear combination of unitaries (LCU) [80].Furthermore, the operator ā † ψ can be implemented with LCU as well, with probability of success It is not easy to determine ρ exactly (more is said on this at the end of this section), but one could use, for example, the techniques described in Ref. [81] to find numerical estimates for it.Then one would have to repeat state preparation O(1/ρ) times to get the desired initial state, or, alternatively, one could apply amplitude amplification [81] to obtain a quadratic speedup at the expense of a larger circuit depth.Typically, we need to prepare two incoming wavepackets.
If we denote by ρ 1 and ρ 2 the probabilities of obtaining the two wavepackets, the total probability is given by ρ = ρ 1 ρ 2 .Generalization to more than two wavepackets is straightforward.
Before moving on, we remark on how to implement the bound-state creation operator (28).The discussion goes along the line described so far, except that now, instead of (33), we have to use the operator up to a shift.We are free to choose this shift in such a way as to eliminate the terms with i = j in (35), which are proportional to 1 1.In this way, instead of α in (34), we obtain Circuit description We want to give a high-level description of a circuit implementing ā † ψ , so we focus only on the dependence on the lattice sites and the time to keep the discussion concise.We take a register of N a = ⌈log(kN S)⌉ ancillary qubits and we label the computational basis as |t i , x j ⟩, with i = 1, . . ., N and j = 1, . . ., S.
We define operators V ψ and V ′ ψ such that Then the circuits in Fig. 3 implement the operator ā † ψ written in (32) up to a normalization factor and when the state ⊗Na is obtained by measurement of the ancillary register.In general we have t 1 < 0, t N > 0 and t i − t i+1 = −δ t ; therefore, the sequence of time evolution operators appearing in Circuit 1 consists of a backward evolution for time |t 1 |, followed by N steps forward, each one of time δ t for a total of t N − t 1 , and by a final backward evolution for time t N .In Circuit 2 we notice that controlling only each φ is equivalent to controlling e −iHti φe iHti because e −iHti e iHti = 1 1.To prepare two wavepackets, we can take two ancillary registers, apply in sequence the circuit of Fig. 3, one for each wavepacket, and measure both ancillary registers at the end.
Complexity We show that the complexity of Circuit 1 is dominated by the sequence of time evolutions.We do not want to discuss here how to implement the time evolution, as this is not in the scope of this work.Moreover, we take k = Õ(1).
The sequence of operators Φ(t 1 ), . . ., Φ(t N ) requires O(kSN ) = Õ(SN ) gates, as well as the operators V and V ′ , which basically provide generic state preparation on N a qubits, and we have S < V. We can estimate the error introduced by discretizing the integral over t, and hence how large N needs to be, in the following way.We take t 0 = t 1 − δ t /2 and T = t N − t 1 + δ t .Then we split the integral from t 0 to t 0 + T into N integrals in the following way: where from the first line to the second line we have shifted the integration variable, t → t+t i .We expand e iHt φ(x)e −iHt using the formula and we expand ψ(t i + t, x) using the Taylor expansion around t = 0 up to order t 2 .Odd orders in t do not contribute because the integration domain is symmetric around zero.The leading order gives us exactly the operators appearing in (32).We use the spectral norm of the next to leading order to estimate the error due to discretization and we apply the triangular inequality: The dominant contribution is given by the term with H, [H, φ(x)] , and the quantity is approximately a constant independent of the lattice and the precision.Finally, given that T = δ t N , we have If we use a first-order Suzuki-Trotter formula to implement the time evolution, we need [82] Trotter steps to keep the error below ϵ.As one may check by looking at Eqs. ( 45) and ( 46), the commutator [H ϕ , H π ] scales with the total size of the lattice, as it involves a summation over all the sites, while H, [H, φ(x)] involves only a few neighbors of x.Thus, N ST clearly dominates over N even if we consider higher-order product formulae.This shows that the complexity is determined by the time evolution.
For the ϕ 4 theory we have On the success probability Here we want to show what to expect regarding the probability of success.We write it again here: Near the continuum limit, ||ā † ψ |Ω⟩ || 2 should approach its continuum value where Z depends on the normalization of the field operator φ(x).The summation in α can be approximated by an integral, so we have Since g(p) f1 (p 0 , p) has support of size δ E in the p 0 direction, and support of size δ p in the other directions, we have In the limit δ p → 0 we can approximate g with a delta function centered in p, and recalling Ē = E(p), we obtain and We now turn our attention to the factor Z/ϕ 2 max . Set where |ψ E ⟩ is a state with maximum energy E and E is the energy of the process being simulated.In Ref. [13], as a consequence of Chebyshev's inequality, it is shown that where ϵ trunc is the error due to truncation of the spectrum of φ(x).In the field representation we have where 1, . . ., V is some labeling of the lattice sites.It is well known [83] that ψ E (ϕ 1 , . . ., ϕ V ) decays rapidly (exponentially or more) outside the classically forbidden region V (ϕ 1 , . . ., ϕ V ) > E, where V (ϕ 1 , . . ., ϕ V ) is the potential associated with the Hamiltonian H. From the discussion in the Appendix it is clear that we can improve the bound (55) to Moreover, from the analysis in Ref. [13] it follows that (at most) σ E = O( √ Eσ 0 ), and thus, ignoring the logarithmic contribution, we have We can use the Kållen-Lehmann representation of the two-point function in the continuum to estimate σ 0 .We have where ρ(µ) ≥ 0 is the so-called spectral function and In renormalizable models such as ϕ 4 , the spectral function is known to decay as 1/µ at all orders of perturbation theory, ensuring the convergence of the integral in (59).Therefore, we may assume that, at coincident points x = y, the behavior of the two-point function is determined by ∆ + (0, m 2 ), which diverges logarithmically for d = 1 and linearly and quadratically (apart from logarithmic factors) for d = 2 and d = 3, respectively.As 1/a is a natural cutoff on the lattice, we may conclude that where we use the notation Õ(•) = O(• polylog(•)), and consequently up to logarithmic factors.As we are interested in the continuum limit, the dominant factor is a d−1 /V.From the discussion in the Appendix we expect this bound to be quite loose.For the perturbative ϕ 4 theory, for instance, we can replace V with √ V.

IV. COMPARISON WITH THE JLP ALGORITHM
In this section we compare the complexity scaling of our strategy with the complexity scaling of the JLP strategy [12,13] in the case when both methods are applicable, that is when we want to prepare elementary particles.We consider only the ϕ 4 theory.The settings in the two cases are very similar and a comparison is immediate.
The JLP algorithm can be summarized in five steps: 1. Free vacuum preparation 2. Creation of free wavepackets 3. Adiabatic transformation of the free wavepackets into interacting wavepackets 4. Time evolution;

Measurements
As stated in the introduction, our strategy requires the preparation of the interacting vacuum.While other techniques, such as variational approaches, might be more suitable in practice, to study the complexity we choose here to do the following: 1. Free vacuum preparation 2. Adiabatic transformation of the free vacuum into the interacting vacuum 3. Creation of interacting wavepackets 4. Time evolution

Measurements
Steps 4 and 5 are the same in the two approaches, and in both cases the bottleneck of complexity is in the initial-state preparation, steps 1,2 and 3, so we focus on these steps only.We consider here the case of two particles in the initial state.As the complexity depends on the success probability, simulating scattering between three or more incoming particles becomes more and more inefficient in our approach, and is of little or no practical interest in general.We assume the reader is familiar with Ref. [13], especially with Secs.3.2 and 4.2.We do not want to discuss the details of how time evolution is implemented, although it may have an impact on complexity.Our goal is to give a rough estimate, certainly not exhaustive, and to compare it with the results obtained by Jordan, Lee, and Preskill.To this end and following their approach, we assume we can implement time evolution with a kth-order Suzuki-Trotter formula with large k.The gate cost for time t on a lattice of V sites is which we indicate as to simplify the presentation.The little-o notation is used also to include logarithmic factors implicit in the Õ notation.The most-time-consuming parts are typically the free vacuum preparation and the adiabatic transformation.The complexity for the free vacuum preparation, using the Kitaev-Webb method [84,85], is O(V 2.376 ).This exponent is determined by the classical computation of the LDL decomposition of the covariance matrix.The quantum circuit has a depth of Õ(V 2 ) and requires O(log V) ancillary qubits [85].
The adiabatic transformation of the wavepackets requires a modification with respect to a traditional treatment to take into account the fact that wavepackets are not eigenstates of the Hamiltonian, as described in Secs.3.1 and 4.2 of [13].In particular, the adiabatic transformation, of total time τ wp , needs to be split into J ∼ √ τ wp steps and interspersed with backward time evolution to suppress the dynamical phases.This causes the adiabatic error ϵ to vanish like J 2 /τ 2 wp ∼ 1/τ wp instead of as 1/τ 2 wp .Our approach has the advantage that the adiabatic transformation is performed on the vacuum, rather than on the wavepackets, and there is no need to suppress the dynamical phases.To continue the discussion, we distinguish between the weak-coupling regime and the strong-coupling regime from now on.In both regimes we take the adiabatic paths chosen in Refs.[12,13].
Weak coupling We consider scaling in the continuum limit.To determine the scaling of the time τ vac required to perform the adiabatic transformation in the vacuum, we can use the analysis of adiabaticity in Sec.4.2 in Ref. [13], taking J = 1.Setting V = a d V, this leads to to be compared with the result obtained by Jordan, Lee and Preskill We next consider the number of gates needed for the free vacuum preparation, G prep , and the two kinds of adiabatic transformation, G wp and G vac .To this end, we take a ∼ √ ϵ and V ∼ log(1/ϵ), as argued in Ref. [13].There, in section 3.2, they also provide G prep = Õ(1/ϵ 1.188d ) and, by (64), Following the same reasoning we find We see that G prep dominates over G vac for all the dimensions.
The advantage obtained by avoiding the adiabatic transformation on wavepackets is partially spoiled by the fact that the wavepacket creation does not succeed with a probability of 1.For simplicity, we consider two similar wavepackets, so, by equation ( 62) (replacing V with √ V in weakly coupled ϕ 4 theory), we find that the probability goes as This means that we need to repeat the state preparation times to obtain the correct initial state.Notice that the LDL decomposition of the covariance matrix does not need to be repeated every time, so the total complexity of the state preparation protocol proposed here is obtained by multiplying the depth of the free vacuum preparation by 1/ρ, which gives with depth given by Õ(V 2 ) = Õ(1/ϵ d ).The total complexity can be improved at the expense of a larger depth by use of amplitude amplification, in which case we have For a direct comparison, the total complexity of the JLP algorithm is given by Strong coupling Because of the triviality issue of the ϕ 4 theory in a three-dimensional space, we can have strong coupling only for d = 1, 2, when we approach the phase transition.For the success probability, we take We limit our discussion to the scaling of complexity with the coupling strength and the momenta of the incoming particle.If the phase transition occurs at the critical value λ c , we can take |λ 0 − λ c | as a measure of the coupling strength.To estimate the scaling of the adiabatic time with this quantity, we use the results in Ref. [86], implying τ vac ∼ 1/m 3 .The temporal size T of the wavepacket ψ in (32), which determines the duration of the time evolution required to create the wavepacket in our approach, grows at most as 1/m, by the uncertainty principle, so τ vac dominates over T .Near the phase transition, the physical mass vanishes as which gives Furthermore, the probability of success (74) does not depend explicitly on the mass gap.However, the volume has to be large enough to contain the wavepackets, which in turn have linear size proportional to 1/m; hence, Considering the volume dependence in equation ( 64), the adiabatic transformation has a stronger scaling than the free vacuum preparation.Taking into account the success probability (74), we find that, at fixed a and incoming momenta, the total complexity scales with the coupling strength as Using amplitude amplification, we have The corresponding result found in Refs.[12,13] is Finally, we consider scaling with incoming momenta, at fixed coupling strength and volume.The free vacuum preparation and the adiabatic transformation in the vacuum do not depend explicitly on the incoming momenta, but the lattice spacing has to be small enough to contain momentum mode p.With a ∼ 1/p, the free vacuum preparation has a cost proportional to p 2d , while the adiabatic transformation has a slower growth.The time evolution required for the wavepacket creation has a cost count of (VT ) 1+o(1) ∼ p d+1+o (1) , since T grows at most linearly with p.The success probability (74) has scaling p 4d+2 .Putting all this together, we find the scaling G strong ∼ p 8+o(1) d = 1, or, using amplitude amplification, The corresponding result in Refs.[12,13] is The scaling with the spread of the wavepackets in momentum space, δ p , can be obtained in a similar way.In this case, we take V ∼ 1/δ d p , while we keep a fixed.The total scaling is given by , which can be improved to with the use of amplitude amplification.

V. CONCLUSIONS
In this paper we provide a quantum algorithm to create single-particle wavepackets of a lattice quantum field theory starting from the vacuum state.The method we propose is quite general and the idea is independent of details of the model.For example, it works equally well for free and interacting theories.The key aspect of our strategy is that it is suitable for preparation of composite particles, which is an important novelty in the context of quantum simulation of relativistic scattering.To our knowledge this is the first work on state preparation of bound states for digital quantum simulation of scattering.
The work is based on the Haag-Ruelle scattering theory in the framework of axiomatic quantum field theory, which is ideal for quantum simulation as it is developed in the operator formalism.In this respect, this work also shows the potential importance that the axiomatic approach might have for quantum computation applied to quantum field theory, as both fields are suited to nonperturbative investigations from first principles.The Haag asymptotic theorem 1 provides a strong limit to obtain scattering states, rather than a weak limit as is the case in the more famous LSZ approach.This feature is what makes the Haag-Ruelle framework particularly suitable for quantum simulation, together with the fact that the convergence rate of the strong limit is very fast, as discussed after Theorem 1.
Our present result shows the potential of the Haag-Ruelle formalism in the context of quantum simulation.Here we decided to use LCU because it seems to us the easiest and most natural route, but there may be other, more efficient techniques in the context of digital quantum computation or in other contexts.Excluding the steps of vacuum preparation, and of time evolution and measurement, our algorithm requires a number of qubits that grows logarithmically with the size of the wavepacket, and a circuit depth equivalent to that of the time evolution.It succeeds with a probability that vanishes polynomially in the continuum limit for highly energetic processes and for narrow wavepackets in momentum space.The size of the mass gap is relevant only because the extension of the wavepacket in the time coordinate is inversely proportional to the mass gap.This work decomposes the problem of state preparation for scattering into more approachable ones.On one hand, efficient techniques to prepare the vacuum state are required.On the other hand, one has to find interpolating operators with the right properties for a given particle in a given theory, and one needs to know the size of the corresponding lower and upper mass gaps.On the first front, much work has already been done in the context of quantum computation, and we propose a way to prepare the interacting vacuum by an adiabatic transformation.We found that a scattering protocol based on creating wavepackets from the adiabatically prepared vacuum is in some cases more efficient than the protocol proposed by Jordan, Lee and Preskill [12,13].For the second front, standard techniques of Euclidean lattice field theory are available.As a next step, the approach of this work needs to be specialized case by case.Also, we need to investigate how gauge invariance and the presence of massless particles affect this approach.
In the last stages of the present work, we became aware of two new preprints [87,88], where similar problems are addressed with different approaches.

FIG. 2 :
FIG.2: Structure of the spectrum of P µ in the theory λ( φ6 − φ4 ).The blue line represents the mass hyperboloid of elementary particles, the cyan line represents the mass hyperboloid of two-particle bound states, and the red region is the multiparticle continuum.

1 .
Axiom IaThe state space H of the system is a separable Hilbert space.It carries a unitary representation U (Λ, x) (Λ is an element of the homogeneous Lorentz group and x is a spacetime coordinate vector) of the proper inhomogeneous Lorentz group (i.e., the Poincaré group).Thus, for all |α⟩ ∈ H, |α⟩ → U (Λ, x) |α⟩, with the U (Λ, x) satisfying the Poincaré algebra 1 e iH(t 1 −t 2 ) (b) Circuit 2: Overview of the operator implementing φψ (t i ).The symbol connecting φ in a squared box to the rounded box containing t i , x j represents the operator φ(x) controlled on the state |t i , x j ⟩.