Resource-efficient simulations of particle scattering on a digital quantum computer - Nature

IntroductionScattering experiments are at the heart of unraveling the internal structure of matter and the interactions between the fundamental particles. Major experimental facilities such as LHC1 and RHIC2 continue to generate valuable experimental data to test the theoretical predictions and drive the search for new physics. Concomitantly, there are significant efforts to develop analytical and numerical methods for improving our understanding of gauge field theories, which provide the theoretical framework for particle physics.In particular, Lattice Field Theory (LFT) provides a powerful tool for exploring non-perturbative regimes from first principles. Discretizing a theory on a Euclidean space-time lattice allows for applying sophisticated Monte Carlo (MC) methods that have been extremely successful for computing properties such as mass spectra, phase diagrams and many other static properties3,4. However, the conventional MC approach to LFT is not suited to directly explore dynamical problems, as it crucially relies on the formulation in Euclidean space-time, and using Minkowski space-time would lead to a sign problem preventing efficient MC sampling. While indirect approaches exist to address scattering problems with Monte Carlo methods, such as extracting scattering phase shifts using Lüscher’s method5,6, these techniques become increasingly challenging to apply at high-energies or for inelastic scattering processes7,8,9. In addition, such indirect approaches do not provide detailed access to the real-time dynamics of the particles during or immediately after a collision. Hence, there is an ongoing effort to find alternative methods allowing for overcoming these limitations.An alternative approach for the simulation of scattering processes is tensor network techniques10. Numerical algorithms based on tensor networks take advantage of the classical simulability of slightly-entangled states11, allowing for an efficient representation of the wave function in terms of low-rank tensors. In addition, numerical algorithms based on tensor networks do not suffer from the sign problem and have been demonstrated to allow for reliable calculations in regimes which are inaccessible with conventional MC methods12,13,14,15,16,17,18,19,20,21. In particular, refs. 22,23,24,25 simulated meson scattering within the Schwinger model in both weak and strong coupling regimes using Matrix Product States (MPS). However, the substantial growth of entanglement following the collision in certain parameter regimes necessitates considerable computational resources to achieve accurate results. This poses a significant challenge for tensor networks, a difficulty that can become even more pronounced when trying to extend the approach to higher-dimensional settings26,27.Quantum computing facilitates the simulation of real-time dynamics by constructing quantum circuits that approximate the time-evolution operator, e.g., through the use of Trotterization28,29,30. Several proposals have been made to extract asymptotic scattering observables from such simulations31,32,33. Quantum computers also have the potential to simulate the full dynamics of the scattering process, allowing for capturing intermediate snapshots at any stage of the evolution34. A central challenge is the preparation of the initial state, which should represent separated particle wave packets in position space with well-defined momenta. Recent works have addressed this and demonstrated simulations and first implementations on quantum processors35,36,37,38,39,40,41,42,43,44,45,46. However, the simulation of dynamics up to long time-scales generally requires deep quantum circuits, making it difficult to obtain accurate results on current and near-future quantum hardware due to the levels of noise. This necessitates the development of novel methods that reduce circuit depth, while preserving fidelity, to fully leverage the potential of current quantum devices.In this work, we seek to make the most of the combined strengths of quantum and classical simulation by leveraging the observation that tensor network methods can simulate the initial state and the dynamics of fermionic scattering processes at short times. In particular, we employ tensor network methods to variationally optimize quantum circuits47, enabling an automated search for shallow-depth unitary blocks that are more compatible with current quantum devices. An increasing body of work has explored the use of tensor network methods to support quantum simulations across a range of applications via variational circuit compilation. These efforts include ground state preparations48,49,50, short-time state evolution51,52,53, and the construction of short-depth circuits that approximate the unitary evolution operator e−iHt governed by a Hamiltonian H54,55,56,57,58,59.Here we demonstrate our scalable tensor network compression strategy on the interacting Thirring model. That is, rather than performing the full real-time evolution on hardware, we perform MPS simulations to compute the early-time dynamics and identify low-entanglement states. These states are then variationally compiled into short-depth circuits and used as starting points for hardware execution. To simulate the remaining real-time evolution, we apply circuit compression to the time-evolution operator as well. Specifically, we approximate short-time evolution segments e−iHt using variational circuits, where the target operator can be expressed as a matrix product operator (MPO) with high accuracy for small t. The variational circuit is optimized to maximize the overlap with the target time-evolution unitary, represented by this MPO. Compared to standard Trotterization, the resulting circuits are shallower in depth, enabling efficient execution of the full scattering process on noisy quantum hardware. Figure 1 provides an overview of our approach.Fig. 1: Overview.Full size imagea The target quantum simulation is a scattering interaction. Initially, spatially separated wave packets only show entanglement (indicated by the brightness of the color) within each wave packet. As they move towards each other during the Hamiltonian evolution and interact, entanglement between the wave packets is generated. This may persist in long-range as they separate, motivating the use of the quantum computer to represent this state with increased entanglement and growing tensor-network cost, despite the weakly entangled initial wave packet configuration. To aid the simulation on near-term quantum hardware, we exploit the low-entanglement early dynamics by simulating these times classically with tensor networks up to their limit; then the state is passed over to the quantum computer for further simulation. b The hybrid use of tensor networks and the quantum computer is enabled by MPS-based circuit optimization, to find efficient and short-depth circuits for performing the simulation. This involves generating two circuits: V1(θ) performs the initial state preparation of \(\left\vert \psi ({t}_{0})\right\rangle\), which has been pre-computed as an MPS; to continue the time evolution, we exploit the low-entangling property of e−iHt for small t, to learn a circuit V2(θ) with a shorter depth than available Trotterizations. c These circuits realize the same simulation as alternate conventional methods with shorter depth circuits, and we find up to a total factor of 3.2 reduction, significantly improving the amenability to near-term quantum hardware. d The combination of tensor network circuit optimization, performant quantum hardware, and modern error mitigation tools, enables a high-fidelity simulation of the full scattering dynamics on 40 qubits.The remainder of this paper is organized as follows. In the Methods section, we briefly review the lattice Thirring model and describe the construction of fermionic wave packets used to simulate scattering processes. The use of MPS pre-computation is described in the Methods section, and the method used for variational circuit optimization for state preparation and time evolution is further detailed in the Supplementary Information (Supplementary Note 1). We motivate the scalability of our MPS-based approach with circuit compilations on up to 160 qubits. The Results section presents our results demonstrating the techniques on quantum hardware. We find our methods enable the target simulation to be realized with a 3.2 times reduction in circuit depth compared to the conventional approach, allowing the full scattering dynamics on 40 qubits to be simulated on quantum hardware with high-accuracy. We further give a hardware demonstration of the TN compressed state at t0 = 10 on 80 qubits. Finally, in the Discussion we summarize our findings and discuss possible directions for extending this approach to more complex scattering processes and to other lattice field theories.ResultsWe execute the resulting tensor network optimized circuits on an IBM quantum computer for different parameter sets to demonstrate that our approach can accurately capture different physical behaviors. The experiments were run on IBM Quantum’s ibm_fez, a Heron r2 processor with 156 fixed-frequency transmon qubits with tunable couplers on a heavy-hex lattice layout (see Fig. 2). Here, we present device properties at the time of experimentation. The mean readout error was 1.2% and the median was 1.0%. A higher mean is indicative of an asymmetry that leads to higher qubit readout errors skewing the overall distribution of errors. The relaxation time, T1, had a mean of 149 μs and a median of 140 μs, and dephasing times, T2, with mean and median times of 103 μs and 98 μs, respectively. The single-qubit gates had a mean error rate of 1.4 × 10−4 and median of 6.4 × 10−5. Lastly, the two-qubit gate errors had a mean of 3.7 × 10−3 and a median of 3.4 × 10−3. In our experiments, a linear chain of qubits was chosen on the heavy-hex lattice trying to exclude individual qubits with low T1 and T2 times as well as large error rates (c.f. Fig. 2).Fig. 2: Layout of ibm_fez.Full size imageThe circles indicate the 156 qubits, lines connecting circles indicate qubit pairs between which a CZ gate can be directly implemented. For all of our experiments, we choose a linear chain of qubits on the heavy-hex lattice depending on the current error rates of the chip. The qubits highlighted in red indicate the linear chain used for the experiment with N = 40 and (m, g) = (0.4, 0.7) at t = 26. The qubits highlighted in blue are those used for the N = 80 initial state preparation experiment with (m, g) = (0.2, 0.4). The qubits highlighted in violet were used in both experiments.Simulating dynamics for 40-qubitsIn this subsection, we demonstrate hardware runs of the scattering process for 40-qubit systems. Three representative cases are considered: (m, g) ∈ {(0.2, 0.4), (0.4, 0.5), (0.4, 0.7)}, corresponding respectively to cases where the fermion and antifermion pass through each other, exhibit partially repulsive interaction and experience strong repulsion. In all cases, we chose the state \(| \psi ({t}_{0})\left.\right\rangle\) shortly before the collision as the starting state for the hardware run, corresponding to t0 ∈ {11, 18, 16}, as mentioned in the Methods section, Overview of Variational Circuit Approximation. These states have relatively low entanglement compared to the states after collision, and can be approximated by shallow, trainable quantum circuits, as outlined in the previous section. The subsequent timesteps are executed on the quantum hardware by applying the compressed quantum circuit for evolving the system for t = 2 to the starting state (see Table 1 for the exact gate numbers).Table 1 Comparison of two-qubit gate depth in terms of distinct CNOT layers and total CNOT gate number using the tensor network optimized circuit versus the conventional approach (in brackets)Full size tableThe Qiskit Estimator primitive is used to evaluate expectation values along the dynamics. This class provides methods to perform an array of error suppression and mitigation techniques in a unified way, making it straightforward to estimate expectation values of observables. We configure an Estimator with dynamical decoupling (DD) using the XY4 sequence to suppress decoherence on idle qubits, while Twirled Readout Error Extinction (TREX)60 is used to mitigate readout errors. In addition, Pauli twirling61 is applied to two-qubit gates to transform coherent errors into stochastic Pauli noise, allowing for more robust extrapolation to the zero-noise limit when estimating physical observables. For each circuit, we generate 200 random twirled instances and perform 1000 measurements per instance to estimate the expectation values of Pauli-Z operators. To further mitigate errors, we use zero-noise extrapolation (ZNE). Specifically, for a given set of parameters, we perform hardware runs at multiple noise factors G, and extrapolate the results to the zero-noise limit G → 0 (see the Supplementary Information (Supplementary Note 4) for details).Figure 3 displays the fermion density extracted from the Pauli-Z expectation after ZNE, where we also subtract the contribution from the vacuum to highlight the wave packet’s distribution$$\Delta {\langle {\xi }_{n}^{\dagger }{\xi }_{n}\rangle }_{t}=\left\langle \right.\psi (t)| {\xi }_{n}^{\dagger }{\xi }_{n}| \psi (t)\left.\right\rangle -\left\langle \right.\Omega | {\xi }_{n}^{\dagger }{\xi }_{n}| \Omega \left.\right\rangle .$$ (1) The fermion density in the vacuum state in the above expression is calculated using the MPS results.Fig. 3: Hardware run results for single time slices.Full size imageFermion densities for (m, g) = (0.4, 0.7) at times a t = 18 and b t = 26, comparing quantum hardware runs with precise MPS simulations. Gray triangles show unmitigated hardware results, orange dots show results after ZNE, and blue diamonds represent MPS data. For clarity, markers are filled for even sites and left empty for odd sites. For both panels, the x-axis represents the site index, and the y-axis represents the fermion density with the vacuum contribution subtracted.As shown in the figure, the unmitigated hardware results deviate significantly from the ideal MPS simulations, especially at t = 26, which requires a deeper circuit and more CZ gates. After applying the ZNE, the corrected results agree well with the ideal MPS data across all sites and both time slices. A more detailed discussion of the hardware difficulty associated with deeper circuits, including a simple fidelity estimate for the conventional implementation and a depth-sensitivity analysis based on folded circuits, is provided in the Supplementary Information (Supplementary Note 6).Finally, the full scattering dynamics for all three cases are shown in Fig. 4.Fig. 4: Dynamics of the scattering process.Full size imageData for times slices above the dotted horizontal line is obtained from simulations on quantum hardware, while those below are computed using MPS simulations. a (m, g) = (0.2, 0.4), data are from MPS simulation for t≤11, and are from hardware run for time t > 11. b (m, g) = (0.4, 0.5), data are from hardware run for time t > 18. c (m, g) = (0.4, 0.7), data are from hardware run for time t > 16. For a comparison with the exact simulation values (computed via tensor network methods) see Fig. S5 (Supplementary Note 5).For time slices before the collision, the results are obtained from the MPS simulation, while for time slices after the collision, the dynamics are executed on the quantum hardware with a timestep of t = 2. Based on the real-time fermion-density evolution, we observe qualitative signatures of different scattering behaviors in the three cases. In Fig. 4(a), the density profiles are consistent with transmission-like behavior, with the fermion and antifermion wave packets predominantly passing through each other after t ~ 11. In Fig. 4(b), the larger mass leads to lower velocities and hence a later collision time; the outgoing density profiles contain both fermion and antifermion components, suggesting qualitatively a mixture of transmission-like and reflection-like behavior. For the stronger interaction in Fig. 4(c), the density evolution is instead more consistent with predominantly reflection-like dynamics after the collision. Detailed fermion-density distributions for each time slice are displayed in Figure S5 (Supplementary Note 5).Our results thus show that using tensor network circuit approximation techniques together with a quantum device allows for obtaining a comprehensive picture of the scattering process. The circuit required to prepare the moderately entangled state at an initial time t0 can be efficiently compiled using classical methods. This enables the preparation of the state and its subsequent evolution on a quantum device, thereby entering a regime in which larger amounts of entanglement are generated and classical methods ultimately encounter limitations.Scaling up to larger system sizes: state preparation for 80 qubitsIn order to demonstrate that our approach can also be scaled up to larger system sizes, we benchmark the performance of the quantum hardware preparing \(| \psi ({t}_{0})\left.\right\rangle\) for a 80-qubit system. Specifically, we consider (m, g) = (0.2, 0.4), and take the MPS-simulated state at t = 10 as the reference. Similarly to the previous section, this state is then approximated by a parametrized circuit with a two-qubit depth of 24, using the tensor network circuit approximation approach from the Methods section, Overview of Variational Circuit Approximation. Executing this circuit on hardware allows us to assess the capabilities of larger-scale scattering simulations on current quantum hardware.The experiment is carried out on the ibm_fez device using Qiskit’s Sampler primitive. Measurement error mitigation via twirling is employed, along with XY4 dynamical decoupling to suppress decoherence. Again, ZNE is employed to mitigate the effects of noise, where we use the noise factors G = {1, 3, 5, 7}, corresponding to depths of two-qubit gates {24, 72, 120, 168} and total two-qubit gate counts {948, 2844, 4740, 6636}.Figure 5 compares the results from quantum hardware with the ideal MPS simulation. Even without error mitigation, the raw data from the quantum device qualitatively reproduces the shape of the wave packet. Applying ZNE, just as before, most of the hardware data aligns more closely with the ideal result, as Fig. 5(a) shows. However, there are a few data points that show no improvement or even an increased deviation after the ZNE. This is likely the effect of qubits with low performance, as for such a large system size, we have to use a significant fraction of the qubits on the chip (see Fig. 2). Alternatively, the deep circuit and large number of two-qubit gates for the largest two noise factors might just introduce too much noise, thus rendering the extrapolation unreliable.Fig. 5: Hardware run results for 80 qubit state preparation.Full size imagea The fermion densities from the MPS simulation (blue diamonds), hardware run result without error mitigation (gray upward triangles), and with ZNE (orange dots). Even sites are shown with filled markers, odd sites with empty markers for visual distinction. b Results after ZNE when averaging data related by CP symmetry, which improves outliers and enhances the agreement with the ideal simulation.To improve our data further, we can take advantage of CP symmetry in the model, which in spin language translates to \(\langle {\sigma }_{n}^{z}\rangle\) and \(\langle {\sigma }_{N-1-n}^{z}\rangle\) being equal. Averaging over these two results effectively allows for increasing the statistics and to mitigate outliers. Fig. 5(b) presents the ZNE results followed by the CP averaging, where the raw data has been removed to enhance visual clarity. In general, using the averaging procedure, a reduction in the deviation of most of the previous outliers and an overall better agreement with the wave packet distribution is observed. While exploiting the symmetry generally yields an improvement, certain data points still deviate from the ideal results and exhibit large error bars, indicating the need for improved hardware or more advanced error mitigation methods to achieve a similar level of precision as for N = 40 at this scale.While we do not perform further time evolution on this 80-qubit state due to the prohibitive circuit depths to simulate the full scattering dynamics on available quantum hardware, this benchmark clearly demonstrates the scalability of our tensor network optimized circuit approach for state preparation, as discussed in the Methods section, Overview of Variational Circuit Approximation. Future work will investigate real-time dynamics on large systems and improve accuracy, utilizing upgraded hardware and more advanced error mitigation techniques, such as ZNE combined with probabilistic error amplification (PEA).DiscussionIn this work, we demonstrated a hybrid quantum-classical strategy for simulating fermion scattering processes using tensor networks and quantum hardware. By using tensor network optimization to prepare low-entangling states and compress short-time evolution circuits, we significantly reduced the circuit depth compared to conventional methods. This enabled successful hardware runs for 40 qubit dynamics and 80 qubit state preparation on IBM superconducting quantum devices; after error mitigation we achieved results consistent with ideal simulations, even for circuits with two-qubit gate depths up to 96 and total two-qubit gate counts up to 1872.Specifically, we use the tensor network simulations to prepare the fermion-antifermion scattering state before the collision, where the entanglement is relatively low, and can be constructed by a shallow quantum circuit. This state serves as the initial configuration for execution on quantum hardware. To simulate the post-collision dynamics, where entanglement between subsystems increases due to interactions, we apply a sequence of time-evolution circuits starting from this state, implemented using Trotterized dynamics. Each Trotter-step time-evolution circuit can also be variationally optimized and compressed to reduce circuit depth and thereby lessen the impact of hardware noise, without significantly compromising accuracy. By executing these circuits for a sufficiently large number of steps, we are able to probe distinct dynamical behaviors following the interaction. In particular, the real-time fermion-density evolution exhibits qualitative signatures consistent with different scattering behaviors, ranging from transmission-like to reflection-like dynamics depending on the interactions encoded in the Hamiltonian.While classical simulations based on MPS remain feasible for the current system sizes and time scales, they already approach the limits of tractability as entanglement grows with time evolution. In this setting, MPS results can still serve as a valuable benchmark for validating quantum hardware runs. However, we anticipate that near-future experiments, leveraging advanced error mitigation techniques such as probabilistic error cancellation (PEC) and probabilistic error amplification (PEA)62,63,64,65, or post processing tensor network approaches (e.g, TEM)66 will enable the simulation of increasingly complex scattering phenomena on larger quantum systems (with 100 qubits and more) and over significantly longer time scales.Our method highlights the potential synergies between classical tensor network circuit compression and state propagation on a quantum computer, offering a practical path for simulating scattering in lattice field theories. While our study focuses on the Thirring model, the approach can be extended straightforwardly to other (1+1)D models, such as the Schwinger model. In addition, there are recent developments about the circuit construction for hadronic wave packets36,40,41. It would therefore be valuable to compare their resource costs with those of the tensor-network optimization strategy adopted here. In the (1+1)D setting considered in this work, we expect the TN-based optimization to offer a more resource-efficient route, whereas such an advantage is much less clear in higher dimensions. Moreover, incorporating finite lattice spacing and volume effects67 will be a necessary step toward connecting quantum simulation results to their continuum counterparts. A further important direction will be to move beyond the present density-based qualitative characterization of the dynamics and connect the real-time simulations more directly to standard scattering observables, such as phase shifts, time delays, or S-matrix elements. Taken together, these advances may eventually push quantum simulations toward regimes that are increasingly challenging for classical methods, opening the door to quantitatively controlled studies of non-perturbative phenomena in quantum field theory on near-term quantum devices.MethodsThe Thirring model and scattering setupLattice formulationIn this work, we use the Thirring model68 to study fermion scattering. While it can be solved exactly in the massless case using bosonization, and its spectrum can be determined using Bethe-ansatz in the massive case, it shares many interesting features with more complicated gauge field theories from the Standard Model. In particular, the Thirring model is renormalizable and can show scale-dependent behavior reminiscent of asymptotic freedom69. Hence, it can serve as a testbed for new lattice techniques.Adopting the Kogut-Susskind staggered formulation70,71, the lattice Hamiltonian of the model reads69$$\begin{array}{rcl}H&=&\sum _{n}\left(\frac{i}{2a}\left({\xi }_{n+1}^{\dagger }{\xi }_{n}-{\xi }_{n}^{\dagger }{\xi }_{n+1}\right)+{(-1)}^{n}m\,{\xi }_{n}^{\dagger }{\xi }_{n}\right)\\ &&+\sum _{n}\frac{g}{a}\,{\xi }_{n}^{\dagger }{\xi }_{n}{\xi }_{n+1}^{\dagger }{\xi }_{n+1}\,,\end{array}$$ (2) where \({\xi }_{n}^{\dagger }\) and ξn are fermion creation and annihilation operators; a is the lattice spacing, m is the fermion mass, and g is coupling strength of the four-fermion interaction term. Without loss of generality, we set a = 1 for the rest of this work.As in ref. 23, in the non-interacting case (g = 0) and with periodic boundary conditions, the Hamiltonian in Eq. (2) can be diagonalized in momentum space, \(H={\sum }_{k}{w}_{k}\left({c}_{k}^{\dagger }{c}_{k}-{d}_{k}{d}_{k}^{\dagger }\right)\), where the momentum dependent operators are given by$$\begin{array}{rcl}{c}_{k}^{\dagger }&=&\frac{1}{\sqrt{N}}\sqrt{\frac{m+{w}_{k}}{{w}_{k}}}\sum _{n}{e}^{ikn}\left({\Pi }_{n0}+{v}_{k}{\Pi }_{n1}\right){\xi }_{n}^{\dagger },\\ {d}_{k}^{\dagger }&=&\frac{1}{\sqrt{N}}\sqrt{\frac{m+{w}_{k}}{{w}_{k}}}\sum _{n}{e}^{ikn}\left({\Pi }_{n1}+{v}_{k}{\Pi }_{n0}\right){\xi }_{n}\,,\end{array}$$ (3) with k ∈ 2π/N × { − ⌊N/4⌋, ⋯ ⌈N/4⌉ − 1} and N the number of sites. The constants vk and wk correspond to$${v}_{k}=\frac{\sin (k)}{m+{w}_{k}},\,\,{w}_{k}=\sqrt{{m}^{2}+{\sin }^{2}(k)},$$ (4) and Πn0(Πn1) are projection operators defined as$$\begin{array}{r}{\Pi }_{nl}=\left\{\begin{array}{ll}1,\quad &n\equiv l\,(\mathrm{mod}\,\,2),\\ 0,\quad &n\not\equiv l\,(\mathrm{mod}\,\,2),\end{array}\right.\qquad l\in \{0,1\}.\end{array}$$ (5) The operators \({c}_{k}^{\dagger }\) (\({d}_{k}^{\dagger }\)) create (annihilate) fermions in the position space, and are therefore referred to as fermion (antifermion) operators.Fermion scattering setupBuilding on the work of refs. 23,35, we investigate the scattering between a fermion and an antifermion wave packet for the interacting case. As shown in the reference, such wave packets can be created on top of the ground state \(| \Omega \left.\right\rangle\), by acting on it with a set of creation operators$$| \psi (t=0)\left.\right\rangle ={D}^{\dagger }{C}^{\dagger }| \Omega \left.\right\rangle .$$ (6) Here C† and D† create, respectively, a Gaussian fermion and antifermion wave packet. These operators can be expressed as the linear combinations$$\begin{array}{rcl}{C}^{\dagger }({\phi }^{c})&=&\sum _{k}{\phi }_{k}^{c}{c}_{k}^{\dagger }=\sum _{n}{\tilde{\phi }}_{n}^{c}{\xi }_{n}^{\dagger },\\ {D}^{\dagger }({\phi }^{d})&=&\sum _{k}{\phi }_{k}^{d}{d}_{k}^{\dagger }=\sum _{n}{\tilde{\phi }}_{n}^{d}{\xi }_{n},\end{array}$$ (7) where the Gaussian coefficients \({\phi }_{k}^{c(d)}\) in momentum space are given by$${\phi }_{k}^{c(d)}=\frac{1}{\sqrt{{{\mathcal{N}}}_{k}^{c(d)}}}{e}^{-ik{\mu }_{n}^{c(d)}}{e}^{-{(k-{\mu }_{k}^{c(d)})}^{2}/4{\sigma }_{k}^{2}}\,.$$ (8) In the expression above \({\mu }_{n}^{c(d)}\) corresponds to the position around which the wave packet is centered, \({\mu }_{k}^{c(d)}\) to the mean momentum, σk represents the width in momentum space, and \(\sqrt{{{\mathcal{N}}}_{k}^{c(d)}}\) is a normalization factor.Throughout our work, we set \(\{{\mu }_{k}^{c},{\mu }_{k}^{d}\}=\{4\times 2\pi /N,-4\times 2\pi /N\},\,{\sigma }_{k}=2\pi /N\), and \(\{{\mu }_{n}^{c},{\mu }_{n}^{d}\}=\{N/4,3N/4-1\}\). The factors \({\tilde{\phi }}_{n}^{c(d)}\) in Eq. (7) represent the coefficients in position space, which can be obtained by Fourier transformation from the ones in momentum space as detailed in Eq. (13) of ref. 35. Although the operators defined above yield exact fermion (antifermion) wave packets only in the noninteracting limit under periodic boundary conditions, they nevertheless provide a good approximation in the interacting regime, provided that the coupling constant g is sufficiently small35.To simulate the fermionic degrees of freedom on a quantum computer, we map them to Pauli matrices using the Jordan-Wigner transformation$${\xi }_{n}^{\dagger }=\prod _{l < n}\,{\sigma }_{l}^{z}{\sigma }_{n}^{-},\quad {\xi }_{n}=\prod _{l < n}\,{\sigma }_{l}^{z}{\sigma }_{n}^{+},$$ (9) where \({\sigma }_{l}^{\pm }=\left({\sigma }_{l}^{x}\pm i{\sigma }_{l}^{y}\right)/2\) and \({\sigma }_{l}^{j},\,j\in \{x,y,z\}\) are the usual Pauli matrices. The resulting qubit Hamiltonian in terms of Pauli operators then reads$$\begin{array}{rcl}H&=&\frac{i}{2}\mathop{\sum }\limits_{n=0}^{N-2}\left({\sigma }_{n+1}^{-}{\sigma }_{n}^{+}-{\sigma }_{n}^{-}{\sigma }_{n+1}^{+}\right)\\ &&+\frac{m}{2}\mathop{\sum }\limits_{n=0}^{N-1}{(-1)}^{n}\left({\mathbb{1}}-{\sigma }_{n}^{z}\right)\\ &&+\frac{g}{4}\mathop{\sum }\limits_{n=0}^{N-2}\left({\mathbb{1}}-{\sigma }_{n}^{z}\right)\left({\mathbb{1}}-{\sigma }_{n+1}^{z}\right),\end{array}$$ (10) where we have chosen open boundary conditions for simplicity. For the spatially localized wave packets we use in our simulations, boundary effects remain negligible, provided the system size is sufficiently large. In both our classical simulations and quantum hardware runs with N = 40, no significant boundary effects were observed.As pointed out in ref. 35, simulating the scattering process generally involves the following steps. First, the ground state of the model has to be prepared. It can be obtained, for example, by a variationally optimized circuit \(| \Omega \left.\right\rangle ={U}_{{\rm{GS}}}| 0\left.\right\rangle\). Second, particle wave packets can be created on top of the vacuum by applying suitable operators as in Eq. (6). A corresponding quantum circuit for these operators can be constructed based upon Givens rotations (see the Supplementary Information (Supplementary Note 3), Figure S3) to prepare this state \(| \psi (0)\left.\right\rangle ={U}_{{\rm{WP}}}| \Omega \left.\right\rangle\), serving as the initial state of the scattering process. Finally, the dynamics can then be realized by performing real-time evolution of the initial state up to a time t. This can be done using standard Trotterization, e.g. approximating the time evolution operator by a series of unitaries implementing a small timestep Δt which can be efficiently realized on quantum hardware, \(| \psi (t)\left.\right\rangle =U{(\Delta t)}^{t/\Delta t}\left\vert \psi (0)\right\rangle\).However, the quantum resource requirements for simulating a complete scattering process using the above protocol remain quite substantial. For a 40-qubit system, the total circuit depth exceeds 300, with over 5000 two-qubit gates (see Table 1 in the “Methods” section, Overall Circuit Depth Reduction). While a total number of around 5000 CNOTs may be feasible in practice, the primary limitation arises from the circuit depth and the associated two-qubit gate error and qubit decoherence. The details of this resource estimation can be found at the Supplementary Information (Supplementary Note 3).Tensor networks can efficiently simulate low-entanglement dynamics. In the scattering process, entanglement typically increases after the collision, making the tensor network simulation more costly. However, before the collision, especially at the early stage, entanglement remains relatively low, allowing for efficient and accurate tensor network simulation. To illustrate this, we take the fermion-antifermion scattering as an example and calculate the entanglement entropy of bipartition. Specifically, we calculate the von Neumann entropy \({S}_{n}(t)=-{\rm{tr}}[{\rho }_{n}(t){\log }_{2}{\rho }_{n}(t)]\), where ρn(t) is the reduced density matrix of the first n qubits at time t. To quantify the entropy generated throughout the scattering process, we subtract the vacuum contribution and obtain the excess entropy ΔSn(t).As shown in Fig. 6(a), before the collision, entanglement is mainly localized within each wave packet. At the collision point, where the wave packets overlap, entanglement peaks at the central sites. After the collision, significant entanglement is generated between the outgoing particles, and long-range entanglement persists as they propagate, leading to increasing simulation costs for tensor networks.Fig. 6: Entanglement growth during scattering dynamics.Full size imageEntanglement entropy and half-chain bond dimension over time in the fermion-antifermion scattering process. a The bipartite entanglement entropy for the case (m, g) = (0.4, 0.5) for H in Eq. (10). The x-axis is the site index, and the y-axis represents time. b The half-chain bond dimension χ versus normalized time t/T. Empty markers represent the time slices simulated by tensor networks and compiled into the state preparation circuit, and the filled markers represent time slices that will be simulated directly on quantum hardware. The time t is normalized with the total time T to visualize different processes in the same plot, T ∈ {22, 29, 27} for (m, g) ∈ {(0.2, 0.4), (0.4, 0.5), (0.4, 0.7)} respectively. All simulations are performed using MPS with an SVD truncation threshold of 10−8.To quantify the associated tensor-network cost during the scattering process, Fig. 6(b) shows the bond dimension χ at the middle bond as a function of time. We observe that χ increases during the time evolution, especially after the collision (marked by filled symbols), indicating the growing tensor-network cost of simulating the post-collision dynamics. This behavior motivates a hybrid strategy: simulate the low-entangled states before collision by tensor networks, and continue the subsequent dynamics on quantum computers. To realize this hybrid approach, we employ a recently developed circuit optimization algorithm based on tensor networks55, which we apply to both state preparation and time-evolution circuits, as demonstrated in Fig. 1. This approach significantly reduces circuit complexity, enabling high-accuracy hardware execution for 40-qubit systems.Tensor network circuit approximationTo implement the hybrid strategy introduced above, we use tensor-network-based circuit optimization to construct compact-depth circuits for the quantum hardware simulation. In this direction, there are two natural circuits to produce using this optimization. The first circuit prepares the initial state followed by a period of time evolution, \({e}^{-iHt}\left\vert \psi (0)\right\rangle\). The initial state, \(\left\vert \psi (0)\right\rangle\), consisting of the two spatially separated wave packets on top of the vacuum, lacks long-range entanglement and therefore can be efficiently represented by an MPS. As the time evolution proceeds and the wave packets move together and interact, entanglement increases and a correspondingly higher bond dimension is required to accurately describe this state. Therefore, we can use a tensor network time evolution algorithm to evolve the state \({e}^{-iHt}\left\vert \psi (0)\right\rangle\) up to the largest time t0 that still permits an accurate representation with a computationally tractable bond dimension. The second circuit exploits the low-entangling nature of e−iHt for small values of t. When the Hamiltonian, H, has a quasi-1D topology with local interactions, then e−iHt can be efficiently represented by a low bond dimension MPO72. Again, this permits a variational optimization to search for circuits approximating this unitary, with shallower circuits compared to those derived from Trotterization.Overview of variational circuit approximationHere we briefly describe the approach taken for the variational circuit approximation, with more details found in the Supplementary Information (Supplementary Note 1) and ref. 55.We focus on minimizing 2-qubit gate count and depth, as these dominate errors on near-term hardware. During the optimization, all variational circuits are constrained to have a 1D brickwork structure of layers of SU(4) gates (the most general expression of a 2-qubit gate as a dimension-4 unitary matrix) with nearest-neighbor connectivity, resulting in layers of gates acting on ‘odd’ and ‘even’ pairs of qubits (as sketched in Fig. S1 in the Supplementary Information (Supplementary Note 1)). All SU(4) gates can be decomposed into at most three 2-qubit gates (c.f. the Supplementary Information (Supplementary Note 3), Eq. (S1)), therefore a circuit composed of l SU(4) layers can be described by at most 3l layers of primitive 2-qubit gates (e.g. CNOT, CZ) when mapped onto the device native gate set. These layers of 2-qubit gates all act on distinct qubits, and therefore can be applied in parallel.All our circuit optimizations can be separated into two categories. First, the generation of a short-depth circuit V1(θ) to approximately prepare a target state$${V}_{1}({\boldsymbol{\theta }})\left\vert 0\right\rangle \approx \left\vert {\psi }_{{\rm{Targ}}}\right\rangle ,$$where \(\left\vert {\psi }_{{\rm{Targ}}}\right\rangle\) is represented by a low-bond dimension MPS. A successful optimization will output a circuit to a target fidelity with a lower gate count/depth compared to the alternate ‘conventional’ approach for wave packet preparation described in ref. 35. This motivates the following cost function for a numerical optimizer to minimize, measuring the infidelity between the target and variational state$${C}_{{\rm{State}}}({\boldsymbol{\theta }})=1-{\left\vert \left\langle {\psi }_{{\rm{Targ}}}\right\vert {V}_{1}({\boldsymbol{\theta }})\left\vert 0\right\rangle \right\vert }^{2}.$$ (11) Our second category of circuits involves simulating real-time evolution under the Hamiltonian. e−iHt, for small t, can be represented by a low-bond dimension MPO and is computed using high accuracy MPO Trotterization. We aim to use the variational optimization to produce a short depth circuit (shorter than a standard Trotter depth circuit) to simulate the time evolution, such that$${V}_{2}({\boldsymbol{\theta }})\approx {e}^{-iHt}.$$ (12) This could be done by variationally searching for a shallower circuit to simulate a given Trotter step with a fixed error. Alternatively, given a fixed depth ansatz for the circuit, we can optimize the ansatz to minimize the error. We here take the latter approach.This motivates the following cost function, measuring a Hilbert-Schmidt inner product between the unitaries$${C}_{{\rm{Uni}}}({\boldsymbol{\theta }})=1-\frac{1}{{2}^{2N}}{\left\vert \text{Tr}({V}_{2}{({\boldsymbol{\theta }})}^{\dagger }{e}^{-iHt})\right\vert }^{2}.$$ (13) Both circuits use the same optimization procedure, where the SU(4) gates defining the circuits are iteratively updated by computing the gate’s environment tensor, after which a Polar Decomposition is applied to find the optimal gate update. This is a standard practice in tensor network optimization73, which has been similarly motivated for circuit optimization49,54,55,74. Further details on the tensor network based optimization can be found in the Supplementary Information (Supplementary Note 1).Initial stateIn this section, we detail how the initial state wave-packet configuration at t = 0 is generated as an MPS. First, the ground state of the Hamiltonian (see Eq. (10)) in the half-filling U(1) symmetry sector is computed as an MPS using the 2-site density matrix renormalization group (DMRG) algorithm75; we use the quantum number preserving MPS functionality provided by the ITensorMPS package76 to perform this calculation. For all Hamiltonian settings considered, the ground state can be represented by an MPS with a maximum bond dimension of 32, after retaining a minimum singular value of 10−12.The ground state MPS is then used to create the state \(\left\vert \psi (0)\right\rangle\) describing the initial configuration of the separated fermion and antifermion wave packets. This state is exactly constructed using (Eqs. (6)–(9)) with MPS arithmetic. The addition of MPS in this fashion substantially increases the bond dimension compared to the ground state. To find a more compact MPS representation, we perform a variational compression75. Namely, we sweep through the site tensors of a lower bond dimension variational MPS, maximizing the fidelity with the uncompressed MPS (Note that here we apply the individual terms of Eq. (7) to the ground state and sum the resulting states. This leads temporarily to a higher bond dimension. However, as shown in ref. 23, Eq. (7) can also be written as an MPO of bond dimension 2, thus only doubling the initial bond dimension. In both cases variational compression can be used to decrease the bond dimension. We perform this compression for increasing bond dimension MPS, finding the smallest with infidelity less than 10−6 with respect to the uncompressed state, which resulted in an MPS with a maximum bond dimension χ ≈ 20.To probe the scalability of MPS-based circuit optimization we first compile circuits for preparing the initial wave packet configurations \(\left\vert \psi (0)\right\rangle\) with the Hamiltonian coefficients (m, g) = (0.2, 0.4) on increasing system sizes N ∈ {40, 80, 120, 160}. We perform the variational optimization for increasing depth circuits, with the value of the infidelity cost function (CState) upon convergence plotted as a function of circuit depth, shown in Fig. 7. When the depth is scaled with respect to the number of qubits in the system, a similar scaling is observed in the achievable infidelity to this target state versus the CNOT layers per qubit in the circuit. The quality of these approximate initial state preparations is further visualized in Figure S2 in the Supplementary Information (Supplementary Note 2), by computing the fermion densities of the initial state. These state preparation results on up to 160 qubits provide strong evidence that MPS circuit optimization can aid quantum simulation on system sizes far beyond the limits of exact state vector simulation.Fig. 7: Scaling of Optimization–Initial State.Full size imageThe variational optimization is performed to generate circuits to prepare the state \(\left\vert {\psi }_{{\rm{Targ}}}\right\rangle =\left\vert \psi (0)\right\rangle\), representing the initial separated wave packets, on increasing size systems, with N ∈ {40, 80, 120, 160} qubits shown by the colored data. Each data point corresponds to the outputted infidelity cost function value after convergence of the optimization. When the depths of the circuits are scaled by the system size, as shown in the main plot where the x-axis plots the number of CNOT layers per qubit, a similar scaling across all qubit numbers is observed in achievable infidelity versus CNOT layers per qubit. The unscaled data (with the x-axis simply the number of CNOT layers) is shown in the inset.For our hardware demonstration, we further push the use of MPS pre-computation by incorporating some early time evolution steps into the target state. Concretely, we simulate the evolution up to where the wave packets begin to interact and residual entanglement entropy begins to be generated (as sketched in Fig. 1a)). We perform this time evolution using the time evolution block decimation (TEBD) method77 with a maximum MPS bond dimension of χ = 150, using a 2nd-order Trotterization of the propagator with a timestep of Δt = 0.25. This computation results in an MPS representing the target state \(\left\vert {\psi }_{{\rm{Targ}}}\right\rangle ={e}^{-iH{t}_{0}}\left\vert \psi (0)\right\rangle\). The variational optimization is performed three times to find circuits to approximately prepare the target state, for Hamiltonians with coefficients (m, g) ∈ {(0.2, 0.4), (0.4, 0.5), (0.4, 0.7)} and corresponding initial evolution times t0 ∈ {11, 18, 16}. The result of these optimizations is presented in Fig. 8. For (m, g) = (0.2, 0.4), a circuit depth of 30 CNOT layers is used for V1(θ), and for (m, g) ∈ {(0.4, 0.5), (0.4, 0.7)} a circuit depth of 36 CNOT layers is used.Fig. 8: Optimization–Time Evolved Initial State.Full size imageData showing final cost function value upon convergence of variational MPS optimization for generating the state \(\left\vert {\psi }_{{\rm{Targ}}}\right\rangle ={e}^{-iH{t}_{0}}\left\vert \psi (0)\right\rangle\). The depth of the variational circuit measured in layers of commuting CNOT gates, the y-axis measures the corresponding value of infidelity (CState) upon convergence of the optimization. The target state is computed using MPS arithmetic, followed by TEBD.The 3 data curves correspond to the Hamiltonian settings considered, (m, g) ∈ {(0.2, 0.4), (0.4, 0.5), (0.4, 0.7)}, with corresponding initial evolution times of t0 ∈ {11, 18, 16} respectively.Time evolution unitaryTo continue the dynamics, we further perform real-time evolution under the Hamiltonian Eq. (10) on the quantum computer. This requires the implementation of the unitary e−iHt. To derive a circuit to approximate this unitary, the standard approach would be to employ a Trotterization. Under this approach, there is always a trade-off between the Trotter error versus the circuit depth. We further take advantage of tensor network circuit optimization to search for both short-depth and low error circuits to approximate this unitary.Given the Trotter errors we can tolerate, and the low depth overhead compared to 1st-order Trotter for 1D nearest neighbor Hamiltonians, the 2nd-order Trotterization gives the shortest depth circuits out of the available Trotter orders. Therefore, the circuit V2(θ) is initialized with the equivalent depth 2nd-order Trotterization, and any reduction in error versus e−iHt is an advantage (i.e., has a lower error than the equivalent-depth Trotter step circuit).The MPO-based optimization updates the SU(4) gates via Polar Decomposition of gate environments to minimize the cost function Eq. (13), as described in ref. 55. Figure 9 shows the results of this optimization.Fig. 9: Optimization–Time Evolution Unitary.Full size imageData comparing the approximation error to the target unitary e−iHt, with timestep t = 2.0, for increasing depth circuits. The data connected by the dotted line labeled `Trot-II' corresponds to increasing depth 2nd-order Trotterizations of the target unitary. The full line labeled `Opt' displays the approximation error, CUni, of increasing depth circuits derived from the MPO-based optimization.The time evolution operator e−iHt with t = 2.0 is constructed using a high-accuracy Trotterization, with the Trotter time step chosen such that the result is converged, and is then represented as an MPO with bond dimension χ = 128. For each of the three Hamiltonian parameter pairs, (m, g) ∈ {(0.2, 0.4), (0.4, 0.5), (0.4, 0.7)}, we find that (excluding the shortest depth circuit of 9 CNOT layers) the optimization is able to reduce the error CUni compared to the equivalent depth 2nd-order Trotterization by 6.60, 8.00 and 5.52 respectively on average across the range of circuit depths tested. For all hardware demonstrations we perform the time evolution using the optimized circuits with a circuit depth of 15 CNOT layers. We highlight this has a lower error than the second-order Trotterization circuits with a depth of 21 CNOT layers, an over 25% reduction in circuit depth.Finally, we note that, unlike the initial state, the short-time evolution operator considered here is generated by a one-dimensional local Hamiltonian with a translation invariant bulk structure. This is expected to make the corresponding compression less sensitive to system size in the regime considered here. This expectation is consistent with recent work using infinite tensor networks, where translation invariance is exploited to avoid an explicit scaling of the optimization complexity with system size and to learn compressed unit-cell circuits for the same Thirring Hamiltonian directly in the infinite-system setting59.Overall circuit depth reductionWe highlight the overall reduction in circuit depth for full simulation of scattering dynamics gained through the use of the tensor network optimization, in comparison to the alternate construction previously described in ref. 35. Here, the circuit suggested can be summarized by \(U{(\Delta t)}^{t/\Delta t}{U}_{{\rm{WP}}}{U}_{{\rm{GS}}}\left\vert 0\right\rangle\), where UGS prepares the vacuum state \(\left\vert \Omega \right\rangle ,\,{U}_{{\rm{WP}}}\) is an operator preparing the initial wave packet configuration on top of the vacuum, and U(Δt) is the Trotter unitary used to simulate real-time evolution. We give circuit depth estimates for these three circuits separately.The depth of UGS to prepare \(\left\vert \Omega \right\rangle\) is negligible compared to the total circuit depth, so its contribution is omitted. For an N-qubit system, the CNOT depth of UWP decomposed by Givens rotations is around 2N − 4. Similarly, the CNOT depth due to time evolution by second-order Trotterization over total time T is given by 6T/Δt + 3, see the Supplementary Information (Supplementary Note 3) for further details. The sum of these two formulas give a lower-bound on the total CNOT depth$${D}_{{\rm{Conv}}}=2N+6T/\Delta t-1.$$ (14) In our hardware experiments of the full scattering process we implement the simulation on N = 40 qubits. We observed in Fig. 9 that an error of CUni ≈ 0.01 for approximating the time-evolution circuits provided results that were sufficiently consistent with high-accuracy MPS simulation. A second-order Trotterization depth with the same (Trotter) error would use Δt = 2/3. The total evolution time is chosen such that at the end of the simulation the scattered wave packets have fully spread to the edge of the boundaries. For the three Hamiltonian settings (m, g) ∈ {(0.2, 0.4), (0.4, 0.5), (0.4, 0.7)}, the corresponding total evolution times are found to be T = 21, 28, and 26, leading to total CNOT depths for the conventional approach of DConv = 268, 331, and 313, respectively.In comparison, the simulation using our MPS-optimized circuits for Hamiltonian settings (m, g) ∈ {(0.2, 0.4), (0.4, 0.5), (0.4, 0.7)} required total circuit depths of DOpt = 90, 96, 96, respectively. By computing the relative reduction in circuit depth, given by DConv/DOpt = 2.98, 3.44, 3.26, we find that on average the total circuit depth of the simulation has been reduced by a factor of 3.23 compared to the conventional approach, a significant reduction which greatly improves the ability to perform these simulations on state-of-the-art noisy processors. For our 80 qubit hardware demonstration of state preparation, the tensor network optimized circuits can realize the target state \({e}^{-iH{t}_{0}}\left\vert \psi (0)\right\rangle\) with a significantly reduced circuit depth of 24 CNOT layers, down from 249 if using the conventional approach, a depth reduction by a factor of 10.0. Table 1 summarizes these improvements.