Scattering of mesons in quantum simulators

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–6]. In the last decade, substantial efforts have been devoted to the implementation of gauge-invariant Hamiltonian dynamics [7–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 [14, 15, 17, 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 figure 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.

Zoom In Zoom Out Reset image size

The paper is organized as follows. In section 2, we introduce the ${\mathbb{Z}}_{2}$-LGT analyzed throughout, and discuss particle confinement and the resulting mesonic spectra and wavefunctions in the regime of large particle mass. In section 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 section 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 gauge-neutral sector onto those of the quantum Ising chain in a tilted magnetic field. In appendices C and D we provide more details on the exact solution of the two- and 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.

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 ${\mathbb{Z}}_{2}$-LGT defined by the following Hamiltonian [25, 26]:

$$\begin{equation}H=m\sum\limits _{j}{c}_{j}^{{\dagger}}{c}_{j}\enspace +\frac{\tau }{2}\sum\limits _{j}{\sigma }_{j+1/2}^{z}+w\sum\limits _{j}\left({c}_{j}^{{\dagger}}-{c}_{j}\right)\enspace {\sigma }_{j+1/2}^{x}\enspace \left({c}_{j+1}+{c}_{j+1}^{{\dagger}}\right).\end{equation} \tag{ 1 }$$

In this equation, ${c}_{j}^{{\dagger}}$ and c_j denote creation and annihilation operators of spinless fermions of mass m > 0 on the sites $j\in \mathbb{Z}$ of a one-dimensional lattice, and ${\sigma }_{b}^{x,y,z}$ denote spin-1/2 operators (Pauli matrices) acting on the bonds $b\in \mathbb{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}={\sigma }_{j-1/2}^{z}\enspace {\sigma }_{j+1/2}^{z}\left(1-2{c}_{j}^{{\dagger}}{c}_{j}\right)$ are conserved, i.e. [G_j , H] = 0. These operators satisfy ${G}_{j}^{2}\equiv 1$ and thus generate local ${\mathbb{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 equation (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 ${\mathcal{E}}_{\ell }\left(k\right)$ and wavefunctions ψ_ℓ,k (j₁, j₂) can be determined exactly by solving the reduced two-body problem, governed by the projection of H in equation (1) onto the two-fermion sector spanned by states $\left\{\left\vert {j}_{1}{< }{j}_{2}\right\rangle \right\}$ (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₂ − j₁). The problem is solved by switching to the center-of-mass and relative variables, s = j₁ + j₂, r = j₂ − j₁ > 0, and using the ansatz ψ(s, r) = e^iks ϕ_k (r). The exact solution [25, 27] yields the quantized mesonic spectra

$$\begin{equation}{\mathcal{E}}_{\ell }\left(k\right)=\tau \enspace {\nu }_{\ell }\left(2{\tilde {w}}_{k}/\tau \right)\end{equation} \tag{ 2 }$$

with ℓ = 1, 2, ..., ${\tilde {w}}_{k}=2w\enspace \mathrm{cos}\enspace k$, k ∈ [−π/2, π/2), and ν_ℓ (x) is the ℓth (real) zero of the map $a{\mapsto}{\mathcal{J}}_{a}\left(x\right)$, where ${\mathcal{J}}_{a}\left(x\right)$ is the Bessel function ⁵ . The associated mesonic wavefunctions read

$$\begin{equation}{\psi }_{\ell ,k}\left(s,r\right)={\mathrm{e}}^{\text{i}ks}\enspace {\mathcal{J}}_{r-{\nu }_{\ell }\left(2{\tilde {w}}_{k}/\tau \right)}\left(2{\tilde {w}}_{k}/\tau \right).\end{equation} \tag{ 3 }$$

As an example, figure 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–30]. For instance, the first correction involves next-nearest-neighbor fermion hopping with amplitude w²/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.

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 non-trivial 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 ${\mathcal{E}}_{\ell }\left(k\right)$, 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 $\left\{\left({\ell }_{1}^{\prime },{k}_{1}^{\prime }\right),\left({\ell }_{2}^{\prime },{k}_{2}^{\prime }\right)\right\}$ compatible with the incoming ones by conservation of total energy and momentum:

$$\begin{equation}\begin{cases}E\equiv {\mathcal{E}}_{{\ell }_{1}}\left({k}_{1}\right)+{\mathcal{E}}_{{\ell }_{2}}\left({k}_{2}\right)={\mathcal{E}}_{{\ell }_{1}^{\prime }}\left({k}_{1}^{\prime }\right)+{\mathcal{E}}_{{\ell }_{2}^{\prime }}\left({k}_{2}^{\prime }\right),\\ K\equiv {k}_{1}+{k}_{2}={k}_{1}^{\prime }+{k}_{2}^{\prime }\quad \mathrm{mod}\enspace \pi .\end{cases}\end{equation} \tag{ 4 }$$

For all choices of incoming states, there always exist two elastic solutions, called transmitted and reflected, having $\left({\ell }_{1}^{\prime },{\ell }_{2}^{\prime }\right)=\left({\ell }_{2},{\ell }_{1}\right)$ and (ℓ₁, ℓ₂) 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}_{{\ell }_{1}^{\prime },{\ell }_{2}^{\prime }}$ of detecting that particular scattering outcome in the asymptotic future [31].

Determining the scattering amplitudes and cross sections requires solving the four-fermion Schrödinger equation. We thus consider the effective Hamiltonian H_4-body for the four-body problem, i.e. equation (1) projected to the four-fermion subspace spanned by the states $\left\{\left\vert {j}_{1}{< }{j}_{2}{< }{j}_{3}{< }{j}_{4}\right\rangle \right\}$. This consists of the hopping terms of amplitude w for the four particles, and the two diagonal confining pairwise potentials τ(j₂ − j₁) and τ(j₄ − j₃). We formulate the ansatz ${\psi }_{{\ell }_{1},q}\left({s}_{1},{r}_{1}\right){\psi }_{{\ell }_{2},K-q}\left({s}_{2},{r}_{2}\right)$, where r_1,2, s_1,2 are the relative distance and the center-of-mass position for the two mesons, the single-meson wavefunctions ψ are defined as in equation (C.7), but, crucially, the momentum $q\in \mathbb{C}$ is allowed to span the complex plane. The ansatz above represents an admissible asymptotic solution, provided $q\in \mathbb{C}$ satisfies the complex energy condition

$$\begin{equation}{\nu }_{{\ell }_{1}}\left(2{\tilde {w}}_{q}/\tau \right)+{\nu }_{{\ell }_{2}}\left(2{\tilde {w}}_{K-q}/\tau \right)\enspace =\enspace E/\tau ,\end{equation} \tag{ 5 }$$

where the total energy E and momentum K are determined by the incoming state {(ℓ₁, k₁ = q_in), (ℓ₂, k₂ = K − q_in)}, with ${q}_{\text{in}}\in \mathbb{R}$, and ν_ℓ (w) is here a complex zero of the analytic function $z{\mapsto}{\mathcal{J}}_{z}\left(w\right)$, labelled by $\ell \in \mathbb{N}$. We index by $\alpha \in \mathbb{N}$ all the triplets $\left({\ell }_{1}^{\alpha },{\ell }_{2}^{\alpha },{q}_{\alpha }\right)$ which satisfy equation (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 $\mathfrak{I}\mathfrak{m}\left(q\right){< }0$) asymptotic solutions:

$$\begin{equation}{\Psi}\left({s}_{1},{r}_{1},{s}_{2},{r}_{2}\right)={\psi }_{{\ell }_{1},{k}_{1}}\left({s}_{1},{r}_{1}\right){\psi }_{{\ell }_{2},{k}_{2}}\left({s}_{2},{r}_{2}\right)+\sum\limits _{\alpha }{A}_{\alpha }\enspace {\psi }_{{\ell }_{1}^{\alpha },{q}_{\alpha }}\left({s}_{1},{r}_{1}\right){\psi }_{{\ell }_{2}^{\alpha },K-{q}_{\alpha }}\left({s}_{2},{r}_{2}\right).\end{equation} \tag{ 6 }$$

The wavefunction Ψ solves the Schrödinger equation in the full region j₃ − j₂ = (s₂ − s₁ − r₁ − r₂)/2 > 0. ⁷ Due to Pauli exclusion at j₂ = j₃, the equation forces the boundary condition ${\Psi}{\vert }_{{s}_{2}-{s}_{1}-{r}_{1}-{r}_{2}=0}\equiv 0$, 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₁, A₂, .... The very nature of confinement, though, provides a natural truncation for this hierarchy: for ${q}_{\alpha }\in \mathbb{R}$, the meson wavefunctions are bound states, and thus fall off rapidly for large distances; for complex solutions ${q}_{\alpha }\notin \mathbb{R}$, the normalizability condition $\mathfrak{I}\mathfrak{m}\left(q\right){< }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.

In figure 2 we plot the cross sections ${P}_{{\ell }_{1}^{\prime },{\ell }_{2}^{\prime }}$ computed as described above, as a function of the incoming momenta k₁, k₂, for the scattering $\left(1,3\right)\to \left({\ell }_{1}^{\prime },{\ell }_{2}^{\prime }\right)$ when w/τ = 0.6.

Zoom In Zoom Out Reset image size

The scattering amplitudes can be readily connected with the products of real-time wavepacket collisions. We verify this by numerically computing instances of the exact time evolution within the four-body subspace. We consider the example in figure 1: in a system with L = 36 fermionic sites, we prepare two Gaussian wavepackets ${\Psi}\left({s}_{1},{s}_{2},{r}_{1},{r}_{2};t=0\right)={\mathrm{e}}^{-\left[{\left({s}_{1}-{s}_{1}^{0}\right)}^{2}+{\left({s}_{2}-{s}_{2}^{0}\right)}^{2}\right]/8{\sigma }^{2}}{\psi }_{{\ell }_{1},{k}_{1}^{0}}\left({s}_{1},{r}_{1}\right){\psi }_{{\ell }_{2},{k}_{2}^{0}}\left({s}_{2},{r}_{2}\right)$ of the meson wavefunctions in equation (C.7) with ℓ₁ = 1, ℓ₂ = 3, centered around momenta ${k}_{1}^{0}=1.3$, ${k}_{2}^{0}=-0.4$ and positions ${s}_{1}^{0}=24$, ${s}_{2}^{0}=48$, with envelopes of width $\sigma =3\sqrt{2}$ lattice sites. Time evolution from this initial state is generated by the four-body Hamiltonian H_4-body with τ = 1, w = 0.6. The final state at time t_f = 50 is examined, when the wavepackets have collided and the products of the collision have not yet reached the boundary of the system. The energy and momentum of the incoming and outgoing states are represented in figure 3(a). In figures 3(b) and (c) we plot the joint probability distribution of the momenta k₁, k₂ at times t = 0 and t = t_f, respectively, obtained via the Fourier transform of Ψ(s₁, s₂, r₁, r₂; t) with respect to the center-of-mass positions s_1,2. While the initial state shows a single density peak at $\left({k}_{1}^{0},{k}_{2}^{0}\right)$, the final state gives three different density peaks, all lying on the line ${k}_{1}+{k}_{2}={k}_{1}^{0}+{k}_{2}^{0}\enspace \mathrm{mod}\enspace \pi$, demonstrating the conservation of total momentum. The three peaks correspond to the channels predicted from the kinematic analysis (the crosses in figure 3(a)): one for the transmitted solution (ℓ₁', ℓ₂') = (3, 1) (with k₁ = −0.4, k₂ = 1.3), one for the reflected solution (1, 3) (k₁ ≃ −1.0, k₂ ≃ −1.2), and one for the inelastic solution (2, 2) (k₁ ≃ −1.5, k₂ ≃ −0.7). The peaks can be better resolved by plotting the distribution of the relative momentum k₁ − k₂, as done in figure 3(e). The relative weights enclosed within the dashed contours in figure 3(c), ${P}_{13}^{\text{num}}\simeq 0.49$, ${P}_{22}^{\text{num}}\simeq 0.1$, ${P}_{31}^{\text{num}}\simeq 0.41$, are compatible with the predicted cross sections P₁₃ ≃ 0.47, P₂₂ ≃ 0.12, P₃₁ ≃ 0.41. ⁸

Zoom In Zoom Out Reset image size

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 equation (1) by exploiting the equivalence with the quantum Ising chain in a tilted magnetic field [25, 32, 33]:

$$\begin{equation}{H}_{\text{qIc}}=\sum\limits _{j}-\frac{m}{2}\enspace {\sigma }_{j-1/2}^{z}{\sigma }_{j+1/2}^{z}+\frac{\tau }{2}\enspace {\sigma }_{j+1/2}^{z}+w\enspace {\sigma }_{j+1/2}^{x}.\end{equation} \tag{ 7 }$$

To obtain this non-trivial identification, one has to exploit the gauge-symmetry constraints to eliminate the fermionic degrees of freedom [17, 25, 26, 34], as reported in appendix B. Within this equivalence, ${\mathbb{Z}}_{2}$-mesons may be pictured as domains of reversed magnetization in a polarized background, arising from domain-wall confinement [35–39]. Recently, confinement in the quantum Ising chain has been shown to underlie a plethora of intriguing non-equilibrium phenomena [25, 40–43]. The mapping above is extremely advantageous for quantum simulations because it implements gauge invariance exactly, similarly to what done in references [15, 17] for the Schwinger model. The dynamics governed by equation (7) can be experimentally realized both with optical lattices [44, 45] and Rydberg atoms trapped in optical tweezers [16, 46].

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 equation (7) to filter meson wavepackets with sharply-defined momenta, at the price of moderately longer chains and a limited amount of post-selection. The numerical simulation in figure 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.

Zoom In Zoom Out Reset image size

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 figure 1), so they can be resolved as spatially separate wavepackets. For implementations based on equation (7), the particle density ${c}_{j}^{{\dagger}}{c}_{j}$ in equation (1) maps to the domain-wall density $\left(1-{\sigma }_{j-1/2}^{z}{\sigma }_{j+1/2}^{z}\right)/2$: thus, it is sufficient to measure the magnetization profile $\langle {\sigma }_{j+1/2}^{z}\left({t}_{\mathrm{f}}\right)\rangle$ in the final state [16, 44–46] 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: reference [16], for example, demonstrated state preparation and single-qubit readout in a chain of 51 ⁸⁷Rb atoms governed by Ising-type dynamics close to equation (7), with excellent coherence control over several tens of time units (2π/w).

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 LGTs (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 [25, 41, 42, 49–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–60]. This would represent a first step towards the ultimate goal of simulating realistic scattering problems in QCD such as heavy-ion collisions [61].

We acknowledge useful discussions and feedback by D A Abanin, A Bastianello, P Calabrese, M Dalmonte, W De Roeck, G Giudici, G Pagano, J Sonner, T Wang, and our co-authors of reference [25]. FMS is partly supported by the ERC under Grant No. 758329 (AGEnTh). AL is supported by the Swiss National Science Foundation.

The data that support the findings of this study are available upon reasonable request from the authors.

While completing the present manuscript, we became aware of a related work [62], appeared simultaneously.

In this work we have focused on the (1 + 1)-dimensional ${\mathbb{Z}}_{2}$-LGT defined by the Hamiltonian in equation (1) of the main text, reported here for convenience:

$$\begin{equation}H=m\sum\limits _{j}{c}_{j}^{{\dagger}}{c}_{j}\enspace +\frac{\tau }{2}\sum\limits _{j}{\sigma }_{j+1/2}^{z}+w\sum\limits _{j}\left({c}_{j}^{{\dagger}}-{c}_{j}\right)\enspace {\sigma }_{j+1/2}^{x}\enspace \left({c}_{j+1}+{c}_{j+1}^{{\dagger}}\right).\end{equation} \tag{ A.1 }$$

Equation (A.1) describes gauge-invariant interactions of fermionic particles mediated by a ${\mathbb{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 ${\sigma }_{j+1/2}^{z}$ on that bond, as illustrated in figure B1(a). More formally, interactions are such that all the local operators ${G}_{j}={\sigma }_{j-1/2}^{z}\enspace {\sigma }_{j+1/2}^{z}\left(1-2{c}_{j}^{{\dagger}}{c}_{j}\right)$ are conserved, i.e. [G_j , H] = 0. These operators satisfy ${G}_{j}^{2}\equiv 1$ and thus generate local ${\mathbb{Z}}_{2}$ symmetries. Accordingly, the complete Hilbert space decomposes into dynamically disconnected subspaces, labelled by the set of eigenvalues $\left\{{\text{e}}^{\text{i}\pi {q}_{j}}={\pm}1\right\}$ of {G_j }, where q_j = 0 or 1 is interpreted as the absence or presence of a static background ${\mathbb{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

$$\begin{equation}{\sigma }_{j-1/2}^{z}\enspace {\sigma }_{j+1/2}^{z}=1-2{c}_{j}^{{\dagger}}{c}_{j}\end{equation} \tag{ A.2 }$$

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 figure B1(a).

The LGT in equation (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 $\left\vert \mathrm{G}\mathrm{S}\right\rangle$ becomes a fermion vacuum with ${c}_{j}^{{\dagger}}{c}_{j}\left\vert \mathrm{G}\mathrm{S}\right\rangle =0$; by gauge-neutrality, (cf equation (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 ${\mathbb{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 ${\mathbb{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.

Here we illustrate the mapping between the ${\mathbb{Z}}_{2}$-LGT in equation (A.1) and the quantum Ising chain in a tilted magnetic field, equation (7). This mapping was proposed in reference [25], and is connected with the one discussed in reference [26].

We consider the neutral sector, defined by the Gauss law (A.2), i.e. ${G}_{j}={\sigma }_{j-1/2}^{z}\enspace {\sigma }_{j+1/2}^{z}\left(1-2{n}_{j}\right)\equiv 1$, with ${n}_{j}={c}_{j}^{{\dagger}}{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 [34], as highlighted in figure B1(b). The latter are redundant due to the infinitely many local constraints in equation (A.2): a classical configuration of the gauge field uniquely fixes the configuration of the fermionic matter via the Gauss law.

Zoom In Zoom Out Reset image size

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 ${\tau }_{j}^{\alpha }$, for α = +, −, z, defined as

$$\begin{equation}{\tau }_{j}^{-}=\prod\limits _{k{< }j}\left(1-2{n}_{k}\right){c}_{j},\hspace{25.0pt}{\tau }_{j}^{+}={\left({\tau }_{j}^{-}\right)}^{{\dagger}},\hspace{25.0pt}{\tau }_{j}^{z}=2{n}_{j}-1.\end{equation} \tag{ B.1 }$$

As can be easily checked, the operators ${\tau }_{j}^{z}$, ${\tau }_{j}^{x}={\tau }_{j}^{+}+{\tau }_{j}^{-}$ and ${\tau }_{j}^{y}=-i\left({\tau }_{j}^{+}-{\tau }_{j}^{-}\right)$ satisfy the usual commutation relations of spin operators. By applying this transformation, we get the Hamiltonian

$$\begin{equation}H=m\sum\limits _{j}\frac{1+{\tau }_{j}^{z}}{2}\enspace +\frac{\tau }{2}\sum\limits _{j}{\sigma }_{j+1/2}^{z}+w\sum\limits _{j}{\tau }_{j}^{x}\enspace {\sigma }_{j+1/2}^{x}\enspace {\tau }_{j+1}^{x},\end{equation} \tag{ B.2 }$$

with the constraint

$$\begin{equation}{G}_{j}=-{\sigma }_{j-1/2}^{z}\enspace {\sigma }_{j+1/2}^{z}{\tau }_{j}^{z}=1.\end{equation} \tag{ B.3 }$$

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 ' = UG_j U^† ≡ 1 from equation (B.3) only depends on the matter degrees of freedom, whereas the transformed Hamiltonian H' = UHU^† only involves the gauge degrees of freedom. This can be accomplished with

$$\begin{align}U=\prod\limits _{j}\mathrm{exp}\left[\frac{\mathrm{i}\pi }{2}\left({\tau }_{j}^{x}-1\right)\frac{1-{\sigma }_{j-1/2}^{z}{\sigma }_{j+1/2}^{z}}{2}\right]=\prod\limits _{j}\left[\frac{1+{\sigma }_{j-1/2}^{z}{\sigma }_{j+1/2}^{z}}{2}+{\tau }_{j}^{x}\frac{1-{\sigma }_{j-1/2}^{z}{\sigma }_{j+1/2}^{z}}{2}\right].\end{align} \tag{ B.4 }$$

This transformation flips the spin ${\tau }_{j}^{z}$ where the neighbouring gauge fields are anti-aligned and does nothing where they are aligned. We get $U{\tau }_{j}^{z}{U}^{{\dagger}}={\tau }_{j}^{z}{\sigma }_{j+1/2}^{z}{\sigma }_{j+1/2}^{z}$ and $U{\tau }_{j}^{x}{\sigma }_{j+1/2}^{x}{\tau }_{j+1}^{x}{U}^{{\dagger}}={\sigma }_{j+1/2}^{x}$. The transformed constraint

$$\begin{equation}{G}_{j}^{\prime }=U{G}_{j}{U}^{{\dagger}}=-{\tau }_{j}^{z}\equiv 1\end{equation} \tag{ B.5 }$$

forces the ${\tau }_{j}^{\alpha }$ spins to be polarized in the $-\hat{z}$ directions. They enter the transformed Hamiltonian only via G_j ':

$$\begin{equation}{H}^{\prime }=m\sum\limits _{j}\frac{1-{G}_{j}^{\prime }{\sigma }_{j-1/2}^{z}{\sigma }_{j+1/2}^{z}}{2}\enspace +\frac{\tau }{2}\sum\limits _{j}{\sigma }_{j+1/2}^{z}+w\sum\limits _{j}{\sigma }_{j+1/2}^{x}.\end{equation} \tag{ B.6 }$$

Thus, in the neutral gauge sector, by equation (B.5) the spins ${\tau }_{j}^{\alpha }$ 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 equation (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

$$\begin{equation}{\tilde {\sigma }}_{j+1/2}^{x}=\left({c}_{j}^{{\dagger}}-{c}_{j}\right){\sigma }_{j+1/2}^{x}\left({c}_{j+1}^{{\dagger}}+{c}_{j+1}\right),\end{equation} \tag{ B.7 }$$

$$\begin{equation}{\tilde {\sigma }}_{j+1/2}^{y}=\left({c}_{j}^{{\dagger}}-{c}_{j}\right){\sigma }_{j+1/2}^{y}\left({c}_{j+1}^{{\dagger}}+{c}_{j+1}\right),\end{equation} \tag{ B.8 }$$

$$\begin{equation}{\tilde {\sigma }}_{j+1/2}^{z}={\sigma }_{j+1/2}^{x}.\end{equation} \tag{ B.9 }$$

These operators satisfy the usual commutation relations of Pauli matrices and are related to the original spins by a local transformation. It is then immediate to write equation (A.1) in terms of the new spin operators, obtaining the quantum Ising chain in the neutral gauge sector:

$$\begin{equation}{H}^{\prime }=m\sum\limits _{j}\frac{1-{\tilde {\sigma }}_{j-1/2}^{z}{\tilde {\sigma }}_{j+1/2}^{z}{G}_{j}}{2}\enspace +\frac{\tau }{2}\sum\limits _{j}{\tilde {\sigma }}_{j+1/2}^{z}+w\sum\limits _{j}{\tilde {\sigma }}_{j+1/2}^{x}.\end{equation} \tag{ B.10 }$$

The two-body Hamiltonian for m → ∞ is obtained by projecting equation (A.1) onto the two-fermion subspace. It can be written in the basis of the fermion positions as

$$\begin{align}{H}_{2-\text{body}}=\sum\limits _{{j}_{1}{< }{j}_{2}}\tau \left({j}_{2}-{j}_{1}\right)\enspace \left\vert {j}_{1},{j}_{2}\right\rangle \left\langle {j}_{1},{j}_{2}\right\vert +w\left(\left\vert {j}_{1}+1,{j}_{2}\right\rangle \left\langle {j}_{1},{j}_{2}\right\vert +\left\vert {j}_{1},{j}_{2}+1\right\rangle \left\langle {j}_{1},{j}_{2}\right\vert +\text{h.c.}\right),\end{align} \tag{ C.1 }$$

where j₁, j₂ label the positions of the two fermions along the chain.

For τ = 0 the eigenstates are (antisymmetric combinations of) plane waves

$$\begin{equation}{{\Psi}}_{{k}_{1},{k}_{2}}={\text{e}}^{\text{i}{k}_{1}{j}_{1}+\mathrm{i}{k}_{2}{j}_{2}}-{\text{e}}^{\text{i}{k}_{2}{j}_{1}+\mathrm{i}{k}_{1}{j}_{2}}\end{equation} \tag{ C.2 }$$

with energy ${\mathcal{E}}_{\text{free}}\left({k}_{1}\right)+{\mathcal{E}}_{\text{free}}\left({k}_{2}\right)$, where ${\mathcal{E}}_{\text{free}}\left(k\right)=2w\enspace \mathrm{cos}\enspace 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 ${\mathcal{E}}_{\ell }\left(k\right)$, where k is the center-of-mass momentum of the bound state.

The meson wavefunctions and dispersion relations can be solved explicitly [25, 27] by switching to the center-of-mass and relative variables, s = j₁ + j₂, r = j₂ − j₁ > 0. Substituting the ansatz

$$\begin{equation}{\psi }_{k}\left(s,r\right)={\mathrm{e}}^{\text{i}ks}{\phi }_{k}\left(r\right)\end{equation} \tag{ C.3 }$$

into the Schrödinger equation, we get a Wannier–Stark equation for the reduced wavefunction ϕ_k (r),

$$\begin{equation}{\tilde {w}}_{k}\left[{\phi }_{k}\left(r+1\right)+{\phi }_{k}\left(r-1\right)\right]+\tau r\enspace {\phi }_{k}\left(r\right)=\mathcal{E}\enspace {\phi }_{k}\left(r\right)\end{equation} \tag{ C.4 }$$

with an effective hopping ${\tilde {w}}_{k}=2w\enspace \mathrm{cos}\enspace 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:

$$\begin{equation}{\phi }_{k}\left(r\right)={\mathcal{J}}_{r-\mathcal{E}/\tau }\left(2{\tilde {w}}_{k}/\tau \right).\end{equation} \tag{ C.5 }$$

The boundary condition ${\mathcal{J}}_{-\mathcal{E}/\tau }\left(2{\tilde {w}}_{k}/\tau \right)=0$ yields the quantization rule

$$\begin{equation}\mathcal{E}={\mathcal{E}}_{\ell }\left(k\right)=\tau \enspace {\nu }_{\ell }\left(2{\tilde {w}}_{k}/\tau \right)\equiv -\tau {\times}\left\{\ell \mathrm{t}\mathrm{h}\;\;\mathrm{z}\mathrm{e}\mathrm{r}\mathrm{o}\;\;\mathrm{o}\mathrm{f}\;\;x{\mapsto}{\mathcal{J}}_{x}\left(2{\tilde {w}}_{k}/{\tau}\right)\right\}\end{equation} \tag{ C.6 }$$

for ℓ = 1, 2, ..., which defines the exact dispersion relations of all mesons. The meson wavefunctions are thus

$$\begin{equation}{\psi }_{\ell ,k}\left(s,r\right)={\mathrm{e}}^{\text{i}ks}\enspace {\mathcal{J}}_{r-{\nu }_{\ell }}\left(2{\tilde {w}}_{k}/\tau \right).\end{equation} \tag{ C.7 }$$

Note that k ∈ [−π/2, π/2), because k and k + π generate the same solution up to a phase: since ${\mathcal{J}}_{\alpha }\left(-z\right)={\text{e}}^{\text{i}\pi \alpha }\enspace {\mathcal{J}}_{\alpha }\left(z\right)$, when k ↦ k + π the wavefunction ψ gets multiplied by ${\left(-\right)}^{s}{\text{e}}^{\text{i}\pi \left(r-{\nu }_{\ell }\right)}={\text{e}}^{-\text{i}\pi {\nu }_{\ell }}{\left(-\right)}^{2{j}_{2}}={\text{e}}^{-\text{i}\pi {\nu }_{\ell }}$, i.e. a global phase.

The most important qualitative aspects of this exact solution are the following. For w → 0, one finds energies ${\mathcal{E}}_{\ell }\left(k\right)=\tau \ell$, 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 [63]

$$\begin{equation}{\nu }_{\ell }\left(x\right)\enspace \underset{x\to 0}{\sim }\enspace \ell \left(1-\frac{{x}^{2\ell }}{{\left(\ell !\right)}^{2}}\right).\end{equation} \tag{ C.8 }$$

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 equation (C.8) [25].

There exist solutions of the Schrödinger equation (C.4) with complex momentum k and energy $\mathcal{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.

The four-body Hamiltonian for m → ∞ is obtained by projecting equation (A.1) onto the four-fermion subspace. It can be written in the basis of the fermion positions as

$$\begin{equation}{H}_{4-\text{body}}=\sum\limits _{{j}_{1}{< }{j}_{2}{< }{j}_{3}{< }{j}_{4}}\left[\tau \left({j}_{2}-{j}_{1}+{j}_{4}-{j}_{3}\right)\enspace \left\vert \overrightarrow {j}\right\rangle \left\langle \overrightarrow {j}\right\vert +w\sum\limits _{n=1}^{4}\left(\left\vert \overrightarrow {j}+{\hat{e}}_{n}\right\rangle \left\langle \overrightarrow {j}\right\vert +\text{h.c.}\right)\right]\end{equation} \tag{ D.1 }$$

where $\overrightarrow {j}=\left({j}_{1},{j}_{2},{j}_{3},{j}_{4}\right)$ and we defined the unit vectors ${\hat{e}}_{1}=\left(1,0,0,0\right)$, ${\hat{e}}_{2}=\left(0,1,0,0\right)$, 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-of-mass positions s₁ = j₁ + j₂, s₂ = j₃ + j₄ and relative distances r₁ = j₂ − j₁, r₂ = j₄ − j₃ for the two mesons.

We are interested in the problem of a scattering event with incoming mesons in states (ℓ₁, k₁), (ℓ₂, k₂). This asymptotic state defines the total energy $E\equiv {\mathcal{E}}_{{\ell }_{1}}\left({k}_{1}\right)+{\mathcal{E}}_{{\ell }_{2}}\left({k}_{2}\right)$ and the total momentum K ≡ k₁ + k₂ mod π of the system. Since meson–meson interactions are local, we formulate an ansatz in terms of the product state

$$\begin{equation}{\chi }_{{\ell }_{1}^{\prime },{k}_{1}^{\prime },{\ell }_{2}^{\prime },{k}_{2}^{\prime }}\left({s}_{1},{r}_{1},{s}_{2},{r}_{2}\right)={\psi }_{{\ell }_{1}^{\prime },{k}_{1}^{\prime }}\left({s}_{1},{r}_{1}\right){\psi }_{{\ell }_{2}^{\prime },{k}_{2}^{\prime }}\left({s}_{2},{r}_{2}\right),\end{equation} \tag{ D.2 }$$

where ${\psi }_{{\ell }_{1,2}^{\prime },{k}_{1,2}^{\prime }}$ are eigenstates of the two-body problem with quantum numbers ℓ_1,2' and generally complex momenta ${k}_{1,2}^{\prime }\in \mathbb{C}$. For $\mathfrak{I}\mathfrak{m}\left({k}_{1}^{\prime }\right){\leqslant}0$, $\mathfrak{I}\mathfrak{m}\left({k}_{2}^{\prime }\right){\geqslant}0$, the ansatz ${\chi }_{{\ell }_{1}^{\prime },{k}_{1}^{\prime },{\ell }_{2}^{\prime },{k}_{2}^{\prime }}$ is an asymptotic solution, and solves the Schrödinger equation in the full domain (s₂ − s₁ − r₁ − r₂)/2 = j₃ − j₂ > 0 away from the scattering region (that is the hyperplane j₂ = j₃).

To obtain the complete solution for given scattering data, we first need to determine the set of parameters ${\left\{\left({\ell }_{1}^{\alpha },{k}_{1}^{\alpha }\right),\left({\ell }_{2}^{\alpha },{k}_{2}^{\alpha }\right)\right\}}_{\alpha =1,2,\dots }$ which simultaneously satisfy the conservation laws of total energy and momentum

$$\begin{equation}E={\mathcal{E}}_{{\ell }_{1}^{\alpha }}\left({k}_{1}^{\alpha }\right)+{\mathcal{E}}_{{\ell }_{2}^{\alpha }}\left({k}_{2}^{\alpha }\right),\end{equation} \tag{ D.3 }$$

$$\begin{equation}K={k}_{1}^{\alpha }+{k}_{2}^{\alpha }\enspace \mathrm{mod}\enspace \pi .\end{equation} \tag{ D.4 }$$

Note that, while E and K are real, ${k}_{1,2}^{\alpha }$ and ${\mathcal{E}}_{{\ell }_{1,2}^{\alpha }}\left({k}_{1,2}^{\alpha }\right)$ 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

$$\begin{equation}{\Psi}\left({s}_{1},{r}_{1},{s}_{2},{r}_{2}\right)={\chi }_{{\ell }_{1},{k}_{1},{\ell }_{2},{k}_{2}}+\sum\limits _{\alpha }{A}_{\alpha }{\chi }_{{\ell }_{1}^{\alpha },{k}_{1}^{\alpha },{\ell }_{2}^{\alpha },{k}_{2}^{\alpha }}.\end{equation} \tag{ D.5 }$$

The sum in equation (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₂ = j₃. We define S = (s₁ + s₂)/2, R = s₂ − s₁, representing the global center-of-mass position and the distance between the two mesons, respectively. Hence, we have

$$\begin{align}{\Psi}\left(S,R,{r}_{1},{r}_{2}\right)={\mathrm{e}}^{\text{i}KS}\left[{\mathrm{e}}^{\text{i}pR}\enspace {\phi }_{{\ell }_{1},K/2-p}\left({r}_{1}\right){\phi }_{{\ell }_{2},K/2+p}\left({r}_{2}\right)+\sum\limits _{\alpha }{A}_{\alpha }\enspace {\mathrm{e}}^{\mathrm{i}{p}_{\alpha }R}\enspace {\phi }_{{\ell }_{1}^{\alpha },K/2-{p}_{\alpha }}\left({r}_{1}\right){\phi }_{{\ell }_{2}^{\alpha },K/2+{p}_{\alpha }}\left({r}_{2}\right)\right]\end{align} \tag{ D.6 }$$

where p = (k₂ − k₁)/2, ${p}_{\alpha }=\left({k}_{2}^{\alpha }-{k}_{1}^{\alpha }\right)/2$, and ϕ_ℓ,k are the solutions to equation (C.4). The boundary condition applied to equation (D.6) reads

$$\begin{equation}\sum\limits _{\alpha }{M}_{\left({r}_{1},{r}_{2}\right),\alpha }{A}_{\alpha }={B}_{\left({r}_{1},{r}_{2}\right)},\end{equation} \tag{ D.7 }$$

where

$$\begin{equation}{M}_{\left({r}_{1},{r}_{2}\right),\alpha }={\mathrm{e}}^{\mathrm{i}{p}_{\alpha }\left({r}_{1}+{r}_{2}\right)}\enspace {\phi }_{{\ell }_{1}^{\alpha },K/2-{p}_{\alpha }}\left({r}_{1}\right){\phi }_{{\ell }_{2}^{\alpha },K/2+{p}_{\alpha }}\left({r}_{2}\right),\end{equation} \tag{ D.8 }$$

$$\begin{equation}{B}_{\left({r}_{1},{r}_{2}\right)}=-{\mathrm{e}}^{\mathrm{i}p\left({r}_{1}+{r}_{2}\right)}\enspace {\phi }_{{\ell }_{1},K/2-p}\left({r}_{1}\right){\phi }_{{\ell }_{2},K/2+p}\left({r}_{2}\right).\end{equation} \tag{ D.9 }$$

The coefficients A_α are then obtained by solving the linear system in equation (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 $\mathfrak{I}\mathfrak{m}\left({p}_{\alpha }\right)$, resulting in negligible contributions. The truncation to r₁, r₂ ⩽ r_max is also legitimate: ${B}_{\left({r}_{1},{r}_{2}\right)}$ decays exponentially fast with r₁ and r₂, thanks to the spatial decay of the mesonic wavefunctions ϕ_ℓ,k ; the coefficients ${M}_{\left({r}_{1},{r}_{2}\right),\alpha }$ decay for the same reason when α represents an outgoing solution with $p\in \mathbb{R}$, whereas the normalizability condition $\mathfrak{I}\mathfrak{m}\left(p\right){ >}0$ guarantees the decay of the prefactor ${\mathrm{e}}^{\mathrm{i}{p}_{\alpha }\left({r}_{1}+{r}_{2}\right)}$ when α represents an evanescent solution. In all the calculation presented in the main text, we have checked convergence with respect to these truncation cutoffs.

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 ith meson at position x

$$\begin{equation}{\chi }_{i}\left(x\right)=\sum\limits _{{j}_{1}{< }{j}_{2}{< }\cdots {< }{j}_{2q}}{\delta }_{{j}_{2i-1}+{j}_{2i},x}\left\vert \overrightarrow {j}\right\rangle \left\langle \overrightarrow {j}\right\vert \end{equation} \tag{ E.1 }$$

with $\overrightarrow {j}=\left({j}_{1},{j}_{2},\dots ,{j}_{2q}\right)$, and the total mesonic density

$$\begin{equation}\rho \left(x\right)=\sum\limits _{i=1}^{q}{\rho }_{i}\left(x\right)\equiv \sum\limits _{i=1}^{q}\left\langle {\Psi}\right\vert {\chi }_{i}\left(x\right)\left\vert {\Psi}\right\rangle .\end{equation} \tag{ E.2 }$$

The mesonic current is defined as $J\left(x\right)={\sum }_{i=1}^{q}{J}_{i}\left(x\right)$, with

$$\begin{equation}{J}_{i}\left(x\right)=-2w\sum\limits _{{j}_{1}{< }{j}_{2}{< }\cdots {< }{j}_{2q}}{\delta }_{{j}_{2i}+{j}_{2i+1},x}\enspace \mathfrak{I}\mathfrak{m}\left\{{{\Psi}}^{{\ast}}\left( \overrightarrow {j}\right)\left[{\Psi}\left( \overrightarrow {j}+{\hat{e}}_{2i}\right)+{\Psi}\left( \overrightarrow {j}+{\hat{e}}_{2i+1}\right)\right]\right\},\end{equation} \tag{ E.3 }$$

where ${\hat{e}}_{n}$ is the unit vector along the direction of j_n . As an example, in figure E1(a) we plot the time-evolving profile of the meson current J(x, t) for the scattering event discussed in figure (3). The red (blue) color is associated to positive (negative) current, i.e. to a meson moving to the right (left).

Zoom In Zoom Out Reset image size

The mesonic density and current satisfy the continuity equation

$$\begin{equation}{\mathrm{\partial }}_{t}\rho \left(x\right)=-J\left(x\right)+J\left(x-1\right).\end{equation} \tag{ E.4 }$$

The proof of this equation is reported in the next section. We now derive the constraints imposed on the scattering solutions of the Schrödinger 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 ${\mathcal{E}}_{\ell }\left(k\right)$ and wavefunction ϕ_ℓ,k (r) as

$$\begin{align}J\left(x\right)& =-2w\enspace \mathfrak{I}\mathfrak{m}\left\{\sum\limits _{r{ >}0}{\phi }_{\ell ,k}^{{\ast}}\left(r\right){\mathrm{e}}^{\mathrm{i}k}\left[{\phi }_{\ell ,k}\left(r+1\right)+{\phi }_{\ell ,k}\left(r-1\right)\right]\right\}\\ & =\sum\limits _{r{ >}0}{\phi }_{\ell ,k}^{{\ast}}\left(r\right)\left(-2w\enspace \mathrm{sin}\enspace k\right)\left[{\phi }_{\ell ,k}\left(r+1\right)+{\phi }_{\ell ,k}\left(r-1\right)\right]\\ & =\left\langle {\phi }_{\ell ,k}\right\vert {\mathrm{\partial }}_{k}{H}_{k}\left\vert {\phi }_{\ell ,k}\right\rangle ={\mathrm{\partial }}_{k}{\mathcal{E}}_{\ell }\left(k\right)={v}_{\ell }\left(k\right)\end{align} \tag{ E.5 }$$

where

$$\begin{equation}{H}_{k}=\sum\limits _{r{ >}0}2w\enspace \mathrm{cos}\enspace k\left(\left\vert r\right\rangle \left\langle r+1\right\vert +\text{h.c.}\right)+\tau r\left\vert r\right\rangle \left\langle r\right\vert \end{equation} \tag{ E.6 }$$

is the reduced Hamiltonian for the internal coordinate of the meson with center-of-mass momentum k (see equation (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₁, r₁, s₂, r₂) 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}_{\alpha }\in \mathbb{R}$) should be taken into account. Therefore, in a scattering process with incoming mesons of quantum numbers ℓ₁ and ℓ₂ and momenta k₁ and k₂, the two currents read

$$\begin{equation}{J}_{\mathrm{L}}=v\left({\ell }_{1},{k}_{1}\right)+\sum\limits _{\alpha \vert {q}_{\alpha }\in \mathbb{R}}\vert {A}_{\alpha }{\vert }^{2}v\left({\ell }_{1}^{\alpha },{k}_{1}^{\alpha }\right)\end{equation} \tag{ E.7 }$$

$$\begin{equation}{J}_{\mathrm{R}}=v\left({\ell }_{2},{k}_{2}\right)+\sum\limits _{\alpha \vert {q}_{\alpha }\in \mathbb{R}}\vert {A}_{\alpha }{\vert }^{2}v\left({\ell }_{2}^{\alpha },{k}_{2}^{\alpha }\right)\end{equation} \tag{ E.8 }$$

where ${k}_{1}^{\alpha }={q}_{\alpha }$, ${k}_{2}^{\alpha }=K-{q}_{\alpha }$ 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

$$\begin{equation}{J}_{\mathrm{i}\mathrm{n}}=v\left({\ell }_{1},{k}_{1}\right)-v\left({\ell }_{2},{k}_{2}\right)\end{equation} \tag{ E.9 }$$

and

$$\begin{equation}{J}_{\text{out}}=\sum\limits _{\alpha }{J}_{\alpha }=\sum\limits _{\alpha }\vert {A}_{\alpha }{\vert }^{2}\left[-v\left({\ell }_{1}^{\alpha },{k}_{1}^{\alpha }\right)+v\left({\ell }_{2}^{\alpha },{k}_{2}^{\alpha }\right)\right].\end{equation} \tag{ E.10 }$$

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 {(ℓ₁, k₁), (ℓ₂, k₂)}; at t = +∞, the outgoing meson states $\left\{\left({\ell }_{1}^{\alpha },{k}_{1}^{\alpha }\right),\left({\ell }_{2}^{\alpha },{k}_{2}^{\alpha }\right)\right\}$ 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_α |² and the group velocity.

We stress that the sign of the total current defines outgoing states, characterized by $-v\left({\ell }_{1}^{\alpha },{k}_{1}^{\alpha }\right)+v\left({\ell }_{2}^{\alpha },{k}_{2}^{\alpha }\right){ >}0$. 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 equation (D.7). In figure E1(b) we plot the relative violation of this conservation law, for the computations involved in figure 2.

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 ${{\Delta}}^{+}={\sum }_{i=1}^{q}{{\Delta}}_{i}^{+}$ where

$$\begin{equation}{{\Delta}}_{i}^{+}=\sum\limits _{{j}_{1}{< }{j}_{2}{< }\cdots {< }{j}_{2q}}\sum\limits _{s=2i,2i+1}w\left(1-{\delta }_{{j}_{s+1},{j}_{s}+1}\right)\left\vert \overrightarrow {j}+{\hat{e}}_{s}\right\rangle \left\langle \overrightarrow {j}\right\vert .\end{equation} \tag{ E.11 }$$

The Hamiltonian can be written as H = Δ⁺ + Δ⁻ + V, where ${{\Delta}}^{-}={\left({{\Delta}}^{+}\right)}^{{\dagger}}$, V = ∑_i V_i and

$$\begin{equation}{V}_{i}=\sum\limits _{{j}_{1}{< }{j}_{2}{< }\cdots {< }{j}_{2q}}\tau \left({j}_{2i+1}-{j}_{2i}\right)\left\vert \overrightarrow {j}\right\rangle \left\langle \overrightarrow {j}\right\vert .\end{equation} \tag{ E.12 }$$

The Heisenberg evolution of the meson density reads

$$\begin{align}{\mathrm{\partial }}_{t}\rho \left(x\right)& =i\sum\limits _{i=1}^{q}\left\langle {\Psi}\right\vert \left[H,{\chi }_{i}\left(x\right)\right]\left\vert {\Psi}\right\rangle \\ & =i\sum\limits _{i=1}^{q}\left\langle {\Psi}\right\vert \left[{{\Delta}}_{i}^{+}+{{\Delta}}_{i}^{-},{\chi }_{i}\left(x\right)\right]\left\vert {\Psi}\right\rangle \\ & =2\sum\limits _{i=1}^{q}\mathfrak{I}\mathfrak{m}\left\langle {\Psi}\right\vert {\chi }_{i}\left(x\right)\left({{\Delta}}_{i}^{+}+{{\Delta}}_{i}^{-}\right)\left\vert {\Psi}\right\rangle .\end{align} \tag{ E.13 }$$

By using the properties ${{\Delta}}_{i}^{-}{\chi }_{i}\left(x\right)={\chi }_{i}\left(x-1\right){{\Delta}}_{i}^{-}$ and $\left\langle {\Psi}\right\vert {{\Delta}}_{i}^{-}{\chi }_{i}\left(x\right)\left\vert {\Psi}\right\rangle ={\left(\left\langle {\Psi}\right\vert {\chi }_{i}\left(x\right){{\Delta}}_{i}^{+}\left\vert {\Psi}\right\rangle \right)}^{{\ast}}$ we get

$$\begin{align}{\mathrm{\partial }}_{t}\rho \left(x\right)& =2\sum\limits _{i=1}^{m}\mathfrak{I}\mathfrak{m}\left\langle {\Psi}\right\vert {\chi }_{i}\left(x\right){{\Delta}}_{i}^{-}+{\chi }_{i}\left(x\right){{\Delta}}_{i}^{+}\left\vert {\Psi}\right\rangle \\ & =2\sum\limits _{i=1}^{m}\mathfrak{I}\mathfrak{m}\left\langle {\Psi}\right\vert {\chi }_{i}\left(x\right){{\Delta}}_{i}^{-}-{{\Delta}}_{i}^{-}{\chi }_{i}\left(x\right)\left\vert {\Psi}\right\rangle \\ & =2\sum\limits _{i=1}^{m}\mathfrak{I}\mathfrak{m}\left\langle {\Psi}\right\vert \left({\chi }_{i}\left(x\right)-{\chi }_{i}\left(x-1\right)\right){{\Delta}}_{i}^{-}\left\vert {\Psi}\right\rangle \\ & =-J\left(x\right)+J\left(x-1\right),\end{align} \tag{ E.14 }$$

i.e. equation (E.4).

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.

F.1. 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–30]. In this scheme, one sequentially solves for unitary transformations ${U}_{n}={\text{e}}^{\text{i}{m}^{-n}{S}_{n}}\dots {\text{e}}^{\text{i}{m}^{-1}{S}_{1}}$, where the nth generator S_n is chosen to exactly cancel all processes violating fermion number conservation, in such a way that the transformed Hamiltonian ${H}_{n}^{\prime }={U}_{n}H{U}_{n}^{{\dagger}}$ at the nth step commutes with ${H}_{0}={\sum }_{j}{c}_{j}^{{\dagger}}{c}_{j}$ up to terms of order m⁻ⁿ :

$$\begin{equation}{H}_{n}^{\prime }=m{H}_{0}+{H}_{1}+{m}^{-1}{H}_{2}+\cdots +{m}^{-n+1}{H}_{n}+\mathcal{O}\left({m}^{-n}\right),\end{equation} \tag{ F.1 }$$

with [H_n , H₀] = 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 nth 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²/2m:

$$\begin{align}{H}_{2q-\text{body}}^{\left(1\right)}=\sum\limits _{{j}_{1}{< }\cdots {< }{j}_{2q}}\left[\tau \sum\limits _{n=1}^{2q}{\left(-\right)}^{n}{j}_{n}\enspace \left\vert \overrightarrow {j}\right\rangle \left\langle \overrightarrow {j}\right\vert +w\sum\limits _{n=1}^{2q}\left(\left\vert \overrightarrow {j}+{\hat{e}}_{n}\right\rangle \left\langle \overrightarrow {j}\right\vert +\text{h.c.}\right)-\frac{{w}^{2}}{2m}\sum\limits _{n=1}^{2q}\left(\left\vert \overrightarrow {j}+\hat{2}{e}_{n}\right\rangle \left\langle \overrightarrow {j}\right\vert +\text{h.c.}\right)\right]\end{align} \tag{ F.2 }$$

with the understanding that

$$\begin{equation}\left\vert {j}_{1},\dots ,{j}_{n}={j}_{n+1},\dots ,{j}_{2q}\right\rangle \equiv 0,\end{equation} \tag{ F.3 }$$

$$\begin{equation}\left\vert {j}_{1},\dots ,{j}_{n}={j}_{n+1}+1,\dots ,{j}_{2q}\right\rangle \equiv -\left\vert {j}_{1},\dots ,{j}_{n+1},{j}_{n},\dots ,{j}_{2q}\right\rangle .\end{equation} \tag{ F.4 }$$

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|.

F.2. 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):

$$\begin{equation}{\mathcal{E}}_{\text{free}}\left(k\right)=m\sqrt{1+\frac{4{w}^{2}}{{m}^{2}}+4\frac{w}{m}\mathrm{cos}\enspace k}.\end{equation} \tag{ F.5 }$$

When m approaches 2|w| (from above), the renormalized mass $\mu \equiv {\mathrm{min}}_{k}\enspace \mathcal{E}\left(k\right)=m-2\vert w\vert$ 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.

F.3. Decay of heavy mesons