Scattering of mesons in quantum simulators

Simulating real-time evolution in theories of fundamental interactions represents one of the central challenges in contemporary theoretical physics. Cold-atom platforms stand as promising candidates to realize quantum simulations of non-perturbative phenomena in gauge theories, such as vacuum decay and hadron collisions, in prohibitive conditions for direct experiments. In this work, we demonstrate that present-day quantum simulators can imitate linear particle accelerators, giving access to S-matrix measurements of elastic and inelastic meson collisions in low-dimensional Abelian gauge theories. Considering for definiteness a $(1+1)$-dimensional $\mathbb{Z}_2$-lattice gauge theory realizable with Rydberg-atom arrays, we present protocols to observe and measure selected meson-meson scattering processes. We provide a benchmark theoretical study of scattering amplitudes in the regime of large fermion mass, including an exact solution valid for arbitrary coupling strength. This allows us to discuss the occurrence of inelastic scattering channels, featuring the production of new mesons with different internal structures. We present numerical simulations of realistic wavepacket collisions, which reproduce the predicted cross section peaks. This work highlights the potential of quantum simulations to give unprecedented access to real-time scattering dynamics.


Introduction
While implementing fault-tolerant quantum computations still requires significant technological advances, highly controllable quantum devices with hundreds of qubits are already being realized in various experimental platforms [1]. The possibility of accessing real-time dynamics and strongly correlated quantum many-body states opens numerous avenues in the theoretical research. One of the most promising directions is represented by the quantum simulation of high-energy physics phenomena [2,3,4,5,6]. In the last decade, substantial efforts have been devoted to the implementation of gauge-invariant Hamiltonian dynamics [7,8,9,10,11,12,13,14]. Recently, certain aspects of the physics of vacuum decay have been explored with trapped ions [15] and Rydberg atom arrays [16,17]. A challenging problem in high-energy physics is simulating collisions of complex composite particles. In quantum chromodynamics (QCD), a first-principle estimation of the distribution of particles produced by hadron scattering would facilitate the search for new physics beyond the Standard Model; moreover, heavy-ion collisions provide fundamental information on the deconfinement transition and on the early Universe evolution [18]. Although simulating higher-dimensional non-Abelian gauge theories is still a far-fetched goal, it is of great interest to understand whether quantum simulators are already capable of studying the scattering of composite particles in a strong coupling regime, at least in simplified settings. Lower-dimensional gauge theories [19,20] exhibit a tractable version of particle confinement leading to an analog of quark-antiquark bound states (mesons). Real-time dynamics of variants of these models witnessed recent developments in both classical [21,22] and quantum [17,15,14,23] simulations, opening the door to investigations of the simplest instances of collisions between complex structured objects arising from confinement.
In this work, we demonstrate that present-day quantum simulators allow to investigate selected meson collisions in 1 + 1-dimensional Abelian lattice gauge theories (LGTs), as sketched in Fig. 1, mimicking scattering experiments with particle accelerators. Quantum simulators offer unprecedented access to full real-time resolution of a complex collision event and to the quantum correlations thereby generated. Here, we particularly focus on the production of new mesonic species, i.e., inelastic events redistributing internal and kinetic energies of mesons emerging from the collision. We propose protocols to experimentally observe this with current facilities, and provide a benchmark theoretical study of scattering amplitudes. While we consider a controlled regime where exact numerical simulations can be pushed and compared with analytical results, quantum simulators may explore conditions inaccessible to traditional methods, including the continuum limit of quantum field theories.
The paper is organized as follows. In Sec. 2, we introduce the Z 2 -LGT analyzed throughout, and discuss particle confinement and the resulting mesonic spectra and wavefunctions in the regime of large particle mass. In Sec. 3, we give a theoretical study of meson-meson scattering amplitudes, based on an exact solution of the Schrödinger equation. We discuss elastic and inelastic processes, and benchmark the results against numerical simulations. Finally, in Sec. 4, we propose concrete protocols to prepare, simulate and observe meson scattering with present-day quantum simulators (e.g., Rydberg-atom arrays). The appendices contain various additional details on the discussion and computations in the main text. In Appendix A we report additional details on gauge invariance and confinement in the model under consideration in the main text. In Appendix B we prove the exact mapping of its dynamics in the gaugeneutral sector onto those of the quantum Ising chain in a tilted magnetic field. In Appendix C and Appendix D we provide more details on the exact solution of the twoand four-fermion problem, i.e., on the spectra of mesons and their scattering amplitudes, in the limit of large fermion mass. In Appendix E we derive the analytic expression of the meson current, we discuss its physical meaning and we prove the associated continuity equation. Finally, in Appendix F we summarize and discuss the effects of having a finite fermion mass.

Confinement and mesons
Particle confinement is a non-perturbative phenomenon arising in certain gauge theories, which consists in the absence of charged asymptotic states: all stable excitations of the theory above the ground state are neutral bound states of fermionic charges [24]. In the context of QCD, confinement underlies the fact that quarks can only be observed in composite structures such as mesons and baryons. Despite the fundamental difference between particle confinement in QCD in (3 + 1) dimensions and in lower-dimensional models [19,20], the emergent composite particles share some basic properties, making the latter convenient settings to gain insights into difficult aspects of the theory. In this work we will be concerned with (1 + 1)-dimensional LGTs of this kind. For the sake of definiteness, we will focus on the Z 2 -LGT defined by the following Hamiltonian [25,26]: In this equation, c † j and c j denote creation and annihilation operators of spinless fermions of mass m > 0 on the sites j ∈ Z of a one-dimensional lattice, and σ x,y,z b denote spin-1/2 operators (Pauli matrices) acting on the bonds b ∈ Z + 1/2 of the lattice, representing a gauge field with string tension τ . Interactions, with coupling strength w, are such that all the local operators G j = σ z j−1/2 σ z j+1/2 (1 − 2c † j c j ) are conserved, i.e., [G j , H] = 0. These operators satisfy G 2 j ≡ 1 and thus generate local Z 2 symmetries. Here we focus on the neutral gauge sector, i.e., the subspace characterized by G j |ψ = |ψ for all j.
(For more details, see Appendix A.) The LGT in Eq. (1) exhibits particle confinement for m > 2|w|, τ = 0. By gauge invariance, a string of excited gauge field extends between two charges created out of the vacuum, inducing a confining potential V (r) ∝ r that grows unbounded at large distances r. Thus, the excitations form a discrete tower of neutral bound states (termed mesons, in analogy with QCD), labelled by their internal quantum number = 1, 2, . . . and their center-of-mass momentum k.
In the large-m limit, mesonic spectra E (k) and wavefunctions ψ ,k (j 1 , j 2 ) can be determined exactly by solving the reduced two-body problem, governed by the projection of H in Eq. (1) onto the two-fermion sector spanned by states {|j 1 < j 2 } (labelled by the positions of the two fermions along the chain). For definiteness, we assume τ > 0 from now on. The projected Hamiltonian H 2-body consists of nearest-neighbor hopping terms of amplitude w for the two particles, plus the diagonal confining potential τ (j 2 −j 1 ). The problem is solved by switching to the center-of-mass and relative variables, s = j 1 + j 2 , r = j 2 − j 1 > 0, and using the ansatz ψ(s, r) = e iks φ k (r). The exact solution [27,25] yields the quantized mesonic spectra with = 1, 2, . . .,w k = 2w cos k, k ∈ [−π/2, π/2), and ν (x) is the -th (real) zero of the map a → J a (x), where J a (x) is the Bessel function ¶. The associated mesonic wavefunctions read As an example, Fig. 3-(a) below reports a plot of the lowest mesonic spectra = 1, 2, 3 for w/τ = 0.6. For the derivation of these results, see Appendix C. While suppressing dynamical fluctuations of the fermion number, the large-m limit encompasses the regime |w| |τ |, where mesons are dominated by strong quantum fluctuations in their spatial extent. A large but finite fermion mass m only produces a perturbative dressing of the vacuum and of mesons, which can be explicitly computed order by order via the so-called Schrieffer-Wolff transformation [28,29,30]. For instance, the first correction involves next-nearest-neighbor fermion hopping with Note that the Z 2 -charge is defined modulo 2, i.e., particle and antiparticle are the same object, so a two-fermion bound state is neutral.
amplitude w 2 /2m. Using this approach, the large-m analysis of meson dynamics can be systematically modified to achieve the desired accuracy for large but finite m, limited only by the practical complexity of high-order computations. Thus, for simplicity, here we will focus on the limit of large fermion mass; for more details on finite-m effects, see Appendix F.

Scattering amplitudes
We first provide a theoretical analysis of meson-meson scattering. We present an exact solution of the problem in the regime of large fermion mass m. The predictions of elastic and inelastic cross section peaks, together with our numerical simulations, provide a nontrivial benchmark for quantum simulations. While our solution is valid for arbitrary coupling strength w/τ and arbitrary incoming states, quantum simulations turn out to be easiest for w/τ ≈ 1 and low mesonic quantum numbers, as discussed below. Armed with the mesonic spectra E (k), we consider the scattering of two incoming mesons with quantum numbers 1,2 , approaching each other with definite momenta k 1,2 . The open elastic and inelastic scattering channels can be found by a kinematic analysis, which consists in determining the set of the outgoing quantum numbers {( 1 , k 1 ), ( 2 , k 2 )} compatible with the incoming ones by conservation of total energy and momentum: For all choices of incoming states, there always exist two elastic solutions, called transmitted and reflected , having ( 1 , 2 ) = ( 2 , 1 ) and ( 1 , 2 ) respectively. The existence of inelastic channels, instead, is not guaranteed for generic incoming states + . The conservation of the number of fermions allows to derive a continuity equation, which defines an associated mesonic current, as derived in Appendix E. The conservation of the total current across the collision yields a constraint on the scattering amplitudes of open channels. The fraction associated with each outgoing asymptotic solution has the physical meaning of a total cross section, as it can be identified with the probability P 1 , 2 of detecting that particular scattering outcome in the asymptotic future [31]. Determining the scattering amplitudes and cross sections requires solving the fourfermion Schrödinger equation. We thus consider the effective Hamiltonian H 4−body for the four-body problem, i.e., Eq. (1) projected to the four-fermion subspace spanned by the states {|j 1 < j 2 < j 3 < j 4 }. This consists of the hopping terms of amplitude w for the four particles, and the two diagonal confining pairwise potentials τ (j 2 − j 1 ) and τ (j 4 − j 3 ). We formulate the ansatz ψ 1 ,q (s 1 , r 1 )ψ 2 ,K−q (s 2 , r 2 ), where r 1,2 , s 1,2 are + For intermediate ratios τ /w (a condition that best suits experiments, see below), it can be seen that inelastic channels are favoured when at least one incoming meson is "heavy", i.e., 2 > 1. The reason is that, for sufficiently small w/τ , the sum of the quantum numbers 1 + 2 is conserved in the scattering. This is a consequence of the conservation of the total energy and of the fact that E( , k) τ in this limit. The example in Fig. 3 comprises an inelastic channel ( 1 , 2 ) = (1, 3) → ( 1 , 2 ) = (2, 2). the relative distance and the center-of-mass position for the two mesons, the singlemeson wavefunctions ψ are defined as in Eq. (C.7), but, crucially, the momentum q ∈ C is allowed to span the complex plane. The ansatz above represents an admissible asymptotic solution, provided q ∈ C satisfies the complex energy condition where the total energy E and momentum K are determined by the incoming state is here a complex zero of the analytic function z → J z (w), labelled by ∈ N. We index by α ∈ N all the triplets ( α 1 , α 2 , q α ) which satisfy Eq. 5. For a given incoming state, the exact solution Ψ of the scattering problem is expressed by a linear superposition including the incoming state and all compatible outgoing (i.e., with outgoing currents) and evanescent (i.e., with Im(q) < 0) asymptotic solutions: The wavefunction Ψ solves the Schrödinger equation in the full region which determines the coefficients A α , including the scattering amplitudes of open channels. In fact, this condition gives rise to an infinite set of inhomogeneous linear equations on varying r 1,2 = 1, 2, . . . for the infinitely many unknowns A 1 , A 2 , . . . The very nature of confinement, though, provides a natural truncation for this hierarchy: For q α ∈ R, the meson wavefunctions are bound states, and thus fall off rapidly for large distances; for complex solutions q α / ∈ R, the normalizability condition Im(q) < 0 guarantees exponential decay. Thus, asymptotic solutions with high mesonic quantum numbers 1,2 have tiny amplitudes, and their contribution is effectively redundant. For more details, see Appendix D.
The scattering amplitudes can be readily connected with the products of realtime wavepacket collisions. We verify this by numerically computing instances of the exact time evolution within the four-body subspace. We consider the example in Fig. 1: In a system with L = 36 fermionic sites, we prepare two Gaussian wavepackets    the Fourier transform of Ψ(s 1 , s 2 , r 1 , r 2 ; t) with respect to the center-of-mass positions s 1,2 . While the initial state shows a single density peak at (k 0 1 , k 0 2 ), the final state gives three different density peaks, all lying on the line demonstrating the conservation of total momentum. The three peaks correspond to the channels predicted from the kinematic analysis [the crosses in Fig. 3-(a)]: one for the trasmitted solution 2), and one for the inelastic solution (2, 2) (k 1 −1.5, k 2 −0.7). The peaks can be better resolved by plotting the distribution of the relative momentum k 1 − k 2 , as done in Fig. 3-(e). The relative weights enclosed within the dashed contours in Fig. 3

Quantum simulation
The analysis above outlines a tractable regime where non-trivial meson scattering phenomena can be accessed and understood. We now discuss how to observe them -and possibly extend their scope -with a quantum simulator, which minimally requires: i) designing the desired Hamiltonian dynamics; ii) preparing the incoming state; iii) detecting the outgoing states.
Crucially, problem i) does not involve any experimental fine-tuning: the basic phenomena only rely on confinement, and are thus robust to any weak perturbation to the model. As a concrete example, we will focus on simulating the Hamiltonian in Eq. (1) by exploiting the equivalence with the quantum Ising chain in a tilted magnetic field [25,32,33]: To obtain this non-trivial identification, one has to exploit the gauge-symmetry constraints to eliminate the fermionic degrees of freedom [34,17,25,26], as reported in Appendix B. Within this equivalence, Z 2 -mesons may be pictured as domains of reversed magnetization in a polarized background, arising from domain-wall confinement [35,36,37,38,39].
Recently, confinement in the quantum Ising chain has been shown to underlie a plethora of intriguing non-equilibrium phenomena [40,41,25,42,43]. The mapping above is extremely advantageous for quantum simulations because it implements gauge invariance exactly, similarly to what done in Refs. [15,17] for the Schwinger model. The dynamics governed by Eq. (7) can be experimentally realized both with optical lattices [44,45] and Rydberg atoms trapped in optical tweezers [46,16].
The preparation of the initial state ii) is subtle, as sharp meson wavepackets involve considerable entanglement between atoms, which is precluded to single-site optical manipulations. We present here an approach exploiting spatially inhomogeneous fields in Eq. (7) to filter meson wavepackets with sharply-defined momenta, at the price A smaller peak can be observed far from the momentum-conserving line k 1 + k 2 = k 0 1 + k 0 2 , and corresponds to the reflection of the second meson on the boundary after scattering in the (3, 1) channel. The missing probability fraction outside the dashed contours in Fig. 3-(c) is due to such effects as well as to the arbitrary cutoff used the define the contours, and amounts to ≈ 20% here.  of moderately longer chains and a limited amount of post-selection. The numerical simulation in Fig. 4-(a) illustrates the core idea: when w/τ 1, a spatially localized spin flip in the left region mostly excites the lowest (and fastest) meson = 1 at all momenta; hence, a sharp spatial variation in the fields τ (j), w(j) (inset) determines a corresponding change in the shape of mesonic bands [from that in panel (b) to that in (c)]; energy conservation (horizontal dashed lines) selects a narrow momentum window k * ± δk (vertical dashed lines) for which rightward propagation is allowed. Panel (d) shows that at time t = 50 the fraction of mesonic wavepacket filtered in the right region is ≈ 20% (the rest is reflected at the interface), and its momentum distribution has support within the selected window. An analogous preparation can be made on the opposite side of the chain for the desired incoming mesonic wavepacket from the right. Similarly, inhomogeneous fields can be used to accelerate mesons.
Finally, detecting the scattering products iii) is conceptually simple, as the mesons involved in the various possible outcomes of a collision have different velocities [cf. Fig.  1], so they can be resolved as spatially separate wavepackets. For implementations based on Eq. (7), the particle density c † j c j in Eq. (1) maps to the domain-wall density (1 − σ z j−1/2 σ z j+1/2 )/2: Thus, it is sufficient to measure the magnetization profile σ z j+1/2 (t f ) in the final state [44,45,46,16] to reconstruct the momenta of the mesons from their positions, the quantum numbers from their extension, and the cross sections from their probabilities. We note that the required time and length scales estimated from the above discussion (50 ÷ 100 lattice sites and units of time) are within reach of present-day experiments: Ref. [16], for example, demonstrated state preparation and single-qubit readout in a chain of 51 87 Rb atoms governed by Ising-type dynamics close to Eq. (7), with excellent coherence control over several tens of time units (2π/w).

Outlook
The analysis of the meson scattering problem and the proposed strategies for quantum simulations presented here can be straightforwardly applied to any one-dimensional model exhibiting confinement, including Abelian and non-Abelian lattice gauge theories (e.g., quantum link models [17,47]). They can also be extended to long-range interacting models, for which confinement effects [48,49] have been recently experimentally investigated with trapped ions [23]. The novel theoretical approach and exact solution to the meson scattering problem presented here will provide the basic building block for understanding the non-equilibrium evolution in quantum spin chains with confinement of excitations [40,48], particularly the recently reported lack of thermalization [42,41,25,50,51,49,52,53]. Compared to real-world scattering experiments, quantum simulations naturally give access to full real-time resolution of the dynamics of a complex collision event, and to the pattern of quantum correlations and entanglement at the level of partons [54,55], for which simplified lower-dimensional models such as the one discussed here could already provide deep insights. In future work, we plan to investigate this, as well as to optimize schemes for cold-atom platforms. Intriguingly, quantum simulators could allow to explore regimes beyond our theoretical analysis such as the continuum limit of quantum field theories [56,57,58,59,60]. This would represent a first step towards the ultimate goal of simulating realistic scattering problems in QCD such as heavy-ion collisions [61]. Note added -While completing the present manuscript, we became aware of a related work [62], appeared simultaneously.

Appendix A. Gauge invariance and confinement
In this work we have focused on the (1 + 1)-dimensional Z 2 -LGT defined by the Hamiltonian in Eq. (1) of the main text, reported here for convenience: Equation (A.1) describes gauge-invariant interactions of fermionic particles mediated by a Z 2 gauge field. The parameter m represents the mass of the fermions, τ quantifies the excitation energy cost per site of the gauge field excitation (string tension), and w is the coupling strength.
Gauge-invariance entails that a fermion hop or pair creation across a bond (j, j + 1) is always accompanied by a flip of the gauge field σ z j+1/2 on that bond, as illustrated in Fig. B1-(a). More formally, interactions are such that all the local operators G j = σ z j−1/2 σ z j+1/2 (1 − 2c † j c j ) are conserved, i.e., [G j , H] = 0. These operators satisfy G 2 j ≡ 1 and thus generate local Z 2 symmetries. Accordingly, the complete Hilbert space decomposes into dynamically disconnected subspaces, labelled by the set of eigenvalues {e iπq j = ±1} of {G j }, where q j = 0 or 1 is interpreted as the absence or presence of a static background Z 2 -charge on site j, respectively. In this work we have focused on the neutral gauge sector with q j ≡ 0, i.e., the subspace characterized by for all j. The meaning of this equation, referred to as the Gauss law, is that the gauge field spatial variations can only take place upon crossing a site on which a fermion is located, as shown in Fig. B1-(a). The LGT in Eq. (A.1) exhibits particle confinement in the regime m > 2|w|, |τ | > 0. The occurrence of fermion confinement in model (A.1) can be understood by considering the limit m → ∞. The ground state |GS becomes a fermion vacuum with c † j c j |GS = 0; by gauge-neutrality, [cf. Eq. (A.2)], the gauge field is uniformly polarized. The interaction term H int acting on the vacuum creates a pair of neighboring fermions. These particles can hop away from each other, with a hopping amplitude w. As they move, gauge-invariance forces a string of excited gauge field with tension τ to extend between them, which costs an energy per unit length equal to the string tension τ . Therefore,the two Z 2 -charges experience a confining potential V (r) = |τ |r that grows unbounded at large distances r whenever τ = 0. As a consequence, the excitations form a discrete tower of neutral (in a Z 2 theory, particle and antiparticle are the same object) bound states, labelled by their internal quantum number = 1, 2, . . . and their center-of-mass momentum k.

Appendix B. Mapping to the quantum Ising chain
Here we illustrate the mapping between the Z 2 -LGT in Eq. (A.1) and the quantum Ising chain in a tilted magnetic field, Eq. (7). This mapping was proposed in Ref. [25], and is connected with the one discussed in Ref. [26]. We consider the neutral sector, defined by the Gauss law (A.2), i.e., G j = σ z j−1/2 σ z j+1/2 (1 − 2n j ) ≡ 1, with n j = c † j c j representing the number of fermions at site j. The essence of the mapping is the exact elimination of the fermion degrees of freedom [63], as highlighted in Fig. B1-(b). The latter are redundant due to the infinitely many local constraints in Eq. (A.2): a classical configuration of the gauge field uniquely fixes the configuration of the fermionic matter via the Gauss law.
We present a formal proof of the equivalence. As a first step, we apply the Jordan-Wigner transformation to turn the fermions into hard-core bosons. To this aim, we introduce the spin-1/2 operators τ α j , for α = +, −, z, defined as As can be easily checked, the operators τ z j , τ x j = τ + j + τ − j and τ y j = −i(τ + j − τ − j ) satisfy the usual commutation relations of spin operators. By applying this transformation, we get the Hamiltonian with the constraint We now define a unitary transformation U that eliminates the matter degrees of freedom. In other words, we seek for a (gauge-variant) unitary U such that the transformed Gauss law G j = U G j U † ≡ 1 from Eq. (B.3) only depends on the matter degrees of freedom, whereas the transformed Hamiltonian H = U HU † only involves the gauge degrees of freedom. This can be accomplished with This transformation flips the spin τ z j where the neighbouring gauge fields are antialigned and does nothing where they are aligned. We get U τ z forces the τ α j spins to be polarized in the −ẑ directions. They enter the transformed Hamiltonian only via G j : Thus, in the neutral gauge sector, by Eq. (B.5) the spins τ α j are eliminated. Equation (B.6), which governs the dynamics within this sector, coincides with the quantum Ising chain in a tilted magnetic field reported in Eq. (7) (up to an irrelevant additive constant).
We note that, while in this derivation we used a non-local transformation to convert the fermions into hard-core boson, it is nevertheless possible to formulate a completely local mapping between the two Hamiltonians: in the neutral gauge sector, the Jordan-Wigner string can be completely reabsorbed using Gauss' law. To show that the mapping is local, it is sufficient to define the transformed spin operators where j 1 , j 2 label the positions of the two fermions along the chain. For τ = 0 the eigenstates are (antisymmetric combinations of) plane waves with energy E free (k 1 )+E free (k 2 ), where E free (k) = 2w cos k, k ∈ [−π, π) is the free-fermion dispersion relation. For τ = 0, however, a linear confining potential emerges between the two fermions, and the spectrum is nonperturbatively modified into a discrete tower of bound states (mesons) labelled by a quantum number = 1, 2, . . .. Each meson has a different dispersion relation E (k), where k is the center-of-mass momentum of the bound state. The meson wavefunctions and dispersion relations can be solved explicitly [27,25] by switching to the center-of-mass and relative variables, s = j 1 + j 2 , r = j 2 − j 1 > 0. Substituting the ansatz into the Schroedinger equation, we get a Wannier-Stark equation for the reduced wavefunction φ k (r), with an effective hoppingw k = 2w cos k and, crucially, the boundary condition φ k (0) ≡ 0 due to Pauli exclusion. This equation is equivalent to the recursion relation of the Bessel functions: The boundary condition J −E/τ (2w k /τ ) = 0 yields the quantization rule Note that k ∈ [−π/2, π/2), because k and k + π generate the same solution up to a phase: Since J α (−z) = e iπα J α (z), when k → k + π the wavefunction ψ gets multiplied by (−) s e iπ(r−ν ) = e −iπν (−) 2j 2 = e −iπν , i.e., a global phase. The most important qualitative aspects of this exact solution are the following. For w → 0, one finds energies E (k) = τ , corresponding to a pair of fermions separated by a string of excited gauge fields of length . In this limit, bound states are dispersionless (flat bands). As w increases, the lightest mesons progressively acquire mobility (band curvature). In particular, one can see that an effective hopping of the -th meson appears at the 2 -th order in perturbation theory in w/τ , which gives rise to a band curvature (and hence a maximal velocity) of this order of magnitude. This can be confirmed by the exact solution above, as [64] ν Interestingly, the flat-band property of heavy mesons is a nonperturbative feature that persists to arbitrarily large values of the ratio w/τ . In fact, for 4w/τ , the band curvature drops to zero faster than exponentially. This phenomenon is due to Wannier-Stark localization of particles in a linear potential: Single, isolated particles can be seen to perform a finite oscillatory motion (Bloch oscillations) with an amplitude of ξ = 2w/τ lattice sites. Correspondingly, their eigenstates are localized around each lattice site, with a localization length ξ. When two particles are initialized at a distance much larger than 2ξ, they perform independent oscillations without touching each other, and the meson is thus immobile and localized. The mobility is provided by the hard-core interaction between the two kinks, which is suppressed as the overlap between the two localized wavefunction tails, corresponding to the estimate in Eq. (C.8) [25].
There exist solutions of the Schroedinger equation (C.4) with complex momentum k and energy E, with the same wavefunction (C.7) and the same (analytically continued) energy-momentum relation (C.6). Such solutions correspond to evanescent waves and are important in the scattering problem that will be analyzed below.

Appendix D. Solution of the four-body problem
The four-body Hamiltonian for m → ∞ is obtained by projecting Eq. (A.1) onto the four-fermion subspace. It can be written in the basis of the fermion positions as where j = (j 1 , j 2 , j 3 , j 4 ) and we defined the unit vectorsê 1 = (1, 0, 0, 0),ê 2 = (0, 1, 0, 0), etc. The sum is constrained by Pauli exclusion. The diagonal term ∝ τ can be viewed as the two confining potentials for the first and last pair of adjacent fermions. These potentials give rise to two bound states (mesons). These mesons experience no residual interactions; they only interact when they bump into each other, due to Pauli exclusion. For later convenience, we define the center-ofmass positions s 1 = j 1 + j 2 , s 2 = j 3 + j 4 and relative distances r 1 = j 2 − j 1 , r 2 = j 4 − j 3 for the two mesons.
We are interested in the problem of a scattering event with incoming mesons in states ( 1 , k 1 ), ( 2 , k 2 ). This asymptotic state defines the total energy E ≡ E 1 (k 1 ) + E 2 (k 2 ) and the total momentum K ≡ k 1 + k 2 mod π of the system. Since meson-meson interactions are local, we formulate an ansatz in terms of the product state where ψ 1,2 ,k 1,2 are eigenstates of the two-body problem with quantum numbers 1,2 and generally complex momenta k 1,2 ∈ C. For Im(k 1 ) ≤ 0, Im(k 2 ) ≥ 0, the ansatz χ 1 ,k 1 , 2 ,k 2 is an asymptotic solution, and solves the Schrödinger equation in the full domain (s 2 − s 1 − r 1 − r 2 )/2 = j 3 − j 2 > 0 away from the scattering region (that is the hyperplane j 2 = j 3 ).
To obtain the complete solution for given scattering data, we first need to determine the set of parameters {( α 1 , k α 1 ), ( α 2 , k α 2 )} α=1,2,... which simultaneously satisfy the conservation laws of total energy and momentum Note that, while E and K are real, k α 1,2 and E α 1,2 (k α 1,2 ) are generally complex. Solutions with real momentum and energy correspond to incoming or outgoing states, are only a finite number. The candidate scattering solution is a linear superposition of the form The sum in Eq. (D.5) has to be restricted to the asymptotic solutions with outgoing current (see below) among those with real energy/momentum, and includes all the evanescent states that decay exponentially with the distance from the scattering region among those with complex energy/momentum. The values of the coefficients A α are then obtained by imposing the continuity of the solution in the scattering region. More explicitly, we have to impose that Ψ ≡ 0 on the hyperplane j 2 = j 3 . We define S = (s 1 + s 2 )/2, R = s 2 − s 1 , representing the global center-of-mass position and the distance between the two mesons, respectively. Hence, we have where M (r 1 ,r 2 ),α = e ipα(r 1 +r 2 ) φ α 1 ,K/2−pα (r 1 )φ α 2 ,K/2+pα (r 2 ), (D.8) The coefficients A α are then obtained by solving the linear system in Eq. (D.7), truncated to a finite set of values of α and r 1,2 ≤ r max . The truncation in α can be safely performed: the open outgoing channels are only a finite number, and the evanescent states with increasingly high quantum numbers have large Im(p α ), resulting in negligible contributions. The truncation to r 1 , r 2 ≤ r max is also legitimate: B (r 1 ,r 2 ) decays exponentially fast with r 1 and r 2 , thanks to the spatial decay of the mesonic wavefunctions φ ,k ; the coefficients M (r 1 ,r 2 ),α decay for the same reason when α represents an outgoing solution with p ∈ R, whereas the normalizability condition Im(p) > 0 guarantees the decay of the prefactor e ipα(r 1 +r 2 ) when α represents an evanescent solution. In all the calculation presented in the main text, we have checked convergence with respect to these truncation cutoffs.

Appendix E. Mesonic current
For sufficiently large m (see Appendix F), in a scattering process, the number of mesons is globally conserved. This conservation law is associated with a continuity equation. We now illustrate this continuity equation for the generic case of q mesons (i.e., in the 2q-fermion sector) in limit m → ∞. We define the density operator for the i-th meson at position x with j = (j 1 , j 2 , . . . , j 2q ), and the total mesonic density The mesonic current is defined as whereê n is the unit vector along the direction of j n . As an example, in Fig. E1-(a) we plot the time-evolving profile of the meson current J(x, t) for the scattering event discussed in Fig. (3). The red (blue) color is associated to positive (negative) current, i.e., to a meson moving to the right (left).
The mesonic density and current satisfy the continuity equation The proof of this equation is reported in the next section. We now derive the constraints imposed on the scattering solutions of the Schroedinger equation by the continuity equation.
Let us first consider the case of a single meson, with internal quantum number and momentum k. The density current associated with it can be written in terms of its dispersion relation E (k) and wavefunction φ ,k (r) as 2w cos k( |r r + 1| + H.c.) + τ r |r r| (E.6) is the reduced Hamiltonian for the internal coordinate of the meson with center-of-mass momentum k [see Eq. (C.4)]. We obtain that the mesonic current corresponds to the group velocity of the meson, in analogy with the case of a structureless quantum particle. We now apply the continuity equation to the solution Ψ(s 1 , r 1 , s 2 , r 2 ) of the stationary Schrödinger equation for the scattering problem discussed in the main text. Equation (E.4) implies that J R = J L , where J R and J L are the currents on the right and on the left, very far from the scattering region. There, the density current can be easily computed as the sum of the currents of the isolated mesons. Since evanescent waves do not contribute far from the scattering region, only propagating waves (i.e., those with q α ∈ R) should be taken into account. Therefore, in a scattering process with incoming mesons of quantum numbers 1 and 2 and momenta k 1 and k 2 , the two currents read where k α 1 = q α , k α 2 = K − q α are the momenta of the outgoing mesons. The condition J R = J L can be equivalently formulated as an equality between the incoming and outgoing currents J in = J out , defined as and The equation J in = J out has an immediate physical interpretation as a conservation of probability: in a scattering event, at t = −∞ the two mesons are with probability 1 in the state {( 1 , k 1 ), ( 2 , k 2 )}; at t = +∞, the outgoing meson states {( α 1 , k α 1 ), ( α 2 , k α 2 )} have fractional probabilities P α = J α /J in . Similarly to the scattering of structureless quantum particles, the probability of finding a certain scattering outcome (or total cross section) is proportional to the width of the wavepacket, which is determined by both the squared amplitude |A α | 2 and the group velocity.
We stress that the sign of the total current defines outgoing states, characterized In computing the amplitudes of a scattering event, one has to select the set of propagating asymptotic solutions according to this criterion, as anticipated in Appendix D above.
We finally note that the conservation law J in = J out represents a consistency check on our results for the coefficients A α obtained from the truncation of the linear system in Eq. (D.7). In Fig. E1-(b) we plot the relative violation of this conservation law, for the computations involved in Fig. 2.

Appendix E.1. Proof of the continuity equation
We prove here the continuity equation (E.4).
In the sector with q mesons, we define the operators The Hamiltonian can be written as The Heisenberg evolution of the meson density reads , (E.14) i.e., Eq. (E.4).

Appendix F. Finite fermion mass
For the sake of simplicity, the discussion above and in the main text focuses on the limit m → ∞. We compactly summarize here the effects of a finite fermion mass.
Perturbative corrections to the exact spectra and scattering solution -The main consequence of the finiteness of the fermion mass m is to produce a perturbative dressing of the vacuum and of the excitations. These effects can be explicitly computed order by order via the so-called Schrieffer-Wolff transformation [28,29,30]. In this scheme, one sequentially solves for unitary transformations U n = e im −n Sn . . . e im −1 S 1 , where the n-th generator S n is chosen to exactly cancel all processes violating fermion number conservation, in such a way that the transformed Hamiltonian H n = U n HU † n at the n-th step commutes with H 0 = j c † j c j up to terms of order m −n : H n = mH 0 + H 1 + m −1 H 2 + . . . + m −n+1 H n + O(m −n ), (F.1) with [H n , H 0 ] = 0 for all n. The approximate Hamiltonian obtained by neglecting the higher-order remainder conserves the total fermion number, and exactly accounts for all perturbative n-th order transitions within each fermion-number sector occurring through up to n virtual transitions involving intermediate states in other sectors. Upon restricting the transformed Hamiltonian to the 2q-fermion sector, one ends up with higher-order corrections to H 2q-body . For example, the first correction involves next-nearest-neighbor particle hopping terms with amplitudes w 2 /2m: Similarly, corrections of order w r /m r−1 introduce new hopping terms of range r and renormalize shorter-range terms. From the perturbatively corrected Hamiltonian, we can in principle derive the mesonic spectra, the scattering amplitudes and the mesonic currents to arbitrarily good accuracy, as long as m |w|.
Particle pair creation in high-energy collisions -In the regime considered in this work, fermionic pair creation is energetically forbidden, because the fermion mass ∼ m exceeds by far the kinetic bandwidth of excitations ∼ w. However, this phenomenon becomes relevant when m 2|w|. This can be inferred from the exact spectrum of the free fermions for τ = 0, obtained from the equivalence with the solvable transverse-field Ising chain (see Appendix B): When m approaches 2|w| (from above), the renormalized mass µ ≡ min k E(k) = m−2|w| of fermionic particles decreases to small values, and the bandwidth is ∼ 2m. Thus, if a weak string tension τ = 0 is considered, the kinetic energy of mesons can reach values much larger than their rest mass ∼ 2µ, and thus high-energy collisions could generate extra mesons. This phenomenon goes beyond the theoretical analysis presented in this work, but could be accessed with quantum simulators.
Decay of heavy mesons -A finite fermion mass may also trigger the instability of heavy mesons, which can decay into two or more lighter mesons when their gauge field string is sufficiently extended (string breaking). The lifetime of unstable mesons is (at least) exponentially long in the ratio m/|w|, as discussed in Ref. [25]; James et al. [42] have argued that it may even be infinite, based on numerical evidence. Thus, this phenomenon is not relevant in the regime studied in this work. The instability threshold is instead relevant when approaching the continuum limit m 2|w|, where the model exhibits an emergent Lorentz invariance (as can be inferred from the exact mapping in Appendix B). In this regime, the lifetimes of mesons with mass M > 4µ is only perturbative ∼ τ 3 , as computed by Rutkevich [57].