State preparation and detection for quantum simulation of particle collisions

1 Introduction

Quantum simulators have emerged as a powerful approach for simulating the dynamics of quantum many-body systems, competing directly with state-of-the-art classical computational techniques [1, 2]. A promising application is the simulation of scattering processes, akin to those occurring in particle colliders [3–12]. In this setting, the theory of interest is modeled as a lattice Hamiltonian, and the particles are the low-energy excitations (quasiparticles) carrying well-defined momentum. Quantum simulators can reproduce the collision of two such particles in real time, probing not only the scattering matrix (which encodes the probabilities of the final asymptotic states), but also the state of the system during the collision or shortly after. This capability opens the door to studying phenomena such as the emergence of hydrodynamics, thermalization, and aspects of heavy-ion collisions in controlled quantum many-body systems [13, 14].

Quantum simulators may become particularly advantageous in the high-energy regime, where the energy is much larger than the mass gap and inelastic channels with many particles become accessible. Since the entanglement entropy generated in a scattering event is expected to increase with the number of particles produced, such processes are challenging to simulate using tensor-network methods [15]. While tensor networks have proven highly successful at low energies [15–23], their applicability is expected to deteriorate as the collision energy and the number of particles increase. Furthermore, although tensor-network methods have achieved remarkable success in one dimension, extending them to two dimensions is substantially more demanding [23], and they have not yet been demonstrated for scattering problems in three dimensions. Quantum simulators therefore offer a potential route toward exploring scattering regimes that are difficult to access with existing classical methods.

Simulating scattering on a quantum simulator requires the following steps [3–5]: (i) preparing an initial state consisting of two particle wave packets traveling toward one another, (ii) evolving this state with the model Hamiltonian until after the particles have collided, (iii) measuring observables to extract relevant properties of the final state. Assuming that the desired Hamiltonian can be implemented, either with analog approaches or with Trotterization in a digital quantum computer, the preparation of the initial state and characterization of the final state remain key challenges. Several methods have been proposed to address these tasks [3, 20, 24–36], particularly for digital quantum simulators. Some of these approaches have been implemented in the first digital quantum simulations of particle collisions [37–41]. These works demonstrated an impressive degree of control over wave packet preparation. However, they also highlighted the difficulty of maintaining high fidelity over the long evolution times required to observe inelastic scattering, due to multiple sources of error. This motivates the development of robust scattering protocols for analog quantum simulators, where the time evolution can be implemented directly, without Trotter errors.

In this work, we propose a method for state preparation and detection that can be used to simulate collisions on several different platforms and is particularly suited for analog simulators. The wave packet is prepared by coupling one boundary site (in 1D) or a corner site (in 2D) of the system to an auxiliary qubit, initialized in the excited state, which then decays by emitting an excitation into the system (figure 1). More generally, the same effect can be achieved by designing a local quench which injects the energy needed for the target excitation, while maintaining a small energy variance. The method we propose has the significant advantage of being robust to imperfections or errors in the calibration of the model’s parameters, since it relies exclusively on the conservation of energy and conservation of energy variance, and does not require fine-tuning or prior knowledge of the energy eigenstate representing the particle. The momenta of the final states are then detected, using only local measurements, by observing the interference pattern when the particles are reflected at the system’s boundary.

Our state preparation protocol shares similarities with the well-studied model of a superconducting qubit or quantum dot (here, the auxiliary qubit) coupled to a waveguide (the system) [42]. In the context of wave packet preparation for particle collisions, related protocols have been suggested in refs. [20, 24], but have not been verified in a quantum many-body system. Here, we demonstrate the preparation of particle wave packets first in a simple single-particle model, using exact time evolution, and then in two interacting many-body systems (a Rydberg atom chain and a mixed-field Ising chain) using tensor-network simulations. By combining this protocol with an adiabatic ramp and a “quantum slide” (in which the Hamiltonian is spatially modulated), we show how to prepare approximately Gaussian wave packets in a regime where inelastic scattering is kinematically allowed.

The paper is structured as follows. In section 2 we focus on the analytically tractable single-particle model for a pedagogical illustration of the state preparation and detection protocols. In section 3 we move on to the case of a generic many-body system, here exemplified by a Rydberg atom chain, demonstrating all the steps needed to simulate elastic scattering. In section 4 we show how to use adiabatic ramps to probe inelastic scattering, focusing on a model with multiple distinct single-particle excitations (a mixed-field quantum Ising chain). In section 5 we demonstrate a method to control the shape of the wave packet using a quantum slide, and discuss its underlying theory. In section 6 we demonstrate the wave packet preparation in a two-dimensional single-particle model and discuss prospects for simulating scattering in more than one spatial dimension.

2 Single-particle model

To illustrate the protocols for state preparation and detection, we start with the textbook model of a single particle hopping in a one-dimensional lattice (tight-binding model). This model captures the general idea of the protocol and is simple enough that we can have analytic understanding, giving us the tools to later approach the case of a full many-body Hamiltonian.

The Hilbert space of a single particle in a one-dimensional lattice is spanned by the states $|j⟩$ with $j\; \in \; \mathbb{Z}$ labeling the position of the particle. We will consider a model where the particle can hop between neighboring sites, with a hopping amplitude $w$:

The Hamiltonian eq. (2.1) is easily diagonalized in momentum space. The eigenvectors $|{\varphi }_k⟩$ are plane waves $|{\varphi }_k⟩={\sum }_j{e}^{ikj}|j⟩$ labeled by the lattice momentum $k\; \in \; [0,2\pi )$ with energy $E(k)\; =\; -2w\cos <!--\; FUNCTION\; APPLICATION\; -->k$.

Our first goal is to prepare a wave packet, i.e., a state with a narrow momentum distribution around a value $k_{\ast }$ of choice, and having a finite width in real space. An example is a Gaussian wave packet of the form $|{\psi }_G⟩={\mathcal{\text{𝒩}}}^{-1}{\sum }_j{e}^{-\alpha {(j-j_0)}^2+ikj}|j⟩$, where ${\alpha }^{-1∕2}$ represents the spatial width of the wave packet and $\mathcal{\text{𝒩}}$ is a normalization constant. Gaussian wave packets are particularly useful because they minimize the uncertainty product between position and momentum, yielding the most localized states for a given momentum uncertainty. We will more generally consider wave packets with non-Gaussian shapes, provided that they are narrow in momentum space and sufficiently localized in real space. We want to prepare such states starting from a simple initial state — a particle localized on the first site of a semi-infinite one-dimensional lattice — by time evolving using a Hamiltonian with locally tuned parameters. As will be discussed in section 3, this protocol will mimic the evolution after a local quench in a quantum simulator. To study a collision in a one-dimensional system, one can use this protocol to “inject” a particle from each edge of a finite chain.

Our second goal is to find a protocol for detecting the momentum of a wave packet. To mimic the measurements that are typically done in an experimental setup such as an array of neutral atoms, of superconducting qubits, or a trapped ion system, we will assume that the observer can measure in a local basis. In our single-particle model, this corresponds to measuring the position of the particle, i.e., having access to the probability distribution $P(j)\; =\; |⟨j|\mathrm{\Psi }⟩{|}^2$ for a state $|\mathrm{\Psi }⟩$, while direct measurements in the momentum basis are not allowed.

2.1 State preparation

Our state preparation scheme for a wave packet moving from the left edge of a chain uses site-dependent hopping amplitudes and on-site potentials in a small region near the edge. As the simplest setup, we now consider a scheme where we tune the Hamiltonian parameters on the first site of a semi-infinite chain. We will later consider a more refined scheme, where the local parameters are modulated in a larger region (section 5). The Hamiltonian for the simple case reads

where we add a potential $V_0$ at the first site $j\; =\; 0$, and we change the hopping parameter between the first two sites to $w^{\prime }$. The system is initialized in the state $|\psi (t\; =\; 0)⟩=|0⟩$, i.e., with the particle localized at site $j\; =\; 0$. In the evolution under $H_{\mathrm{\text{prep}}}$, the particle will hop and eventually move away from the edge. The crucial ingredient for our protocol is that the energy and energy variance are conserved in the evolution, i.e.,

Conservation of energy and energy variance strongly constrains the modes that can be excited in the bulk of the system, once the particle has moved away from the edge: we expect that only modes in an energy window $\approx V_0\pm w^{\prime }$ will be excited. As shown in figure 2(a), by choosing $V_0=\; E(k_{\ast })\; =\; -2w\cos <!--\; FUNCTION\; APPLICATION\; -->k_{\ast }$ and $w^{\prime }\ll \; w$, this small energy window corresponds to a small range of momenta around $k_{\ast }$.

2.1.1 Spectrum

We now make the above statement more precise by solving for the spectrum of $H_{\mathrm{\text{prep}}}$, and showing that only eigenstates with $k\; \approx k_{\ast }$ have large overlap with the initial state $|\psi (t\; =\; 0)⟩=|0⟩$.

The eigenstates of $H_{\mathrm{\text{prep}}}$ in a semi-infinite chain are superpositions of an incoming plane wave with momentum $-\; k$ and an outgoing plane wave with momentum $+\; k$ and energy $E(k)\; =\; -2w\cos <!--\; FUNCTION\; APPLICATION\; -->k$. Their exact form (including the phase shift) is found by imposing the boundary condition at the edge $j\; =\; 0$ (see appendix A)

where $k\; \in \; [0,\pi ]$ and

The eigenstates for some values of $k$ are plotted in figure 2(b). From eq. (2.5) we immediately find that the eigenstate’s overlap $|{\psi }_k(j\; =\; 0)|$ with the initial state is largest when $z_k=\; \pi ∕2$, i.e., so that $\cos <!--\; FUNCTION\; APPLICATION\; -->k\; =\; w{V}_0∕[{(w^{\prime })}^2-\; 2w^2]$. In the limit $w^{\prime }\ll \; w$, this condition corresponds to the one predicted by energy conservation $E(k_{\ast })\; =V_0$ (figure 2(c)).

2.1.2 Time evolution

We now study the time-evolved state $|\psi (t)⟩=e^{-i{H}_{\mathrm{\text{prep}}}t}|0⟩$, demonstrating the formation of a wave packet with momentum $k_{\ast }$. The real-space probability distribution $P(j,t)\; =\; |⟨j|\psi (t)⟩{|}^2$ is plotted in figure 3.

A localized pulse is emitted from the boundary site $j\; =\; 0$ at time $t\; =\; 0$ and propagates into the bulk with the characteristic group velocity $v(k_{\ast })\; =\; (dE∕dk){|}_{k_{\ast }}=\; 2w\sin <!--\; FUNCTION\; APPLICATION\; -->(k_{\ast })$, corresponding to the target momentum $k_{\ast }$, indicated by the white dashed line. The resulting wave packet is clearly non-Gaussian, exhibiting pronounced oscillations as well as a long asymmetric tail. Nevertheless, it is sharply peaked in momentum space, as anticipated from our spectral analysis and as we will further illustrate in the next subsection. In the present simulation we consider a finite system of fixed size $L\; =\; 400$; consequently, once the wave packet reaches the opposite boundary it undergoes reflection.

2.2 State detection

As shown in the previous subsection, the group velocity of the wave packet can be inferred from the local probability distribution, which in turn allows one to determine the momentum, assuming the dispersion relation is known. While this approach is useful in many cases, it becomes less practical when the velocity varies little with $k$, as, for example, in the relativistic regime. Here, we illustrate an alternative method that, like the previous one, relies solely on measurements of the local probability distribution, and does not require previous knowledge about the dispersion relation.

As shown in figure 3, when the wave packet is reflected at the boundary, the real-space probability distribution exhibits pronounced oscillations. These oscillations arise from the interference between the incoming wave packet, with momentum $k_{\ast }$, and the reflected wave packet, with momentum $-k_{\ast }$. The resulting interference pattern is $|e^{ik\ast j}+e^{-ik\ast j}{e}^{i\varphi }{|}^2=\; 2[1\; +\; cos(2k\; \ast j\; -\varphi )]$, where $\varphi$ is the phase shift acquired upon reflection. The momentum $k_{\ast }$ can then be extracted from the Fourier transform of the probability distribution $P(j,t)$ with respect to the space coordinate $j$:

Since we are interested in the interference fringes on the right end of the chain, in the definition of $\mathcal{\text{𝒫}}(q,t)$ we omit the first site ($j\; =\; 0$), thereby obtaining a cleaner signal. Note that $\mathcal{\text{𝒫}}(q,t)$ differs from the actual momentum-space probability distribution, which requires first Fourier transforming $\psi (j,t)$ and then taking the modulus squared. This distinction is important: the true momentum-space distribution cannot be accessed through single-site measurements alone, whereas $\mathcal{\text{𝒫}}(q,t)$ can be reconstructed directly from measurements of the local probability $P(j,t)$.

The function $|\mathcal{\text{𝒫}}(q,t){|}^2$, obtained from the Fourier transform of the probability distribution in figure 3, is plotted in figure 4. The plot shows a clear peak at $q\; \approx \; 2k_{\ast }$, consistent with the interference argument. When the wave packet approaches the boundary, the peak initially appears at slightly larger momenta; as the reflection develops, it shifts toward $2k_{\ast }$ and increases in magnitude. The early appearance of the peak at larger $q$ can be attributed to interference involving the high-momentum components of the wave packet, which propagate with a larger group velocity. Their contribution, however, is weaker, since they originate from the small high-$k$ tail of the wave packet.

3 Many-body Hamiltonian

In this section, we show how the protocols for state preparation and detection that we introduced in the single-particle model can be generalized to a many-body Hamiltonian. We will consider the model of a one-dimensional array of Rydberg atoms, where each atom of the chain can be in its internal ground state $|g⟩$ or in a Rydberg state $|r⟩$. The model is described by the Hamiltonian

Here we defined ${\sigma }_j^x={(|r⟩⟨g|+|g⟩⟨r|)}_j$ and $n_j=|r⟩{⟨r|}_j$ is the occupation of the Rydberg state in the $j$-th atom. The Rydberg-Rydberg interaction is of the form $V_{jk}=\; C∕|j\; -k{|}^6$.

The model (3.1) supports several phases, including a disordered phase and various ordered phases with spontaneous breaking of translational symmetry. Here we work in the disordered phase and make the following choice of parameters for the numerical simulations:¹ $\mathrm{\Delta }=\; -0.5\mathrm{\Omega },C\; =\; 4\mathrm{\Omega }$. This point in the disordered phase can be adiabatically connected to the limit $\mathrm{\Delta }\to \; -\infty$ where the ground state is the empty state $|ggg\dots ⟩$ and quasiparticles correspond to Rydberg excitations. The approximate identification of quasiparticles with Rydberg excitations will be useful to make connections with the single-particle model in section 2.

3.1 State preparation and detection

To prepare a propagating wave packet, we introduce an auxiliary atom (labeled $j\; =\; 0$) positioned at the edge of the chain, at a distance $d$ from site $j\; =\; 1$. We choose $d$ to be larger than the spacing $a$ between the other atoms, such that the auxiliary atom is only weakly coupled with the rest of the system. This setup mirrors the role of the $j\; =\; 0$ site in the single-particle model (figure 1). We control the state of the auxiliary atom by modulating its time-dependent detuning ${\mathrm{\Delta }}_0(t)$, leading to the single-site Hamiltonian

This local control over ${\mathrm{\Delta }}_0(t)$ allows us to engineer the emission of a quasiparticle from the edge, analogously to the single-particle case. The interaction between the auxiliary atom and the chain is described by

where the dimensionless parameter $\kappa \; =\; (d∕a)\; -\; 1$ effectively tunes the coupling strength between the auxiliary atom and the chain. For large $\kappa$, the auxiliary atom is nearly decoupled from the chain and behaves as an isolated two-level system with an energy gap $V_0=\sqrt{{\mathrm{\Omega }}^2+{\mathrm{\Delta }}_0^2}$.

To prepare the wave packet, we initialize the auxiliary qubit in its (approximate) excited state, while the chain remains in its ground state. The detuning ${\mathrm{\Delta }}_0={\mathrm{\Delta }}_{\ast }$ is selected so that $V_0$ matches the energy of the desired quasiparticle excitation. As the system evolves, the weak coupling between the auxiliary site and the chain enables a gradual energy transfer from the qubit to the chain, exciting a propagating mode with energy $V_0$.

As in the single-particle case, the dynamics can be understood in terms of the conserved energy and energy variance. The initial state carries energy $V_0$ (relative to the ground state) and has a narrow energy variance ${(\delta E)}^2$ determined by the weak interaction $H_{\mathrm{\text{int}}}$. Consequently, only system eigenmodes with energy within $V_0\pm \; \delta E$ are significantly populated, leading to the formation of a localized wave packet with a well-defined momentum peak at $k_{\ast }$ and a small momentum variance.

An important caveat in this protocol is that the single-particle modes must be well separated from the continuum of multi-particle excitations. In other words, the mass gap should be sufficiently large to ensure that only a single quasiparticle is excited during the process. If this condition is not satisfied, the energy injected into the system may instead populate multi-particle states. In section 4, we will discuss how to prepare a wave packet even when the mass gap is small and the single-particle band is not well separated from the rest of the spectrum.

3.1.1 Time modulation of the Hamiltonian parameters

We now describe how to prepare the qubit in its excited state by modulating the (global) Rabi frequency $\mathrm{\Omega }(t)$ and the local detuning on the auxiliary spin ${\mathrm{\Delta }}_0(t)$. The Rabi frequency is varied only globally, and the time-dependent modulation of the local detuning that we propose is compatible with current experimental capabilities in Rydberg atom array platforms [43–47].

The system is first initialized in the global ground state with ${\mathrm{\Delta }}_0=\; -{\mathrm{\Delta }}_{\ast }$, using an adiabatic ramp of the Rabi frequency from $0$ to $\mathrm{\Omega }$, starting from the fully polarized state $|ggg\dots ⟩$. Since we operate in the disordered phase, this ramp does not encounter any phase transition, allowing the ground state to be prepared reliably.

Next, we apply a $\pi$-rotation around the ${\sigma }^z$ axis to the auxiliary site by evolving the system for a short time $\delta t$ with a large detuning ${\mathrm{\Delta }}_0=\; \pi ∕\delta t$. The evolution of the qubit’s state on the Bloch sphere is illustrated in figure 5(a). This pulse flips the sign of the qubit’s $x$-polarization, effectively preparing the ground state of the Hamiltonian $(-\mathrm{\Omega }{\sigma }_0^x+{\mathrm{\Delta }}_{\ast }{\sigma }_0^z)∕2$. This state is, by construction, the excited state of the target Hamiltonian $H_0=\; (\mathrm{\Omega }{\sigma }_0^x-{\mathrm{\Delta }}_{\ast }{\sigma }_0^z)∕2$.

3.1.2 Numerical results

We now verify, using matrix-product state (MPS) numerical simulations,² that the protocol we propose is preparing the desired wave packet, and that the momentum can be efficiently reconstructed with our detection method.

In our simulations we choose ${\mathrm{\Delta }}_{\ast }=\; -0.5\mathrm{\Omega }$, $\kappa \; =\; 0.5$, and $\delta t\; =\; 0.1{\mathrm{\Omega }}^{-1}=\; 0.1{(2\pi )}^{-1}$ cycles. We define $t\; =\; 0$ as the time at which the ${\pi }_z$ rotation is completed. The energy of the auxiliary qubit is $V_0=\sqrt{{\mathrm{\Omega }}^2+{\mathrm{\Delta }}_{\ast }^2}\approx \; 1.12\mathrm{\Omega }$, which corresponds to a target momentum $k_{\ast }\approx \; 2.3$ (see dispersion relation in figure 5(b), obtained as discussed in appendix B). We choose a relatively small value of $\kappa$, which results in a large energy variance and, consequently, in a sizeable momentum spread of the emitted wave packet. This choice is instrumental in generating a wave packet that extends over only $\mathcal{\text{𝒪}}(10)$ sites, suited for a proof-of-principle experimental implementation on small systems that are currently within reach. We will demonstrate this protocol for narrower momentum spreads in section 4.

As shown in figure 5(c), the auxiliary qubit, initially prepared in its excited state, gradually decays to the ground state. The probability $P_{\mathrm{\text{aux}}}$ of finding it in the excited state decays exponentially in time. As the auxiliary qubit decays, it releases an excitation into the system (figure 5(d)–(h)), which propagates from the left edge of the chain. A simple way to observe this propagation experimentally is to monitor the variation of a local observable, for instance the local Rydberg occupation $n_j$, between time $t$ and time $0$, when the system is locally in the ground state. As indicated by the white dashed line, the resulting wave packet propagates with the expected group velocity $v(k_{\ast })$. Upon reaching the right boundary, the wave packet is reflected, leading to pronounced spatial oscillations. As in the single-particle case, these oscillations provide a useful probe of the wave packet momentum. In figure 6 we therefore plot the absolute value squared of ${\mathcal{\text{F}}}_n(q,t)$, defined as the Fourier transform, with respect to the spatial coordinate, of $n_j(t)\; -n_j(0)$ (we again omit the site $j\; =\; 0$):

Similarly to the single-particle case, we observe that a peak develops when the wave packet is reflected at the boundary (at time $t\; \approx \; 26$). Because the momentum uncertainty is relatively large, this peak is not very sharp in momentum space and extends over a broad region around $q\; =\; 2k_{\ast }\approx \; 4.6$. Note that in figure 6 we plot only momenta in the interval $[0,\pi ]$; as a result, the peak appears at the equivalent momentum $q\; =\; 2\pi \; -\; 2k_{\ast }$. We also observe a similar temporal trend from higher-velocity components to lower-velocity ones, as in the single-particle case. Here, however, the faster components correspond to smaller momenta, as can be inferred from the dispersion relation in figure 5(b). Consequently, the peak around $q\; =\; 2(\pi \; -k_{\ast })$ moves from larger to smaller values of $q$ as a function of time.

In addition to this main feature, we identify a distinct signal at $q\; \approx k_{\ast }$, which is already visible before the wave packet reaches the boundary and persists after its reflection (see, for example, the profiles at $t\; =\; 18$ and $t\; =\; 34$). This peak is significantly smaller in amplitude than the one at $2(\pi \; -k_{\ast })$ observed at $t\; =\; 26$, but is observed systematically over a wide range of parameters ${\mathrm{\Delta }}_{\ast }$ and $\kappa$. This feature can be understood as follows. Suppose that the preparation of the wave packet is not perfect, so that the time-evolved state retains a small overlap with the ground state of the system. We may then write, schematically,

where $|g⟩$ denotes the ground state and $|e_k⟩$ is the eigenstate corresponding to a single-particle excitation with momentum $k$ and energy $E_k$. For the purpose of this argument, we treat the system as translationally invariant with periodic boundary conditions, which is justified since we are interested in the propagation of the wave packet in the bulk, far from the boundaries. We also assume that the population of higher-energy eigenstates involving two or more particles can be neglected.

Using translational invariance, the local operator $n_j$ can be related to the operator at a reference site $j\; =\; 0$ via the translation operator $T$ as

Together with the relations $T|g⟩=|g⟩$ and $T|e_k⟩=e^{ik}|e_k⟩$, this leads to

Since the coefficients $c_k$ are sizeable only within a narrow window around $k_{\ast }$, the first two terms contribute to $F_n(q,t)$ only at small values of $q$. The third term, instead, gives rise to a peak at $q\; \approx k_{\ast }$, which oscillates in time with a frequency set by $E_{k_{\ast }}$. This argument fully accounts for the features observed in figure 6, including both the low-$q$ signal and the oscillating peak at $q\; \approx k_{\ast }$.

Although this latter peak originates from an undesired effect — namely the imperfect preparation of the wave packet and its residual overlap with the ground state — it can nevertheless be exploited as an alternative probe of the momentum. Finally, we note that this discussion is not specific to the choice of the operator $n_j$, but is quite general and applies to generic local observables.

3.2 Scattering

Having established methods for both the preparation of wave packets and the detection of their momenta, we now turn to the simulation of a scattering process. An auxiliary qubit is placed at each end of the chain and prepared in its excited state, following the same protocol used for the preparation of a single wave packet. The results for a chain of length $L\; =\; 50$ are shown in figure 7, using the same set of parameters as in the previous subsection. The two wave packets propagate inward from opposite ends of the chain, collide near the center, and are subsequently reflected at the boundaries (figure 7(a)–(e)).

During the collision, conservation of energy and momentum restricts the dynamics to elastic scattering. As a result, two outgoing wave packets emerge from the scattering region with momenta $\pm k_{\ast }$, identical to those of the incoming wave packets. The outgoing momenta are detected by analyzing the Fourier transform of $n_j(t)\; -n_j(0)$, as described in the previous subsection (figure 7(f)–(i)). When the wave packets are reflected at the boundaries, a peak develops near $q\; =\; 2k_{\ast }$ (and, equivalently, at $q\; =\; 2\pi \; -\; 2k_{\ast }$), in agreement with the expectation for elastic scattering.

For elastic scattering of this type, the scattering matrix reduces to a single parameter, namely a phase shift. This phase shift can be extracted by measuring the time delay accumulated during the collision, for different values of the momentum $k_{\ast }$. As shown in figure 8, the time delay $\tau$ is obtained by comparing the free propagation of a single wave packet with the propagation of a wave packet that undergoes a collision with another wave packet.

4 Probing inelastic scattering at high energies

In the previous section, we demonstrated our methods for state preparation and detection, and showed how to extract the time delay associated with elastic scattering. A key requirement for applying this protocol was the presence of a sufficiently isolated mode, so that an energy filter (implemented by fixing the average energy and its variance) selectively targets a single-particle excitation with well-defined momentum. However, many physically interesting scattering regimes do not satisfy this condition. In particular, one often aims to probe scattering at relatively high energies, approaching a relativistic regime in which the bandwidth greatly exceeds the mass gap and strongly inelastic processes can occur.

In this section, we show how to address this limitation by employing adiabatic ramps, and demonstrate that this approach enables the study of inelastic scattering at high energies. To this end, we consider a model different from the Rydberg-atom chain discussed above, and instead focus on a system that supports multiple bands, or “species”, of single-particle excitations: the mixed-field quantum Ising chain. This model enables inelastic scattering between two-particle states with discrete allowed momenta, providing a clean setting in which such processes can be resolved. Furthermore, the mixed-field quantum Ising chain can be realized in several experimental platforms, including ultracold atoms in optical lattices [48], Rydberg atom arrays [49, 50] and, most notably, trapped-ion systems [51], where multiple quasiparticle excitations have already been observed [52]. The protocol proposed here is directly compatible with modern analog trapped-ion architectures featuring local control [53, 54]. Scattering in the same model has also been investigated on digital platforms based on superconducting qubits [40], providing a natural benchmark for comparison with analog implementations.

We note that, for a different choice of parameters, the Rydberg-atom Hamiltonian considered in the previous section can also support inelastic scattering processes in which two incoming particles scatter into three outgoing particles. However, in that case the allowed final-state momenta form a continuum, which significantly complicates their detection. For this reason, we defer a detailed study of such processes to future work.

In this section, we also introduce a slightly different preparation scheme, in which the single-site terms are tuned in a site-dependent manner while the interaction terms remain uniform across the chain, without the need for weak-interaction links. This variation may be advantageous for a range of experimental platforms.

4.1 Ramp protocol and inelastic scattering at the boundary

To demonstrate the ramp protocol, we consider here the mixed-field quantum Ising chain, defined by the Hamiltonian

For $J\; >\; |g|$ and nonzero $h$ (which we here set to $h\; =\; 0.1J$), the low-energy spectrum of this model consists of a tower of single-particle excitations. These excitations are known as mesons, as they correspond to confined states of two domain walls (which are the elementary excitations of the model for $h\; =\; 0$) [55–57]. The bands corresponding to different meson species are well separated for small $g$, and increasingly overlap as $g$ approaches $J$ from below (figure 9(a),(b)).

If the energy ranges of two bands $\ell$ and ${\ell }^{\prime }$ overlap, inelastic scattering processes can occur in which two $\ell$-mesons scatter into two ${\ell }^{\prime }$-mesons while conserving the total energy and total momentum $k_{\mathrm{\text{tot}}}=\; 0$. For example, for a specific choice of parameters, shown in figure 9(b), a $1$-meson with momentum $k_{\ast }$ has the same energy as a $2$-meson with momentum $k_s$. As a consequence, an incoming state consisting of two $1$-mesons with momenta $\pm k_{\ast }$ can scatter inelastically into a pair of $2$-mesons with momenta $\pm k_s$. Another possible inelastic channel, which may be present even if the two bands do not overlap in energy, is the conversion of two $1$-mesons into a final state consisting of a $1$-meson and a $2$-meson. These various inelastic processes are the phenomena we here aim to investigate.

Because the protocol relies on energy conservation, a scheme such as the one used above would generically prepare a superposition of states at the target energy (marked by a grey horizontal line in figure 9(b)) belonging to different bands ($\ell \; =\; 1$ and $\ell \; =\; 2$) and carrying different momenta ($k_{\ast }$ and $k_s$). To achieve a faithful preparation of an excitation within a single band and with a well-defined momentum, we therefore employ an adiabatic ramp. Specifically, we first prepare the wave packet at a smaller value of $g$, where the $1$-meson with momentum $k_{\ast }$ is well separated in energy from other bands (figure 9(a)). Once the wave packet has been completely emitted from the first site, we adiabatically ramp $g$ to its final value (figure 9(c)). Since $g$ is ramped uniformly throughout the system, momentum is conserved during the evolution, provided the wave packet remains far from the boundaries. Moreover, if the ramp is sufficiently slow, adiabaticity ensures that transitions between different bands are suppressed. As a result, we obtain the desired wave packet: an excitation confined to a single band with a narrow momentum distribution.

We demonstrate this protocol for the Ising chain defined in eq. (4.1) by coupling the system to an auxiliary site, following the approach introduced in the previous sections. The full Hamiltonian used for state preparation reads

The Ising coupling between the auxiliary site and the first site of the chain is taken to be equal to $J$, matching the couplings within the bulk of the chain, while the on-site fields $g_0$ and $h_0$ acting on the auxiliary site are distinct and can be tuned independently. A model of this type, with local time-dependent control of the transverse and longitudinal fields, has been demonstrated, for example, in recent trapped ion experiments [53, 54].

In this setup, the auxiliary qubit is not weakly coupled to the system. Nevertheless, the state-preparation protocol can still be applied, provided that one can prepare a state with the desired energy and a sufficiently small energy variance. One possible approach is to first prepare the ground state of the system with $g_0=\; 0$ and a large negative $h_0$, which initializes the auxiliary qubit in the $\uparrow$ state. Subsequently, one performs a sudden quench of the on-site fields, setting $g_0$ to a small but nonzero value and $h_0$ to a positive value: after this quench, the energy variance is $g_0^2$, while the energy can be controlled by tuning the value of $h_0$.

In our numerical simulation, this first step of the wave packet preparation is carried out at $g∕J\; =\; 0.5$. The parameters $g_0=\; 0.09$ and $h_0=\; 0.95$ are chosen so as to match the energy of a $1$-meson with momentum $k_{\ast }\approx \; 1.28$, while ensuring a small momentum spread of the prepared wave packet (figure 9(a)).

This state preparation protocol is demonstrated in figure 9(d), where we plot the excess energy density after the quench. The latter is defined as the expectation value of the local energy-density operator

from which we subtract the expectation value of the same operator evaluated in the ground state. As shown in figure 9(d), the wave packet propagates with the expected group velocity.

Once the wave packet has been fully emitted from site $j\; =\; 0$, we slowly ramp the transverse field $g$ from $0.5$ to $0.9$ (figure 9(c)), thereby entering a regime in which inelastic scattering processes are allowed. As a consequence of the ramp, the group velocity of the wave packet changes, matching the group velocity at the target momentum in the final Hamiltonian. We then let the wave packet scatter off the boundary to probe the inelastic process converting a 1-meson into a 2-meson. Importantly, no symmetry-based selection rules forbid such transitions. To characterize the resulting dynamics, we analyze the Fourier transform of a local observable. In this case, we focus on the operator ${\sigma }_j^x$, as we expect it to have sizeable matrix elements between different energy eigenstates. Such off-diagonal matrix elements play an important role in this detection scheme, as already implied by eq. (3.7).

In this signal, shown in figure 10, we observe three distinct peaks: one located at $q\; =\; 2k_{\ast }$, as expected for elastic scattering, and two peaks at $q\; =k_{\ast }\pm k_s$, demonstrating that in the collision with the boundary the $1$-meson with momentum $+k_{\ast }$ was partially scattered into a $2$-meson with momentum $-k_s$. Similar peaks at these values of $q$ are observed for other observables, such as ${\mathcal{\text{E}}}_j$ or ${\sigma }_j^z$. However, in these cases the signal is noticeably less clean than for ${\sigma }_j^x$.

4.2 Inelastic two-particle scattering

Having demonstrated the state-preparation protocol based on adiabatic ramps, and after testing the inelastic scattering from a collision with a boundary, we now analyze the case of a two-particle inelastic scattering process.

We employ the same ramp-based state preparation scheme as in the previous subsection, but we apply it at both ends of the chain to generate two incoming wave packets. The ramp brings the system into a regime where several inelastic channels are kinematically allowed (see appendix B for a detailed analysis of the kinematics). The outcome of the collision is shown in figure 11. In addition to the elastic channel with two $1$-mesons in the final state, we observe several outgoing traces consistent with inelastic processes: (i) a final state with a $1$-meson and a $2$-meson, (ii) a final state with two $2$-mesons, (iii) a continuum of states corresponding to three $1$-mesons. The contributions of channels (i) and (ii) can be distinguished from the elastic channel by their characteristic group velocities, which can be predicted from a simple kinematic analysis (see appendix B). We note, however, that the $(1\; +\; 3)$ channel is also kinematically allowed. For the system sizes considered here, its contribution cannot be resolved because its group velocity is too similar to those of the other inelastic channels, and we therefore cannot exclude its presence. The kinematically allowed three-$1$-meson states involve mesons with very small momenta and, consequently, very small group velocities. As a result, these excitations produce nearly vertical traces in the spacetime plot in figure 11.

5 Controlled shaping of the wave packet

While having a sharp momentum distribution, the wave packets produced with our method are not Gaussian. They tend to have a rather irregular shape and long tails. This is a drawback that adversely affects the sensitivity of detection. A wave packet with long tails corresponds to smaller local amplitudes in the relevant observables, like the excess energy density. Moreover, since in our ramp scheme we need to wait until the wave packet has been emitted completely before applying the ramp, we may have to wait for a very long time if the wave packet has a long trailing tail. In this section we show how to address this limitation and improve the shape of the wave packet.

To achieve this goal, we take inspiration, once again, from the simple single-particle model presented in section 2. Even in that case, our protocol produced wave packets with long tails and irregular shapes with wiggles, rather than Gaussian wave packets. A protocol that has previously been shown to produce nearly Gaussian wave packets is the so-called “quantum slide”. The idea is to modulate the parameters of the model in a region near the edge of the chain: while in section 2 the hopping amplitude and local potential have a sudden jump between $j\; =\; 0$ and $j\; =\; 1$, in a slide the parameters vary more smoothly in a region near the edge. A careful choice of the parameters along the slide is guaranteed to produce an approximate Gaussian wave packet (see appendix C and refs. [58, 59]).

We here would like to employ a similar scheme to control the shape of the wave packets in the many-body case, as exemplified by the mixed-field quantum Ising chain. In contrast with the single-particle case, where the shape of the slide can be justified analytically, in the many-body case we do not have a similarly rigorous approach. We therefore choose a simple linear slide of size $R$, where the (site-dependent) transverse and longitudinal fields take the form³ (see inset in figure 12(b))

Here $g$ and $h$ are chosen to match the values in the interior of the system, while $h_0$ and $g_0$ are the (post-quench) values of the longitudinal field and transverse field on the auxiliary site in eq. (4.2).⁴

The results of the wave packet preparation using the previous protocol (without the slide) and the modified protocol with a slide of size $R\; =\; 5$ are compared in figure 12(a),(b). In the slide protocol, the wave packet is significantly narrower, with a larger local amplitude and a more regular shape. As shown in figure 12(b), in the slide protocol the wave packet is fully emitted from the first site at a time when, in the simple protocol, it has not yet been emitted. This represents a major practical advantage, as it substantially reduces the waiting time before applying the adiabatic ramp.

5.1 Analysis of the quantum slide method

To gain a theoretical understanding of these numerical observations, we analyze the time required to eject the wave packet from the auxiliary site. This timescale, multiplied by the group velocity, sets the spatial extent of the wave packet. It cannot be made arbitrarily short, since the minimal spatial extent is constrained by the inverse of the momentum uncertainty. Our goal is therefore to minimize the decay time of the auxiliary site, subject to a fixed energy variance (which determines the momentum uncertainty). This optimization leads to a wave packet that is more localized in real space without increasing its momentum spread, thereby yielding a profile closer to a Gaussian wave packet.

We now show how to compute the decay rate. As before, the system is initialized in the ground state with a large negative $h_0$: the ground state $|{\psi }_0⟩$ of the combined system and auxiliary site has the latter in the $\uparrow$ state, while the system is in the ground state of $H\; -\; J{\sigma }_1^z$, i.e., the ground state subject to the boundary condition imposed by the auxiliary site. We denote this state by $|G{S}_{\uparrow }⟩$:

The combined system is thus in the initial state

and evolves after the local quench under $H_{\text{prep}}=H^{\prime }+\; V$, where

and

The initial state $|{\psi }_0⟩$ is an eigenstate of $H^{\prime }$ with energy $E_0=h_0+{𝜖}_{\uparrow }^{(GS)}$ (although not its ground state), and $V$ can be treated as a weak perturbation. The spectrum of $H^{\prime }$ separates into two sectors, depending on whether the auxiliary site is in the $\uparrow$ or $\downarrow$ state. In the post-quench Hamiltonian, the parameters are chosen such that the ground state lies in the $\downarrow$ sector. For ease of notation, we set the zero of the energy as the ground state energy of $H^{\prime }$, i.e., $-h_0+{𝜖}_{\downarrow }^{(GS)}=\; 0$ (note that this choice needs to be accounted for in the definition of ${𝜖}_{\uparrow }^{(GS)}$). Here $-h_0$ is the energy of the auxiliary site in the $|\downarrow ⟩$ state, ${𝜖}_{\downarrow }^{(GS)}$ is the ground state energy of $H\; +\; J{\sigma }_1^z$ (i.e., the Hamiltonian of the system subject to the boundary condition corresponding to having the auxiliary site in the $\downarrow$ state).

As before, a key requirement of our protocol is that the energy $E_0$ of $|{\psi }_0⟩$ matches the energy of the single-particle excitation at the desired momentum. The decay rate of the auxiliary excitation into a bulk excitation can then be estimated using Fermi’s golden rule:

where $|E⟩$ are eigenstates of $H^{\prime }$ with energy $E$, and $\rho (E)$ is the density of states.

On the other hand, the energy variance of the initial state is

where we have used that $⟨{\psi }_0|V\; |{\psi }_0⟩=\; 0$. Comparing eqs. (5.7) and (5.8), we see that enhancing the decay rate $\mathrm{\Gamma }$ while keeping the energy variance fixed requires $V$ to couple $|{\psi }_0⟩$ preferentially to states at the target energy. In the standard language of a small system (here, the auxiliary qubit) coupled to a large bath (the chain), this amounts to requiring that the spectral density of the system-bath coupling be peaked at the energy of the target excitation. Although this argument was presented in the context of the Ising model, the same considerations apply more generally to this class of protocols.

We now examine the decay rate explicitly in this model. The relevant matrix elements are those between $|{\psi }_0⟩$ and states in the opposite sector, $|E⟩=|\downarrow ⟩\otimes |{𝜖}_{\downarrow }⟩$, where $|{𝜖}_{\downarrow }⟩$ is an eigenstate of $H\; +\; J{\sigma }_1^z$ with eigenvalue ${𝜖}_{\downarrow }(E)\; =\; E\; +h_0$ (such that $H^{\prime }|E⟩=\; E|E⟩$):

The corresponding decay rate is therefore

where ${\sigma }^2=⟨{\psi }_0|V^2|{\psi }_0⟩=\; |g_0{|}^2$ is the energy variance after the local quench, and

Here ${\rho }_{\downarrow }$ is the density of states of the Hamiltonian $H\; +\; J{\sigma }_1^z$, and ${𝜖}_{\downarrow ,0}={𝜖}_{\downarrow }(E_0)\; =E_0+h_0$. Note that $\gamma (E_0)$ does not depend on $g_0$, and therefore the decay factor $\gamma (E_0)$ can be tuned independently from the energy variance.

The slide protocol modifies the factor $\gamma (E)$, effectively enhancing the spectral density at the target energy $E_0$, as shown in figure 12(c) (see appendix D for details on the calculations). As a result, for a fixed energy variance ${\sigma }^2=\; |g_0{|}^2$, the decay rate of the auxiliary qubit can be significantly increased, as illustrated in figure 12(d).

We now comment on the generality of this result and on how to design quantum slides more broadly. The key intuition is that the parameters should vary smoothly in space, rather than exhibiting the discontinuities present in the absence of a slide. When the variation is sufficiently smooth, a local density approximation can be employed to analyze how the dispersion relation depends on position. Within this framework, the slide should be designed so that the energy after the local quench is never in a gapped region of the (local) spectrum. This ensures that an excitation can propagate resonantly across the slide. A detailed investigation of this approach and of the optimal slide design is left for future work.

5.2 Scattering

We now apply the quantum slide method for wave packet shaping, as demonstrated in figure 12, followed by the adiabatic ramp, and simulate the two-particle scattering in the same regime as in figure 11. The results are shown in figure 13. The scattering channels remain visible in this case; however, the inelastic channels are less prominent. The origin of this suppression is not yet understood and is left for future investigation.

6 Higher dimensions

Our protocol can be straightforwardly extended to higher dimensions. For instance, auxiliary sites can be placed at opposite corners of a square lattice. As in the one-dimensional case, we use conservation of energy and energy variance to constrain the momentum of the propagating wave packet.

A key difference from the one-dimensional case is that energy conservation no longer selects a single momentum $\overrightarrow{k}$, but rather a manifold of momenta. In systems with emergent rotational symmetry, this constant energy contour reduces approximately to a fixed-$|\overrightarrow{k}|$ arc at low energy. This is illustrated in figure 14 for the two-dimensional tight-binding model with hopping amplitude $w$ along both the $x$ and $y$ directions.

In the initial state, the particle is localized at a corner site. The hopping amplitudes from this site to its neighbors, denoted $w^{\prime }$, can be tuned to control the momentum spread of the resulting wave packet. The onsite energy of the corner site determines the energy of the state, and hence the target value of $|\overrightarrow{k}|$.

As shown in figure 14(a), the wave packet propagates from the corner toward the center of the system. During this evolution, its momentum distribution becomes peaked along an arc with fixed $|\overrightarrow{k}|\; \approx \; 1.27$ and arbitrary angle (figure 14(b)). These momentum values correspond to those satisfying $E(\overrightarrow{k})\; \approx V_0$. As shown in appendix E, the momentum components can be measured using a method analogous to that employed in the one-dimensional case, based on reflection from a boundary.

This protocol could be optimized to achieve a better control of the momentum distribution and target a specific angle. Possible strategies include modifying the lattice geometry, by attaching a one-dimensional system that acts as a waveguide, or by introducing a funnel-like structure at the corner. Alternatively, a site-dependent potential could be used to bias the particle toward a given direction. A demonstration of the latter scheme is shown in figure 15(a), where the focusing potential is illustrated via its equipotential lines. The potential near the corner is chosen to be flat along a diagonal region and to increase away from it, thereby energetically favoring propagation in the diagonal direction. This funnel-like profile gradually opens up, so that the emitted wave packet eventually evolves in a homogeneous two-dimensional region. The resulting wave packet is less dispersed in the transverse direction and exhibits a more sharply defined propagation angle in momentum space (figure 15(b)), leading to improved momentum resolution.

7 Conclusions

We propose and numerically validate a protocol for the preparation and detection of particle collisions. The protocol is compatible with currently available quantum simulation platforms and can be implemented on both analog simulators and digital simulators via Trotterization.

A key consideration is the resource cost associated with the implementation. In our scheme, the time required to generate a wave packet is fundamentally limited by the desired momentum resolution. Specifically, achieving a small momentum uncertainty $\delta k$ requires a wave packet of spatial width $W\; \gtrsim \; \delta k^{-1}$. Since the packet has to be emitted from the boundary, the corresponding preparation time scales as $\mathrm{\Delta }t\; \sim \; W∕v\; \gtrsim {(v\delta k)}^{-1}$, where $v$ denotes the group velocity at the target momentum. In general, the width $W$ is larger than the optimal value $\sim \; \delta k^{-1}$, but we have shown that the use of a quantum slide enables the preparation of significantly narrower wave packets, approaching the optimal width.

Another contribution to the preparation cost is the adiabatic ramp, which is used to reach the regime of inelastic scattering, where the mass gap is much smaller than the bandwidth. Accessing inelastic channels with $n$ particles in the final state requires a bandwidth-to-mass-gap ratio of order $n$, making small mass gaps necessary for probing highly inelastic processes. However, as the mass gap decreases, maintaining adiabaticity becomes increasingly challenging. The adiabatic condition requires the rate of change of the control parameter to scale no faster than the square of the instantaneous energy gap. Consequently, preparing states in the small-gap regime demands progressively slower ramps, making the adiabatic evolution a potentially significant contribution to the overall state-preparation time. (Note, however, that this limitation does not arise specifically from the wave packet preparation itself, but rather from the more general problem of ground-state preparation.) Moreover, this strategy is only applicable provided the target regime is adiabatically connected to a parameter regime with an isolated single-particle band. In particular, the interpolation path must avoid quantum phase transitions or other gap-closing points that would invalidate adiabatic state preparation.

Another important resource is the system size. While elastic scattering processes can be investigated in relatively small systems, the observation of inelastic scattering is considerably more demanding. To unambiguously identify the different scattering channels, the outgoing particles must propagate sufficiently far apart after the scattering event so that the final states can be resolved. This limitation is not specific to our protocol; rather, it reflects a fundamental requirement of scattering theory, namely that asymptotic particle states must be well separated in order to define and characterize scattering outcomes.

These considerations suggest that generating wave packets at the two ends of a large system may not always be the most efficient strategy. Since the packets must travel a substantial distance before colliding, the protocol requires a correspondingly long evolution time. A potentially more efficient approach is to create the wave packets closer to the center of the chain by introducing temporary, fictitious boundaries. Once the desired wave packets have been prepared, these boundaries can be removed adiabatically, allowing the outgoing particles to evolve in the full system after the collision. Such a scheme could substantially reduce the overall state-preparation time.

Additional reductions in the required preparation time may be achieved through further optimization of the protocol. In particular, one could improve the design of the adiabatic ramp used to generate the wave packets or refine the implementation of the quantum slide to approach the optimal packet shape more efficiently. Exploring such optimizations is an interesting direction for future work and could further enhance the practicality of scattering experiments on quantum simulators.

Beyond the study of asymptotic states and the scattering matrix, investigating the transient dynamics immediately after a collision is another interesting direction. In particular, the emergence of hydrodynamic behavior as a localized lump of energy spreads and cools raises important questions in the high-energy regime, already before the lump fragments into well-separated particles. Such physics can be accessed at relatively short times after the collision and represents a promising avenue for near-term investigation on quantum simulation platforms.

Finally, testing our protocols in a two-dimensional many-body system is another interesting direction for future work. While tensor-network simulations in higher dimensions are considerably more computationally demanding than in one dimension, recent studies have demonstrated the simulation of wave packet dynamics on lattices as large as $24\; \times \; 24$ [23]. Such system sizes may already be sufficient to verify the preparation of wave packets that are relatively narrow in momentum space.

Acknowledgments

We acknowledge useful discussions with Zohreh Davoudi, Manuel Endres, Roland Farrell, and Alessio Lerose. This work was supported by RIT (Research IT, Trinity College Dublin). FMS and JP acknowledge support provided by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, (DE-SC0025572); DOE QuantISED program through the theory consortium “Intersections of QIS and Theoretical Particle Physics” at Fermilab; DOE National Quantum Information Science Research Centers, Quantum Systems Accelerator; and by the Institute for Quantum Information and Matter, an NSF Physics Frontiers Center (PHY-2317110). FMS acknowledges support from Amazon Web Services, AWS Quantum Program. This research was supported in part by grant NSF PHY-2309135 to the Kavli Institute for Theoretical Physics (KITP). SB acknowledges support from Caltech’s Summer Undergraduate Research Fellowship (SURF) program. Calculations were performed using the TeNPy Library (version 1.0.0) [60].

Appendix

A Exact solution for the single-particle model

The eigenstates of the Hamiltonian in eq. (2.2) can be obtained using the ansatz

where $k\; \in \; [0,\pi ]$. For $j\; >\; 1$ this state satisfies the Schroedinger equation

with $E_k=\; -2w\cos <!--\; FUNCTION\; APPLICATION\; -->k$. We then have to impose that it also satisfies the Schroedinger equation for $j\; =\; 0,1$:

from which we get

To solve this system of equations, we note that the conservation of probability implies $|A_k{|}^2=\; |B_k{|}^2$, since $B_k$ and $A_k$ represent the amplitudes of an incoming and a reflected wave, respectively. We can then set $|A_k|\; =\; |B_k|\; =\; 1$ (fixing the normalization of the propagating wave) and define $z_k$ such that $A_k=\; -i{e}^{i{z}_k}$ and $B_k=A_k^{\ast }$ (fixing the global phase of the wavefunction). From eq. (A.6) we obtain

The condition in eq. (A.5) then results in the following equation for $z_k$:

Additionally, for $|V_0|\; \le \; 2|w|-\frac{{(w^{\prime })}^2}{|w|}$ there are no bound states (i.e., no eigenstates that decay as $e^{-\kappa j}$ for large $j$), so all energy eigenstates are of the form in eq. (A.1).

We can use the fact that $|{\psi }_k⟩$ form a complete orthogonal set and $⟨{\psi }_{k^{\prime }}|{\psi }_k⟩=\; 2\pi \delta (k\; -k^{\prime })$ to write the initial state $|\psi (t\; =\; 0)⟩=|0⟩$ as

The distribution $|C_k{|}^2=\; 4{\left(\frac{w}{w^{\prime }}\right)}^2{\sin <!--\; FUNCTION\; APPLICATION\; -->}^2{z}_k$ then corresponds to the expected momentum distribution at long times.

B Dispersion relations and kinematic analysis

We here show the numerical results about the dispersion relations in the Rydberg atom chain and the mixed-field quantum Ising chain.

B.1 Rydberg atom chain

We used exact diagonalization to compute the momentum-resolved low-energy spectrum of the Hamiltonian in eq. (3.1) with periodic boundary conditions. The results are shown in figure 16(a) for different system sizes $L$, and appear to be well converged with $L$, exhibiting no visible finite size effects. We used the single-particle band obtained with this method to further compute the 2-particle and 3-particle continua (figure 16(b)): these data show that, for the parameters chosen in the numerical simulation, the 2-particle and 3-particle continua do not overlap, so no inelastic channel is available in the collision of two particles.

B.2 Mixed-field quantum Ising chain

Together with exact diagonalization, the dispersion relations of the meson bands in the mixed-field quantum Ising chain can be computed with high accuracy using a two-fermion approximation. This approach is expected to be reliable in the regime of weak longitudinal field $h$. Here we briefly summarize the method and refer to refs. [61–65] for a detailed derivation.

The key observation is that the Hamiltonian is exactly solvable when $h\; =\; 0$. Applying a Jordan-Wigner transformation followed by a Bogoliubov transformation maps the model to free fermions,

where $\gamma (𝜃)$ are fermionic quasiparticle operators and $\omega (𝜃)$ is the single-particle dispersion relation. We then introduce the real-space Bogoliubov fermions

Within the two-fermion approximation, the $\ell$-th meson state with momentum $k$ is described by the ansatz

where $|0⟩$ denotes the ground state of $H_0$. The function $f_{\ell }(n)$ describes the relative wavefunction of the two fermions, while the index $\ell$ labels the different meson bands. It is determined by requiring that $|{\mathrm{\Psi }}_{k,\ell }⟩$ be an eigenstate of the full Hamiltonian projected onto the two-fermion sector.

This condition leads to the effective eigenvalue equation

where

and $M\; ={(1\; -g^2)}^{1∕8}$ is the spontaneous magnetization of the transverse-field Ising model in the limit $h\; \to \; 0$. The term proportional to $hMn$ acts as a linear confining potential between the two fermions, leading to the formation of mesonic bound states. Solving the above eigenvalue problem yields the meson dispersion relations $E_{\ell }(k)$.

The results of the two-fermion approximation are compared with exact diagonalization in figure 17(a),(b). In the regions where the exact diagonalization spectrum appears converged with system size, the agreement between the two methods is excellent. This confirms that the low-energy excitations are accurately captured by the two-fermion description.

At $k\; =\; 0$, the exact diagonalization spectrum contains an additional state whose energy lies within the range spanned by the single-meson bands. This state can be interpreted as a false vacuum state [61, 63, 66], whose energy (approximately $2hML$) grows proportionally to the system size, indicating that it does not belong to the low-energy spectrum in the thermodynamic limit.

Additional finite-size effects are visible in regions where the meson bands overlap with the continuum. This behavior is expected, since mesons are no longer stable in this regime and can hybridize with multi-particle states. As a result, the finite-size spectrum exhibits stronger size dependence, and the assumptions underlying the two-fermion approximation become less accurate.

C Quantum slide in a single-particle model

We now illustrate how to prepare a nearly ideal Gaussian wave packet in the single-particle regime using a quantum slide. The construction follows the approach of refs. [58, 59]. Both the hopping amplitudes and the local potential are made site dependent, with profiles $w_j=\sqrt{(j\; +\; 1)(2R\; -\; j\; -\; 4)}$ and $V_j=V_0\left(1\; -\frac{j}{R}\right)$ (figure 18(c)).

The hopping profile is inspired by perfect-state-transfer protocols [67], in which the hopping Hamiltonian can be identified with the operator $H\; \propto S_x$ of a fictitious spin $S\; =\; (L\; -\; 1)∕2$. In such protocols, the state evolves into an approximately Gaussian distribution centered at the middle of the chain at half of the state-transfer time. The quantum-slide construction adapts this idea by engineering a hopping matrix that resembles one half of the spin Hamiltonian. As shown in refs. [58, 59], this modification still produces wave packets with an approximately Gaussian spatial profile.

We test this protocol using the parameters of figure 3. To obtain a wave packet with the same variance as the one considered there, the slide must contain $R\; =\; 200$ sites. The resulting wave packet dynamics, with and without the quantum slide, are shown in figure 18(a),(b). Without the slide, the wave packet displays large oscillatory tails and significant deviations from a Gaussian profile. In contrast, the slide produces a wave packet that remains approximately Gaussian throughout its evolution, undergoing only a gradual broadening as it moves along the chain.

D Numerical estimation of the decay rate

We now discuss how the decay factor $\gamma (E)$ introduced in section 5 is estimated numerically.

Using exact diagonalization, we compute the low-energy spectrum of the Hamiltonian $H_{\downarrow }=\; H\; +\; J{\sigma }_1^z$. We then evaluate the overlaps between $|G{S}_{\uparrow }⟩$ — the ground state of $H\; -\; J{\sigma }_1^z$ — and the eigenstates $|{𝜖}_{\downarrow }⟩$ of $H_{\downarrow }$. The low-energy excited states exhibiting large overlap with $|G{S}_{\uparrow }⟩$ (shown as white dots in figure 19) correspond to the single-particle excitations populated during the state-preparation protocol.

The energies of these states, displayed in the insets, are used to estimate the density of states, $\rho \; ={(dE∕dn)}^{-1}$, through a linear fit of the excitation energies as a function of their index. This procedure assumes that the density of single-particle states is approximately constant within the relevant energy window. The validity of this assumption is supported by the quality of the linear fit. Under this approximation, $\rho$ depends only on the system size $L$.

We then plot the quantity $\rho |⟨{𝜖}_{\downarrow }|G{S}_{\uparrow }⟩{|}^2$ as a function of the excitation energy $E\; ={𝜖}_{\downarrow }-h_0$ (figure 19). The resulting data are fitted with a parabolic function for $R\; =\; 0$ and with a Gaussian function for $R\; =\; 5$. The fitted curves exhibit only weak dependence on system size, suggesting that the extracted estimates are representative of the thermodynamic-limit behavior.

E Momentum detection in 2D

We now present numerical results for the state-detection protocol in two dimensions. As in the one-dimensional case, the method relies on observing the interference between the incoming and reflected wave packet after reflection from a boundary. In two dimensions, however, an additional complication arises: depending on the orientation of the boundary, the interference pattern probes different projections of the momentum vector $\overrightarrow{k}$.

We consider a wave packet prepared as described in figure 14(a) of the main text. Figure 20(a) shows its reflection from a diagonal wall. Upon reflection, the component of the momentum perpendicular to the wall changes sign, generating a characteristic interference pattern. Taking the Fourier transform of the probability distribution $P(x,y)$ along both spatial directions reveals a peak at $\overrightarrow{q}\approx \; 2k_{\ast }\widehat{n}$, where $\widehat{n}$ is the unit vector normal to the wall (figure 20(b)). More generally, this approach provides a measurement of the momentum component $k_{\perp }=\overrightarrow{k}\cdot \widehat{n}$ perpendicular to the reflecting boundary. In the present example, the wave packet is emitted predominantly along the diagonal direction, so the observed peak provides a good estimate of the momentum magnitude $k_{\ast }$.

Other boundary orientations can be used to probe different momentum components. As an illustration, figure 20(c) shows the interference pattern produced by reflection from a wall normal to the $y$ direction. In this case, the Fourier-space signature is less pronounced, but a peak can still be identified near $\overrightarrow{q}\approx \; 2({\overrightarrow{k}}_{\ast }\cdot ŷ)ŷ$ (figure 20(d)). This behavior is consistent with the fact that the interference pattern is sensitive only to the component of the momentum perpendicular to the reflecting boundary.

References

Notes

^1. In our plots, the time is measured in units of a Rabi cycle, i.e. $2\pi ∕\mathrm{\Omega }$.

^2. Using the TenPy library, we employ, in particular, the two-site density matrix renormalization group (DMRG) method to prepare the initial ground state and the two-site time-dependent variational principle (TDVP) for the time evolution. In all the simulations reported in this paper, the scattering involves only few particles, and the entanglement growth is bounded, such that relatively small bond dimensions can be used, even for long-time evolutions. Where not explicitly stated, the bond dimension used in the simulations was $\chi \; =\; 200$.

^3. Explicitly, the Hamiltonian is $H\; =\; -J{\sum }_{j=1}^{L-1}{\sigma }_j^z{\sigma }_{j+1}^z+{\sum }_{j=1}^L(g_j{\sigma }_j^x+h_j{\sigma }_j^z)$.

^4. It may seem natural to choose a linear interpolation from $h_0$ to $h$ for the longitudinal field. However, one must take into account that the auxiliary site, unlike the sites in the chain, has only a single neighbor rather than two. Ensuring a smooth evolution of the excitation energy along the interpolation therefore requires accounting for this asymmetry, modeling such effect as an additional longitudinal field. To obtain a smooth potential for the meson, we then design a ramp that linearly interpolates between $h_0-\; 1$ and $h$. We find empirically that the choice of interpolation reported here performs well for our chosen value of $h_0$, although a different choice may be required for other values.