Real-time quantum calculations of phase shifts using wave packet time delays

We present a method to extract the phase shift of a scattering process using the real-time evolution in the early and intermediate stages of the collision in order to estimate the time delay of a wave packet. This procedure is convenient when using noisy quantum computers for which the asymptotic out-state behavior is unreachable. We demonstrate that the challenging Fourier transforms involved in the state preparation and measurements can be implemented in $1+1$ dimensions with current trapped ion devices and IBM quantum computers. We compare quantum computation of the time delays obtained in the one-particle quantum mechanics limit and the scalable quantum field theory formulation with accurate numerical results. We discuss the finite volume effects in the Wigner formula connecting time delays to phase shifts. The results reported involve two- and four-qubit calculations, and we discuss the possibility of larger scale computations in the near future.

1. Introduction. The idea of simulating quantum field theory with quantum computers has gained considerable interest recently [1]. In the context of high-energy and nuclear physics, a long term motivation is to develop quantum computing methods that perform real-time evolution for lattice quantum chromodynamics (QCD). This is an ab-initio, ultraviolet complete, theory of strong interactions which has been very successful in describing the static properties of hadrons and nuclei [2]. Importance sampling methods, used successfully today for lattice QCD at Euclidean time, are not effective for dealing with the rapid oscillations of real-time unitary operators acting on large Hilbert spaces. Currently, physicists resort to semi-empirical algorithms such as Pythia and Herwig [3,4]) to interpret hadronic collider data. Doing ab-initio calculations for jet physics would represent a major accomplishment and remains a long-term goal for the high-energy community and strategies to deal with parton distributions are outlined in Refs. [5,6]. Related methods for out-of-equilibrium processes in many-electron systems and information paradoxes in quantum gravity [7] would also have a large potential impact.
Quantum computers offer an alternative solution to the sign problem plaguing current efforts. A sensible approach would be to follow the sequence of models that has been successful for the development of lattice QCD at Euclidean time on classical computers [8,9]. The first step in this sequence is to study and understand the properties and behavior of the quantum Ising model (QIM) on today's Noisy Intermediate-Scale Quantum (NISQ) hardware.
Real time evolution involving a limited number of sites for the QIM has already been attempted using a few qubits on gate based quantum computers [10][11][12][13][14][15][16][17][18][19][20][21][22][23], as well as developments in progress for more complicated models [24][25][26][27][28]. Today researchers must grapple with discrepancies between the measured qubit-state occupations on quantum platforms and the exact evolution after only a few Trotter steps. Currently, error-mitigation methods such as zero-point extrapolation [27,29] are necessary to improve the computational results. It has been shown [13] that by modelling four qubits on an IBM Q quantum computing hardware platform these mitigation methods provide a reasonable extrapolation for times of the order of the approximate periodicity of the problems considered. It is possible to use Trotter steps that are more than 20 times larger than rigorous bounds in (δt) 2 or (δt) 3 would naively suggest, without encountering major discretization errors [13,30]. This allows us to reach larger time scales for studying scattering processes.
In this Letter, we demonstrate that it is possible to use state-of-the-art NISQ devices to prepare and evolve suitable wave packets for the QIM. We show that it is possible to project the wave-function in the early stages of a collision process onto momentum states and to pinpoint a time that corresponds to the middle of the collision. By introducing an extra interaction, we obtain a time delay ∆t illustrated in Fig. 1 that is half of the time delay ∆t W invoked in Wigner formula [31] given in Eq. (6) to estimate the derivative of the phase shift with respect to the momentum.
Phase shifts are a key measurement in the scattering process and represent the total change of phase due to interactions. Significant progress has been  Figure 1. Illustration of the measurement of the time delay between the free and interacting wave packets. The renormalized reflection probability R−(t) is defined in Eq. (13) made in calculating them from lattice QCD in Euclidean time [32][33][34][35][36][37]. In standard textbooks, phase shifts and scattering amplitudes are estimated from asymptotic data long after the collision processes have occurred. However, for NISQ devices with limited coherence time or gate-depth, using the information from the early stages of the collision is advantageous. We show that this idea can be implemented on both a quantum computer using superconducting transmon qubits and a trapped ion system operating at the University of Maryland [38].
The article is organized as follows. We introduce the quantum Ising model with an extra interaction and its Hilbert space. We show that it is possible to implement the three steps of the calculation of the S-matrix elements: 1) preparation of the initial state, 2) real-time evolution, and 3) measurement of the probability for a particular final state. We then extract the phase shift by comparing the cases with and without an external potential. Steps 1) and 3) involve Fourier transforms and are very challenging with NISQ devices. This is why we first restrict ourselves to the quantum mechanics limit where oneparticle states interact with an external potential localized at one site. We then show that it is possible to extend the computations to the case of the quantum field theory formulations [12] that require more qubits but are guaranteed to scale efficiently for larger volume [39].
2. The model. We consider the transverse-field Ising model in one spatial dimension, This model is very well understood [9] and discussed for NISQ devices [10][11][12][13]. It is equivalent to a theory of free fermions with subtle effects from the boundary. Non-trivial interactions can be introduced with an extra termĤ In order to perform Fourier transforms with a reduced number of qubits we first consider the quantum mechanics limit J h T , where the model consists of energy bands that can be assigned a particle number. This amounts to neglecting terms of the formσ + iσ + i+1 and their conjugates. For the reduced one-particle problem, the interaction introduces an effective potential at the end of the chain. This reduces the size of the Hilbert space from 2 N to N and also allows approximate analytic calculations [12].
Specializing to the case N = 4, we have, up to an unimportant additive constant, the following effective Hamiltonian matrix: This allows us to reduce the Hilbert space from four qubits needed for the N = 4 field theory problem to two qubits for the one particle limit. The remapping is shown in Eq.( 4) with the correspondence illustrated in Fig. 2: The Hamiltonian in Eq. (3) can now be written aŝ The subscripted roman numerals are used to indicate the use of our two-qubit decomposition.
3. Real-time scattering. The time delay ∆t W of a wave packet with a sharply defined momentum k is related to the derivative of the phase shift by the Wigner formula [31]: where ∂E/∂k is the group velocity, which in our case is 2J sin(k). We will explain that ∆t W can be estimated from the first half of the real-time evolution of the scattering process and that it is actually twice the time delay illustrated in Fig. 1. U 00 01 10 11 Figure 2. Visual depiction of the qubit states (black dashes), potential for the interacting (blue) and noninteracting case (dark red), and initial wave packet (light red).
We now report quantum computations for J = 0.02 and U = 0.03. In the 2-qubit Hilbert space, the momentum states for k = ±π/2 read It is necessary for the initial wave packet to have some localization in space so that a distinct scattering event is visible, i.e., there is a point in time when the particle reaches the interaction region. As a side effect, the wave packet will have some momentum distribution because it no longer is a plane wave. As depicted in Fig. 2, we chose our initial wave packet to be |π/2 but restricted it to be non-zero only in the middle: We construct this wave packet on two sites with the following quantum circuit with all qubits initialized in the |0 state: (9) The time evolution operator can be written as a combination of XX, YY, X, rotations and a controlled phase operation (R φ ): where ρ = Jδt, θ = U δt, and δt = 12.5. We use standard notations [40] for the gates. The very slow growth of the one-step error for large δt [30] allows us to reach t = 75 with only six Trotter steps [13].
We then perform a quantum Fourier transform (QFTr) on these two qubits to take this state into momentum space: After applying the QFTr, the qubit states |10 and |11 correspond to the momentum states |k and | − k respectively, with k = π/2. We define the probabilities to be in the |±k state and their renormalized versions which by design satisfy The real-time evolution provides the time t necessary to reach the symmetric situation where P + (t ) = P − (t ). This corresponds to the time where a classical particle would hit the barrier. We can then compare t in the case where U = 0 and U = 0.03 and we call these times t f ree and t int. respectively. In order to determine the difference The relation to the time delay ∆t W used in Eq. (6) is supported by numerical calculations provided later. This relation can also be justified from the timereversal argument that after t int. only half of the phase shift, δ(k), has built up while the other half builds after t . This is why historically, the total phase shift is denoted 2δ(k). The time t is determined by the symmetric condition R − (t ) = R + (t ) = 0.5. This requires some continuous interpolation among the discrete values at the Trotter steps. We find that the sigmoid parametrization provides very good fits of the numerical data and satisfies the sum rule of Eq. (14) because of its reflection property. This also implies that the time reversal with respect to t amounts to the momentum reversal. The parameter w describes the width of the sigmoid. The data for R − obtained from both quantum computers are shown in Fig. 3 together with the sigmoid fits where the values of t f ree and t int. are indicated by vertical lines crossing the fits at the 0.5 horizontal line. This provides the numerical values of ∆t given in Table I. For comparison, we give the values obtained by doing sigmoid fits of the continuous-time evolution (first column) and the Trotter steps (second column) calculated numerically at the same discrete times as the experimental data. The systematic errors are expected to be larger than the statistical errors and very difficult to estimate. Details of the fits are given in the Supplementary Material (SM).
We see that both the IBM and trapped ion estimates provide larger absolute values of ∆t * than the target values. This can be in part explained by the fact that the fits for the free process tend to lag below the Trotter steps for t > 50 indicating a loss of coherence.
Measurements from the IBMQ Bogota machine contain both the noisy data with just readout corrections and a mitigated version obtained using methods discussed in Ref. [13] and which account for some slightly negative occupations at low t. The trapped ion simulations include only readout corrections without noise mitigation (see SM).
We see that the quantum mechanics approximation allows us to perform the QFTr and get reasonable estimates of ∆t * , (Table I). We expect to improve the accuracy of these estimates in the near future. The extension of this procedure for more than four sites requires an all-to-all connectivity and a CNOT depth increasing with the number of sites. In contrast, the field theory calculation discussed in the next section, and which is our ultimate goal, requires more qubits but remains local [39] with a constant CNOT depth.
4. Towards scalable field theory calculations. We have performed field theory computations with the same devices using four sites and four qubits. This allows for shallower local circuits but requires a more expensive Fourier transform [41] as discussed in the SM. In Fig. 4    ion trap machines using the full Hamiltonian with a Trotter step of δt = 15 with six Trotter steps. In both cases, unmitigated data is used for the analysis. Fig. 4 indicates that the discrepancies between the individual fits and the Trotter exact results are more pronounced than in the previous calculation. However the effects somehow compensate when we calculate the differences and we obtain time delays closer to the Trotter exact ones as shown in Table I. As discussed in the SM, the fit methods are similar to the previous ones. As improved quantum computing hardware platforms become available, we plan to use these upgraded facilities to get more accurate values of ∆t and extend these calculations to 8 and 16 sites with one qubit per site. 5. Extracting the phase shifts from time delays. We are now in position to discuss how to use Eq. (6) to estimate the phase shifts. First, we need to say that the wave packet in Eq. (8) is not narrow and volume effects should be significant for N = 4. In order to estimate these effects we consider the one particle quantum mechanical problem at larger N with narrow Gaussian wave packets. The details are given in the SM, where we also show that in the infinite volume limit .
The results at finite and infinite volume are compared in Fig. 5 illustrate the magnitude of the volume effects for N = 32, 64 and 128. Some integration procedure is needed in order to extract δ(k). In the quantum mechanics limit, we noticed that for small enough k, δ(k) could be approximated very well by 2E(k)∆t . It is shown in the SM that this feature can be explained from the Taylor expansion of the exact formula. This provides initial data that is reliable for k ≤ 0.5. For larger k an integration procedure based on Eq. (6) needs to be used.

Conclusions.
We have proposed a method to estimate the time delay using the early steps of the real-time evolution. We have given a proof of principle that actual computations of the time delays can be implemented on both the IBM superconducting transmon and trapped ion hardware platforms for a quantum mechanics limit with two qubits and the field theory formulation with four qubits. There is plenty of room for optimization for both devices and we do not claim that our results allow a systematic comparison between the devices. We expect that the field theory calculations should feasible for a larger number of qubits in the near future. A detailed comparison with existing real-time methods in one spatial dimension [42][43][44][45] would be of great interest. Quantum computations for quantum Ising models in two spatial dimensions could offer the possibility to reach quantum advantage.

Quantum Computing Hardware Platforms
This project used both the IBM Quantum Network superconducting transmon hardware platforms and a trapped ion machine. For the IBM quantum computers, the project used IBM Q Bogota (5 qubit) and IBM Q Casablanca (7 qubit) machines. Both IBM platforms are rated with a quantum volume of 32. Access to these machines was through the IBM cloud via the NC State IBM Quantum Hub.
The ion trap quantum computer used in this study consists of a chain of 171 Yb + ions confined in a linear Paul trap and laser-cooled to their motional ground states [38,46]. The qubits are defined in the hyperfine-split 2 S 1/2 manifold as |0 = |F = 0, m F = 0 , |1 = |F = 1, m F = 0 , with a splitting of 12.642821 GHz and are insensitive to magnetic field fluctuations to first order. Qubits are initialized to the |0 state by optical pumping [47]. Coherent operations are achieved by illuminating the ions with a pair of counter-propagating Raman beams at 355-nm which have a beat-note that is set resonant with the qubits or near the motional sidebands of the chain [48]. One of the Raman beams illuminates the entire chain. The other one is split into individual beams that are each matched onto one channel on a multi-channel acoustic-optic modulator (AOM). The amplitude, frequency, and phase of each beam can be modulated independently by the correspond-ing RF signal. Individual addressing is achieved by focusing each of these beams onto one single ion. Each ion is also matched onto a distinct channel of a photo-multiplier tube (PMT) that collects statedependent fluorescence for individual state readout [47]. There are two mechanisms of quantum control. Single qubit rotations are performed by driving Rabi rotations of the target ion on resonance. Twoqubit entanglement is implemented by creating an effective Ising spin-to-spin interaction via transient entanglement between the qubits and the motional modes with the Raman beat-note tuned near the motional sidebands [49][50][51]. Only qubit-to-qubit entanglement remains after the motional modes are dis-entangled at the end of the scheme [52]. Entanglement can be created between any pair of ions owing to the long-range Coulomb interaction [53]. Typical gate fidelities are around 99.5(2)% for the single-qubit gates and 98(1)% for the two-qubit entangling gates. The main sources of two-qubit gate errors are residual entanglement between the motional modes and the qubits due to laser intensity fluctuations and motional heating. State preparation and measurement (SPAM) errors are accounted for by applying the inverse of an independently measured state-to-state error matrix when analyzing the data.

Determining time delay with sigmoids
Using the sigmoid for the fit of the quantities R − (t) we find the time delay. We have separated out each of the various data sets into lines. The continuous sampling uses 100 equally spaced time steps between t = 0 and t = 75 with a corresponding uncertainty of 0.001. Similarly we also used 6 equally spaced data points sampled from the continuous distribution at time which are multiples of t = 12.5 with the same flat error rate of 0.001. The Trotter-exact was sampled at times steps of δt = 12.5 and had an uncertainty of δP ± (t) = P ± − P 2 ± (t)/ √ 1000. The remaining three lines were sampled by fitting to the first five data points for the interacting case (int) and six data points for the free case. This was done to make sure that the value of R − (t) is at least 0.5 and that the system has not reached an asymptotic state. Stated in a different way, in the interacting case, the signal flattens below 1 after 5 Trotter steps due to the reflection on the left wall. The last data point was not used for the fit and is always below the fit.
The systematic errors are significant and very difficult to estimate. It is possible to get an idea of their magnitude by looking at the first time step where a detailed examination shows that the measured values of R − are significantly larger than the Trotter exact values and the estimated statistical errors. In absence of a reliable model for the systematic errors, we estimated the error on the parameters using the jackknife method.
The fitting methods for the field theory results are similar, the only difference being a slightly larger Trotter step.

QFT circuit
Initially we want to scale up from the quantum mechanics picture to the field theory perspective, in other words the full Ising model. As in the quantum mechanical case, we will need to break up this circuit into three parts: 1) state Preparation, 2) Trotter Evolution, 3) momentum Projection. The various circuits and the relative costs are discussed below.
We want to prepare the same initial state that we proposed in the main text mapped on to the Ising Model. This initial state that we prepared: can be remapped to the Ising model by using the inverse mapping originally described in the text: This remapping will give us the initial state: This state can easily be prepared with 1 XX gate  (2) 59(1) -13(2) 18(1) 14.2(9) Trapped Ions 69(1) 57(1) -12 (2) 17(1) 13(1) and one single qubit rotation. The time evolution of the system can be easily written in terms of a circuit requiring at most 3 XX gates per Trotter step when using a four site system and this will be 2 XX gates deep. This circuit will end up being an application of entangling rotations followed by a set of single qubit phasing operations.
The final portion of the quantum simulation is a Fourier transformation to take the circuit into mo-mentum space [54,55]. This transformation can easily be done using four two-qubit operations.
Where F is given by the matrix: This can be decomposed into at worst 4 CNOT gates. It should be noted this is different from the traditional quantum Fourier transform which uses binary encoding for the qubits. The measurement is done then in the z-basis and we select out the states |0010 and |0001 for the |k and | − k state respectively.

Phase Shift Derivation
The discrete Schroedinger equation obtained in the limit of small J, with U = 0, and no boundaries admits plane wave solutions e ±ikx with energy 2J(1 − cos(k)). The additive constant has been adjusted in order to have zero energy at zero momentum. If we introduce a right wall and impose ψ(N + 1) = 0, we need a relative phase −e i2k(N +1) for the reflected wave e −ikx . If, in addition, we also introduce a repulsive potential U > 0 at the rightmost site N , the problem with a right wall can be solved by introducing the phase shift. A short calculation shows that: We are now in position to make two independent calculations of ∆t . In the first one, we use the gen-eral Eqs. (6) and (15) together with the exact expression for δ (k) from Eq. (S8) to get: ∆t * = (−1 + J U cos(k) + J U 2 + J 2 + 2JU cos(k) )/2J sin(k).
(S9) The other method consists in using the real-time evolution to calculate R − (t) and determine the time t where R − (t ) = 0.5. This second method depends on the volume and the details of the initial wavepacket. It is important to prepare a wave packet that is broad enough in space to have sharply defined momentum but not broad enough to "bounce back on the left wall" too early. In the numerical calculations, we used an initial wavepacket ψ(x) ∝ e ikx−(x−N/2) 2 /(N/4) 2 . (S10) The comparison between these results and Wigner formula are shown in the text. A few words of explanation regarding the low k approximate formula for δ(k). From Eq. (S8), we can expand δ(k) in power of k. For our numerical values J = 0.02 and U = 0.03, we obtain δ(k) −0.6k + 0.008k 3 + . . . (S11) As we see, for 0 ≤ k ≤ 0.5 the correction to the linear approximation are less then one percent and we can in good approximation use δ(k) δ (k)k. (S12) [1] J. Preskill, in Proceedings of Science, Lattice Conference 2018 (2018), URL https://pos.sissa.it/ 334/024/pdf.