Scattering Amplitude from Quantum Computing with Reduction Formula

Utilizing the Lehmann-Symanzik-Zimmermann reduction formula, we present a new general framework for computing scattering amplitudes in quantum field theory with quantum computers in a fully nonperturbative way. In this framework, one only has to construct one-particle states of zero momentum, and no wave packets of incoming particles are needed. The framework is able to incorporate scatterings of bound states, and is ideal for scatterings involving a small number of particles. We expect this framework to have particular advantages when applied to exclusive hadron scatterings. As a proof of concept, by simulations on classical hardware, we demonstrate that in the one-flavor Gross-Neveu model, the fermion propagator, the connected fermion four-point function, and the propagator of a fermion-antifermion bound state obtained from our proposed quantum algorithm have the desired pole structure crucial to the implementation of the Lehmann-Symanzik-Zimmermann reduction formula.


I. INTRODUCTION
The calculation of scattering amplitudes in quantum field theory has long been a core topic in theoretical particle physics [1][2][3][4][5].All tests of theories against experiments in particle accelerators entail theoretical predictions of scattering amplitudes.Despite the huge success of the perturbative approach to the calculation of scattering amplitudes [6][7][8], there are still circumstances in which the perturbative framework does not work, namely the cases where the coupling constants are large, as is the case for quantum chromodynamics at low energies for instance.To date, first-principle nonperturbative calculations of scattering amplitudes in quantum field theory are not available.The main obstacle is that real-time dynamics cannot be simulated in traditional path-integral lattice quantum field theory [9], while simulating real-time Hamiltonian evolutions in quantum field theory requires unbearable computational cost on a classical computer.In a series of papers, Jordan, Lee, and Preskill (JLP) proposed that, with the help of quantum computers, simulations of Hamiltonian evolutions of scattering processes in quantum field theory can be achieved with affordable computational cost on the lattice, making nonperturbative evaluations of scattering amplitudes possible [10][11][12].
Although the quantum-computational framework developed by JLP is fully general, it may encounter difficulties in practice.One major difficulty is that one has to prepare spatially well-separated wave packets of incoming particles in the initial state.The central values of the 4-momenta of these wave packets p i , as well as their Lorentz-invariant products p i • p j , set a constraint on the lattice spacing a and the lattice size L; while the spatial separation distances d ij between the initial-state wave packets set a constraint on the lattice size.In addition, the separations d ij have to be wide, meaning that d ij ≫ 1/∆p µ i , where ∆p µ i is the uncertainty of the wave packet of the ith incoming particle in momentum space, and we require ∆p µ i ≪ |p µ i | so that the wave packets are narrow enough to mimic a scattering process of definite incoming momenta.Constraint Eq. ( 1) is required for reliable simulations of incoming particles with definite 4-momenta on the lattice, and can be potentially improved by introducing factorization theorems using the method of effective field theory [29].Constraint Eq. ( 2) is due to the introduction of wave packets, which generally implies a larger lattice size than required by Eq. (1).Another feature of the JLP formalism is that the wave packets are first prepared with the coupling constant turned off.Therefore, the incoming particles cannot be bound states.The coupling constant is subsequently adiabatically turned on before the scattering occurs and adiabatically turned off after the scattering occurs.To ensure adiabaticity, a long time span of evolution is required.In addition, one has to insert backward evolutions in order to eliminate unwanted broadening of wave packets during the adiabatic turn on and turn off of the coupling constant, thus increasing the time complexity.In fact, in the strong coupling regime, most theoretical uncertainties come from the adiabatic turn on and turn off of the coupling constant [10].
We note that, in the conventional perturbative approach, scattering amplitudes are computed using the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula [37], which relates scattering amplitudes to n-point correlation functions, which in turn can be expanded as power series in the coupling constant using the Feynman diagram technique.The LSZ reduction formula, being a nonperturbative relation, is a natural alternative starting point for the evaluation of scattering amplitudes with quantum computers in a fully nonperturbative way.In this work, we propose a quantum computational framework for calculating scattering amplitudes using the LSZ reduction formula.In this framework, in order to evaluate scattering amplitudes, one calculate n-point correlation functions on a quantum computer.We will see that, this approach is ideal for scattering processes involving a small number of particles, and will have potential applications in exclusive strong-interaction processes such as 2 → 2 scatterings of pions or nucleons.
In the following, we first give a short review of the LSZ reduction formula.We then propose a quantum algorithm which utilizes the LSZ reduction formula to compute scattering amplitudes, and discuss its features and advantages.After that, as a proof of concept, in a simple model, the one-flavor Gross-Neveu model, we simulate the fermion propagator, the connected fermion four-point function, and the propagator of a fermion-antifermion bound state with our proposed quantum algorithm on classical hardware.We give a conclusion at the end.

II. LSZ REDUCTION FORMULA
The LSZ reduction formula relates the scattering amplitude of a given scattering process to correlation functions of fields in the vacuum [37].For instance, consider the scattering process h(k , where h is some spin-0 particle with mass m annihilated by a scalar field ϕ.Using the LSZ reduction formula, the scattering amplitude M can be written as where n = n in + n out .The G({p i }, {k j }) is the connected n-point function in momentum space, given by where T denotes time ordering, the subscript "con" denotes the connected part, and |Ω⟩ is the interacting vacuum, i.e. the ground state.The K(p) is the two-point function in momentum space, also called the propagator, given by The factor R is the field normalization, defined by where |h(p = 0)⟩ denotes the state with a single particle h with zero spatial momentum.The generalization of Eq. ( 3) to cases which involve multiple types of massive particles with arbitrary spin is trivial, with suitable inclusions of polarization tensors and spinors on the righthand side of Eq. ( 3).In essence, the LSZ reduction formula Eq. ( 3) says that the scattering amplitude is simply a connected n-point function in momentum space with momenta put on-shell, with external-leg propagators amputated.The field normalization factors √ R on the righthand side of Eq. ( 3) ensure that the scattering amplitude, as a physical observable, is independent of the normalization of the field operators which create or annihilate the external particles.It should be noted that the connected n-point function G({p i }, {k j }) has simple poles at p 2 i , k 2 j = m 2 , and so is divergent when the momenta are put on shell.On the other hand, the propagator K(p) also has a simple pole at p 2 = m 2 , namely Therefore, in Eq. ( 3), the pole singularities in G({p i }, {k j }) cancel with those in the K(p) factors, giving a finite scattering amplitude.In practice, when the continuum theory is approximated by a theory on the lattice, these singularities are tamed and the pole structure is approximated by some bounded function of p 2 that approaches it in the continuum and infinite-volume limits.A study of finite-volume effects on Minkowski correlation functions has been done in Ref. [38].
According to Eq. ( 3), the computation of the scattering amplitude is broken down to the computation of three objects: the connected n-point function G({p i }, {k j }), the propagator K(p), and the field normalization R. Implementing the computation of these objects on a quantum computer will involve three steps: (1) the spatial dimensions are discretized into a lattice, (2) the field degrees of freedom are mapped to qubits, (3) a suitable quantum algorithm is constructed to evaluate the three objects individually.For gauge theories, step (1) can be achieved in the standard way under the Kogut-Susskind Hamiltonian formalism [39,40], and alternative approaches have been proposed [41][42][43].Step (2) can be done straightforwardly for fermionic degrees of freedom [44][45][46], while for bosonic degrees of freedom and in particular gauge bosons considerable progress has been made [42,43,[47][48][49][50][51][52][53][54][55].In this work, we will focus on step (3), assuming that steps (1) and ( 2) have been achieved.It should be remarked that, in step (3), the three objects to be calculated could be ultraviolet divergent, meaning that their individual values blow up in the continuum limit.However, the scattering amplitude, as a physical observable, remains a finite constant when the continuum limit is taken.The large cancellation in the continuum limit among the components in the LSZ reduction formula could potentially cause problems on numerical stability in practical calculations.We leave the detailed study of the approach to the continuum limit of the LSZ reduction formula for the future.

III. THE QUANTUM ALGORITHM
Here we propose a quantum algorithm to compute the three objects involved in the LSZ reduction formula Eq. ( 3), namely the connected n-point function, the propagator, and the field normalization.Accordingly to Eq. ( 6), the field normalization R involves the field operator ϕ evaluated at x = 0 sandwiched between the vacuum and a single-particle state with zero spatial momentum (p = 0).Since no time evolution of the field operator is involved, the value of R can be readily determined once the vacuum and the single-particle state are obtained.To obtain the vacuum and the single-particle state, one can employ the quantum algorithm proposed in Ref. [31], which shows that both the vacuum and the single-particle state can be obtained efficiently with the quantum alternating operator ansatz (QAOA) and the quantumnumber-resolving variational quantum eigensolver [56][57][58][59].It should be noted that, only states with zero spatial momentum are involved in our formalism. 1Since these states are translation invariant, the QAOA can be applied easily: one simply uses input reference states and multidimensional quantum Fourier transform with complexity O( [62], which is negligible compared to O(nn q N n+1 T n ).To get the connected part of the n-point function, we also have to calculate its disconnected part, the evaluation of which according to what we have just derived has complexity Similarly, evaluating the n propagators in Eq. (3) costs O(nn q N 3 T 2 ) operations.Consequently, the overall complexity in our approach is O(2 n nn q N n+1 T n ), which is exponential in n.In the JLP formalism, it was shown that the complexity scales with n polynomially [10].Therefore, our approach is ideal only when the number of external particles is small, e.g. in 2 → 2 scatterings.
With the estimate of complexity above, we can see how the complexity scales with the largest energy scale in the scattering process, denoted by Λ max .In general, in order to satisfy constraint Eq. ( 1), both N and T are taken to be proportional to Λ max .For fermionic degrees of freedom, n q is fixed and does not depend on Λ max .For bosonic degrees of freedom, n q is chosen as a finite number which is large enough so that modes of energy or momentum of order Λ max are properly represented.It is found that n q can be taken as O(log Λ max ) for scalar fields [47] and gauge fields [48], if we estimate the maximum value of the field strength tensor as ∼ Λ 2 max .Therefore, the complexity scales at most as ∼ Λ 2n+1 max log Λ max , which is polynomial in Λ max .A polynomial dependence of the complexity on energy was also obtained in the JLP formalism [10,12].
An important feature of our approach using the LSZ reduction formula is that bound states are allowed as incoming or outgoing particles.This is because the interaction coupling constant is never turned off in our approach, as opposed to the JLP formalism.In Eqs. ( 4)-( 6), the field operator ϕ is not necessarily a fundamental field of the theory.In fact, any operator which has the same quantum numbers as the external particle h can be used.For instance, in a theory with only a spin-1/2 fundamental field ψ, there might exist a spin-0 scalar bound state h made of a fundamental fermion and its antiparticle.One can then simply take the composite operator ϕ = ψψ as the operator which annihilates h in the LSZ reduction formula for scattering processes involving h as external particles.This is an ideal feature of our formalism, since in the most interesting potential application of quantum computing in particle physics, namely quantum chromodynamics, all incoming and outgoing particles are bound states owing to quark confinement.Our framework is therefore most useful for scattering processes involving a small number of composite particles in a strongly coupled theory, such as 2 → 2 scatterings of pions or nucleons in quantum chromodynamics, for which first-principle calculations are currently only possible below the three-hadron thresholds [63][64][65].

IV. POLOLOGY IN THE ONE-FLAVOR GROSS-NEVEU MODEL
With simulations using the proposed quantum algorithm on classical hardware, we can demonstrate the feature of poles of the propagator and the connected four-point function in a simple model, the (1+1)dimensional one-flavor Nambu-Jona-Lasinio model [66,67], also known as the one-flavor Gross-Neveu model [68].The Lagrangian of this model is given by where ψ is a Dirac field in 1+1 dimensions, which we will refer to as the quark field, and m q and g are the bare quark mass and the bare coupling constant, respectively.Both m q and g are free parameters.We choose the values of them in such a way that the particle states h we consider below have masses m h satisfying π L < m h < π a , where a and L are the lattice spacing and the lattice size, respectively.In the following, we take m q a = 0.84 and g = 0.40.
Consider the propagator of the quark field, Similar to Eq. ( 7), near the mass shell of a particle state h which has the same quantum numbers as the quark field, we have where R h > 0 is the field normalization defined by where u(p) is the positive-energy solution to the free Dirac equation ( / p − m h )u(p) = 0, normalized to ū(p)u(p) = 2m h .The field normalization R h can be calculated using our proposed method by taking p = 0 in Eq. (11).
According to Eq. ( 10), when we take p 1 = 0 and treat K ψ (p) as a function of p 0 , K ψ (p) should have poles at p 0 = ±m h .We evaluate K ψ (p) as a function p 0 , with p 1 taken to be zero, with the proposed quantum algorithm using classical hardware.The calculation is performed on a desktop workstation with 16 cores, using opensource packages QuSpin [69] and projectQ [70], with 14 qubits (seven lattice sites), using seven temporal sites, with the temporal spacing taken to be the same as the lattice spacing.We follow the mapping of the Gross-Neveu model onto qubits and the method to evaluate the correlation function with a quantum algorithm as discussed in Ref. [31].To obtain K ψ (p), we first calculate the integrand in Eq. ( 9) in position space and then implement a discrete Fourier transform to perform the Fourier integral.In order to show the main features of our result, we present our result for TrK ψ (p), which according to Eq. ( 10) takes the following form near a particle mass shell: Equation (12) shows that Re [TrK ψ (p)] → +∞ and Im [TrK ψ (p)] → 0 as p 2 → m 2 h .Figure 1 shows our results for the real part (solid line) and the imaginary part (dashed line) of TrK ψ (p) as a function of p 0 a.The peaks of the real part at p 0 a = ±1.14(solid blobs) correspond to the poles from the lowest-lying state with the same quantum numbers as the quark field, as is verified by solving for the mass spectrum with direct numerical diagonalization of the discretized Hamiltonian, which gives m h a = 1.18 (±m h a shown by dotted vertical lines).This state can be interpreted as a quark. 2 In the continuum limit, a pole corresponds to a peak of infinite height, while in the discretized model we consider here the peaks have finite height.
According to Eq. ( 12), when viewed as a function of p 0 with p 1 = 0, Im [TrK ψ (p)] has zeros at p 0 = m h and p 0 = −m h , at which points the derivative of Im [TrK ψ (p)] with respect to p 0 is positive and negative, respectively.These features are exhibited in our result for the imaginary part of TrK ψ (p), shown by the dashed line in Fig. 1.In Fig. 1, the dashed line has positive and negative derivative at the zeros at p 0 a = 1.28 and p 0 a = −1.28(hollow circles), respectively.These zeros are close to ±m h a, and are expected to be there by virtue of Eq. ( 12) in the continuum theory.
For the 2 → 2 scattering of a quark and an antiquark, q(k 1 )q(k 2 ) → q(p 1 )q(p 2 ), one has to calculate the following connected four-point function: Similar to K ψ (p), the connected four-point function h +iϵ when k 1 is close to the mass shell of particle h which has the same quantum numbers as the quark field.To demonstrate this feature, with the same method as we evaluate K ψ (p), we evaluate G αβγδ ψ (p 1 , p 2 , k 1 ) as a function of k 0 1 with the external-leg momenta set to k 1 = (k 0 1 , 0), p 1 = (0, 0), p 2 = (k 0 1 , π/a).This setting of external-leg momenta makes sure that k 2 , p 1 , p 2 are off shell. Figure 2 shows our results for the real part (solid line) and the imaginary part (dashed line) of G αβαβ ψ (p 1 , p 2 , k 1 ) as a function of k 0 1 a.Similar to the case of the propagator, the peaks of the real part 2 Note that in this model there is no quark confinement.2, has positive and negative derivative, respectively, at the zeros at k 0 1 a = 0.87 and k 0 1 a = −0.87(hollow circles).These zeros, being close to ±m h a, are expected by the pole structure of the connected four-point function in the continuum theory.
In order to demonstrate the power of the LSZ reduction formula in handling scatterings of bound-state particles, we also simulate the propagator of the composite operator O(x) = ψ(x)ψ(x), given by Figure 3 shows our results for the real part (solid line) and the imaginary part (dashed line) of TrK O (p) as a function of p 0 a, with p 1 = 0.The peaks of the real part at p 0 a = ±2.02(solid blobs) correspond to the poles from the second lowest-lying state h O with the same quantum numbers as the vacuum, as is verified by solving for the mass spectrum with direct numerical diagonalization, which gives m h O a = 1.96 (±m h O a shown by dotted vertical lines).This state can be interpreted as a quark-antiquark bound state.As expected again from the pole structure of K O (p) in the continuum theory, the imaginary part has positive and negative derivative, respectively, at the zeros at p 0 a = 1.89 and p 0 a = −1.89(hollow circles), which are close to ±m h O a.Note that the peak of the real part at p 0 a = 0 does not correspond to any single-particle pole, since the pole structure due to such a massless state would imply that the imaginary part has a zero at p 0 a = 0 which is a local minimum, in contrary to the local maximum we found.
This simple example shows that our proposed quantum algorithm succeeds in recovering the expected pole structure of both the propagator and the connected npoint function, which is crucial to the implementation of the LSZ reduction formula.
It should be noted that, although we find that both the quark propagator and the connected four-point function exhibit convergent behavior when the number of qubits is increased to 14 from a smaller number, the speeds of their convergence are quite different, with the propagator converging faster.For 2 → 2 quark-antiquark scattering, since the scattering amplitude involves four powers of the ratio of peak values of the connected four-point function and the quark propagator, its convergence requires that both the quark propagator and the connected four-point function reach a similar speed of convergence.In our initial attempt to calculate this scattering amplitude, owing to the limited number of qubits one can use in simulations on classical hardware with reasonable computational time cost, a satisfactory convergence has not been observed.However, with the quantum advantage, in a quantum computer with more than a hundred qubits in near future, we believe that a convergent result for the scattering amplitude can be expected for a (1+1)-dimensional model.

V. CONCLUSIONS
In this work, we proposed a new framework for evaluating scattering amplitudes in quantum field theory on quantum computers in a fully nonperturbative way.The framework was based on the LSZ reduction formula, which relates scattering amplitudes to correlation functions.In this framework, as opposed to a direct Hamiltonian simulation of the scattering process, no preparation of wave packets of incoming particles is required, and one only has to prepare one-particle states of zero momentum.The framework is capable of incorporating scatterings of bound-state particles, and is ideal for scatterings that involve a small number of particles.This framework is expected to have potential applications in exclusive processes in a strongly coupled theory, such as 2 → 2 scatterings of pions or nucleons.As a proof of concept, in a simple model, the one-flavor Gross-Neveu model, we demonstrated by simulations on classical hardware that the propagator and the connected four-point function obtained from the quantum algorithm has the desired pole structure crucial to the implementation of the LSZ reduction formula.

FIG. 1 .
FIG. 1. Real part (solid line) and imaginary part (dashed line) of TrK ψ (p) in the one-flavor Gross-Neveu model as a function of p 0 a with p 1 = 0, simulated with the proposed quantum algorithm on classical hardware.The solid blobs on the solid line and the hollow circles on the dashed line, respectively, show the peaks of the real part and the zeros of the imaginary part of TrK ψ (p) corresponding to the pole structure due to the quark.The dotted vertical lines show the locations p 0 a = ±m h a with the quark mass m h obtained from exact diagonalization of the discretized Hamiltonian.

(p 1
, p 2 , k 1 ) at k 0 1 a = ±1.01(solid blobs) correspond to the poles from the lowest-lying state with the same quantum numbers as the quark field, namely the quark.The imaginary part of G αβαβ ψ (p 1 , p 2 , k 1 ), shown by the dashed line in Fig.

FIG. 2 .
FIG. 2. Real part (solid line) and imaginary part (dashed line) of G αβαβ ψ (p1, p2, k1) in the one-flavor Gross-Neveu model as a function of k 0 1 a, with k1 = (k 0 1 , 0), p1 = (0, 0), p2 = (k 0 1 , π/a), simulated with the proposed quantum algorithm on classical hardware.The solid blobs on the solid line and the hollow circles on the dashed line, respectively, show the peaks of the real part and the zeros of the imaginary part of G αβαβ ψ (p1, p2, k1) corresponding to the pole structure due to the quark.The dotted vertical lines show the locations p 0 a = ±m h a with the quark mass m h obtained from exact diagonalization of the discretized Hamiltonian.

FIG. 3 .
FIG. 3. Real part (solid line) and imaginary part (dashed line) of TrKO(p) in the one-flavor Gross-Neveu model as a function of p 0 a with p 1 = 0, simulated with the proposed quantum algorithm on classical hardware.The solid blobs on the solid line and the hollow circles on the dashed line, respectively, show the peaks of the real part and the zeros of the imaginary part of TrKO(p) corresponding to the pole structure due to the quark-antiquark bound state.The dotted vertical lines show the locations p 0 a = ±m h O a with the bound-state mass m h O obtained from exact diagonalization of the discretized Hamiltonian.