Quasi-soliton scattering in quantum spin chains

The quantum scattering of magnon bound states in the anisotropic Heisenberg spin chain is shown to display features similar to the scattering of solitons in classical exactly solvable models. Localized colliding Gaussian wave packets of bound magnons are constructed from string solutions of the Bethe equations and subsequently evolved in time, relying on an algebraic Bethe ansatz based framework for the computation of local expectation values in real space-time. The local magnetization profile shows the trajectories of colliding wave packets of bound magnons, which obtain a spatial displacement upon scattering. Analytic predictions on the displacements for various values of anisotropy and string lengths are derived from scattering theory and Bethe ansatz phase shifts, matching time evolution fits on the displacements. The time evolved block decimation (TEBD) algorithm allows for the study of scattering displacements from spin-block states, showing similar scattering displacement features.


I. INTRODUCTION
The study of classical dynamics in nonlinear media has proven to be a source of astonishing surprises over the last century. Two observations, both based on numerical simulations, have challenged prejudices and fundamentally altered traditonal ways of thinking. First, the famous observation of Fermi, Pasta, Ulam, and Tsingou 1 of a simple nonlinearly-coupled set of oscillators showing nontrivial recurrences, has shattered the long-held assumption that all nonlinear dynamical systems ergodically explore their full phase space. Second, the pioneering numerical analysis of Zabusky and Kruskal 2 on the Korteweg-de Vries equation 3 demonstrated that this equation supports excitations, which they coined 'solitons', displaying a number of surprising fundamental features. The solitons are localized in space, with a form remaining stable under time evolution which sees them moving uniformly at a speed linearly proportional to their amplitude. Additionally the astounding characteristic was observed of solitons emerging intact from mutual scattering processes, during which they simply '"pass through" one other without losing their identity', 2 the only effect of the collision being a relative spatial displacement as compared to their free propagation. The proper understanding of solitons in nonlinear media ultimately led to the development of the classical inverse scattering method, 4,5 which is the overarching framework for classical integrable models.
Nonlinear classical systems find their quantum mechanical analogue in the shape of interacting many-body systems. In parallel to the classical case, some quantum models have been shown to be special, in the sense of being exactly solvable using the quantum version of the inverse scattering method 6 (one might also say integrable, though the quantum notion of integrability is not as well-defined as its classical counterpart 7 ). Fundamental representatives of this family are the Heisenberg spin chain, solved by Bethe using what is now known as the Bethe ansatz, 8 along with the Lieb-Liniger model of δ-interacting bosons on a line. 9 For the latter, particlelike excitations called Lieb Type I modes 10 exist due to the interparticle repulsion, along with Type II hole-like excitations visualized as holes in an effective Fermi sea. It is possible to distinguish the presence of Type I and II modes in correlated bosonic gases in optical lattices using Bragg spectroscopy. 11,12 Quantum magnets such as the Heisenberg spin chain similarly carry particle-like magnon modes when the magnetization is close to saturation. In the limit of small magnetization, hole-like modes again appear, which in zero field are known as spinons. 13 Their dynamics can be experimentally observed using inelastic neutron scattering. 14,15 The Heisenberg chain supports distinct bound states of magnons, whose dynamics has been investigated theoretically 16 and has recently been observed experimentally. 17 One could view such excitations as the quantum equivalents of classical solitons. This equivalence is however only partial: on the one hand, these quantum mechanical modes represent exact eigenstates and are stable under time evolution; on the other hand, being exact eigenstates of translationally-invariant systems, they are not spatially localized. That said, as is usually the case in quantum mechanics, it is possible to adopt a 'complementary' picture and create spatially localized excitations by forming wave packets of fundamental excitations by linearly combining states over a range of differing momenta. Locality however comes at a price: the wave packet, mixing together states at different energies, will disperse and is thus not stable over long timescales, unlike its classical counterpart. We use the term quasisoliton for such a wave packet construction, an example of which was recently studied in the context of the Lieb-Liniger model, 18 while their mutual scattering has been studied in quantum spin chains. 19,20 The spectroscopic methods traditionally employed to experimentally study condensed-matter systems typically provide momentum-and energy-resolved measurements. However, current experimental developments provide motivation to obtain a better theoretical understanding of spatially localized dynamics. Time-resolved experiments are now able to track quantum many-body systems at timescales smaller than relaxation timescales, particularly in experiments with ultracold atoms, 17,21,22 but also in NMR setups, 23 and potentially in pump-probe spectroscopy experiments. 24 As a result, non-dissipative dynamics of many-body systems out of equilibrium is now a rapidly growing field of experimental as well as theoretical investigation. 25 Novel cold-atom experimental techniques for spatially resolved manipulation and observation at the single-site level 17,22,26 have opened the door to explicit high-resolution tracking of spatial propagation phenomena. Moreover, an experiment with interacting bosonic atoms has highlighted the interactioninduced longevity of repulsive pairs. 27 This has motivated increased theoretical attention to the (anti-)binding of localized excitations and interactions of these bound clusters, both in itinerant systems 19,28 and in spin chains. 16,19,20,29 Very recently, the spatial dynamics of itinerant clusters has been studied in experiment. 30 Propagation of quantum solitons in Bose-Hubbard chains has also been studied numerically. 31 In Ref. 19, the scattering of a magnon wave packet on approximate bound eigenstates of n particles was studied numerically in the Heisenberg chain. In the present work, we present an algebraic framework and exact calculations based on Bethe ansatz. We therefore consider quantum scattering of localized excitations over a ferromagnetic background in the Bethe ansatz solvable anisotropic spin-1/2 Heisenberg chain (XXZ model) 8,32,33 (1) The XXZ model is experimentally realizable, for example in the setups of Refs. 17 and 22, which use two hyperfine states of bosonic atoms in the Mott phase to experimentally realize the spin-up and spin-down states. The effective model is a nearly isotropic (∆ ≈ 1) Heisenberg chain, while an experimental setup with variable anisotropy ∆ is under development. 34 In addition, the XXZ model has been shown to describe Josephson junction arrays of the flux qubit type, 35 and may also be realizable in optical lattices 36 or with polaritons in coupled arrays of cavities. 37 It is conceivable that these or similar experimental setups may provide time-resolved observations of propagating and interacting localized excitations.
The parameter J in Hamiltonian (1) is given by the exchange interaction of two neighboring electrons or, in the aforementioned experimental setup, by the exchange interaction between two neighboring bosons in an optical lattice at unit filling. Its sign does not matter for our purposes here; 19,38,39 for definiteness, we set it to J > 0. We distinguish two regions for the anisotropy parameter ∆, namely the planar xy (|∆| < 1) or axial z (|∆| > 1) cases. Due to the interaction term (∝ ∆) the XXZ chain displays a whole zoology of fundamental excitations: isolated down spins can form spatially bound states (as was understood by Bethe already in his original publication 8 ), whose bond size decreases as ∆ increases. These excitations are often referred to as 'string states', as the rapidities describing Bethe states containing such bound multi-magnons appear as approximately equallyspaced vertical strings in the complex plane (the precise set of available bound states depends on the value of ∆; at the isotropic point ∆ = 1 and in the axial regime, all string lengths are allowed).
We construct spatially localized wave packets of n bound magnons using linear combinations of these string states of length n with Gaussian-distributed momenta. We call them 'n-string wave packets'. We further demonstrate that exact methods based on the algebraic Bethe ansatz 6 provide a framework to evaluate the timedependent expectation value of the local magnetization S z j (t) algebraically, which can be used to track those localized magnon-like wave packets. We investigate their stability and mutual scattering using a combination of scattering theory, Bethe ansatz, and numerically exact calculations.
At large anisotropy ∆ 1, magnon bound states resemble having downturned spins on neighboring sites. An n-string wave packet is thus closely approximated by a consecutive block of n downturned spins. While this correspondence breaks down at smaller ∆, this provides motivation to study the evolution of states with downturned spins on a consecutive block of sites. In addition, this is exactly the type of initial state prepared in experiments. 17,22 Scattering of such blocks has been explored in Ref. 19, where a spatial displacement of two sites was observed for the block at several ∆, and explained at large ∆ in terms of energy conservation. In the present work, we connect scattering phase shifts with trajectory displacements in order to provide a Bethe ansatz derivation of the observed displacements.
The paper is organized as follows. In Sec. II we will introduce the concepts of string solutions and the scattering phase shifts associated with n-strings, along with details on the algebraic Bethe ansatz 6 based evaluation of the time-dependent expectation value of the local magnetization S z j (t) , which can be used to track localized magnon-like excitations of the spin chain. In Sec. III, we elaborate on the construction, stability and time evolution of scattering of n-string wave packets. We also compare with the stability of consecutive-site spin blocks. In Sec. IV we consider the scattering trajectory displacements, derive analytical results for them and compare with numerical measurements. The ∆ → ∞ limit is treated analytically, comparing with numerical data for spin block initial states. The appendices provide details involving scattering theory and details on obtaining the phase shift directly from the phase of the time-evolving wave function.

II. TIME EVOLUTION FROM BETHE ANSATZ
In this section we present the basic formulas of the Bethe ansatz for the spin-1/2 XXZ model and explain the concept of string solutions and the scattering phase shift associated with n-strings. In Subsec. II C, expressions are given for the time evolved expectation value of the local magnetization by using results from algebraic Bethe ansatz. The method used here relies on the availability of determinant expressions for matrix elements between Bethe states. [40][41][42] This last subsection is relatively technical and could be skipped on first reading.

A. Coordinate Bethe ansatz for the XXZ model
The eigenstates of the XXZ spin chain (1) can be constructed via the Bethe ansatz, 8,33 and have the form where M denotes the number of downturned spins and therefore fixes the magnetization. The sum over Q is a sum over all permutations of M objects and the amplitudes A Q are related to the scattering phases. The set of M complex rapidities {λ} ≡ {λ j } M j=1 completely determines the Bethe state and is simply related to the physical energy and momentum. By imposing periodic boundary conditions, i.e. S α N +1 = S α 1 for α = x, y, z in the Hamiltonian (1), each set of rapidities {λ} corresponding to an M -magnon eigenstate must obey Bethe equations, The different definitions of φ n (λ), θ n (λ) and ζ for various regions in anisotropy ∆ are given in Tab. I. Within Bethe ansatz, the momenta of single downturned spins can be parametrized in terms of rapidities, By invoking Schrödinger's equation H|{λ} = E|{λ} for a Bethe state consisting of a single downturned spin, the magnon dispersion relation is easily derived. In the case of just two single magnons, their scattering phase shift χ can be obtained from the permutation of two magnons in the Schrödinger equation, where Q is the identity and Q interchanges the two indices 1 and 2. Furthermore, the magnon momenta p 1 = p(λ 1 ) and p 2 = p(λ 2 ) as well as the scattering phase shift χ = θ 2 (λ 1 − λ 2 ) are parametrized by the two rapidities λ 1 and λ 2 .

B. Strings and magnons
The sets of rapidities solving Bethe Eqs (3) are selfconjugate and arrange themselves in patterns of string solutions, where the string center λ (n) j ∈ R and a = 1, . . . , n is the internal label of a rapidity within a string of length n and parity ν j . In the planar regime |∆| < 1, periodicity of the trigonometric functions also allows for string centers to be located on the line iπ/2, resulting in negative parity strings (ν j = −1). This type of strings will be left out of consideration for the analysis of scattering magnons, restricting to ν j = 1.
At finite size, solutions are not exactly given by strings, but rather contain string deformations δ (n) j,a ∈ C, under the constraint that the full set of rapidities {λ} remains self-conjugate. In the cases considered here the deviations are exponentially small in system size and it is therefore sufficient to take the limit of vanishing deviations. In this limit, the product of the Bethe equations of all rapidities within a string reduces to the Bethe-Gaudin-Takahashi equations, 43 which are similar to the Bethe equations but given in terms of the n-string centers λ (n) j . In logarithmic form they read where M n denotes the number of n-strings present, satisfying n nM n = M . The logarithmic scattering kernels θ n (λ) are defined in Tab. I and the scattering phase shift between two individual strings of arbitrary length is given by (for |∆| < 1 we consider only strings with positive parity) The logarithmic form of the Bethe-Gaudin-Takahashi equations allows for the introduction of string quantum numbers I (n) j , obeying an exclusion principle for all Bethe states, meaning that every Bethe state is characterized by a unique set of string quantum numbers. By considering the limit of sending a string center to infinity, the maximum allowed string quantum number can be derived. These limiting quantum numbers define the dimensions of sub-sectors of the Hilbert space containing a specific string content.
Moreover, in the planar case |∆| < 1, the existence of strings with a specific length n in the spectrum is determined by the anisotropy. 43 Therefore restrictions on the availability of n-string wave packets as well as on their momenta 44 are present in the planar case |∆| < 1, depending on the value of ∆.
We get all Bethe states with the desired string content by solving the Bethe-Gaudin-Takahashi equations (8) using an iterative algorithm for all combinations of allowed string quantum numbers. After obtaining the rapidities, the energy of a Bethe state containing strings is easily computed as where θ n is the derivative of θ n . The energy contribution E (n) (p) of a string of length n to the energy E {λ} is where n = cos(nζ), 1, cosh(nζ) for |∆| < 1, ∆ = 1, and ∆ > 1, respectively. The momentum p = p (n) (λ) of an n-string with string center λ is given by In the following we use the convention −π < p (n) ≤ π.
Note that for a single n-string the minimum of the energy dispersion is always at λ (n) = 0, i.e. at momentum p (n) = π. The total momentum of a Bethe state can be extracted from its string quantum numbers, In the thermodynamic limit and for a finite number M of rapidities, each separate bound magnon represented by a string quantum number I (n) j can be associated to a single particle momentum p To summarise, the rapiditites belonging to each eigenstate are obtained by iteratively solving the Bethe-Gaudin-Takahashi equations (8). They can be used to evaluate the determinant expressions (18) for the normalised matrix elements of the following section.

C. Magnetization expectation value
Time evolution of the expectation value of the local magnetization S z j (t) is performed by making use of the algebraic Bethe ansatz. 6 The time dependent wave function is computed using unitary time evolution in a basis of Bethe states |{λ} , where the coefficients C {λ} = {λ}|Ψ(0) are determined by the initial state |Ψ(0) , which is given in Sec. III, Eq. (23), for the construction of n-string wave packets.
The expectation value of the local magnetization at site j is given by (16) As the states |{λ} , |{µ} are Bethe states determined respectively by sets of rapidities {λ j } M j=1 , {µ j } M j=1 that obey Bethe Eqs (3) (here we have M rapidities in both sets since the operator S z j does not change the magnetization), the matrix elements are given by the normalised expressions obtained from algebraic Bethe ansatz Here we make use of the determinant expressions obtained in Ref. 42, The normalization N ({λ}) is computed from the Gaudin determinant, 45,46 where the Gaudin matrix is given by the Jacobian of the Bethe equations, The time-dependent expectation value of the local magnetization is then obtained by evaluating the double sum over matrix elements in Eq. (16). By construction, the double sum only includes eigenstates with the same particle content, which is not large for the few-magnon states we will consider. As a result, the double summation is still tractable at lattice sites N ∼ O(10 2 ). In the case of dealing with string solutions for the magnon bound states, reduced determinants for strings described in Ref. 47 must be used.

III. BOUND MAGNON WAVE PACKETS
The strings described in the previous section do not correspond one to one to localized bound states of downturned spins, but rather are translationally invariant constituents of Bethe eigenstates. In order to create states of n bound magnons with localized magnetization features, we construct Gaussian wave packets by summing over single n-string states (labeled by the string center λ (n) ) with momenta distributed around p, where N 0 is a normalization constant. Unitary time evolution under Hamiltonian (1) implies an expression of the velocities of the wave packets by expanding the dispersion relation (11) around p = p to first order (see The prefactor can be expanded for large anisotropy, implying that wave packets constructed from higher strings have a lower velocity in real space.
In the remainder of this section we will analyze the stability of n-string wave packets, along with the stability of another form of a localized multi-magnon, a consecutivesite spin block. The two constructions are closely related to each other for large ∆. Subsequently, scattering processes of n-string wave packets are visualised by computing the time evolution of the local magnetization. Furthermore, we describe a method for direct observation of the phase accumulated by the wave function of a finite spin chain during a scattering process.

A. Stability of n-string wave packets
The center of the wave packet x(t) and its width ∆x(t) are respectively given by the expectation value and the variance of the position operator In the continuum limit, the sum over all possible n-string momenta p = p (n) can be approximated by an integral. For the width of the Gaussian wave packet in real space we obtain String solutions, being associated with bound states of magnons with exponentially decaying wave functions, will add exponential terms to the shape of the n-string wave packets in real space. By inserting the complex momenta of the individual constituents of the string solutions to the Bethe wave function, it can be shown that for example a 2-string state with momentum π/2 contains exponentials in the Bethe wave function reading e −(x2−x1)/ξ , with an effective binding length ξ = 2/ ln(2∆ 2 ). This average distance between the constituent particles in the string states provides a lower bound to how localized the wave packets constructed out of such string states can be.
For α < ξ(∆), the wave packet will start to lose its Gaussian shape in real space in favour of a simple exponential decay around the wave packet center. This issue can be circumvented by choosing large enough Gaussian widths, but will require much higher system sizes.
The effective binding length becomes smaller at higher anisotropy, making the problem of the extra exponentially decaying shape to the Gaussian wave packet construction relevant only at |∆| 1.
Induced by the nonlinear dispersion relation of the magnons, the width of a wave packet in real space will furthermore increase in the course of time (see Eq. (A2)), obtained by expanding the dispersion relation (Eq. (11)) to second order around the average momentum, yield- The broadening of n-string wave packets in real space is therefore described by the initial width α and To first order, all strings of arbitrary length are stable at momenta p = ±π/2, but possess a non-trivial dependence on the anisotropy for other momenta. Moreover, the broadening of 1-string wave packets is not influenced by the anisotropy. The anisotropy dependent factor can again be approximated in the large anisotropy limit given by Eq. (25), showing that the stability of the wave packets increases with increasing anisotropy and string length. In Fig. 1 the magnetization profile of a diffusing 2string wave packet with zero group velocity computed from algebraic Bethe ansatz is shown. We used as average momentum p = π since the energy dispersion relation has its minimum there. Furthermore, fitted parameters on the time-dependent wave packet width are compared to the theoretical values of δ 2 n for 2-and 3-strings. The diverging behaviour at low anisotropy of both δ 2 n and the effective binding length of strings constrain the applicability of scattering theory results for the planar regime |∆| < 1.

B. Stability of spin blocks
A similar or even more drastic dispersing behaviour can be observed for blocks of n adjacent sites. Fig. 2 shows the time evolution of a block of 20 upturned spins in a ferromagnetic chain of downturned spins for different anisotropy parameters ∆. The results shown in Figs. 2 and 6 were obtained using the time evolving block decimation (TEBD) algorithm. 16,19,48 In the axial regime ∆ > 1, the initial state mostly projects onto many-magnon bound states, namely nstrings with n ≤ 20, which for such anisotropy values are tightly bound and thus have a large overlap with the initial spin block. A quantitative analysis of the overlap between the initial spin block and large strings is provided by Refs. 49 and 50 for a comparable situation involving Since the spin blocks mostly contain large strings for large anisotropy, they display slow dispersion (see Eq. (11)), meaning that the initial spin block remains more or less intact in time over long time scales. However, as the isotropic point ∆ = 1 is approached, the nature of the overlaps drastically changes. Eigenstates with combinations of smaller strings start carrying a larger fraction of the total overlap with the initial state. Under time evolution, one thus sees spacetime propagation lines corresponding to shorter strings, which disperse more rapidly. By the time one has entered the planar regime, 0 < ∆ < 1, the initial spin block decomposes into all available string lengths including the most rapidly-dispersing 1-string states, leading to a rapid dispersion of the magnetization throughout the 'light cone' defined by the maximal group velocity of the 1-strings.
A comparison of Figs 1 and 2 (both at p = π) shows that spin blocks decay faster than n-string wave packets of the previous subsection, which has a twofold explanation. First, an n-string wave packet is a superposition of n-string Bethe states which have no decay channel into strings of smaller lengths. In contrast, in the superposition of the initial spin block state all smaller string lengths are allowed and actually present. Second, in the momentum distribution of the spin block, momenta belonging to high velocities are not Gaussian-like suppressed. Therefore, states with high velocities (p ≈ π/2) are more dominant in the spin block state than in the n-string wave packet both with p = π. However, the spin-block state and the n-string wave packets become similar at higher ∆, such that the analytic predictions on scattering displacement becomes applicable to both cases.

C. Scattering n-string wave packets
The pre-scattering bound magnon initial states |Ψ(0) are composed of two Gaussian wave packets, where the construction relies on allocation of individual momenta to distinct strings within a single Bethe state according to Eq. (14). We localize two Gaussian wave packets labeled by j = 1, 2 with average momenta p j around two wellseparated lattice sites where For simplicity, we label Bethe states here only by the two string centers λ (n) and λ (m) instead of the whole The two centers are respectively uniquely determined by the momenta p 1 and p 2 via Eq. (12). Due to the energy dispersion (11) of the bound magnons, the relative velocity of the wave packets is maximized at p 1 = −p 2 = −π/2. Fig. 3 shows time-dependent magnetization profiles for scattering of 1-and 2-string wave packets respectively, computed from the algebraic Bethe ansatz matrix elements described in Sec. II C. A close examination of the profiles shows distinctive features akin to soliton scattering, namely that the wave packets emerge out of a collision intact, but spatially displaced. This displacement will be quantified in the next section.

D. Direct phase shift measurements
The real-time scattering trajectories in this work are analyzed using the idea that the scattering between two localized wave packets is well-described by the scattering phase shift corresponding to the average momenta of the two wave packets. This idea can also be verified through direct observation of the phase accumulated by the wave function of a finite chain during a scattering process. We therefore compare the evolution of an interacting chain (∆ = 0) with the evolution of a noninteracting chain (∆ = 0), each containing two localized 1-string wave packets. The overlap between the two wave functions, gives the phase acquired due to the interaction between the magnons. The value of ∆ is here indicated in the superscript.
In Appendix C we show how the phase of the quantity Γ ∆ after a single scattering event, calculated using numerical exact diagonalization, matches the Bethe ansatz phase shift of Eq. (6), even for wave packets that are spatially well localized. Such overlaps between time-evolved wave functions are currently not directly accessible by Bethe ansatz.

IV. SCATTERING DISPLACEMENT
The initial state constructed from Bethe states is prepared as two wave packets of bound states of an arbitrary finite number of magnons with initial average positions and momenta (x j , p j ). The m-and n-string wave packets are constructed separately at large separation, such that their motion before and after scattering can be considered to be free. In particular, the motion of the center of each wave packet is (see Appendix A) x j (t) = x j + v j t before scattering, x j + v j t − χ j (p 1 , p 2 ) after scattering, (33) given in units of lattice distance, where the velocity v j was defined in Eq. (24). Note that, in the case of single magnon scattering for example, a negative average momentum p j yields a positive velocity v j and vice versa (see Eq. (5) with J positive).
The displacement χ j (p 1 , p 2 ) can be obtained by expanding the scattering phase shift around the average momenta, see also Eq. (B11), and is therefore given as where we introduced the notation with subscript j to refer to the momentum derivative χ j = ∂χ ∂pj of the scattering phase shift χ.
Eq. (33) assumes that all scatterings occur without particle production, which is the case for the integrable model we are dealing with.

A. Displacements from Bethe ansatz
An analytic expression for the displacement as a function of anisotropy and incoming momenta can be extracted from the Bethe ansatz scattering phase. The phase shift of two bound magnons of arbitrary length is obtained from the scattering kernel Θ nm of the Bethe-Gaudin-Takahashi Eqs (8), which consists of a sum over the functions θ s , s = |n−m|, . . . , n+m, defined in Tab. I, see also Eq. (9).
For the scattering of an n-string wave packet labeled by 1 with an m-string wave packet labeled by 2, the displacements on the trajectories of the n-and m-string wave packets (j = 1, 2 respectively) are given as a function of momenta as The n-and m-strings with centers λ (n) and λ (m) carry momenta p 1 and p 2 , respectively. The wave packets are located such that x 1 x 2 and v 1 > v 2 . We first discuss the planar case |∆| = | cos(ζ)| < 1. By inverting Eq. (12), we express the string center of an n-string in terms of its momentum as λ (n) (p) = atanh tan nζ 2 tan π − p 2 .
Due to Eq. (9), expression (35) for the displacement consists of a sum of the momentum derivative of the functions θ s of which the individual terms are computed as sin(nζ) cos(nζ) − cos(p 1 ) .
For the scattering phase shift of two 1-string wave packets with p 1 = −p 2 = −π/2, the displacement is given as a function of anisotropy as The latter result could have been also obtained directly by taking the derivative of Eq. (6). Similarly, the scattering displacement of two 2-string wave packets becomes The validity of the latter equation only extends to the region where ∆ > 1/ √ 2, as 2-strings with momentum p = ±π/2 do not exist for lower anisotropy 16,44 , which can be shown from the anisotropy dependent maximum string quantum numbers.
For the regime ∆ > 1, Eqs (38) and (39) hold as well, since the θ s (λ) for both regimes are just rotated in the complex rapidity plane with respect to each other. Starting from θ s (λ) for ∆ > 1 with ζ = acosh(∆) therefore yields identical results for the scattering displacements.
The displacements for the scattering of an 1-string wave packet at a 3-string wave packet is given by

B. Comparison of scattering theory and time evolution
The displacement in the trajectories induced by scattering effects is easily deduced from the time evolution data of Fig. 3. For symmetric cases with identical particles, the average position of the magnetization of a single wave packet can be computed on one half of the system as a function of time, The result for the scattering between two 2-string wave packets is plotted in Fig. 4. Note that the plotted average location of the wave packet in Fig. 4 Fig. 4 are specifically given by Eqs (38) and (39) as function of anisotropy. The results are plotted in Fig. 4 as well, showing agreement between both approaches. The explicit time evolution relying on algebraic Bethe ansatz matrix elements provides confirmation of the analytical predictions for the displacements.
Besides symmetric scattering situations, colliding distinct n-string and m-string wave packets can be constructed and traced in the time-evolved magnetization profile. Fig. 5 shows the scattering between an 1-and a 3-string wave packet. Two situations are distinguished, the former with both wave packets at maximal velocity at p 1 = −p 2 = −π/2, where the larger string moves much slower because of its effective mass, see Eq. (25). The latter situation consists of an incoming 1-string wave packet scattering on a stationary wall of a 3-string wave packet at p 2 = π. The corresponding analytic scattering displacements, Eqs (40)- (43), are shown adjacent to the time evolution plots in Fig. 5. The scattering displacements are measured from the time evolution by imposing a Gaussian fit on the wave packet after scattering and comparing the average location to the time evolution of a single wave packet without scattering. In the lower right panel (p 2 = π), the fitting procedure of the average location of the 3-string wave packet becomes less accurate for decreasing ∆. In the planar regime where |∆| < 1, we encounter substantial limitations (as described in Sec. III B) on both the construction of the scattering wave packets, as well as on the comparison to results of scattering theory. Due to the effective binding length of the individual constituents of the string states, the tails of the magnetization profile of the wave packets start overlapping with each other significantly at low ∆, invalidating important assumptions of scattering theory which include asymptotic separation of the wave packets before and after scattering. Comparison of the measured displacements to scattering theory is therefore not meaningful for higher strings in the planar regime. Only going to much larger system sizes would resolve the aforementioned issue.

C. Displacements from TEBD calculations for scattering of spin blocks compared to the Ising limit
The scattering displacement turns out to have a particularly simple form at large anisotropy. In particular, it was found in Ref. 19 that when a propagating n-particle cluster is incident on a larger block, the block is displaced by 2n sites. This can be explained from the Bethe ansatz results presented in previous subsections, by taking the Ising limit ∆ → ∞ of Eq. (35) for the displacement of an n-string wave packet scattering at an m-string wave packet.
First, we obtain for all s which is independent of n, m, p 1 , and p 2 . Using Eq. (9) for the phase shift between an n-string and an m-string eventually yields lim ∆→∞ χ (n,m) 1 (p 1 , p 2 ) = −2 min(n, m) + δ nm . Thus, the scattering displacement for unequal wave packets (n = m) is equal to twice the number of particles in the smaller wave packet.
The leading order term of the Ising limit in Eq. (46) can be given for the case where cos(p 1 − p 2 ) = 0 as yielding an error estimate for Eq. (46). A systematic expansion for all momenta becomes cumbersome due to the summation in Eq. (9) and is left out of consideration. If cos(p 1 − p 2 ) = 0, the leading contribution is given by the third line of (47). If further cos(p 1 ) = 0, the leading order terms will be formed by higher powers like ∆ −2n for n = m or ∆ −2|n−m| for m = n. Fig. 6 shows TEBD results of the scattering of a twospin block excitation with a block consisting of 10 adjacent spins at ∆ = 5.0. The initial 2-spin block was created by upturning two neighboring spins at lattice sites 2 and 3. The first site (with open boundary conditions) is energetically inaccessible for large ∆; hence the 2-spin block travels to the right and is incident on the 10-spin block. The displacement is clearly by 4 sites, as predicted by Eq. (46). A similar displacement by 2 sites in the case of a single incident particle was highlighted in Ref. 19. Although the 2-spin block and the 10-spin block are not explicity prepared as wave packets in this case, at large ∆, string states are tightly bound, and therefore these blocks may be understood intuitively to be close to 2-string wave packets and 10-string wave packets. The Ising limit Eq. (46) thus provides a satisfactory explanation to the shift observed in these numerical experiments.
Even more, it was observed 19 that the behaviour of the scattering displacement of a large spin-block state close to the isotropic point at ∆ = 1.1 still resembles the scattering behaviour of the large ∆ limit by shifting the spin-block by two lattice sites upon scattering. In order to explain this from Bethe ansatz, we take the limit of Eq. (35) for the scattering of a large string (m 1) with center λ (m) (p 2 ) and an 1-string with center λ (1) (p 1 ) at all values of ∆ > 1, where s can only take on the values s = m − 1 or s = m + 1, due to Eq. (9). The latter equation finally gives for the displacement of the large m-string at ∆ > 1, implying that already at ∆ = 1.1 (ζ = 0.4436), the scattering displacement of a large spin block should still be close to two sites, as is in agreement with the aforementioned observation.

V. CONCLUSIONS
In this work, we have studied the quantum analogue of soliton-like scattering phenomena in the anisotropic spin-1/2 Heisenberg chain, by utilising the algebraic Bethe ansatz. We considered quantum scattering of localized excitations, created from linear combinations of Bethe states with Gaussian-distributed momenta, constructing wave packets of n bound magnons. This construction allows to study scattering phenomena of wave packets containing an arbitrary number of bound magnons.
Exact methods based on the algebraic Bethe ansatz provide a framework to evaluate the time-dependent expectation value of the local magnetization profile, which allows for a spatial tracking of the localized excitations. This explicit unitary time evolution of the initial state relies on the availability of determinant expressions for matrix elements of local spin operators.
The algebraic Bethe ansatz time evolution of colliding wave packets of bound magnons displays a spatial displacement in the trajectories of the wave packets under scattering, consistent with scattering theory results. For different values of anisotropy, fits on the displacements of the time evolved trajectories are in agreement with analytical results on the displacement from the derivative of the Bethe ansatz scattering phase shifts, for several combinations of string lengths.
The scattering phase shift can also be measured directly as well for the scattering between two localized single-magnon wave packets, again matching phase shift expressions provided by Bethe ansatz. Using TEBD, scattering displacements from spin-block states have been studied, showing similar scattering features and validating the analytic predictions of the Ising limit for the scattering displacement.
The experimental realizability of real time tracking of localized excitations in the Heisenberg spin chain 17,22 might provide an opportunity to study dynamical scattering phenomena of (bound) magnons. A possible manifestation of such phenomena might be provided by the soliton-like scattering effects analysed in this work.
The results on the scattering displacements can be extended to other Bethe ansatz solvable models. Finally, the time evolution method relying on matrix element expressions from algebraic Bethe ansatz can be used to construct other initial states in spin chains as well and to study their respective relaxation phenomena.
In this Appendix, we review some general results from quantum scattering theory, emphasizing the direct connection between physically observable quantities and (derivatives of) the scattering phase shift. 51,52 We consider an initial state with two well-separated quasiparticles (e.g. magnons and magnon bound states) with almost well-defined positions (x 1 , x 2 ) and momenta (p 1 , p 2 ). In the asymptotic region, i.e. when the distance between the two wave packets is much larger than the radius of the interaction, the time evolution is free. Hence, the centers of the wave packets x j (t) = x j (t) translate rigidly. The scattering between the two particles however introduces a displacement proportional to the derivative of the scattering phase χ(p 1 , p 2 ). More precisely, in the asymptotic regions, the motion of the center of each wave packet is after scattering, (A1) where χ j = ∂χ ∂pj is the displacement, while the group velocity v j is given by the derivative of the dispersion relation, i.e. v j = ∂Ej ∂pj | pj =p j . Similarly, the scattering has an effect on the width of each wave packet ∆x 2 For Gaussian wave packets we have to first order . An analogous formula holds for ∆x 2 2 (t) as well. Scattering also builds up correlations between the (initially uncorrelated) Gaussian wave packets, as can be seen from the time evolution of the correlator All the aforementioned quantities carry information about the derivatives of the scattering phase and can be in principle measured in a scattering experiment. These results hold for any one dimensional theory with stable particles. For the XXZ spin chain the stability of magnon bound states above the binding energy threshold is preserved by the integrability of the theory. In this Appendix, we review some general results for the scattering of two particles (magnons, magnon bound states, etc.), 51,52 and derive Eqs (A1)-(A3). For simplicity, we consider the scattering of two distinguishable particles in a continuum integrable model. The same results can be obtained for identical particles.
So, let us consider two particles with asymptotic momenta p 1 and p 2 and different dispersion relations E i (p i ), i = 1, 2 in an infinite volume (zero density). The statement of the coordinate Bethe ansatz, see Eq. (2), is that in the asymptotic region where the two particles are very far apart the eigenfunctions of the system are plane waves, where S(p 1 , p 2 ) = −e iχ(p1,p2) is the scattering matrix and χ(p 1 , p 2 ) the scattering phase shift. At zero density, the energy is simply Let us briefly comment on the structure of the wave function (B1). First of all, as we discussed in the main text for the XXZ model, bound states are characterized by complex conjugate rapidities, which leads to exponentially decaying terms in the Bethe wave function (2) with respect to the relative coordinate. This is a feature of bound states that are characterized by a center of mass coordinate. In what follows, we do not denote these exponentially decaying terms and label the bound states only with the position of the center of mass. For elementary particles, Eq. (B1) is a consequence of the conservation of energy and momentum in one dimension, and as such it is valid for any model with a sufficiently shortrange potential. Instead, if one of these particles is not elementary but it is a bound state, then for a general theory the previous simple form of the wave function is not true anymore. The bound state can decay and scattering be diffractive. However, there exist models (such as the XXZ spin chain) for which the scattering is always non diffractive. Hence, the coordinate Bethe ansatz describes a complete set of asymptotic eigenfunctions, as thoughtfully discussed in Sutherland's book. 44 For such theories, bound states cannot decay, but are protected by integrability. Let us consider the scattering problem. At time t = 0 the two particles are far apart and have (almost) welldefined positions x j and momenta p j , j = 1, 2. Without loss of generality, we may assume that x 1 When they are far apart, the two wave packets move with velocities v j . If v 1 < v 2 , the evolution is always free, while for v 1 > v 2 at some time the two particles become close and the interaction plays a role. The time evolution of a two-body wave function is given by ϕ p1,p2 is a complete set of eigenfunctions, and ψ(x 1 , x 2 ) is the initial wave function. Now, we consider an initial state which is factorized, Since we are interested in computing also the spreading in time of the wave packet, we assume that the functions ψ j (x j ) are Gaussian wave packets with minimal indetermination, where ∆x 2 j = dx j (x j − x j ) 2 |ψ j (x j )| 2 = α 2 j /4. This assumption may be dropped if we are interested only in the displacement, Eq. (A1).
Up to now, everything is exact. However, we can make the following approximations: • Since the initial state ψ(x 1 , x 2 ) is sharply peaked around x 1 and x 2 with x 1 x 2 , the only relevant contributions to (B5) come from the regions around these two points. We can then use the asymptotic formula (B1) for ϕ p1,p2 (x 1 , x 2 ) with x 1 x 2 . Inserting Eq. (B1) into (B5), and using definition (B7) yields • Since the functionsψ j (p j ), j = 1, 2, in (B9) are peaked around p j ≈ p j , we can expand the dispersion relations E j (p j ) around these momenta in the integrand of Eq. (B4). Moreover, for x 1 x 2 we may substitute the asymptotic expression (B1) for ϕ p1,p2 into (B4) and similarly expand the scattering phase shift χ(p 1 , p 2 ). Up to second order these expansions read with δ j = ∂ 2 Ej ∂p 2 j pj =p j , and χ(p 1 , p 2 ) = χ(p 1 , p 2 ) + χ i (p 1 , p 2 )(p i − p i ) with χ i = ∂χ ∂pi and χ ij = ∂ 2 χ ∂pi∂pj , and we sum over repeated indices. These expansions up to the second order are physically meaningful only if we can ignore the distortion of the wave packet, and so they are no more valid for times long enough.
Taking advantage of these approximations, we are now in the position to derive Eqs (A1)-(A3). Before scattering, the two wave packets propagate freely. Since they start around x 1 and x 2 with x 1 x 2 , they are centered around x j (t) with x 1 (t) x 2 (t) for small times. Thus, we can use the asymptotic formula of ϕ p1,p2 (x 1 , x 2 ) valid for x 1 x 2 in Eq. (B1), and hence perform the Gaussian integrations thus obtaining the pre-scattering results, Eqs (A1)-(A3). Similarly, long after the scattering, we have x 1 (t) x 2 (t). Thus, taking advantage of the proper asymptotic formula (B1) for the eigenfunctions and the expansion (B11) of the scattering phase shift we also obtain the post-scattering formulas (A1)-(A3).
Appendix C: Obtaining the phase shift from the phase of the wave function of the full chain A prominent theme of this work is that real-time scattering between two localized wave packets is welldescribed by the phase shift corresponding to the average momenta of the two wave packets. In this Appendix we connect real-time scattering data for two single-magnon (1-string) wave packets in a finite chain directly to the definition of the scattering phase shift, namely, that the scattering phase shift is the phase picked up by the system wave function during the scattering process. To demonstrate the robustness of this idea, we show results from a stringent situation of rather small wave packets (width of a few sites) in a rather small chain (≈ 25 sites), far from the usual idealized limit of infinitely extended excitations.
Here the Hamiltonian of the XXZ model is used in the following form: (C1) The extra terms compared to Eq. (1) are convenient for considering the phase of the time-evolving wave function in the two-magnon sector. When evolving with the Hamiltonian (1), there is a constant accumulation of ∆dependent phase in the time evolution even when the magnons are spatially separated; this is avoided with the above form. With this form, the interaction affects the phase of the chain wave function only when the magnons meet each other.
In the initial state, each magnon wave packet is prepared as a Gaussian, localized approximately around L/4 and 3L/4 respectively, with opposite momenta ±k. The initial state is thus at some appropriately chosen 'final' time t = T . Evolution with the ∆ = 0 (non-interacting) Hamiltonian gives a 'reference' state. We are interested in the phase accumulated in the evolution with the ∆ = 0 Hamiltonian in comparison with the reference state, namely the phase of Γ ∆ . The final time T is chosen such that the particles have completed their scattering, but have not reached the boundaries of the chain. So, there are no edge effects. For example, in the process shown in Fig. 7(left), it would be reasonable to compare overlaps at T ∼13J −1 .
(C6) Fig. 7(right) compares the phases obtained from the time-evolving wave function of a finite chain containing two relatively narrow wave packets, with the Bethe ansatz phase shift expressions which are strictly valid for delocalized excitations. Even with wave packets as narrow as α ≈ 3.5, the Bethe ansatz expressions match extremely well the phase acquired in real-time evolution.