Real-time scattering of interacting quasiparticles in quantum spin chains

We develop a method based on tensor networks to create localized single particle excitations on top of strongly-correlated quantum spin chains. In analogy to the problem of creating localized Wannier modes, this is achieved by optimizing the gauge freedom of momentum excitations on top of matrix product states. The corresponding wavepackets propagate almost dispersionless. The time-dependent variational principle is used to scatter two such wavepackets, and we extract the phase shift from the collision data. We also study reflection and transmission coefficients of a wavepacket impinging on an impurity.

In quantum many-body systems, it is often possible to understand various phenomena such as transport of charges or correlation spreading in terms of their effective quasiparticles 1,2 . In integrable spin chains this idea has been made precise by the recent development of generalized hydrodynamics [3][4][5] , which fully describes the largescale dynamics of a system only in terms of the semiclassical motion of stable quasiparticles. This idea can also in principle work for non-integrable systems, where local quasiparticles are usually thought of as almost stable excitations with long lifetimes that can be treated semi-classically 6,7 . However, there is much less information available about the dynamics of local excitations in generic interacting systems and how to extract them from microscopics.
As a first step towards a generic framework for capturing quasiparticles in one-dimensional quantum system, a variational quasiparticle ansatz was introduced 8 based on the framework of matrix product states [9][10][11] . Inspired by the single mode approximation 12,13 , the ansatz effectively boils down to a plane wave of a local perturbation on top of a strongly-correlated ground state. This approach is very suggestive of a quasiparticle running through the system, but differs from the traditional quasiparticle concept since these states are exact eigenstates (up to numerical errors) with infinite lifetimes with respect to the fully-interacting hamiltonian 14 . In a natural next step, two-particle wavefunctions were constructed and a generic definition of the two-particle S matrix was proposed 15,16 .
In this paper we present a method that allows us to gain insight in how these quasiparticle excitations propagate and interact with each other in real space and time. While in some past works this was done for integrable systems on top of product states 17,18 , our method is used to study particle-particle interactions in both integrable and non-integrable systems on the correlated ground state, and is therefore truly capturing the low-energy degrees of freedom for strongly-interacting systems. We also study the scattering of a particle off an impurity and how this impedes the transport of energy. Methodology. As a first step, we approximate the ground state in the thermodynamic limit as a uniform matrix product state 19,20 (MPS), which we can represent as where every node corresponds to a tensor and every connected node is a contraction over the corresponding index. This class of states has a number of desirable properties: the ground state of generic quantum spin chains can be faithfully represented as an MPS 21,22 , an MPS is numerically efficient to work with 9-11 and because the tensor is identical on every site, A defines a manifold of translation-invariant states. Approximating the ground state can be done using gradient descent or more advanced methods 20 such as the vumps 23 algorithm. Excitations on this ground state can be found using the MPS quasiparticle ansatz 8 Because the tensor B acts on the virtual level of the MPS, it can capture the deformation of the ground state over an extended region, and therefore represents a dressed object on top of a correlated ground state. For a given momentum k the problem of minimizing the energy can be reduced to an eigenvalue problem for the tensor B, which can be efficiently solved using iterative methods. It was shown 8,24,25 that this approach yields quasi-exact results for the dispersion relation of generic spin chains. We can localize a wavepacket in real space using the momentum-resolved quasiparticle states via a procedure reminiscent to the construction of Wannier orbitals 26 . A real-space localized wavepacket is given by where the real-space tensor B n is arXiv:1907.02474v1 [cond-mat.str-el] 4 Jul 2019 The determination of B n is ambiguous in three ways. First, since the tensors B k are found as a solution of an eigenvalue problem, they are only defined up to a phase. In addition, there are additional k-dependent gauge redundancies in the tensor B k , since the transformation gives rise to the same wavefunction |Φ k (B k ) . A third difficulty concerns the fact that, in a practical computer simulation, we only have access to a finite set of momenta, such that the real-space tensor contains an infinite number of shifted copies of the wavepacket localized in real space. First, we impose a condition on the gauge freedom for the B k such that it is as local and symmetric as possible; a natural choice is to minimize the matrix norm of tensors in the center gauge. Next, we find the sampling function f (k) such that B n is confined to the smallest possible region in space, i.e. we minimize the region Λ such that we can find an f (k) for which |B n | < , ∀n / ∈ Λ.
Because we have shifted copies, we make sure this confinement region is smaller than the distance between the copies. This allows us to find an optimal set of B n , for some unknown function f (k). Another way of looking at this is as an interpolation method, trying to find a smooth function f (k)B k through the exactly known points. Because of the localization within a region Λ, we can safely truncate the shifted copies by convoluting with a filter function.
We can now prepare wavepackets centered around a given momentum k 0 by multiplying the sampling function f with an additional gaussian distribution, In this way we can prepare a state with two wavepackets centered around different momenta at different locations in real space. In order to actually simulate the collision in real time, we represent this state as a finite string of tensors V i surrounded by the uniform MPS tensors, and simulate the real-time evolution on this finite window using the time-dependent variational principle (TDVP) for site-dependent MPS 27 . The time evolution can be fully captured by changing V i , so long as the wavepackets propagate within the window 28,29 . The Ising chain. We will illustrate our method by demonstrating how to simulate the scattering process between two wavepackets in an infinite Ising spin chain defined by the hamiltonian In the fully polarized case (λ → ∞) the excitations are spin flips, and at finite λ these become dressed by the quantum fluctuations in the ground state. We first illustrate the construction of the wavepacket. In Fig. 1 we have plotted the norm of the tensor B n that was optimized as described before, showing that we find an infinite number of shifted localized wavepackets within a certain confinement region. In the same figure, we plot the energy density of the wave function, showing that we have indeed found nicely localized packets of energy -i.e., particles -on top of the correlated ground state. This form of B n is then truncated with a filter function to keep a single wavepacket.
In Fig. 2 we display the time evolution of two wavepackets that were prepared around momenta k 0 = ±1.5 with a spreading of v = 0.5. From this figure we can see that the two wavepackets nicely propagate through the chain -there is a small dispersion effect, because of the finite width v of the gaussian -until they collide and interact. After the collision we find two outgoing wavepackets. We observe that no position shift in the wavepackets' trajectory compared to the freely propagating case, which reflects the fact that particles in the Ising chain behave as free fermions. To confirm this quantita- tively, we subtracted away what we would expect during free propagation and find an energy density difference of at most 5 * 10 −4 , in agreement with the expected TDVP error of order dt 3 . Spinon-spinon scattering. As a first application, we consider the XXZ hamiltonian which is integrable 30,31 , and for ∆ > 1 the elementary excitations are massive particles with fractional spin s = 1/2 32,33 . The quasiparticle ansatz readily generalizes to this case 8 , and we find a spinon dispersion relation that agrees with the exact result to arbitrary precision. Similarly as before, we create two wavepackets with parameters v = 0.3 and k 0 = ±0.7. Similarly as in the Ising case, we find two well-defined outgoing modes, a consequence of the fact that particle scattering is purely elastic in integrable systems. In Fig. 3 we plot the location of the particles as a function of time, where the approximate position of the wavepacket was found by where i is the energy density at site i. This can then be compared to the freely propagating case (see Fig. 3) and in contrast to the Ising case, we find that the collided particles have undergone a displacement. As one can see in the inset, we find a displacement that is constant in time with a magnitude d ≈ 0.476.
Using the exact (dressed) scattering phase of excitations on top of the ground state, provided by the integrability framework, the predicted displacement at momenta for which it is well known 34,35 that the elementary excitation is a gapped magnon with spin s = 1. The magnonmagnon scattering properties are qualitatively described by the S matrix of the non-linear sigma model, and this was confirmed quantitatively in Refs. 15 and 36. Because this model is not integrable, two-particle scattering is not purely elastic, and an initial state with two particles can potentially scatter to three-or four-particle modes with the same total energy.
A two-particle state can have total spin s = 0, s = 1 or s = 2, and we can prepare two-wavepacket states within these three sectors separately. In figure 4 we show the approximate position of the wavepacket throughout time. While we don't see a well defined phase shift, probably due to non-integrable effects, we still find that the spin-0 sector leads the freely propagating wavefront, whereas the wavepackets lag in the spin-1 and 2 sector. This is in agreement with previously known results for the signs for the scattering lengths in the different sectors 15,36 . Impurity scattering. Our technique is general enough to consider the case of a wavepacket impinging on an impurity in an infinite spin chain. Consider again the spin-1 Heisenberg but with a single spin-2 site at x = 0 coupling  Figure 4. Comparison between the freely propagating (blue) and colliding wavepackets (red) trajectories in the spin−1 Heisenberg chain. The different total spin sectors spin-0 (red) spin-1 (blue) spin-2 (green) are compared to a freely propagating wavepacket (purple).

to its neighbors as
where T x,y,z are the spin-2 operators acting at the impurity site and h is a small magnetic field polarizing the impurity spin. We find the ground state by optimizing the variational energy within a finite window embedded in a uniform MPS. We have prepared a wavepacket polarized along the z direction with k 0 = 2.65 and let it impinge on the impurity. It is interesting to note that while there may exist local excitations on this impurity spin, scattering states will always be orthogonal to those bound states. Therefore, the wavepacket can only reflect on or tunnel through the impurity. By comparing the total energy density left and right from the impurity, we can calculate the transmission and reflection coefficients to be 0.4206 and 0.5794 respectively (see Fig.5). Conclusions. We demonstrated how to obtain wavepackets of quasiparticles on top of correlated ground states of generic spin chains, and simulate the collision processes in real space and time. Being quasi-stable real-space localized particles, these wavepackets are the quantum version of solitons in classical non-linear field equations. Our results for the scattering displacements are in excellent agreement with known analytical results in integrable spin chains. The absence of integrability does not change the fundamental picture as demonstrated by magnon magnong scattering in the spin-1 Heisenberg model. Finally, we can simulate the scattering of particles off an impurity in the spin chain.
Our results open up the possibility of simulating the low-energy degrees of freedom of generic spin chains and one-dimensional electrons directly in real space and time. In particular, it would be interesting to look at bound- state formation and particle confinement in spin chains, or hole dynamics in electronic chains and ladders.
In the context of integrability, we can study the effect of violations of the Yang-Baxter equation and threeparticle processes on the real-time dynamics. These results could be used to add extra dissipative terms to the generalized-hydrodynamic equations responsible for quasiparticles decay processes in perturbed integrable systems 37 and to compute transport coefficients, related to the dressed scattering shifts 38,39 , in non-integrable models.
Finally, we believe that our approach can be extended to capture stable localized modes on a representative finite temperature state, and, starting from the higherdimensional quasiparticle ansatz 40 , can be applied to scatter quasiparticles in spin and electron systems in two dimension.

I. SUPPLEMENTARY MATERIAL
We discus how to fix the gauge freedom in the quasiparticle ansatz and how to localize the wavepacket by optimizing over the prefactors in T j=1 f (k j )B kj e ixkj .

A. Fixing the gauge freedom
It is usually convenient to first gauge transform the groundstate as (1) this gauge is described in (ref vumps). A L and A R are left and right isometries, which means that the leading respectively left and right eigenvector is the identity matrix. With this gauge transform we can see that the norm of our state reduces to the norm of A C .
The excitation ansatz then becomes : We can see that |Φ k (B ) with B = B + A L X − e ik XA R would represent the same physical state as |Φ k (B) . This can be used to greatly simplify the overlap between two momentum states. Consider picking X such that The overlap between Φ k1 (B 1 )|Φ k2 (B 2 ) then reduces to T r(B † 1 B 2 ). It is this gauge transform that is used when minimizing the energy because instead of solving a generalized eigenvalue problem, the effective norm becomes the identity and we get a regular eigenvalue problem. This gauge is also asymmetric, consider the wavepacket we obtain in spin-1 heisenberg when we follow the procedure from this paper without picking a new gauge, as shown in fig.6. A simple way of obtaining a more symmetric gauge is by minimizing the matrix norm of with respect to X. This gauge is chosen in such a way that the influence of B on the virtual space is as local as possible. When the original groundstate is reflection invariant, we would therefore also expect the optimal B n to be symmetric.

B. Localizing the wavepacket
We have the signal T j=1 f (k j )B kj e ixkj (shown in fig.7) and because we uniformly sampled k j out j(2L+1)

2π
, it repeats with periodicity 2L + 1. In the slice [−L, L] we want the signal to be mostly zero except near x = 0, so that we can truncate away the copies.
To do this we can create one big vector F j = f (k j ) and write |x−L|<s | T j=1 f (k j )B kj e ixkj | 2 as F † N s F . We can then minimize this function for every s efficiently as this is an eigenvalue problem. We then picked the solution F for which s is as large as possible while the cost-function is still less then some threshold (we used 1e −5 ).

C. Spinon scattering shift
The lowest lying excitation on top of the ground state of the XXZ spin 1/2 chain at half-filling is given by the one spinon excitation 32,33 . This excitation can be derived exactly by Bethe ansatz. The logic is to construct the ground state of the chain with a sea of Bethe ansatz excitations (which are magnons) and then perform a single particle excitation on top of that 41 . The interacting nature of the system is such that all the others particles in the ground state are shifted by a fraction of the system size, causing a so-called dressing for the particle dispersion relation and scattering shift. Such a dressing effect can be computed exactly by thermodynamic Bethe ansatz 42 . We here report the main results. The (dressed) energy ε and momentum k of a single spinon is parametrised by the quasi-momentum λ ∈ [−π/2, π/2] and it is given by ε(λ) = π sinh η s(λ), k(λ) = 2π λ −π/2 s(λ), with s(λ) = ∞ n=−∞ e 2inλ cosh nη and where the parameter η is related to the XXZ anisotropy as ∆ = cosh η. therefore the group (effective) velocity of the spinon excitation as function of momentum q is given by v eff (q) = ε (λ) k (λ) λ=λ(q) .
The scattering shift K dr of a spinon excitation is similarly dressed by the non-trivial background of the ground state sea of particles. This is given by the action of the dressing linear operator on the bare scattering shift where the bare scattering shift of Bethe ansatz magnons is K(λ) = sinh η/(π(cosh η − cos(2λ))) and the distributions of Bethe roots in the ground state is given by ρ(λ) = s(λ)/π. The above equation should be intended as the action of a linear operator on a function, where 1 denote the identity operator (Dirac delta function). This can be written as a linear integral equation whose solution can be obtained via Fourier transform. We then obtain Therefore the scattering of two wave-packets with momentaq 1 andq 2 produces a displacement in the trajectory χ j of the j−th packet (j = 1, 2) given by d j = K dr (λ(q 1 ) − λ((q 2 )) k (λ(q j )) .