Dynamical hadron formation in long-range interacting quantum spin chains

The study of confinement in quantum spin chains has seen a large surge of interest in recent years. It is not only important for understanding a range of effective one-dimensional condensed matter realizations, but also shares some of the non-perturbative physics with quantum chromodynamics (QCD) which makes it a prime target for current quantum simulation efforts. In analogy with QCD, the confinement-induced two-particle boundstates that appear in these models are dubbed mesons. Here, we study scattering events due to meson collisions in a quantum spin chain with long-range interactions such that two mesons have an extended interaction. We show how novel hadronic boundstates, e.g. with four constituent particles akin to tetraquarks, may form dynamically in fusion events. In a natural collision their signal is weak as elastic meson scattering dominates. However, we propose two controllable protocols which allow for a clear observation of dynamical hadron formation. We discuss how this physics can be simulated in trapped ion or Rydberg atom set-ups.


I. INTRODUCTION
The formation of hadrons in nature has its origin in the peculiar potential energy of the constituent elementary particles, i.e. the pairwise interaction energy of quarks grows indefinitely for increasing separation. This phenomenon, called confinement, is the key feature of quantum chromodynamics (QCD) [1 and 2]; and because of its non-perturbative nature, the QCD quantum many-body problem is one of the prime targets of recent quantum simulation efforts [3]. A timehonoured way of studying interacting quantum particles are controlled scattering experiments between different types of constituents of matter, most famously at the large hadron colliders [4 and 5].
Confinement is also important in condensed matter physics [6][7][8] and can be observed, for example, in quantum spin chains [9][10][11][12][13][14], whose elementary excitations are mesonlike boundstates of domain-walls. Of course, these onedimensional (1D) models are much simpler than the full SU(3) gauge theory of QCD, but they do share some of the underlying physics. Thus, basic 1D quantum spin chains may be a first step in simulating a full treatment of the strong force. Furthermore, the simulation of quantum spin chain models is within the capabilities of current quantum simulators, with promising first results [15][16][17][18][19][20][21][22].
Confinement has been shown to lead to exotic nonequilibrium dynamics manifest in a suppression of transport [11, 23, and 24] and entanglement spreading [25], a long lifetime of the metastable false vacuum [26][27][28][29] and exotic prethermal phases [30][31][32]. Recent efforts have been made in understanding scattering events among mesons [33][34][35] in the short-range Ising model with both transverse and longitudinal fields. Despite the rich scattering phenomenology, this setup does not host composite excitations beyond the twoquark mesons. Therefore, observing deep inelastic scattering with exotic particle formation needs a richer microscopic dynamics, as it can be realized, for example, by introducing heavy impurities in the short-ranged Ising chain [36].
In this work, we investigate another mechanism that can FIG. 1. A schematic of a scattering event between two wavepackets of confinement-induced mesonic (two-domain-wall) states. Apart from the elastic meson reflection, the long-range interaction between the composite particles may also lead to the formation of multiparticle boundstates, here in the form of a metastable four-domainwall state akin to a tetraquark.
lead to dynamical hadron formation by addressing the transverse Ising chain with long-range interactions which, in addition to confinement [12], naturally induces interactions among mesons and hosts multi-meson boundstates, akin to hadrons. We will show that the presence of long-range interactions between mesons can result in the formation of long-lived hadronic metastable boundstates, which can form dynamically in scattering events between fundamental mesons as depicted in Fig. 1.
The resulting metastable boundstates are different in nature from infinitely long-lived hadronic boundstates, which cannot be excited in real time scattering processes. Hence, natural metastable hadron formation is usually weak and the scattering events are dominated by reflection processes. To enhance the effect, we propose two different modifications of the scattering protocol. Firstly, by an abrupt dynamical change of an external field one can modify the mesons' kinetic energy and induce a non-trivial overlap between the initial asymptotic states and the hadronic bound-states of the post-quench Hamiltonian. As a result, infinitely long-lived boundstates are created and clearly observable. Secondly, we consider time-independent Hamiltonians with a modified long-range spin interaction, which results in a non-monotonic mesonmeson interaction. As a result, scattering mesons can tunnel to a local minimum of the relative interaction and dynamically form a long-lived metastable boundstate, which again becomes clearly observable in the collision.
Our work is structured as follows. In section II, we consider the low-energy sector of the Hamiltonian by projecting on the four domain-wall subspace. The reliability of such an approximation is related to the extremely long-lifetime of domainwall excitations, which scales exponentially in the weak transverse field akin to the short range Ising chain [24]. In section III, we present simulations of the real time dynamics of a collision between two large mesons and find signatures of metastable tetraquark formation. In section IV, we then show how a simple abrupt change of the kinetic energy can be utilized to allow tetraquarks to form dynamically in a more controllable setting. In section V, we further show a second method of inducing fusion by performing a local alteration of the long-range interactions in the spin chain Hamiltonian. We find that tetraquarks with a long lifetime can form dynamically without the need for a time dependent Hamiltonian. In section VI, we discuss the experimental feasibility of dynamical hadron formation in current quantum simulator platforms. In particular, we discuss initial state preparation as well as how to implement the modified long-range spin interaction within Rydberg atom and trapped ion systems. Finally, we conclude with a summary and discussion.

II. HADRONS IN THE LONG-RANGE ISING MODEL
In this work, we consider the one dimensional Ising spin chain with long-range interactions where i labels lattice sites, J gives the overall energy scale (we set J = 1 throughout this work without loss of generality), h is the transverse field strength and σ i are the Pauli matrices. The rich dynamics of the long-ranged Ising model has been recently probed in a trapped ion quantum simulator [15], showing clear signatures of confined excitations. We are interested in the regime where the transverse field is weak h J. The nature of the excitation is best understood by starting with the classical model obtained for h = 0: here, the natural excitations are domain-walls in the z−magnetization and the exponent α crucially determines their interactions. In shown for L = 20, α = 2.6 and h = 0.1. The energies strongly depend on the average size of the constituent mesons resulting in large gaps between levels of different meson sizes. A schematic picture of the local spin configurations within each subspace are depicted above the corresponding energy levels. Importantly, due to the long-range interaction the energy within each subspace also depends on the distance between the two mesons, i.e. the closer the mesons are together the lower their energy. The inset shows the energies of the 1-meson subspace as a function of momentum. The lowest energy levels in this subspace are boundstates of multiple constituent domain-walls. the limit of large α, the short-ranged Ising chain is recovered. Here, domain-walls are non-interacting objects. As α is reduced, longer ranged domain-wall interactions are induced. In particular, two domain-walls n sites apart feel an attractive potential V (n) of the form [12] and ζ(α) is the Riemann zeta function. For arbitrary α, the potential V (n) is a monotonically increasing function of n with a finite maximum lim n→∞ V (n) = V max < +∞ for α > 2.
In contrast, in the regime 1 < α < 2, V (n) is unbounded and domain-walls cannot be pulled infinitely far apart without paying infinite energy, thus showing confinement. Finally, for α < 1 the potential V (n) diverges, signaling a transition to an infinitely-long-ranged model without well-defined local excitations, therefore we will always consider the regime α > 1. The activation of a small transverse field h has two effects: i) domain-walls become dynamical objects and can hop along the chain, and ii) domain-walls in principle are no longer conserved. Nonetheless, resonances between different particle numbers states are greatly suppressed for small transverse field and can be safely neglected [24]. Thus, we can project the dynamics in subspaces with conserved number of domain-walls. For example, in the two domain-wall sector one can con-sider the basis |j, n = |↑ ... ↑↓ j ... ↓↑ j+n ... ↑ . In this notation, j labels the position of the first domain-wall from the left and n is the relative distance of the second domain-wall with respect to the first. The projected Hamiltonian takes the form [12] H = j,n V (n) |j, n j, n| − h |j + 1, n − 1 + |j − 1, n + 1 + |j, n + 1 + |j, n − 1 j, n| , and can be understood as a kinetic term for each domain-wall with hopping strength h and a potential V (n). As mentioned above, confinement in a strict sense requires 1 < α < 2 such that domain-walls cannot be separated. Nonetheless, also for larger α the two-kink subspace shows the presence of deep boundstates and asymptotic states of freely propagating domain-walls can be very high in energy and extremely difficult to excite. Therefore, as long as asymptotically propagating domain-wall states can be neglected, the dynamics of these deep two-kink boundstates closely resemble those observed in the strictly confined regime, and with a slight abuse of jargon we refer to both as "mesons". Due to the fact that larger values of α mitigate finite-size corrections caused by the long-range potential, we will mainly focus on the regime α ∼ 2.5 but stress that similar physics appears for 1 < α < 2.

(4)
Here, I can be seen as the meson interaction such that [28] I(j 1 , n 1 , j 2 , n 2 ) = −4 j1<r≤j1+n1 j2<s≤j2+n2 1 (s − r) α (5) and "hopping terms" refer to those of the form |j 1 ± 1, n 1 ∓, j 2 , n 2 j 1 , n 1 , j 2 , n 2 | or equivalently for the second meson. Note, the form of the interaction, I, depends on the choice of boundary conditions. Here, we consider open boundaries. This effective interaction between mesons falls of as d −α in which d is the meson separation ∼ (j 2 − j 1 ). Thus, similar to domain-walls that are spatially distant from each other, individual mesons which are far apart from each other interact only weakly.
In Fig. 2 we show some of the lower energy levels of the four domain-wall subspace. We clearly see large energy gaps in correspondence with internal quantum numbers, labeling the energy levels of the two kink-mesons. Additional structure is then provided by boundstates of fundamental mesons and asymptotic scattering states (see inset). In the case of deep boundstates where the binding energy is much larger than the transverse field h, the two-kink mesonic wavefunction is very peaked on integer values of the relative distance, hence the meson has approximately fixed length which is in one-to-one correspondence with the internal energy levels. In this regime, we can pictorially use the meson's size as a good quantum number. Nonetheless, this correspondence is blurred as boundstates become shallower and domain walls can oscillate with internal dynamics.
Before turning to real-time numerical simulations of scattering events, it is useful to further comment on the mesonmeson boundstate structure, i.e. our cartoon picture of hadrons. As shown in Fig. 2, energy levels corresponding to hadrons are clearly visible in the spectrum, their number and gap with respect to the asymptotic scattering states depends on the "size" (i.e. the energy) of the binding mesons. These boundstates are clearly orthogonal to asymptotic scattering states, hence scattering events cannot couple to them. If we wish to observe dynamical hadron formation, we should aim for metastable states arising from asymptotic scattering states whose wavefunction is large when mesons are close, which means that they are qualitatively close to true boundstates. Heuristically, these states are most likely to be present where the spectrum shows a smooth transition between boundstates and scattering states, i.e. when the energy gap between the two is small. As it is clear from Fig. 2, this is the case for large mesons. This picture is also suggested by semiclassical arguments, since for large mesons the relative position of the kinks can oscillate. Therefore, part of the scattering energy of the two incoming mesons can be converted to internal energy of the mesons upon scattering and a boundstate may form. This is not possible for tightly bound mesons, where the relative position of the kinks cannot be changed. These consideration motivate the use of large mesons in the following natural scattering protocol.

III. INELASTIC COLLISIONS OF LARGE MESONS
We now consider the Hamiltonian projected in the fourkink subspace and numerically explore scattering events between mesons. As initial states, we choose Gaussian meson wavepackets with a well-defined momentum. Furthermore, we fix the average length of each meson to cover a few lattice sites ∼ 4 − 10. Details on the wavepacket wavefunction are relegated to the Supplementary Material (SM) 37 . Simulations in the four kink subspace are challenging for large system sizes. Therefore, we take advantage of global momentum conservation and focus on the sector with zero total momentum. Using translational invariance we measure the kink coordinates with respect to the leftmost domain-wall which can reduce the Hilbert space dimension by a factor L allowing us to access much larger systems.
In the zero-momentum sector we can then label the Hilbert space by only three variables |n 1 , j 2 , n 2 as the first kink can be pinned to j 1 = 0. As a consequence of this convention, in our real-time simulation the leftmost meson will appear sta- The probability distribution of m d = j2 + n2 at times tJ = 0, the initial wavepacket, and tJ = 5000, after the collision event. After the collision, we clearly see that the majority of the wavepacket is reflected, however, there is a subtle peak at site ∼ 25 that corresponds to a bound tetraquark state. tionary.
In Fig. 3(a) we show a collision event between a left meson of initial width n 1 = 10 with a right meson with initial width n 2 = 4. We employ the observable ) |ψ | 2 which can be seen as the probability that a domain-wall is located on site i.
In Fig. 3 we observe that scattering is dominated by elastic reflection events. However, by plotting the logarithm of the data we are able to highlight subtle details, e.g. a small portion of the wavepacket, ∼ 10 −3 remains close to the stationary meson for long time after the collision which can be seen as a vertical line of intensity at site ∼ 25 Fig. 3(a). Next, we study a basic measure of the 'distance' between the two mesons m d = j 2 + n 2 − j 1 . In Fig. 3(b) we plot the probability distribution of m d before (tJ = 0) and after the collision (tJ = 5000). After the collision the distribution has two components. First, the overwhelming probability accounts for reflection events seen as the large hump for sites > 40. However, there is a second, small peak around site ∼ 25 which is a signature of our sought-after metastable tetraquark state.
For all parameter and initial state choices we have explored the signatures of dynamical tetraquark formation are weak as elastic scattering dominates. Therefore, in the following, we propose two simple extensions of the natural meson scattering protocol which allow for an unambiguous observation of dynamical tetraquark formation.

IV. ABRUPT CHANGE IN THE TRANSVERSE FIELD
In the previous section we have only observed weak signals of dynamically formed tetraquark state. However, true boundstates are clearly present in the spectrum but, as we have already mentioned, these are orthogonal to the scattering states of our wavepacket. In the following we explore the possibility to artificially induce non-trivial overlaps between these two classes of states by dynamical changes in the Hamiltonian. As the most basic example, we consider abrupt changes in the transverse field h at the time of collision t * . We note that the very nature of the so-obtained hadrons is very different from the previously considered metastable states because now true boundstates are excited with an infinitely-long lifetime. Another advantage is that within this scenario, we do not necessarily need to target large mesons, which are challenging for experiments. Therefore, we focus on small mesons whose energy levels are approximately in one-to-one correspondence with their length. Hence, we refer to a meson of length n as a n−meson.
In the following we focus on the simplest case of 1mesons. However, as one could argue that these are just simple magnons, we confirm similar physics for nontrivial 2mesons in the SM [37] For a state consisting of two 1−mesons, a single spin flip naturally leaves the 1−meson subspace. Thus, in order to get an intuitive feeling about the dynamics within a restricted 1−meson Hamiltonian we must consider higher order processes. The derivation can be obtained through perturbation theory as detailed in the SM [37]. The 1−meson subspace has the simplified basis |j 1 , j 2 = |↑ ... ↑↓ j1 ↑ ... ↑↓ j2 ↑ ... ↑ and the Hamiltonian is We obtain two contributions from second order perturbation theory, an inhomogeneous hopping term, h j2−j1 , and an FIG. 4. The dynamical formation of a tetraquark with an abrupt change of the transverse field (calculated in the full four domainwall subspace). Here, the Hamiltonian parameters are L = 50 and for α = 2.6. The full details of the initial wavepacket can be found in the SM [37]. In section A (up to the yellow dashed line) two 1meson wavepackets move toward each other with a large kinetic energy (transverse field h = 0.2). Then at the point of meson collision (tJ = 125) a abrupt change to a smaller transverse field (h = 0.1) is performed reducing the kinetic energy in section B. This reduction of kinetic energy can be seen by comparing the fast velocity of the meson wavepackets in section A to the slow velocity in section B. We see that this process induces the formation of a tetraquark with exotic internal dynamics. additional effective interaction, U j2−j1 . We can take advantage of translational invariance and focus on the sector with global momentum k, where we define the relative distance j = j 2 − j 1 and the momentum-dependent Hamiltonian H k acting on the states |k, j is given by which we can use to calculate the energies as a function of momentum. We note that the energies of the low energy 1−meson sector are in very good quantitative agreement with those of the full four domain-wall subspace as well as full exact diagnonalization (ED), see SM [37].
In the inset of Fig. 2 we show the spectrum of this 1−meson subspace. One can clearly see that the lowest energy levels are discrete with large gaps between them, which correspond to boundstates of two mesons that reside only a few sites away from each other, i.e. they describe tetraquarks. In addition, there is a continuum of levels at higher energy of far-apart mesons which interact only weakly, i.e they are free to propagate.
In Fig. 4 we implement the proposed protocol. We initialize the 1−meson Gaussian wavepackets with well-defined momentum and suddenly reduce the transverse field at a time t * when the scattering takes place (dotted horizontal line in the figure). The main part of the signal remains indeed trapped in a boundstate at the center of the chain, thus realising the de-  Fig. 4 showing the internal dynamics of the dynamically formed tetraquark. An inset shows the spectrum of the oscillation frequencies observed. Clear agreement is seen between the dominant frequencies observed and differences in the energy levels Ei of the 1−meson subspace after the change of field to h = 0.1, these are shown by the vertical red lines. sired dynamical hadron-forming event.
In Fig. 5 we show the dispersion relations of the pre/postquench Hamiltonians as a function of the total momentum k. Boundstates appear as well-separated bands below a highenergy continuum of asymptotic states. Notice that the boundstates for h = 0.1 have an almost flat dispersion, hence very small velocities, which explains the trapped stationary signal observed in Fig. 4. Albeit hardly moving, the trapped boundstate is oscillating in time. In Fig. 5(b) we show that the oscillating frequencies are indeed compatible with the energy gaps between the boundstate energies of the post-quench Hamiltonian.

V. MODIFIED LONG-RANGE INTERACTIONS
While the abrupt change of the transverse field protocol is very efficient in exciting boundstates, it can arguably be regarded as an artificial shortcut to our goal of observing multimeson binding in scattering events. Therefore, we now consider possible modifications of the interactions which can enhance the phenomenon without time-dependent changes.
It is again useful to focus on the 1−meson subspace and consider the relative interaction between two mesons shown in Fig. 6(a) which is a monotonically decreasing function of the relative distance (dashed line). An enhancement of bound- In addition to the main elastic scattering channel, a small portion of the wavepacket quantum tunnels into the tetraquark well forming a long-lived tetraquark state. In the left hand plot we present the real time results for PK (i) of a meson collision. Here we see that while the majority of the incoming meson wavepackets are reflected, some tunnel into the tetraquark well and persist for long times. In the right hand plot we explicitly show the probability of a tetraquark forming by presenting PT Q(i) for the same collision event.
state trapping can then be achieved by modifying the spin interactions to form a potential well. For example, this can be achieved with the minor modification of the Hamiltonian inducing a potential well in the relative potential at separation d, see solid curve in Fig. 6(a). Two colliding mesons will now see a potential barrier of finite width and a quantum tunneling event may occur in the collision. In the following, we again focus on 1-mesons for simplicity, but similar boundstate formation can be observed for scattering of larger mesons, see SM [37]. Physically, this can be understood as increasing the energy of states that have up spins d sites away from each other. In terms of the 1−meson Hamiltonian this is equivalent to removing the interaction between mesons that are d sites away from each other. Fig. 6(b) shows the effect of the modified interaction for allowing a long-lived tetraquark to form dynamically. Here, we initialized two meson wavepackets with opposite momentum. We observe, as expected, that in the presence of the effective potential well there is a small probability that the mesons tunnels through the potential barrier forming a metastable tetraquark decaying with a lifetime τ . We also show the results of measuring P T Q (i) = C |ψ(j 1 , n 1 , j 2 , n 2 )| 2 where C is the set of basis states such that there is a domain-wall at the i th site and j 2 − (j 1 + n 1 ) ≤ d. This observable can be seen as the probability that a tetraquark is located around site i. From this we indeed see that, at the point of collision, a tetraquark forms and has a long lifetime (right panel).
The decay of a metastable boundstate via quantum tunnelling is reminiscent of one of the first paradigms of quantum physics -α-particle decay as observed about a century ago and explained by Gamow's famous semi-classical the-ory [38]. We can go one step further with our spin chain analogy and compare the decay times of bound tetraquarks for varying energies, which can be achieved through varying the depth of the potential well, J d α → J d α − . In Fig. 7(a) we plot i P T Q (i) to observe the decay of states initialized as a tetraquark wavepacket for varying potential well depths. We observe that tetraquarks decay exponentially in time. The lifetime is longer for tetraquarks with lower energy which is consistent with the larger potential barrier between the bound tetraquark and free meson states.
Similar to Gamow's theory of α-particle decay, we can utilize a WKB approximation to calculate the transmission probability of a tetraquark through the potential barrier. We note that our lattice calculation with a rapidly changing potential is beyond the strict range of applicability of any semi-classical approximation. Nonetheless, we will show decent agreement with our numerical calculations which corroborates the use of spin chain simulators of small size for probing confinement physics.
We only summarize the results of the WKB calculation with full details given in the SM [37]. The transmission coefficient of a meson wavepacket is where x 0 and x 1 are the classical turning points of the potential well and the momentum is One can estimate the probability of emission at any given time by 2x1 v T in which v is the average velocity of the boundstate. Due to the exponential nature of T it gives the dominating contribution to the lifetime τ , which can be predicted within the WKB theory as with a constant C [38]. In Fig. 7(b) we show that, by fitting the constant C, we see a remarkably good agreement between the WKB theory calculations of the lifetime of tetraquarks and those extracted from the numerical simulations up to the limit of a deep tetraquark well, E J < 10.045.

VI. EXPERIMENTAL FEASIBILITY
A physical realization of our meson scattering protocols requires the implementation of the one dimensional transverse field Ising model with power-law interactions as described by Eq. 1. Although there are several platforms available for the experimental realization of long-range Ising models with hundreds of spins, e.g. with polar molecules [39 and 40], trapped ions [41] and Rydberg atoms [42][43][44][45], they are often realized in higher dimensional lattices. However for one dimensional systems, the number of spins available for quantum simulation range from 3-53 for ions trapped in a Paul trap [46][47][48][49] or 20-51 in the case of trapped array of Rydberg atoms in optical tweezers [50 and 51], thus reaching the necessary system sizes for our scattering protocols.
In these setups, the many-body state with all spins down can be naturally realised. Then the initial state preparation of a pair of mesons (domains) located far apart from each other can be achieved through local quenches with help of lasers by performing flips on one or more spins at appropriate sites. The lattice spacing and the width of the laser beam determines the size of the domains which has been demonstrated both for trapped ions [20, 52, and 53] and Rydberg atoms [50, 51, and 54]. For initial states to propagate along the spin chain with well defined momentum requires controlled nearest neighbour spin exchange. For trapped ions, spin exchange between sites is engineered using Mølmer-Sørenson type protocols [52]. Effective transitions between |↑ and |↓ are driven with the help of bichromatic laser with frequencies ω ± = ω 0 ± ∆ that shines upon the ion array, where ∆ is the detuning of the laser field from the transition frequency ω 0 of the atomic transition |↑ ↔ |↓ . While for Rydberg systems, the spin exchanges are achieved either by resonant dipole-dipole interaction or by off-diagonal van der Waals flip-flop interaction between Rydberg states [55]. In some case, spin-exchange can also be induced via the standard van der Waals interaction [56].
Apart from the natural scattering dynamics of the mesons, we have proposed two other effective routes for the dynamical formation of tetraquarks: (i) by abruptly changing the kinetic energy as depicted in Fig. 4, and by (ii) fusion of mesons induced via a potential well as shown in Fig. 6(a). Method (i) can be implemented by abruptly changing the time dependent transverse field. Such type of quench dynamics are commonly executed in both systems, trapped ions [15 and 57] as well as Rydberg atoms [58 and 59]. An advantage of trapped ions is that one can easily tune the strength of the long-range interactions (1 < α < 3) in the Ising model with the help of Raman lasers as recently demonstrated in the investigation of domain-wall dynamics [15].
However, the approach (ii) requires a non-trivial modification of the spin-spin interactions to a non-monotonic longrange form. Such exotic potentials may not be obvious for trapped ions, but in the following we sketch how they can be implemented using the highly tunable effective interactions of laser-dressed Rydberg atoms [60][61][62][63]. A dressed Rydberg atom is primarily a ground state atom that is weakly superposed with an excited state corresponding to a large principal quantum number (n ∼ 20 − 100) [64]. The amount and type of the Rydberg character in the dressed superposition state is controlled using dressing lasers that eventually determine the strength and shape of the effective potential [65 and 66]. Recent experimental validations of Rydberg-dressed interactions include the measurement of pairwise interaction between two atoms [67], many-body Ising interactions [68 and 69], as well as distance selective interactions [70].
The key idea for our modified long-range interaction is to construct a spin-1/2 state using two-long lived states (|g ± ) of an atom which can either be a pair of hyperfine ground states or a ground state and a meta-stable state found in alkalineearth atoms [66]. A pair of lasers (left and right circularly polarized with phases set to zero without loss of generality) drive the atoms from |g ± states to Rydberg states |e ± offresonantly with detunings ∆ ± and Rabi frequencies Ω ± . The full Hamiltonian in the atomic basis is where the individual terms are the following where V (r ij ) = C 6 (θ, φ)/r 6 ij is the van der Waals interaction between Rydberg atoms located at sites i and j and α, α ∈ {±}.The effective spin-spin interactions between two atoms can then be obtained with respect to their ground states in the weak coupling limit Ω ± E α,α , where E α,α are the eigenenergies of the relevant atomic systems, Using perturbation theory, one derives the effective interactions between the ground states which are expressed in the two-atom basis as follows, (15) The second order terms in the perturbation theory correspond to light shifts which serve as longitudinal fields for the spin Hamiltonian, while the fourth order terms provide the Ising interaction and transverse field term. The anisotropy in the van der Waals interaction depends on the magnetic quantum numbers and the relative angles (θ, φ) between the atoms with respect to the laser.
Overall, our main point is that there is sufficient tunability of the spatial interaction profile in order to engineer a local minimum of the effective interaction similar to the one of Fig. 6. However, while this effective modified potential between mesons can be achieved through the van der Waals interactions, its asymptotic decay 1/r 6 ij is not enough to ensure the stability of individual meson states. In order to stabilise them one can exploit the long-range dipole-dipole interactions which is also present. For this purpose, one then needs to consider both levels of the two-level system as Rydberg states for which bound states have indeed been shown to exist [71]. Alternatively, switching on a weak longitudinal field would also ensure the formation of individual meson states while there mutual interaction would be governed by the longer range potential discussed above.
A detailed quantitative modeling of experimental protocols is beyond the scope of this work, however we provide here some estimate for coherence times. In case of Rydbergdressed setups, the typical values of Rabi frequencies are in the range of tens of kHz to few MHz while the detunings are an order of magnitude larger. Rydberg states ranging from n = 50 − 70 for Rb atoms will have bare lifetimes around 80 − 130µs [72]. However as a result of the weak coupling to Rydberg states, the lifetime of the effective two-level system is extended to milliseconds [70]. For such Rydberg states, the interactions are in the order of GHz [72] and thus for lattice spacings of 0.5-1.5 µm, one can simulate effective spin-spin interactions which are in the order of few kHz [68] to hundred of kHz [67]. Assuming J = 800 kHz and the time of scattering occurring at tJ = 100 (after optimization of the protocol dynamics), coherence times of 0.125 ms are well within the reach of many-body Rydberg-dressed state lifetimes [70]. Although the dynamical timescales shown in this theoretical work are beyond what has been observed in experiments previously, the use of quantum optimal control theory [73][74][75] to reduce these timescales is promising especially with the access to a variety of control parameters such as maximising transverse field h while maintaining the domain wall approximation, minimising the initial meson separation, tuning the initial meson size and exponent α. Indeed achieving the dynamical timescales and coherence times needed for this work is a challenge for future experiments.

VII. CONCLUSIONS
We have shown that confinement-induced boundstates of many elementary domain-wall excitations exist in the longrange Ising model akin to composite hadronic particles in QCD. We studied the collision of simple two-particle mesonic bound states and found signatures of the dynamical formation of metastable four-domain-wall states in analogy to tetraquarks. However, natural meson collisions are mainly elastic and the signal for natural hadron formation is weak. Therefore, we proposed two alternative protocols to induce dynamical fusion events in the long-range Ising model that are much more controllable and allow for an unambiguous detection of dynamical hadron formation.
First, we show that an abrupt change in the transverse field at the time of collision results in a strong signal of tetraquark formation with interesting internal dynamics. Second, we find that a modification of the long-range interaction leads to a tetraquark potential well. Again, we obtain a clear signal of hadron formation and subsequent decay which can be understood via a semi-classical WKB approximation in analogy the to famous example of α-particle decay.
Finally, we argue that all three of these methods, while challenging, are in principle realisable with current quantum simulator set-ups. In particular, we sketch the experimental requirements for initial state preparation as well as the use of laser-dressed Rydberg atoms for engineering the modified long-range interactions.
Our work motivates a number of future research questions. For example, we have focused on the limit in which the elementary low energy excitations are well approximated by domain-walls. However, it is well known that for the shortrange transverse field Ising part of our Hamiltonion domainwalls continuously evolve into fermionic excitations beyond the small transverse field regime. It would be interesting to gain insight departing from the limit of extremely weak transverse field where, most likely, sharp domain walls will be deformed similarly to what happens in the short-range Ising upon activation of a finite transverse field. A next step toward the long-time goal of understanding fusion events in full QCD would be the quantum simulation of scattering events with dynamical hadron formation in Hamiltonians that lead to non-mesonic hadrons, such as the q-state Potts model with q > 2 (q = 2 corresponds to the Ising model) [76][77][78][79]. Furthermore, one could also consider simplified lattice gauge theories (LGTs). For example, 1D versions of U (1) LGTs have been studied intensively in the past years with many connections to quantum simulation architectures [19 and 20]. A prime candidate would be the 1D Schwinger model in which confinement dynamics can be studied efficiently with time-dependent DMRG methods [80]. In the longer run, similar albeit richer physics is expected in higher dimensional confining LGTs. We hope that emulating particle physics scattering experiments in toy models, whose realisation is feasible for current quantum simulators, can provide a first step towards a better understanding of the fascinating physics of large hadron colliders. As quantum simulators are slowly but surely increasing in size and quality, the simulation of dynamical hadron formation in 1D and beyond would also provide a crucial benchmark towards achieving a genuine quantum advantage. To first order in h, projecting the long-range Ising model given in Eq. 4 to the 1−meson subspace does not include dynamics of the mesons as any spin flip process will leave this subspace. This is depicted in Fig. S1. However, if we consider a second order effect in h we then find terms that allow for these 1−mesons to hop. The simplest of these are two spin flips which would act similarly to |j 1 , n 1 , j 2 , n 2 → |j 1 ± 1, n 1 ∓ 1, j 2 , n 2 (or the equivalent expression for the other meson). Thus, we must use perturbation theory to accurately describe these 1−mesons.
FIG. S1. A schematic of the second order hopping process required for a 1−meson to move.
Our starting point for this projection is the Hamiltonian H given in Eq. 4 which may be factored into a diagonal part H 0 and an off-diagonal part V proportional to what we will consider the perturbative parameter h, that is H = H 0 + hV . Our goal is to find a unitary transformation generated by an operator S such that the expansion H ef f = e S He −S ∼ H 0 + h 2 M + O(h 3 ) i.e. the lowest order hopping events are proportional to h 2 . V contains the hopping terms of the Hamiltonian; the resulting effective Hamiltonian thus contains only on-site potentials and two-site hopping interactions (plus higher order terms which are absorbed into the O(h 3 ) factor). This rotation is known as the Schrieffer-Wolff transformation. To proceed, we write down the Baker-Campbell-Hausdorff expansion of e S He −S up to second order in h: We can eliminate the first order terms by choosing S such that [S, H 0 ] = −hV . The simplest way to do this is element-byelement, enforcing that where |p , |q are eigenvectors of H 0 with eigenvalues E p , E q . With this enforced, the effective Hamiltonian becomes We project this effective Hamiltonian into the 1−meson subspace and obtain two contributions to the Hamiltonian, a hopping coefficient, h n,n+1 in which the 1−meson moves one site. There is also a second order effect in which a 1−meson can 'hop to itself' through a process in which a neighbouring spin flips and then flips again to leave the initial state unchanged, h n,n . In the 1−meson subspace this can be seen as an effective on site potential. These contributions are given by where j = j 2 − j 1 . Given these contributions the 1−meson subspace can be written as in Eq. 6. In Fig. S2 we compare the energies of the 1−meson subspace with the corresponding energies of the four kink subspace as well as the full Hilbert space.
Here we see that not only does the 1−meson subspace have good agreement with the four kink subspace as expected but also with the full Hilbert space.

MESON WAVEPACKETS USED AS INITIAL STATES IN SIMULATIONS
The initial state used in the collision of mesons in Fig. 3 is given by where q is the average momentum of the mobile meson with a standard deviation σ k , c is the average centre of mass of the mobile meson with standard deviation σ c and N i is the average width of the i th meson with standard deviation σ i . In particular we use the parameters q = −1, σ k = 0.1, σ c = 2, N 1 = 10, N 2 = 4, σ i = 1 for i ∈ 1, 2 and c = 60.
The initial states used in the collision of mesons in Fig. 4 and Fig. 6 are given by where k is the initial momentum of each meson, c i is the centre of the wavepacket of the i th meson and σ is the standard deviation of the wavepacket in real space. In particular the parameters used in Fig. 4  Finally, the tetraquark wavepacket used in the simulations presented in Fig. 7 is given by where k is the tetraquark momentum (which we set to zero), c is the average position of the first meson and σ is the standard deviation of the wavepacket. In particular the parameters used are d = 4, k = 0, c = 18 and σ = 5.

TRUNCATION OF MESON WIDTH IN SIMULATIONS
In order to simulate large system sizes we turned to the four domain-wall subspace described by the Hamiltonian given by Eq. 4. The size of this subspace grows as L 4 , here L is system size. This still poses a problem when considering the system sizes required to observe collision events. In order to reduce the computational cost of simulations to quadratic growth, we used a truncation in the meson width, n 1 , n 2 ≤ N for some constant N . When considering small mesons, varying the width of a meson vastly changes the energy of the state. Thus, in order to conserve energy, the meson width will only experience small fluctuations. As a result, accurate approximations of the full four domain-wall subspace are achieved by only considering n 1 and n 2 values below some constant that is at least larger than the meson widths of the initial state. Fig. S3 compares time dynamics of mesons for initial states of meson width 1 and 2 simulated in the four domain-wall subspace with varying truncations. We clearly see that both 1−mesons and 2−mesons are accurately simulated with a truncation of N ∼ 5. This approximation will naturally become less accurate as the initial state contains larger mesons.

DYNAMICAL HADRON FORMATION FOR LARGER MESONS
In this section we would like to confirm that our dynamical protocols enabling fusion events in meson collisions is not exclusive for 1−mesons. One could argue that 1−mesons are special because they do not have any internal dynamics and are simply equivalent to magnon spin flips. Here, we repeat both the abrupt change in the transverse field as well as the local alteration of the long-range interactions in the Hamiltonian for initial wavepackets composed of 2−meson. Fig. S4 shows the result of the collisions which are very similar as the 1-meson collision shown in the manuscript. In these we clearly see that the protocols outlined in the this work are not unique to 1−mesons. The case of the abrupt change in the transverse field also results in continued oscillations localized for a long time around the scattering region for larger mesons. In the case of the induced tetraquark potential well we again see a long lived metastable tetraquark formed during the collision.  Fig. 4 are observable except that the oscillations of the tetraquark formed have a larger amplitude as well as slower frequencies, which is expected for a tetraquark formed from mesons with larger masses and non-trivial internal dynamics. (b) Performing a local alteration of the long-range interactions consisting of two terms of the form 8 with d = 4 and d = 5. Here, the Hamiltonian parameters are α = 2.6, h = 0.1. The wavepacket parameters are k = −1, c1 = 10, c2 = 30 and σ = 3. We clearly see that in the collision a long lived metastable tetraquark is formed for 2−mesons.

WKB APPROXIMATION FOR TETRAQUARK TRANSMISSION PROBABILITY
In this section we derive transmission probability of a meson through a potential barrier of V (n) for a Hamiltonian of the form of Eq. 7 with k = 0 given by H = n J n |n + 1 n| + |n n + 1| + V (n) |n n| (S9) where, J n is the site dependent hopping strength. We use the WKB approximation, ψ(n) = N e φ(n) , in which φ is a complex function and N is a normalization constant. This leads to the eigenvalue equation Ee φ(n) = J n e φ(n+1) + J n−1 e φ(n−1) + V (n)e φ(n) (S10) We make the approximation that e iφ(n+dx) ∼ e φ(n)+dnφ (n) such that dn = 1 leading to E = J n e φ (n) + J n−1 e −φ (n) + V (n). (S11) Given that φ(n) is smooth, we let φ(n) = n q( )d . Now (V (n) − E) + J n e q(n) + J n−1 e −q(n) = 0. (S12) We can solve this quadratic equation in e q(n) to obtain It turns out that q + ∼ −q − away from the potential barrier, thus, we can think of q + as a particle travelling to the right and q− as a particle travelling to the left. Thus, given that ψ 1 is the wavefunction in the potential well, ψ 2 is wavefunction in the classically forbidden region and ψ 3 is the wavefunction after the potential wall we set up the problem as where R is the reflection coefficient, T is the transmission coefficient and n = a and n = b are the classical turning points.