Collisions of false-vacuum bubble walls in a quantum spin chain

We simulate, using nonperturbative methods, the real-time dynamics of small bubbles of"false vacuum"in a quantum spin chain near criticality, where the low-energy physics is described by a relativistic (1+1)-dimensional quantum field theory. We consider bubbles whose walls are kink and antikink quasiparticle excitations, so that wall collisions are kink-antikink scattering events. To construct these bubbles in the presence of strong correlations, we extend a recently proposed matrix product state (MPS) ansatz for quasiparticle wavepackets. We simulate dynamics within a window of about 1000 spins embedded in an infinite chain at energies of up to about 5 times the mass gap. By choosing the wavepacket width and the bubble size appropriately, we avoid strong lattice effects and observe relativistic kink-antikink collisions. We use the MPS quasiparticle ansatz to detect scattering outcomes. (i) In the Ising model, with transverse and longitudinal fields, we do not observe particle production despite nonintegrability (supporting recent observations of nonthermalizing states in this model). (ii) Switching on an additional interaction, we see production of confined and unconfined particle pairs. We characterize the amount of entanglement generated as a function of energy and time and conclude that our classical simulation methods will ultimately fail as these increase. We anticipate that kink-antikink scattering in 1+1 dimensions will be an instructive benchmark problem for future quantum computers and analog quantum simulators.

It is possible that the known universe is built on top of a metastable, or "false" vacuum state.In this scenario, there is a tiny but nonzero probability of a small bubble of "true" vacuum forming via tunneling.The bubble interior has a lower energy density than its surroundings and so it expands, its walls accelerating, bulldozing everything in their path.If multiple bubbles of true vacuum form far apart, their walls will rush toward each other and eventually collide, producing showers of entangled particles.It is also possible that such events have already occurred, and are thus relevant for the evolution of the early universe [1][2][3][4][5].
Vacuum bubble walls are topological excitations and can be modelled as domain walls or kinks.Particle production is known to be suppressed in some weaklycoupled models of relativistic kink-antikink collision (see e.g.[6,7]), but much less is known about scattering of topological excitations under strong interactions, where a lot of entanglement is typically generated and semiclassical methods break down.Whether or not these phenomena are important for simulations of the early universe is an open question.Simulations that can handle such non-classical dynamics could provide an important window into these high-energy, strongly-coupled dynamical processes.
However, one does not easily simulate the dynamics of a strongly interacting quantum field theory (QFT), at least using classical computers.Quantum Monte Carlo -the workhorse for simulations of equilibrium phenomena in lattice systems (such as lattice QCD [8]) -is hard to apply efficiently to real-time dynamics, although recently some progress has been made (e.g.[9]).Tensornetwork methods also show promise, and have been used to simulate nontrivial dynamical phenomena in (1+1)D systems, such as string breaking in lattice gauge theory [10][11][12][13][14][15].Nevertheless, although the computational cost is typically linear in spatial volume, it increases exponentially with time in the general case (due to linear scaling of entanglement entropy), the exponent increasing with the number of spatial dimensions.
In principle, quantum computers (both analog and digital) can simulate dynamics of quantum field theories at long timescales with polynomial costs [16,17], a topic which has recently attracted great interest, with many digital [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32] and analog [33][34][35][36][37][38] proposals (see [39] for a review).However, existing or near-term digital quantum devices are noisy, such that only shallow quantum circuits can avoid being overwhelmed by errors.Analog quantum simulators are currently more capable, supporting longer coherence times, but also have no general mechanism for correcting errors.As such, for the time being, the physical systems that can practically be simulated are limited to simple models and we expect classical simulations to perform better, with quantum devices catching up as the hardware improves, ultimately beating classical computers by an exponential margin with the arrival of large-scale, fault-tolerant dig-ital quantum computers, permitting more complex simulations.Thus dynamical simulations of phenomena like false-vacuum collapse are physically-motivated applications for quantum computers and analog quantum simulators, and simple models can be used as benchmark problems for present and near-term quantum devices.To understand which problems make the most suitable quantum benchmarks, it is important to explore the limits of classical methods.
We address this question by developing a framework for classically simulating the full quantum dynamics of relativistic false-vacuum bubble-wall collisions in (1+1) dimensions on the lattice, using Matrix Product States [40,41] to represent the evolving state.This is computationally feasible as long as the state does not become too entangled.We demonstrate our methods on a simple lattice model with emergent relativistic, strongly-coupled false-vacuum physics, showing that they can work away from perturbative approximations.By highlighting the limitations of such classical methods, our work clarifies where quantum advantage might potentially arise in relatively near-term quantum simulators.In the following paragraphs, we outline the structure of this paper.
To improve both the interpretability and the computational efficiency of our simulations, we choose our initial states to be false-vacuum bubbles whose walls are single topological particles: a kink and an antikink.We motivate this choice in Sec.I.
In Sec.II we explain our selection of lattice model: Rather than a lattice-regularized QFT Hamiltonian (see e.g.[10-15, 17, 42, 43] for MPS simulations of these), which in the case of bosonic fields requires a truncation of an otherwise infinite Hilbert space for each lattice site, we consider a quantum spin chain, chosen and tuned so that its low-energy physics is governed by an emergent relativistic QFT.This is known to occur in the vicinity of many continuous phase transitions, where the emergent QFT is often a (by definition relativistic) conformal field theory (CFT) (see e.g.[44]).To avoid strong lattice effects, we ensure that the kink and antikink lattice velocities remain below their maximum values.
The initial state preparation for our dynamical simulations is subtle, because the kink-antikink pair is not an energy eigenstate.As discussed in Sec.III, we take care to prepare initial kink and antikink wavepackets (efficiently represented as MPS) which are broad compared to the lattice spacing, and are not contaminated by additional unwanted excitations.
We successfully simulate inelastic kink-antikink collisions in our spin-chain model, including particle production, over ∼ 1000 lattice sites at energies of up to ∼ 5m µ , where m µ is the mass of the lightest quasiparticle, as we detail in Sec.IV.We develop tools for analysis of the outgoing particles produced in inelastic kink-antikink collisions and use them to show that particle production is strongly suppressed in the Ising model with intrinsic Z 2 symmetry breaking, even though the model is noninte-grable in that case.We further show that copious particle production occurs once a Z 2 -symmetric three-site local interaction turns on.We also quantitatively track the growth of entanglement entropy during repeated kinkantikink collisions, thus inferring how large a bond dimension is needed to provide an accurate approximation to the evolving quantum state.
Sec. V contains concluding remarks, and further details of our methods and results are provided in the appendices.
Before proceeding, we note that although our work is partially motivated by a potential connection with early universe cosmology, we would need to reach energies several orders of magnitude higher than those of our present simulations, as well as increase the number of spatial dimensions, in order to study models of direct relevance to the early universe.Both of these are challenging and our work can represent, at most, a small step toward performing such simulations at strong coupling.On the other hand, we feel that our studies of inelastic scattering events in a strongly coupled relativistic quantum field theory, and of the entropy generated in such events, are of intrinsic interest apart from any cosmological motivation.

I. BUBBLE STRUCTURE
In (1+1) dimensions, a bubble of false vacuum separates two regions of true vacuum.Such a bubble may also be viewed as a confined pair of topological excitations -a kink and an antikink -with the false vacuum playing the role of a low-energy string that provides the confining force.The metastability of the false vacuum corresponds to suppression of string breaking [45].We simulate the relaxation of bubbles in which the kink and antikink are initially single, localized, spatially separated topological particles of low mass.We henceforth use "kink" and "antikink" to refer exclusively to such topological particles, and "bubble" to refer to a false-vacuum bubbles so composed.Although a general false-vacuum bubble could have a more complicated wall structure (whose dynamics one could also simulate using MPS) we restrict our initial states in this way in order to make the simulation outcome easier to interpret: If the bubble walls are topological quasiparticles, then no particles will be produced until the walls collide and the collision can be thought of as a kink-antikink scattering event with a single input channel.Such scattering events may be free (no interaction), elastic (interaction, but no particle production), or inelastic (particle production), as illustrated in Fig. 1.Particles produced may include non-topological particles, which we call mesons (following [46]) because they can be thought of as bound kinkantikink pairs in our model.These are particularly easy to observe, since outgoing mesons are not subject to a confining force and will propagate ballistically, quickly Cartoon illustrating the relaxation of a falsevacuum bubble in a spin chain.The magnetization ⟨Z⟩ is positive in the true vacuum, but negative in the false vacuum.A bubble-wall collision is a scattering process, which may be (a) free (no interaction), (b) elastic (no particle production), or (c) inelastic (particle production).Note that in free and elastic scattering of topological particles the left (right) particle always remains a kink (antikink).separating from any confined particles.Previous work on related phenomena includes MPS simulations of string breaking in the Schwinger model [11][12][13][14][15]47], where the initial state is typically prepared by applying a bare string operator to the vacuum.This generally creates excitations involving multiple particles of different energies, highly localized at the string edges, leading to relatively complex dynamics.Although such dynamics can still be simulated using standard MPS techniques, the rapid resulting entanglement growth can make it difficult to reach long times and to treat large systems.Such strings can be smeared out into wavepackets (as considered, for example in some recent work on confined excitations in spin chains [35,48]), which focuses the wavepacket momenta, significantly reducing the energy and entanglement growth.Nevertheless, the smeared string will still generally create multiple species of topological excitation, as shown in Fig. 2.
Recently, techniques have been developed [49] to con-  [51,52], both without and with a small longitudinal field h that breaks the Z2 symmetry.At h = 0, g = 1 there is a continuous (symmetry-breaking) phase transition described by the Ising CFT for λ ≲ 0.428, and by the Tri-Critical Ising (TCI) CFT at λ ≈ 0.428.By studying the spin chain near to these transitions (g → 1, h → 0), we can access emergent, relativistic quantum field theories with confined kinks [53,54].Points (i), (ii), and (iii) correspond to the data shown in the figures below.
struct wavepacket states with selective particle content in generic (1+1)-dimensional systems using MPS 1 .A main result of our paper is that we can extend those techniques to build initial bubbles whose walls are kink and antikink quasiparticle wavepackets.Apart from improving the interpretability of results, as discussed above, this further reduces the energy and the rate of entanglement growth compared to smeared string excitations, which enables us to treat larger systems and simulate for longer times.

II. SELECTING A SPIN CHAIN
We seek a spin chain whose IR physics is described by a relativistic emergent field theory supporting confined kinks.In principle, there are many suitable models: An emergent field theory of confined kink-antikink pairs can be engineered by starting with a spontaneouslybroken discrete symmetry, which provides multiple vacua and topological excitations.We then tune close to a symmetry-breaking phase transition, typically described by a CFT, and finally add a weak symmetry-breaking field to lift the vacuum degeneracy and confine the kinks.We must take care, however, since emergent field theories of such spin chains are sometimes integrable [55][56][57][58], in which case scattering, including kink collision, is always elastic.If we want to observe particle production in the emergent field theory, we must avoid integrability.
We choose an extension [51,52] of the transverse-field Ising chain where X, Z are Pauli matrices.At h = 0, this model has a Z 2 symmetry (Z j → −Z j , X j → X j ) that is spontaneously-broken when g < 1 for a large range of λ: see Fig. 3 for a phase diagram.For λ = 0, we have the transverse-field Ising chain, which already supports confined kinks for g < 1 and 0 < |h| ≪ 1 [53,59].However, for small |h| it is very close to being integrable [53,56] (both the emergent field theory and the spin chain itself are noninteracting for h = 0).Previous work has shown that kink-antikink pair states in this model can have extremely long lifetimes [60][61][62][63][64], with recent numerical studies suggesting that these long-lived states can have energies well above the threshold for inelastic scattering [45,65,66].This is despite the lack of any exact conservation law protecting these excited meson states from decay.
Turning on λ allows us to go beyond this "almostintegrable" regime, since both the spin chain itself and the emergent field theory are nonintegrable for λ > 0, g < 1, even at h = 0 [54,58,67,68].We present simulations at a point along the Ising line λ = 0, labelled (i) in Fig. 3, as well as at two points, labelled (ii) and (iii), closer to the Tri-Critical Ising (TCI) point at λ → 0.428, In the vicinity of the TCI point, (1) exhibits the same universal properties [67,68] as the Landau-Ginzburg theory with action with the parameters corresponding as g ∼ a 2 , λ ∼ a 4 , and the Z 2 -breaking field controlled by h mapping to odd powers of ϕ.

III. METHODS
We first describe how to construct the states relevant to our simulations.This includes the bubble states, consisting of a localized kink and antikink pair, which we use to initialize our simulations.We choose the kink and the antikink to be topological quasiparticles of our spin chain model.
For simplicity, we begin with the construction of states in the bare setting λ = 0, g = 0, where quantum fluctuations vanish, before moving onto the dressed setting, where we use MPS to capture, non-perturbatively, the fluctuations that appear.In both cases we first describe the true and false vacua, then the kink and antikink states, before explaining how to combine them into a bubble.We work in the Z-basis throughout: Z|↑⟩ = |↑⟩, Z|↓⟩ = −|↓⟩.

A. Bare states
Let us first consider λ = 0, g = 0, for which the terms in (1) commute and there are no quantum fluctuations (eigenstates of H can always be chosen to have definite spin orientations in the Z basis).In this bare case, the These highly-localized excitations have maximal momentum uncertainty.By smearing them out into wavepackets, e.g.j f j |κ bare j ⟩, we can make them quasilocal in both position and momentum space.We consider Gaussian packets centered at position x and momentum p, with spatial width σ.In the maximally delocalized limit σ → ∞ we obtain a momentum eigenstate with momentum p.By combining kink and antikink wavepackets we can construct a false-vacuum bubble with quasilocalized walls at positions x L and x R where x R − x L determines the size of the bubble, p L and p R specify the expected momenta of the bubble walls, and we define the localized kink-antikink pair states Note that the restriction j < k (the kink must be to the left of the antikink) means that the Gaussian packets (3) are truncated.In practice, one can ensure that this truncation is negligible by choosing x L , x R , and σ so that the coefficients are very small when j ∼ k.
In addition to kinks and kink-antikink pairs, the bare "meson" states |µ bare

B. Dressed states as Matrix Product States
While the bare states introduced above illustrate many relevant properties of the states we wish to construct for H(g > 0, λ ≥ 0), they are all eigenstates of the bare Hamiltonian H(g = 0, λ = 0), implying that kinks and antikinks do not propagate2 .They are also product states, devoid of entanglement.To obtain interesting dynamics, we need g > 0, for which all the bare states have counterparts, non-perturbatively dressed by fluctuations, possessing exponentially decaying correlations and entanglement between lattice sites.In addition to this entanglement, for g > 0 evolution by H can generate new entanglement, in contrast to the bare case.
The dressed states and their dynamics under H can be represented, non-perturbatively, using Matrix Product States (MPS) [40,41,[69][70][71], a variational class of states with the form where N is the number of spins (lattice sites), s j = ↑, ↓ for our model and each A (s) j is a D j−1 ×D j matrix, making A j a rank-3 tensor.At the ends of the chain we have MPS can also be illustrated using tensor network diagrams, which we will use in the following for convenience.For example, we can rewrite (6) as where represents a rank-3 tensor.Computations with MPS typically scale as O(D 3 ), with linear scaling in the number of spins involved in the computation.The dimensions D j , called bond dimensions, limit the amount of entanglement that can be represented.We define the cut entropy at location j S(ρ >j ) = − tr(ρ >j log 2 (ρ >j )) (8) to be the von Neumann entropy of the subsystem consisting of all sites > j.In an MPS, the cut-entropy at j is upper-bounded as S(ρ >j ) ≤ log 2 D j .
With appropriate choice of bond dimensions D j , any state of a quantum spin chain with N spins can be represented exactly as an MPS [69].Furthermore, low-energy eigenstates of locally-interacting one-dimensional lattice systems with a spectral gap, such as our model (1) at off-critical parameters, can be represented exponentially accurately in the bond dimensions D j , independently of the system size N [72].We therefore expect an MPS to accurately represent both the initial and evolved states of our simulations, provided that we use sufficiently large bond dimensions D j .If the bond dimensions are too small, this will introduce errors.
In our constructions, to avoid boundary effects, we work directly in the thermodynamic limit N → ∞, using infinite MPS (see App. A), which we illustrate as using trailing legs on the left and right to indicate that the chain of tensors extends from −∞ to ∞.The numbered sites 1 . . .N w indicate a finite window inside the bulk of the infinite chain.In Fig. 4, we use similar diagrams to indicate how the our dressed states are constructed as infinite MPS.
Vacuum and false vacuum.-TheMPS |Ω⟩ and |Ω⟩ are approximations to the dressed true and false vacua.They are uniform, infinite MPS built from tensors and , respectively.We optimize using variational methods [10,[73][74][75][76] to minimize the energy of |Ω⟩.To find the metastable false vacuum |Ω⟩, we first apply a global spin-flip to |Ω⟩, resulting in another uniform MPS, whose energy we then minimize.For sufficiently small g ≪ 1, we observe that the energy-minimization procedure does not find a path to the true ground state, resulting in a metastable false vacuum state |Ω⟩, with MPS tensor , that behaves as an energy eigenstate for all practical purposes (see App. B 2).
Kinks and antikinks.-TheMPS |κ j ⟩ approximates a dressed, localized kink state.It is constructed by introducing a new tensor that sits at position j, between two semi-infinite chains, one consisting of on the left and one of on the right.The tensor parameterizes the spatial transition between the true and false vacuum regions.Unlike in the bare case, the transition region may encompass many lattice sites, as illustrated by the ⟨Z⟩ plots in Fig. 4. The antikink |κ k ⟩ is similarly constructed by introducing a tensor between a chain of on the left and on the right.We select and using an MPS Bloch-state approach [41,74,75,77] so that |κ j ⟩ and |κ k ⟩ states can be thought of as "position bases" for the kink and antikink quasiparticles of lowest energy.We may use these states to construct topological quasiparticle wavepackets j f j |κ j ⟩ and k f k |κ k ⟩ [49].If D is the bond dimension of the vacuum MPS, such wavepackets have MPS representations with bond dimension 2D (see App. C for details).
The Bloch-state approach for finding and is conceptually simpler when h = 0, so that there is no confining force acting on the kinks and antikinks.We consider the h ̸ = 0 case further below.For h = 0, we solve an effective Hamiltonian for and such that the momentum eigenstates, j e ipj |κ j ⟩ and k e ipk |κ k ⟩, approximate the lowest-energy topological eigenstates of H with momentum p [75].The resulting tensors , (and hence |κ j ⟩, |κ k ⟩) generally depend on the momentum p, but we ignore this dependence when building wavepackets, aside from choosing the p used to solve for and to match the expectation value of momentum in the wavepacket state.This is justified for Gaussian wavepackets with large σ, and hence small momentum variance, if the tensors vary sufficiently slowly with p.
There is also an important physical reason for choosing σ to be large: In the presence of fluctuations, localized packets can no longer be truly static, since they are not eigenstates of H even for h = 0. Instead, they will spread out as time passes, at a rate dependent on σ.
Wavepackets can be made to spread slowly relative to other processes, such as the collapse of a false-vacuum bubble, by choosing σ ≫ ξ, where ξ is the correlation length in lattice units.It is desirable for the kink and antikink wavepackets comprising a bubble to spread only minimally prior to collision, since then the wavepackets of outgoing quasiparticles also tend to be well localized, which makes them easier to characterize.We now explain how to find and in case |h| > 0, where the confining force on kinks and antikinks means that j e ipj |κ j ⟩ and k e ipk |κ k ⟩ can no longer be eigenstates of H.In this case, we find and by optimizing modified energy functions that subtract away the false-vacuum contributions, which in |κ j ⟩ and |κ k ⟩ depend on the positions j and k (thus providing an accelerating force).We explain this for the case of |κ j ⟩ and , since the procedure is completely analogous for |κ k ⟩ and .We first note that it is possible to choose , by exploiting a redundancy in the representation of momentum eigenstates, to achieve ⟨κ j |κ k ⟩ = δ jk (see App. C).After making this choice, we minimize where ∆E j := ∞ j e false captures the infinite bulk contributions to the energy present in |κ j ⟩, coming from the true and false vacua (e true and e false are the energy densities of the true and false vacua).Subtracting them in this position-dependent way makes the contribution of each |κ j ⟩ term to Ẽ finite and independent of j.Note that the ∆E j correction does not affect off-diagonal terms j ̸ = k in the sum of (10).The energy minimization procedure is easily carried out by slightly adapting the methods of [74] (see App. C for details).
False vacuum bubbles.-Tobuild dressed bubble states |Ψ⟩ as MPS, we proceed analogously to the bare case by combining a kink and an antikink wavepacket where, in our simulations, we choose the momenta p L = p R = 0 and set x R − x L ≫ σ so that f j (x L )f k (x R ) is small for small k − j.The kink-antikink pair states |κκ jk ⟩ are constructed by combining the tensors , , and (already optimized to represent the vacua, kinks, and antikinks) without further modification, as illustrated in Fig. 4 (see App. C for further details).With this scheme, |κκ jk ⟩ accurately describes a kink-antikink pair at asymptotically large separations k − j.However, at small separations, corrections would generally be needed due to interaction effects3 .We again rely on x R − x L ≫ σ here, which ensures that terms with small separation are strongly suppressed, so that the error incurred by ignoring interactions is small.
Mesons and meson pairs.-Inaddition to the topological excitations |κ j ⟩ and |κ j ⟩, we can construct MPS representations |µ j ⟩ of topologically trivial "meson" states.The meson states are built from the vacuum tensor and the meson tensor as illustrated in Fig. 4, where we optimize to represent the topologically trivial particle of lowest energy using the same Bloch state approach that we use to find the kink and antikink tensors (except that no special accommodation is needed for h ̸ = 0 as mesons remain unconfined).As in the kinkantikink case, we can construct meson-pair states |µµ jk ⟩ by combining two instances of , separated by vacuum tensors (see Fig. 4).As discussed below, we can use these meson-pair states to detect the presence of outgoing mesons (even before they are visible as particle tracks distinct from an outgoing false-vacuum bubble).

C. Time evolution
To classically evolve an initial MPS |Ψ(t = 0)⟩ in time, we apply the time-dependent variational principle (TDVP) [73] within a finite window of the infinite chain surrounding the initial bubble [76].The TDVP provides effective equations of motion for the MPS tensors so that the evolution of the MPS approximates evolution by the Hamiltonian H. Assuming these equations are integrated accurately, any systematic errors come from restricting the bond dimension D, and hence the entanglement, in the MPS.If D were allowed to grow arbitrarily large, the evolution of the state could be computed exactly.In our simulations, we allow D to grow [78,79] up to a predetermined maximum value in each simulation, running the evolution several times with different limits on D in order to detect any systematic errors due to the resulting entanglement restriction.Note that we do not restrict the MPS tensors within the window surrounding the initial bubble in any other way: during the evolution the state is not constrained to the forms illustrated in Fig. 4, but is allowed to be a general MPS.
For the numerical integration of the TDVP equations of motion, we primarily use the Runge-Kutta 4/5 algorithm, which we find provides a good balance of speed and accuracy except at very early times, where we use the better-conditioned, but more computationally intensive, "split-step" integrator of [79].These methods are implemented in the evoMPS python package [80].
As the state evolves, we monitor its spin and energy expectation values as well as its entanglement properties.This allows us to draw conclusions about collision (scattering) outcomes.For instance, elastic and inelastic scattering are easily distinguished from the trivial case, as interaction generically results in entanglement between any outgoing kinks or particles 4 , whereas trivial scattering never does.We can also easily distinguish elastic and inelastic scattering in many cases.For example, if a collision produces a pair of mesons, their wavepackets will spread ballistically since two mesons are not subject to a confining force.Indeed, any sustained ballistic spread of energy implies particle production.Importantly, the converse does not always hold, since confined topological particles different from those of the initial state may also be produced.

D. Particle detection
Aside from constructing the initial state |Ψ⟩, we can also use MPS such as |κκ jk ⟩, representing kink-antikink pairs, as a kind of particle detector, the inner product ⟨κκ jk |Ψ(t)⟩ corresponding approximately to the amplitude of a kink-antikink pair with position j, k at time t.Conveniently, these states can be made to fulfill ⟨κκ jk |κκ lm ⟩ = δ jl δ kl (see App. C).We can treat the subspace spanned by these basis states as an approximate kink-antikink pair "sector", which we denote κκ.One reason for its approximate nature should be familiar from the discussion above: The basis captures a kink-antikink pair most accurately if the kink and antikink are smeared out into wavepackets that are sufficiently broad, so that the wavepacket momenta are focused around the momentum p used to compute and .Additionally, the kink and antikink must be sufficiently separated so that interaction effects are insignificant.Fortunately, these two properties can be checked after projecting the wavefunction into the |κκ jk ⟩ subspace.The simplest way to deal with terms in which the kink and antikink are too close together is to simply exclude them from the projection subspace.Inaccuracies due to the momentumdependence of and can be mitigated in a few ways: Assuming the momentum dependence is not too strong, the simplest strategy is to tune and to match the expected momentum of the projected wavefunction.A more precise result can be had via a Fourier analysis, in which the detection subspace is further restricted to a range of momenta that match and .See App.D for a more detailed, technical discussion.
The κκ subspace is already sufficient to detect inelastic scattering: if the portion of the wavefunction within the subspace drops significantly during evolution, particle production has likely occurred.Going further, we can construct quasiparticle position bases for other quasiparticle types, both topological and nontopological.To this end, we define the MPS |κ  (a) , and (a) , to be approximate position bases for the a th kink, antikink, and meson quasiparticles, with a = 0, 1, . . . in ascending order of energy.We sometimes suppress the superscript a when considering the lowest-energy quasiparticles of each type a = 0. We compute (a) , (a) , and (a) by simply solving for multiple energies in the Bloch-state approach used above to generate the lowest-energy kink and antikink tensors , [74,77] 5 .This procedure can deliver accurate quasiparticle states for quasiparticles with energy E a below the two-particle threshold [75].
Above that threshold, these tensors may correspond to unstable excitations.The procedure also guarantees that The kink, antikink and meson single-particle bases are mutually orthogonal by construction, due to the orthogonality of the true and false vacua in the thermodynamic limit.
We can construct pair states |κκ jk ⟩ from this extended set of single-quasiparticle states following Fig. 4.These extended bases are not orthonormal at small separations k − j due to interaction effects.Nevertheless, we can compute a minimum separation d for each set of Hamiltonian parameters g, ∆, h so that the bases are approximately orthonormal when k − j ≥ d (see App. C).These restricted bases give us access to extended kink-antikink κκ (a,b) and meson-pair µµ (a,b) "sectors", allowing a much finer analysis of particle content.
Note that, while it is possible to construct basis states containing k particles, including k ≥ 3, the cost of computing inner products of these basis states with the evolved state |Ψ(t)⟩ scales as O(N k w ), where N w is the size of the lattice window where the state is allowed to differ from the vacuum.For N w ∼ 1000 this makes accessing sectors with k > 2 more challenging.That said, the presence of outgoing scattering channels can still be inferred indirectly by looking for a probability deficit after accounting for all the kinematically allowed k = 2 sectors.We further note that this limitation does not affect our dynamical simulations, which represent the state as a general MPS rather than using the quasiparticle basis states, and as such do not affect the suitability of such simulations as benchmark problems for quantum dynamical simulations.

E. Sources of error
Vacuum and false vacuum.-Thequality of our MPS approximations, |Ω⟩ and |Ω⟩, to the true and false vacua, is dependent on the MPS bond dimension D and on the success of the optimization procedure used to find the tensors and .We choose D high enough to ensure that after optimization, the smallest Schmidt coefficient under a cut is O(10 −6 ) or smaller (in amplitude -this corresponds to probabilities O(10 −12 )).We optimize the vacua until the norm of the energy gradient vector is < 10 −11 , indicating that we have, to very good approximation, an energy eigenstate.It is worth noting that, since the vacua are forced to be translation invariant, in simulations any inaccuracy leads to spatially uniform (global quench) dynamics that are easily distinguished from localized quasiparticle wavepackets.

relations.
Kinks, antikinks, and mesons.-Theaccuracy of our kink, antikink, and meson Fourier modes is limited by the vacuum bond dimension.The quality of the localized wavepackets constructed from the tensors , , which we use to initialize our simulations, can be evaluated by projecting the wavepacket states onto the basis of Fourier modes.We do this for individual kinks in App.C, finding errors of O(10 −6 ).Although we do not directly check the accuracy of the Fourier modes themselves, the observation that our kinks and antikinks propagate, to very good accuracy, in a stable fashion (see the results section below) until they collide demonstrates that we are successfully capturing the targeted topological quasiparticles.
Bubble states.-Beyondany errors in the individual kink and antikink wavepackets, we ensure that, in our bubble states, the kink and antikink are sufficiently well separated such that the initial kink and antikink do not interact significantly.What can be considered a sufficiently large separation is evaluated in App.D. Time evolution.-Wesimulate the full quantum dynamics of our initial states under our spin chain Hamiltonian.The only sources of error, aside from the storage and manipulation of complex numbers using a floating point representation (we use 128-bit, "double" precision for complex numbers), are (i) the numerical integration method used to step through time and (ii) the restriction of the allowed MPS bond dimensions D j to a chosen maximum value (for computational efficiency).In our simulations, (ii) is by far the dominant source of error (see App. E).For this reason, we carry out our simulations using several different values for the maximum allowed bond dimension, in order to probe the sensitivity of each result to this parameter.
Particle detection.-Asfor the bubble states, the particle-pair states we use to detect outgoing particles are dependent, for accuracy, on the quality of the kink, antikink, and meson tensors, as well as the validity of neglecting interaction effects.Again, a minimum separation between excitation tensors in the pair states is needed to avoid the latter source of error: see App.D. In some cases, components of the evolving state may not neatly fall into particle sectors due to lack of separation.For example, confinement may prevent all terms in an outgoing bubble state from being sufficiently well separated.We note that, even in such cases, our main result -that we can accurately simulate relativistic kinkantikink collisions -is not affected by these limitations of the particle detection scheme.

IV. RESULTS
In the following text and figures, positions are given in lattice sites (relative to the leftmost point in the simulation window) and times are scaled so that the maximum kink velocity, determined from dynamical simulation, is That there is a good match shows that the kink-quasiparticle ansatz accurately captures the confined quasiparticles present in the Z2-broken Ising model. 1 lattice site per unit time.Note that we only scale times in this way: energies are given, where not specified as ratios, in unscaled lattice units of (1).Momenta p are given in lattice units −π < p ≤ π.

A. Kink dynamics
In the following we consider kinks, but the discussion applies equally to antikinks.The evolution of a kinkquasiparticle wavepacket will generically involve propagation and spreading (delocalization).We wish to con-struct wavepackets that are sufficiently broad so that they spread slowly, relative to propagation.Broader spatial wavepackets lead to slower spread because they have narrower momentum support; furthermore, spreading is reduced for wavepackets with higher momentum, because the relevant part of the kink-quasiparticle dispersion relation E κ (p) looks increasingly linear.
For h = 0, we observe that our kink wavepackets indeed spread slowly as they propagate at their initial set momentum (see App. H).In the presence of a confining force from a symmetry-breaking field h > 0, kinks undergo acceleration, as expected.A stationary kink is initially accelerated in the direction of the false vacuum, as the energy of the false vacuum is converted into kinetic energy of the kink, as would also be expected in a relativistic QFT, but the long term behavior is strongly influenced by the lattice.The lattice momentum p is bounded −π < p ≤ π, and the momentum expectation value of the kink wavepacket precesses around the unit circle with ṗ = constant.To understand how the position of the kink evolves as this happens, we must consider the wavepacket group velocity v(p) := ∂E κ (p)/∂p.With an emergent relativistic QFT governing the IR physics, the dispersion relation is approximately relativistic (E κ (p) ∼ p 2 + m 2 κ for a kink of mass m κ ) for small |p|, becoming almost linear as p increases.However, due to the bounded nature of p on the lattice, E κ (p) must deviate from relativistic behavior as |p| continues to increase.Indeed, assuming E κ (p) is smooth, including at the boundary value , it is also bounded from above and below.As such, a wavepacket will typically reach a maximum group velocity for some p(v max ), after which it will begin to slow down.Assuming E κ (p) = E κ (−p) it will ultimately reverse and retrace its path back to its original position and momentum (with some wavepacket spread), performing Bloch oscillations.These effects are demonstrated in the single-kink simulation of Figs. 5  and 6.

B. Bubble dynamics
Instead of Bloch oscillations of individual kinks, we wish to study the emergent relativistic dynamics of falsevacuum bubbles comprised of a kink wavepacket and an antikink wavepacket.In particular, we want to simulate kink-antikink collisions at large kinetic energies (to increase particle-production amplitudes).Since the kink and antikink accelerate toward each other under the confining force, we can increase the kinetic energy at the time of collision by increasing the initial bubble size x R − x L (and hence the amount of energy stored in the false vacuum).However, if we allow the kink and antikink to evolve for too long prior to collision, their momenta will exceed |p(v max )| and they will begin to undergo Bloch oscillations, deviating from their relativis- tic behavior.We can ensure that this does not occur by limiting the initial bubble size, with the maximum size depending on the Hamiltonian parameters g, λ, h.In general, a smaller mass gap (since h ̸ = 0, this is the meson mass m µ ), measured in lattice units, increases the maximum bubble size, measured in physical units (multiples of the lattice correlation length ξ).Moving closer to criticality thus allows us to reach higher collision energies relative to the mass gap while keeping |p| below We simulated bubble dynamics for the Ising model (λ = 0) as well as near to the Tri-Critical Ising point of the extended model (λ > 0) for a range of parameters.We first focus on the two points marked (i) and (ii) in Fig. 3.

C. The Ising model
In the Ising case (λ = 0) with h = 0, known to be a theory of free kinks, our simulations reproduce the expected trivial scattering: kinks given an initial nonzero momentum collide without generating any additional entanglement.With explicit symmetry-breaking 0 100 200 300 400 500 600 700 time For Ising (i), the small probability after the first collision of t ≈ 140 indicates elastic scattering of kinks, in stark contrast with the TCI case (ii), where the probability remains high after the first collision at t ≈ 250.In (i), the growth of the post-collision probability with subsequent collisions is consistent with delocalization of the wavepackets, since contributions from kink-antikink pairs with k − j < 60 are not counted.Hence the larger values for (i) after the second collision at t ≈ 430 should not be taken as evidence of inelastic scattering.
0 < h ≪ 1 we find nontrivial scattering, as evidenced by entanglement between the post-collision kink wavepackets.However, even when the energy is significantly above the meson pair-production threshold E > 2m µ there is no obvious ballistic spread to indicate production of unconfined particles.In Fig. 7, we show results for g = 0.8, h = 0.007, where we observe the model to have a mass gap (meson mass) of m µ ≈ 1.04 and a maximum kink velocity of v max = 1.6, both in unscaled units of the lattice Hamiltonian (1).The vacuum correlation length (the length scale of the exponentially-decaying vacuum correlations) is ξ ≈ 1.64 lattice sites.To confirm that no particle production (even of confined particles) is occurring, we compute overlaps ⟨κκ jk |Ψ(t)⟩.We compute the probability of being in the κκ (0,0) "sector" by approximating the integral where and the minimum separation d is chosen to avoid interaction effects between the kink and antikink.See App.D for details of the approximation.In this case d = 60 lattice sites.We plot the probability 1 − P κκ of not being in the κκ (0,0) sector as a function of time for g = 0.8 in Fig. 8.The results are consistent with purely elastic scattering of kinks: The probability is estimated to be around O(10 −5 ) before the first collision, which occurs at t ≈ 140.This is the same order of magnitude as our numerical estimate of the accuracy of the kink-antikink quasiparticle basis states, as detailed in App. C.After the first collision, 1 − P κκ again drops to O(10 −5 ), at a value slightly higher than the pre-collision value.The difference closes as the maximum allowed MPS bond dimension is increased, leaving little room for any inelastic scattering process.
During the first collision, components of the state leave the space of well-separated localized quasiparticles as the kink and antikink approach each other and begin to interact.We explicitly exclude these interacting states from the κκ subspace via the minimum separation d = 60 in (13), hence 1 − P κκ increases significantly until the kink and and antikink wavefronts separate again.
After the second collision, at t ≈ 430, 1 − P κκ does not drop as low as O(10 −5 ).However, this should not be interpreted as evidence of inelastic scattering.As is apparent in both Fig. 7 and Fig. 8, the kink and antikink wavepackets broaden significantly with time and successive collisions.This broadening leads to more terms in the wavefunction in which the separation of the kink and antikink remains less than d even between collisions.We cannot unambiguously count these terms toward the κκ sector due to interaction effects.
While we cannot entirely rule out inelastic scattering using our data, the very small value of 1 − P κκ after the first collision tells us that any inelastic process would have to be extremely unlikely to be consistent with these results.This observation is surprising given that the spin chain and its emergent field theory are not integrable, but consistent with recent observations of nonthermalizing states in the Ising model [45,65,66].We further find that elastic scattering persists even if we allow the kink lattice momentum to exceed p(v max ), as it does in simulation (i) of Figs. 7 and 8 (see App. G), so that the emergent relativistic field theory is no longer a good description of the physics.This is strong evidence that, in the Ising chain with a weak longitudinal field, bubbles are stable up to arbitrarily high energies: When a bubble is large enough, its walls will not meet due to Bloch oscillations, so no scattering can occur while it remains localized.When bubbles are small enough for the kinks to collide, our evidence suggests they do so elastically with extremely high probability.

D. Near the Tri-Critical Ising point
Going to nonzero λ = 0.41, with g = 0.98, h = 0.001, we find a mass gap of m µ ≈ 0.43 and a maximum kink velocity v max ≈ 1.58, both in lattice units of (1).The vacuum correlation length is ξ ≈ 3.6 lattice sites.We choose the initial bubble size so that the energy, shown For mesons, energies are shown with and without a weak longitudinal field.Individual kinks do not have a finite energy for h > 0. Threshold energies for pair production are shown (computed assuming h = 0 for kinks and h = 0.001 for mesons), as is the energy (labelled Ψ) of the simulation shown in Fig. 7 for parameter-set (ii).The dispersion relations reach their maximum gradient at momentum p(vmax) ≈ 0.6.Figure 10.Spin expectation values, after projecting into selected quasiparticle subspaces and normalizing, for simulation (ii) of Fig. 7, bond dimension D ≤ 128, after the first collision (t ≈ 426).The amount of wavefunction captured by each (approximately orthogonal) subspace is given as a probability P (see App. D).Included subspaces are µµ (0,0) , a pair of mesons of lowest energy, and κκ (a,b) , a bubble made of a kink of type a and an antikink of type b (where 0 is the lowest-energy kink quasiparticle, and 1 is the next highestsee Fig. 9).
in Fig. 9, is well above the pair-production threshold, but still low enough to keep the kink velocity ≪ v max at all times.Here we find clear evidence that unconfined particles are produced.Most apparently, Fig. 7 (ii) clearly shows ballistic spread of wavepackets emanating from the first collision event.To further resolve the scattering outcomes, we project onto meson-pair and kink-antikink-pair quasiparticle bases, finding four dominant "sectors", illustrated in Fig. 10, where we tune the quasiparticle basis MPS to match the momentum expectation value of the outgoing quasiparticle wavepackets (as estimated from the Fourier transform of the projected wavefunctions) and compute the spin expectation values of the projected wavefunction for each sector.We also compute the scattering outcome probabilities (the norms of the projected wavefunctions) 6 .We find the most likely outgoing configurations to be: a bubble made of type-0 kinks κκ (0,0) (elastic channel) with probability P = 62%, then a type-0 meson pair µµ (0,0) with P = 19%, and finally a bubble made either of a type-0 kink paired with a type-1 antikink (higher energy) κκ (0,1) , or a type-1 kink paired with a type-0 antikink κκ (1,0) , each with P ≈ 7% (reflection symmetry).These outcomes are all kinematically allowed, according to the energetic thresholds shown in Fig. 9.
We note that the (rounded) projection probabilities in Fig. 10 only add to 95%.This may indicate the presence of other sectors we have not accounted for, such as a µ (0) paired with a small κκ (0,0) bubble, or a κκ (0,0) bubble containing one or more quasiparticle-excitations of the false vacuum.Unfortunately, since these "sectors" each involve at least three quasiparticles, the corresponding position bases have many more terms (O(N 3 w ) versus O(N 2 w ) for pairs), making it difficult to compute these projections 7 .We emphasize, however, that this limitation only hinders this particular form of analysis of the evolved state.The dynamical simulation itself is not affected as it is not restricted to these multi-particle basis states.
It is also possible that various sources of error have affected results: (i) when excitation tensors in a 2quasiparticle MPS are close together, so that interactions are relevant, the state may not accurately represent quasiparticles, (ii) the MPS representations of the quasiparticle position states are variational approximations subject to some error (which also affects the initial state of the simulation), and (iii) although we allow the MPS bond dimension to increase up to some maximum during simulations (D ≤ 128 in this case), errors can still accumulate if that maximum is insufficient to capture all entanglement, as well as due to errors in the numerical integration steps.We did not explicitly characterize the effects of (ii), but we have indirect evidence that they are small: see Sec.III E. By varying the minimum quasiparticle separation used in the projection, as well as the maximum bond dimension of the simulation, we were able to characterize effects (i) and (iii), finding them to amount to changes in the outcome probabilities of ≪ 0.01, except in the case of κκ (0,1) and κκ (1,0) , in which one of the quasiparticles is heavier than the other, leading to a smaller separation between the kink and antikink.In this case, our analysis suggests the error here amounts to a change of around ±0.01 in the outcome probability, possibly more (see App. D).This outcome might be better resolved at higher energies, at which the kink-antikink separations would increase.
In case of the µµ (0,0) outcome, we cross-check the computed outcome probability by comparing it with the excess energy (relative to the vacuum) E pkts of the regions containing the ballistic wavepackets, visible in Fig. 7 (ii).If these wavepackets belong to a two-meson "branch" of the wavefunction, that branch (the portion of the wavefunction in the µµ (0,0) subspace) must contribute EP to the energy, where E is the total energy and P is the probability of the µµ (0,0) scattering outcome.We can therefore estimate P as E pkts /E.This gives us a P within the range 19% to 20% at t ≈ 757 (after separation), depending on the precise extent of the region we sum over (e.g. from site 0 to site 250 for the left packet), compatible with the projected µµ (0,0) wavefunction.

E. Entropy and computational cost
As evidenced by Fig. 8, the bond dimension of the MPS representing the evolving state must continue to grow as time goes on, in order to maintain accuracy.The cut entropy at location j is a proxy for the required bond dimension D j .Fig. 11 shows the evolution of the maximum cut entropy for the simulations of Fig. 7.At early times, we observe that the maximum cut entropy jumps dramatically during scattering events, whether elastic or inelastic, remaining almost constant in between.This is consistent with a model of interacting quasiparticle wavepackets: Separated wavepackets undergo stable propagation until they collide, at which point interactions generate entanglement corresponding to the different possible scattering outcomes.At late times, we observe a temporal broadening of the jumps, consistent with spatial broadening of the wavepackets involved.Fig. 11 also shows that, although there is cut entropy associated with the wavepackets themselves, this is quickly surpassed by the cut entropy in the center of the chain, associated with entanglement between the left and right outgoing packets.It is this entanglement between outgoing quasiparticles that is responsible for the post-collision plateaus visible in the maximum cut entropy.
The entropy jumps clearly make MPS simulations of long-time dynamics demanding.However, for the purposes of studying the quasiparticle content of scattering outcomes, with the incoming quasiparticles chosen via the initial state, it is enough to accurately simulate a single collision and then wait until the outgoing wavepackets have separated sufficiently so that interactions between outgoing quasiparticles may be neglected (we assume the simulation parameters are chosen such that wavepackets remain localized for sufficiently long times) 8 . 8Simulating for only relatively short times also increases the like-We expect the entropy generated in a collision of localized quasiparticle wavepackets to depend on the collision energy relative to the masses of quasiparticles: the more scattering outcomes there are, and the greater the probability of those outcomes, the larger the post-collision entanglement entropy can be.In Fig. 12, we explore the maximum cut entropy occurring during the first collision as a function of energy, controlled via the initial kink separation, for two sets of Hamiltonian parameters, one (ii) closer and one (iii) further from criticality.As each datapoint on this Figure requires a full dynamical simulation, we limited the bond dimension to D ≤ 96 (versus D ≤ 128 of Fig. 7) to save some computational time.We find that the trend with increasing bond dimension is nevertheless clearly visible.We find that the entropy indeed grows with energy, smoothly increasing even as thresholds are crossed, e.g. the 2m µ (1) , 3m µ (0) and 4m µ (0) thresholds in case of (ii) (see also Fig. 9).
The entropy continues to increase at least until the energy is sufficient for the kinks to approach the maximum possible kink velocity prior to collision, at which point we expect deviations from the emergent relativistic dynamics to become apparent as Bloch oscillations emerge.In the case of (iii), Fig. 12 shows that the postcollision entropy eventually decreases as lattice effects kick in, coincident with deceleration of the kinks prior to collision.Note that we are able to reach much higher relative energies with parameters (ii) before encountering obvious lattice effects.This illustrates the general principle that more of the emergent relativistic QFT is revealed as one approaches criticality: the relative energies accessible by quasiparticles, while avoiding Bloch oscillations (momenta |p| < p(v max )), grows as the lattice meson mass drops.
We also observe that much more entropy is generated in the first collision for parameters (ii) than for parameters (iii), even when the relative energy is similar 9 .A significant part of this difference likely comes from a much higher probability of meson pair production, as well as the availability of the κκ (0,1) outcomes, in case (ii): the probability of particle production is < 10% in case (iii) at energy E/m µ ≈ 2.52, in contrast with ∼ 38% at energy E/m µ ≈ 2.62 in case (ii), according to κκ (0,0) basis overlaps.This is possible, since these two parameter sets were not chosen to be part of an RG trajectory, so that their emergent QFTs need not be the same.

V. DISCUSSION
Building on recent innovations in the classical simulation of quasiparticle dynamics using Matrix Prodlihood that these tasks can be carried out on NISQ-era quantum devices. 9The difference in post-collision entropy is not attributable to different vacuum entropies, since these are < 0.1 in both cases.
uct States [49], we proposed a framework for simulating and characterizing the full (nonperturbative) quantum dynamics of false-vacuum bubbles in relativistic QFTs that govern the IR physics of one-dimensional lattice systems.While we chose to simulate a quantum spin chain, the methods we use are general and could also be applied directly to, for instance, a spatially-discretized QFT such as the Schwinger model or λϕ 4 theory [10-15, 17, 42, 43].We also demonstrated that the MPS quasiparticle ansatz, with which we initialized our simulations, can be used to detect quasiparticles that are produced as time evolves.This allowed us to verify quasiparticle pair-production in the modified Ising model we studied, including production of different species of confined kink that were not obvious from examining energy density and spin expectation values alone.We used the same kind of analysis to confirm a lack of particleproduction in the unmodified Ising model (with transverse field and small longitudinal field), supporting other recent studies that suggest particle production is very strongly suppressed [45,65,66].
We were able to significantly improve the efficiency and interpretability of our simulations by carefully choosing our initial states in two different ways: Firstly, constructing spatially broad wavepackets allowed us to access the dynamics of the emergent IR QFT without the spoiling effects of UV, high-momentum components that are strongly influenced by the lattice.Broad wavepackets also lead to localization of quasiparticles over long times, making it easier to characterize scattering outcomes, and improve the numerical conditioning of the dynamical simulation (see App. F).Secondly, by precisely tuning the quasiparticle content of the initial wavepackets [49], we were able to study individual scattering events in isolation, while further reducing the computational demands of the simulation by lowering entanglement.
Entanglement growth is the most significant barrier to dynamical simulations with MPS, as the computational cost of each time step scales exponentially with the cut entropy.By choosing broad quasiparticle wavepackets, we reduce entanglement growth at the expense of growing the number of lattice sites involved in the simulation.This is a good tradeoff for MPS simulations, as the computational cost scales only linearly in the number of lattice sites in our simulation window.Even with this tradeoff, we found that the large jumps in cut entanglement with each collision (of confined quasiparticles in the system) preclude simulating more than a handful of successive collisions.Furthermore, we found clear evidence of entanglement growth with the collision energy, although the onset of Bloch oscillations prevented us from drawing strong conclusions about how this growth continues in the emergent IR QFT.Nevertheless, in the absence of lattice effects that obscure the IR QFT, it seems reasonable to expect the entropy to continue to grow with energy, which would eventually preclude ac-curate simulation using MPS.
We note that our current simulations are already well beyond the reach of present digital NISQ devices, as they would require at least ∼ 1000 qubits and circuit depths that are larger still.
A natural next step would be to perform simulations along RG trajectories in the Hamiltonian parameter space, so that results can be extrapolated to the continuum.This is equivalent to finding paths toward criticality of the lattice model, along which the low-energy spectrum remains consistent with a particular emergent (IR) QFT.Moving closer to the continuum would also allow us to reach higher (relative) collision energies while avoiding lattice effects, such as Bloch oscillations.In turn, this would permit a more thorough exploration of the energy-dependence of the entanglement generated in collisions.
As we approach criticality, the bond dimension of the MPS vacua must grow to maintain accuracy, as must the size of the simulation window, since the wavepacket width in lattice units would have to increase with the lattice correlation length in order to maintain localization.Getting closer to criticality seems feasible: The simulations featured in the main text, with maximum MPS bond dimension 128, took between one and two weeks to complete on 8 cores each and this time could likely be reduced significantly with further work to optimize the code10 and the use of a numerical integrator with an adaptive time-step size.
Increasing the number of spatial dimensions presents a significantly greater challenge for classical algorithms: while the computational cost of MPS simulations scales with the bond dimension D as O(D 3 ), the scaling for tensor-networks capable of handling large (2+1)dimensional systems, such as PEPS [82], is much worse (albeit still polynomial) 11 .As an intermediate step, one could consider systems with a small, compactified second dimension of space, which are often within reach of MPS methods.By performing a Fourier transform of the Hamiltonian in the compactified direction only [84], one could study scattering of quasiparticles that are spatially localized in one direction, while being momentum eigenstates of the other.Compared to the purely (1+1)dimensional case, the additional "Kaluza-Klein" excitations associated with the Fourier modes of the compactified dimension would already open up a much greater range of scattering outcomes.
Compared to the simulations we performed, increas-ing the variety of scattering outcomes, whether by raising the relative energy in a given model, choosing a lattice model with a richer set of low-energy excitations (e.g.near a phase transition described by a CFT with larger central charge), or adding spatial dimensions, seems necessary in order to find problems that exhaust tensornetwork methods due to the additional entanglement generated.Such problems appear more amenable to simulation on quantum hardware, which is not a priori limited in the amount of entanglement it can deal with.However, raising the energy may be problematic for near-term quantum devices, which are limited both in their size and coherence times, since avoiding lattice effects (such as Bloch oscillations) at higher relative energies requires moving closer to criticality while increasing system sizes and evolution times.
is a D j−1 × D j matrix assigned to site j in basis-state s and v L and v R are appropriately-sized boundary vectors.In a uniform (translation invariant) iMPS, we use the same tensor A everywhere: A (s) j = A (s) ∀j ∈ Z.Such a state has a well defined norm for generic choices of v L and v R if the D 2 × D 2 "transfer matrix" where * indicates the complex conjugate, has a nondegenerate eigenvalue of largest magnitude, with A normalized so that this eigenvalue is equal to 1 [40].This condition implies exponential decay of correlations with distance.By additionally normalizing v L and v R appropriately, we can achieve ⟨ψ|ψ⟩ = 1.The precise choice of boundary vectors does not affect bulk expectation values due to the aforementioned exponential decay of correlators.

Nonuniform windows
To build the bubble states of the main text, and to simulate their evolution in time, we allow the tensors of an otherwise uniform iMPS to vary within a "window", consisting of N w contiguous sites.These states have the form where A L and A R parameterize the semi-infinite left and right bulk parts of the chain and A 1 . . .A Nw parameterize the nonuniform window.The above transfermatrix conditions for a well-defined uniform iMPS must be satisfied for both the A L and A R tensors.The norm of the state is then determined by the content of the window tensors A 1 . . .A Nw .For the bubble states, we let , where is the tensor optimized for the uniform ground state (true vacuum) of the spin chain.We then choose A 1 . . .A Nw to represent a falsevacuum bubble, as described below in App. C. For example, a fully localized bubble state (the kink-antikink state of Fig. 4) has A 1 = (representing a kink), Appendix B: Finding the true and false vacua

Finding the true vacuum
The tensor A defining a uniform iMPS (A1) can be optimized to represent low-energy, translation-invariant states of gapped quantum spin chains using various algorithms.We use the nonlinear conjugate-gradient method described in [10] and implemented in the evoMPS package [80] to find a uniform iMPS that approximately describe the ground states of gapped quantum spin chains.We denote the optimized iMPS tensor .

Finding the false vacuum
We explain how we obtain an iMPS representation of the false vacuum in practice in App.B 3 below.In this section, we consider the nature of the false vacuum more generally.
For the broken Z 2 symmetry of the Ising-like chain in the main text, the false vacuum |Ω⟩ is a state that has opposite spin orientation to the true vacuum |Ω⟩.It should also be like a vacuum, in that it should be a spatially uniform, approximately static state near a local energetic minimum with respect to some constraint, such as locality.
A candidate state is the "flipped" vacuum It is spatially uniform and typically close to an energetic minimum in the following sense: If we apply a finite string L j=1 X j of length L, the change in energy will be positive for small values of L, becoming negative only after the O(2hL) energy lost by replacing false vacuum with true vacuum on L sites is larger than the O(1) energy penalty of spin anti-alignment at the boundaries.However, the flipped vacuum is generally not close to being an eigenstate in case of a nonzero symmetry-breaking field parameter h and therefore is not suitably static.Nevertheless, one might begin with the flipped vacuum and attempt to bring it closer to a false-vacuum eigenstate by lowering the energy, for example via imaginary time evolution: A problem with this approach is that, since we are not in a true energetic minimum, imaginary time evolution will ultimately take us back to the true vacuum |Ω⟩.At finite system sizes, this corresponds to a nonzero inner product between the flipped vacuum and the vacuum.Let us consider the Z 2 -broken Ising Hamiltonian (λ = 0, h > 0) at finite system size N .Although in our simulations we work directly in the thermodynamic limit N → ∞ using iMPS, finite N is more convenient for the following calculation.We will see that the key result is independent of N .If we perturb around the bare theory of g = 0, we find where g ≪ 1. Overlaps ⟨E i | j X j |Ω⟩ with energy eigenstates |E i ⟩ that are close to the true vacuum (e.g.low-energy excitations) are also exponentially suppressed.
A simplified model allows us to estimate the timescale for "decay" to the true vacuum under imaginary time evolution.Take |ψ⟩ := ψ Ω |Ω⟩ + i ψ i |E i ⟩, where ⟨ψ|ψ⟩ = 1 and |E i ⟩ represents an eigenstate in the false-vacuum "sector", i.e. with a flipped spin orientation versus |Ω⟩.This will be our model for the flipped vacuum j X j |Ω⟩.From our perturbative calculation, we take Now we take E i − E Ω ∼ 2N h, since |E i ⟩ are flipped states which suffer an extensive energy penalty compared to |Ω⟩.The relative contribution of the vacuum after a time τ is then which goes to 1 at τ Ω , independently of N : Hence, for small h, one must evolve for a "long" time to see a significant vacuum contribution.For sufficiently large τ still satisfying τ ≪ τ Ω , assuming initial occupancy and energetic separation of the |E 0 ⟩ state, we end up with: where |E 0 ⟩ is a hypothetical lowest-energy contribution from the false-vacuum "sector".This picture is supported by numerical observations in which performing some imaginary-time evolution on j X j |Ω⟩ quickly results in something that is (numerically) approximately an eigenstate.

Finding an iMPS for the false vacuum
To find an iMPS for the false vacuum, we begin with an iMPS approximation of the flipped vacuum (B1), obtained from the iMPS approximation of the true vacuum.
Instead of using imaginary time evolution to reduce the energy of this state, as considered in the previous section, we use the same conjugate-gradient optimization method used to find the true vacuum [80].Like imaginary time evolution, such variational methods should eventually take the flipped state to the true vacuum state.In practice, however, we observe that for small symmetry-breaking fields |h| ≪ 1 this does not happen.Instead, the state converges to a false-vacuum iMPS (parameterized by a tensor we denote ) that is numerically indistinguishable from an energy eigenstate.This may be because of the limited available numerical precision 12 , which could preclude accurate representation of the gradient components that would lead to the true vacuum.

Appendix C: MPS quasiparticle states
We use a Bloch-state approach to represent low-energy excitations [41,74,75].A localized quasiparticle state is constructed from vacuum tensors A L and A R , which remain constant, together with an "excitation tensor" B that can be chosen to represent different excitations:

This ansatz can represent topological excitations, in case
A L and A R refer to different vacua, as well as nontopological excitations, in case A L and A R represent the same vacuum state.We use the symbol ϕ to denote a generic excitation, and κ, κ, or µ to refer to kinks, antikinks, or mesons specifically.For example, the kink states |κ j ⟩ of the main text have A L = , B = and A R = , while the meson states |µ j ⟩ have A L = , B = and For tensor-network diagrams showing the parts of the tensor networks surrounding B (for |κ j ⟩ and |µ j ⟩) see Fig. 4. Due to exponential decay of correlations in the vacua represented by A L and A R , the tensor B represents a quasilocal excitation and may affect expectation values across many lattice sites.In the following, we assume for simplicity that the MPS |ϕ j ⟩ have uniform bond dimension D.
Momentum eigenstates can be constructed as Fourier modes of the spatially localized excitations: These momentum eigenstates enjoy a "gauge" freedom (parameter redundancy): the B tensor may be transformed as where x is a D × D matrix, without affecting the Fourier mode |ϕ(A L , B, A R , p)⟩.This freedom can be fixed in many ways.For example, the "left orthogonality"13 conditions [74,85] are where ⟨l L | is the dominant left eigenvector of the MPS transfer matrix of the left uniform bulk: Similarly, the "right orthogonality" conditions are where |r R ⟩ is the dominant right eigenvector of the MPS transfer matrix of the right uniform bulk: These conditions can be achieved for any initial tensor B by transforming it with an appropriate choice of x in (C3).Imposing either the left or right conditions, (C4) or (C6), implies orthogonality of the position states: This is particularly convenient for working with the momentum eigenstates, as it greatly simplifies the computation of their inner products and expectation values [74].

Optimizing the excitation tensor B
To find a tensor B that accurately represents a particular quasiparticle excitation, we use the methods of [74] with some modifications for dealing with the case in which one of A L and A R represents a false vacuum (with a different energy density compared to the true vacuum).The basic idea is to project the Hamiltonian onto the ansatz space of momentum eigenstates |ϕ(A L , B, A R , p)⟩, resulting in an effective Hamiltonian for the tensor B that can be solved using a standard sparse eigenvalue solver.

Note that, since ⟨ϕ(A
, it is natural to do this for a particular, chosen value of p.By solving for multiple eigenvalue-eigenvector pairs, a set of orthogonal tensors B (a) (the eigenvectors) can be found that accurately approximate several different low-energy excitations (labeled by the index a), as long as they are all below the two-particle threshold [75].The eigenvalues are the energies of these excitations.By computing them for a range of p, one can obtain an approximate dispersion relation E(p) for the quasiparticles in the system.Importantly, not only the energies, but also the optimized tensors B (a) (p), and hence the position states |ϕ j (B)⟩, generally depend nontrivially on the value of p.This is illustrated in Fig. 13, which shows the error made in using a tensor B (0) (p = 0) optimized for p = 0, to represent the lowest-lying excitations at other momenta.

a. Broken symmetry and kinks
In the presence of explicit symmetry breaking (h ̸ = 0), topological excitations such as kinks and antikinks involve the false vacuum.The energy of a localized kink (or antikink) depends on its position, since different positions lead to different extensive contributions from the false vacuum 14 .As such, there are no energy-momentum eigenstates (of the Hamiltonian and momentum operators) corresponding to these excitations and we cannot find them by solving the effective Hamiltonian for the B tensors considered above.Nevertheless, we expect there to be excitations that behave as quasiparticles subject to a confining force (which makes them accelerate).If one could somehow cancel the confining force, as if by accelerating at the same rate as the quasiparticle, the latter would appear to propagate freely.
With this picture in mind, we define a modified energy function where ∆E j := ∞ j+1 e R and e L , e R are the energy densities of the vacua parameterized by A L and A R , respectively.Here, it is assumed that the momentumeigenstate gauge-freedom on B, (C3), has been fixed so that ⟨ϕ j |ϕ k ⟩ = δ jk .With orthogonality of the position states, the identity term simply shifts the energy of the excitation in a position-dependent way, cancelling the position-dependent contribution due to the differing bulk energy densities.One can also write down a modified Hamiltonian where P j is a projector onto the space of states spanned by |ϕ j (B)⟩ for all B satisfying the chosen orthogonality conditions: for such B, we thus have We can thus optimize B by computing eigenvalueeigenvector pairs of H, after pushing it into the ansatz space, analogously to the symmetric case above.
In this formulation it is manifest that the optimization procedure depends on the conditions used to achieve position-state orthogonality, since different conditions will lead to different P j in (C9).In practice, we find that the difference this makes to the resulting optimized states |ϕ(A L , B, A R , p)⟩ is small: We consider the infidelity per site where B LF and B RF are optimized for the same quasiparticle (and momentum p) using the left and right orthogonality conditions, (C4) and (C6), respectively, and |Z| is the cardinality of the integers (accounting for the infinite norm of the momentum eigenstates).We find empirically the infidelity scales as h 2 , where h ≪ 1 is the symmetry-breaking parameter of the Hamiltonian, for the lowest-energy kink states and the Hamiltonian parameters considered in this paper.This dependence is shown in Fig. 14.
It may be possible to design improved optimization techniques that avoid this ambiguity.A prerequisite would be a cost function other than ⟨ H⟩, presumably related to the stability of the quasiparticle wavepackets, that distinguishes usefully between the different choices that can be made in parameterizing B.
Another important observation is that, in (C9), we assume that the location j of the B tensor reliably indicates the position of the kink (or antikink).In fact, since the excitation described by B is quasilocal, the location of the kink (defined as the point in space at which the spin expectation value crosses zero) may differ from j. Indeed, as shown below in Fig. 15, there may be a relative shift of several lattice sites, depending on the choice of left or right orthogonality conditions on B. The shift will generally also depend on the momentum p, such that a dispersion relation E(p) computed from the eigenvalues of H should really be interpreted as a function of the B-tensor momentum, derived from the position j of B, considered distinct from the kink momentum, derived from the kink position.
The momentum-dependent energy-shift due to these position shifts can be calculated by first computing the actual kink positions, relative to j, for each |ϕ j (p)⟩, as a function of p.These shifts can then be multiplied by the false-vacuum excess energy density to compute the energy shift, which can in turn be used to "correct" the dispersion relation.This provides a more intuitive definition of the kink dispersion in the symmetry-broken setting.Note also, however, that since kink quasiparticles do not have a well-defined energy gap with respect to the vacuum state, these dispersion relations still cannot be used to compute particle-production thresholds.They could, however, be used to estimate the kink velocity dE(p)/dp, this being independent of energy shifts

Wavepackets
Analogously to the momentum eigenstates of (C2), we can construct wavepackets from the localized quasipar- ticle states as where in our simulations we choose f i to be a Gaussian centered at position x with width σ.Importantly for our purposes, it is straightforward to turn such a wavepacket state into a single MPS: where is a 2D × 2D matrix, given that A s L , A s R , and B s are all D × D matrices.We set the boundary conditions to be for some generic choice of v L and v R .If |f i | falls below some numerical threshold for all i less than some i L and for all i greater than some i R > i L , we can truncate it to zero and reduce the bond dimension to D in those regions without introducing significant errors.This allows us to represent a truncated wavepacket in the thermodynamic limit using the nonuniform window ansatz (A3).
As discussed in the main text, this wavepacket construction ignores any momentum dependence of the tensor B. While this introduces errors in the form of contributions from other excitations, as shown in Fig. 13, these become small for large σ, as indicated in Fig. 17 below.

Localized states and Bloch-state parameter redundancy
The parameter redundancy, or "gauge freedom", on the B-tensors of the momentum eigenstates (C2) does not leave the localized states |ϕ j (B)⟩ of (C1) unchanged.These states, as well as the wavepacket states (C12) built from them, depend on how these degrees of freedom are fixed.However, the procedure we use for choosing optimal B tensors is based on momentum eigenstates and does not tell us how to optimally fix the remaining freedom.
That said, the impact of this choice on Gaussian wavepacket states must vanish in the limit σ → ∞, where packets become momentum eigenstates.Hence we can reasonably expect the impact on wavepackets with finite width to become small as σ increases.Since we already have a physical reason to choose broad wavepackets in our simulations (for slow wavepacket spread), this issue is not as severe as it may at first appear.
Nevertheless, we choose to use the reflectionsymmetric conditions of [49], slightly adapted for the topologically nontrivial setting, to fix the remaining freedom on the B tensors used to construct our initial bubble states.To be precise, we fix B by choosing x in (C3) as In terms of tensor networks, we can rewrite this as Unlike the left and right orthogonality conditions, (C4) and (C6), these conditions are manifestly symmetric under spatial reflections.The reflection-symmetric conditions can be formulated as an overdetermined linearleast-squares optimization problem and then solved using standard techniques.See Fig. 15 for a comparison of the reflectionsymmetric conditions to the left and right orthogonality conditions, in terms of the spin expectation values of localized topological states.We plot these results for A L and A R representing the two vacua of the Ising chain (λ = 0, g = 0.8) for both zero and nonzero longitudinal field strength h.The tensor B is variationally optimized, as described above, so that the momentum eigenstate |ϕ(A L , B, A R , p = 0)⟩ approximates the lowest-lying topological excitation.We note that, although the symmetrized states exhibit "smoother" spin expectation values in both cases, they are not perfectly symmetric in the case h > 0. The spatial asymmetry in the spin likely reflects the energetic asymmetry of topological states in this case.In both cases, there is certainly an aesthetic improvement to be had by imposing the symmetrization conditions, but it remains to be seen whether the symmetrized states are better representations of localized quasiparticles.To see that they are, we consider how well wavepackets built from them fit into the corresponding quasiparticle subspace.
In Fig. 17, we plot the portion of a kink wavepacket state (by probability) outside of the targeted kinkquasiparticle subspace for both the orthogonal and reflection-symmetric conditions.We explain how to carry out this kind of projection in App.D. We observe that the symmetrized states result in a much more accurate wavepacket than the orthonormal states, by almost two orders of magnitude, confirming that the symmetrized choices are more than just aesthetically pleasing.In the symmetric Ising model (h = 0), we can also compare the energy of wavepackets.In Fig. 16, we see

E E orthonormal symmetrized
Figure 16.Single-kink wavepacket energy as a function of the width σ, for the Z2-symmetric Ising model g = 0.8, λ = 0, h = 0 (in the SSB phase).We plot the energy (relative to the energy of the kink-quasiparticle momentum eigenstate) for two different ways of fixing the momentum-eigenstate freedom on the B tensors used to construct the wavepacket state: an orthogonal choice ⟨κj|κ k ⟩ = δ jk (left and right conditions are equivalent when h = 0) and the reflection-symmetrized choice (C17).In both cases the MPS tensors used to construct the state are tuned to the wavepacket momentum p = 0. Figure 17.Portion of a single-kink wavepacket state outside of the single-kink subspace κ (0) as a function of the wavepacket width σ, for the Z2-broken Ising model g = 0.8, λ = 0, h = 0.007.We plot the error for two different ways of fixing the momentum-eigenstate freedom on the B tensors used to construct the wavepacket: the left orthogonal choice ⟨κj|κ k ⟩ = δ jk and the reflection-symmetrized choice (with B optimized using the left orthogonal conditions).In both cases the MPS tensors used to construct the state are tuned to the wavepacket momentum, which is p = 0.The projection into the κ (0) subspace uses B tensors optimized using the left orthogonal conditions and fully accounts for momentum dependence via a Fourier analysis.
that the kink wavepackets created with the symmetrized states have consistently lower energy, which indicates improved accuracy, since we are targeting the lowest-lying kink quasiparticle.

Two-particle states
To create false-vacuum bubbles, we need to combine two quasiparticles, namely a kink and an antikink.Twoparticle states have the form where, compared to (C1), we now have a central (false) vacuum tensor A C , as well as two excitation tensors, B L and B R , instead of one.Analogously to the one-particle states, we use ϕϕ to denote a generic pair of quasiparticles, specifying κκ or µµ when we are discussing a kinkantikink pair or a meson pair specifically.Such states are illustrated in Fig. 4. We will assume that the quasiparticles are sufficiently well separated, so that interactions may be neglected and B L and B R can be held constant irrespective of the separation d := k−j.A sufficient condition for this to be justified, is that the reduced state for sites i in between the two excitations j < i < k reverts to that of the central vacuum MPS parameterized by A C for some range of i.If this happens for the reduced state on at least r contiguous sites, where r is the range of interactions in the Hamiltonian (for our model, r = 2 when λ = 0 and r = 3 otherwise), this implies that the energy of the two-particle state, as a function of the separation d, is not affected by interaction between the particles at that location (only by differences in the "vacuum" energies).This means there can be no interaction energy.This condition can be made more precise: we define the strength of interaction effects at location i (with j < i < k) as the deviation of the left and right environment tensors of the 2-particle state from the corresponding environment tensors of the central (false) vacuum MPS parameterized by A C .Here, the environment at site i is the tensor network for the reduced state on site i, excluding the tensors assigned to site i itself.It naturally splits into left and right components, consisting of the tensors to the left and to the right of i, respectively.We compute the deviation for the left and right parts separately, as the norm of the difference between the central (false) vacuum environment and the 2-particlestate environments.For the left environment, we define where the ellipsis indicates that the center transfer matrix should be repeated as many times as is necessary to reach site i from the position of B L .Similarly, for the right environment we define ϵ R (i) := . . .
We define the magnitude of interaction effects at separation d to be This value is plotted, for the Hamiltonian parameters (ii) of the main text, in Fig. 19 for a selection of low-energy quasiparticles.
Wavepackets can be constructed from the two-particle position states (C19) analogously to the single-particle case (C12): where, for our simulations, we choose the packet functions f j and g k to be Gaussians centered at x L and x R , with momenta p and −p, respectively.As in the singleparticle case (C13), these wavepackets can be rewritten as a single MPS with bond dimension 2D, where D is the bond dimension of the vacua A L , A C , A R .Note that, because position states are not defined for k ≥ j, the wavepacket functions f and g are effectively truncated, leaving only the terms j < k, in this ansatz.Of course, if the wavepacket functions have negligible overlap, the effects of this truncation can themselves be neglected.

Appendix D: Particle detection via quasiparticle basis states
As discussed in the main text, it is possible to use the optimized quasiparticle states and their two-particle combinations to estimate the particle content of a wavefunction |ψ(t)⟩ as it evolves during simulation.For example, the inner product ⟨κκ The single-particle quasiparticle position states |ϕ j ⟩ of (C1) can be made orthogonal by imposing either the left or right orthogonality conditions, (C4) or (C6).For the two-particle states |ϕϕ A L and A R in the two-particle states (C19), it is then straightforward to (approximately) project |ψ(t)⟩ into the subspace spanned by these states.
However, we must take care when interpreting the overlaps of a wavefunction |ψ(t)⟩ with quasiparticle position states such as |ϕ j (B)⟩, as should be clear from the discussion of excitation tensors and quasiparticle position states above.In particular, the momentumdependence of the B tensors optimized to represent each quasiparticle, as well as the ambiguity in fixing the degrees of freedom (C3) in B that are not fixed by the optimization procedure (see App. C 1), make the interpretation of the overlap unclear unless the quasiparticle content of |ψ(t)⟩ consists of broad spatial wavepackets, whose momentum support is focused around the momentum p used to optimize the B tensor.We next discuss two methods for avoiding these issues.

Checking consistency in the projected wavefunction
As described in the main text, one way to avoid issues with momentum-dependence and wavepacket breadth is to choose some momentum p, optimize a tensor B(p) at that momentum for the quasiparticle being targeted, then examine the overlaps ψ j := ⟨ϕ j (B(p))|ψ(t)⟩ (or ψ jk := ⟨ϕϕ jk (B(p), B ′ (p ′ ))|ψ(t)⟩ for a two-particle basis).If the wavepacket width of the projected wavefunction ψ j is sufficiently large, and the momentum support (computed via Fourier transformation) sufficiently close to the chosen value of p, we know that the error made is small and can trust that the projection is accurately telling us about the quasiparticle content.We can quantify how broad the wavepacket must be, and how close the wavepacket momenta should be to p, via analyses such as those of Fig. 13 and Fig. 17.
If the distribution of momenta in the wavepacket indicates a large error due to the choice of p made while constructing the basis, it may be possible to iteratively tune p to achieve a better match.This procedure fails, of course, if the wavepackets in |ψ(t)⟩ are too narrow in position space (and hence too broad in momentum space) for the error to be kept small.
We use this procedure to compute the spin expectations values within each particle "sector" in Fig. 10, tun-ing the basis momenta p L and p R for two-particle bases |ϕϕ jk (B L (p L ), B R (p R )⟩ to match the observed momenta.

Fourier analysis
To fully account for the momentum-dependence of the quasiparticle states, and to eliminate any issues due to the momentum-eigenstate "gauge-freedom" (C3) on B tensors, we can simply take overlaps with momentum eigenstates instead of with the position states.These states are exactly invariant under (C3), and the B tensors used to construct them can be optimized for the momentum of the eigenstate to avoid momentum mismatch.
For example, to project onto a single-particle subspace at momentum p, one can compute ⟨ϕ(p)|ψ(t)⟩, with |ϕ(p)⟩ from (C2).We can expand this overlap in terms of position states: where B(p) is an excitation tensor optimized to represent the quasiparticle ϕ at momentum p.This overlap can be computed in practice, despite the infinite sum over j, because the position of excitations in |ψ(t)⟩ is limited to the nonuniform window of (A3) in which the initial quasiparticles (comprising the bubble) were placed, so that there are only ∼ N w nonzero position terms in this overlap.
To compute the projection onto the entire quasiparticle subspace, we must evaluate the integral This can be done approximately by sampling, for example using a numerical integration scheme.In practice, we use the Fast Fourier-Transform (FFT) algorithm to transform the spatial components ⟨ϕ j (B(p))|ψ(t)⟩ into a fixed sampling of momentum components at a resolution determined by the number of lattice sites N w summed over in (D3).We use this method to compute the projected singlekink wavefunction described in Fig. 17, to compute the kink-antikink scattering outcome probability in Fig. 8 of the main text, and to compute the various scattering outcome probabilities reported in Fig. 10.

a. Efficient computation of the quasiparticle Fourier analysis
The projection of a wavefunction |ψ(t)⟩ onto many momentum modes is relatively computationally intensive, since for each momentum mode we must first compute a tuned excitation tensor B(p), followed by a full set of position overlaps ⟨ϕ j (B(p))|ψ(t)⟩.In the case of singleparticle states, the cost is O(N w M ), where M is the number of momentum samples and N w is the nonuniform window size of |ψ(t)⟩.For two-particle states, since we must consider cases in which the two particles have different momenta/positions, making the cost O(M 2 N 2 w ).If M ∼ N w , and N w ∼ 1000, this may be prohibitive!To reduce the cost, we can use the observation that the optimized excitation tensor B(p), for a given quasiparticle, usually varies only slowly with the quasiparticle momentum p by introducing a small momentum mismatch in a controlled way: We project the excitation tensors B(p), optimized for each mode of momentum p that we wish to sum over, onto a small basis of excitation tensors that capture the momentum dependence accurately across a wide range of p: for appropriately chosen coefficients b α (p) and suitably chosen basis tensors B α .A suitable basis can be built by orthonormalizing (via a Gram-Schmidt procedure) a set of B(p α ) obtained at a selection of momenta (say, p = −2, −1, 0, 1, 2).We find that < 10 basis tensors is sufficient, for our chosen Hamiltonian parameters, to achieve an accuracy of ∼ 10 −8 in (D5).Given such a basis, we then compute position-state overlaps only for the basis tensors, from which we can compute the projection of the wavefunction onto an arbitrary momentum mode while making only a small error.For example, in the case of a two-particle basis, we first compute ψ jk;αβ := ⟨ϕϕ jk (B L,α , B R,β )|ψ(t)⟩.Then, using the coefficients of (D5), we approximate the overlap with the momentum-tuned position states as (D6) From here, the momentum-mode overlaps are but an FFT away.

Other sources of error
Even for the Fourier analysis, there are at least two sources of inaccuracy beyond the choice of iMPS bond dimensions in the vacua and the quality of the optimization procedures used to find the vacua and the excitations tensors B. First, there is the further ambiguity (see App. C 1 a) in the case of topological excitations in the symmetry-broken setting (h > 0) owing to the dependence of the B-tensor optimization on an arbitrary choice of orthogonality conditions made during optimization.We do not currently know of a way to avoid this source of error.Fortunately, as illustrated in Fig 14, we have good evidence that it is small.
The other source of error comes from interaction effects, which are not captured properly by the twoparticle states.As discussed above, one can choose the minimum separation of the two particles to avoid interaction, by throwing out position states where the "interaction strength" ϵ, defined in (C22), rises above some threshold.In choosing the threshold, there is a tradeoff between capturing (potentially large) components of the wavefunction that have smaller separation, but likely incur some (possibly small) error due to interaction, and the magnitude of ϵ, which is a conservative estimate of that error and is exponential in the separation.How to make this tradeoff optimally depends on the target wavefunction and the required precision of the projected wavefunction.
For our computations, we examined the dependence of the norm of the projected wavefunctions on the minimum separation d min allowed in the two-particle basis states.Exemplary results are shown in Figs.18 and 19, which also show the dependence of the norms on the simulation bond-dimension limit.In Fig. 18, we observe that, for the largest bond-dimension, the norm of the projected wavefunction is essentially constant for d < 70, despite the rising magnitude of interaction effects.This suggests that the wavefunction has negligible support at small separations (which we confirm via a Fourier analysis).We also observe that a minimum separation of d min = 60 is sufficient to keep ϵ < 10 −6 .Noting that the interaction strength is a property of the quasiparticle basis, and hence independent of time, we make the choice d min = 60 to avoid interaction effects when computing the data shown in Fig. 8 of the main text.
Fig. 19 shows similar data for three quasiparticle-pair "sectors" for the Hamiltonian parameters (ii) of the main text.This data was used to estimate the scattering outcome probabilities shown in Fig. 10 of the main text.We see that, for this target wavefunction |ψ(t)⟩, interaction effects are not important for the κκ (0,0) and µµ (0,0) "sectors".However, they may influence the result at the level of ∼ 0.01, possibly more, for κκ (0,1) (and, by reflection symmetry of the initial bubble state, κκ (1,0) ).The κκ (0,1) and κκ (1,0) results are likely to be more sensitive to interaction than the κκ (0,0) result because the former two "sectors" describe bubble states in which one of the walls (the kink or the antikink) is heavier than in the κκ (0,0) case.These "lopsided" bubbles will be smaller than a κκ (0,0) bubble at the same energy, leading to larger components of the wavefunction at small separations, where interaction effects are stronger.

Appendix E: Evolving through time
To evolve an initial bubble iMPS in time, we define a window of N w lattice sites surrounding the bubble and allow the MPS tensors belonging to those sites to vary during the evolution, while keeping the rest of the MPS tensors fixed.In other words, we use the nonuniform window ansatz (A3) with fixed bulk tensors A L , A R .To compute the evolved state, we use methods based on the Time-Dependent Variational Principle (TDVP), which is set out for this class of states in [76] and implemented in the evoMPS package [80].
The TDVP provides flow equations that describe the evolution of the MPS tensors needed to optimally approximate the evolution of the state by the Hamiltonian, given the constraint that the MPS bond dimension must remain fixed.Various schemes can be used to integrate the flow equations: we use the popular Runge-Kutta 4/5 (RK4) numerical integrator (to directly integrate the global flow equations) as well as the "projector-splitting" (PS) integrator of [79].Since we want the MPS bond dimension to grow as needed (up to some maximum) as the entanglement of the state increases, we combine TDVP flow with techniques for increasing the bond dimension.In particular, we use the "dynamical expansion" scheme described in [73] together with the RK4 integrator, as well as the two-site projector-splitting method of [79].
The PS method and the RK4 integrator (with dynamical expansion), despite having similar theoretical error rates for a given time-step size, behave differently in important ways.For a given step size, the PS method has a larger computational overhead per step, but has better numerical stability and precision since, unlike the "traditional" TDVP scheme of [73], it does not require the inversion of matrices with small eigenvalues.
We find that the RK4 scheme with dynamical expansion is too unstable to use reliably during the initial timesteps of our simulations, which begin with an MPS of relatively small bond dimension (at most twice the vacuum bond dimension).However, we find RK4 can be used successfully after performing a small number of initial steps using the two-site PS scheme.During these initial PS steps, the bond dimension increases significantly.Later in the evolution, once the bond dimension has stabilized, we find the much faster RK4 scheme is able to take over without significant impact on the results.During the evolution, we do allow the MPS bond dimension to grow beyond a predefined maximum value.
To better understand the effects of the integration scheme on our simulations, as well as the impact of the bond-dimension limit, we compute two quantities indicative of numerical error: the energy drift and the truncation error.Although the exact evolution of the quantum state conserves the energy, the imperfect integration of the TDVP flow equations, combined with the limited bond dimension, leads to a small energy drift.This drift is an indicator of error incurred more generally during the evolution.We estimate the truncation error -the portion of the state by norm that is lost due to the bond dimension limit -as the maximum value, taken over position, of the minimum Schmidt coefficient (the Dth-largest) for the left-right bipartition at that position.Since there are rarely large jumps in the Schmidt spectrum, this value provides a good estimate of the magnitude of the terms that cannot be represented due to the bond-dimension limit.
As shown for simulation (ii) of the main text in Fig. 20, the energy drift is particularly sensitive to the integration scheme and step size, whereas the truncation error is, unsurprisingly, most sensitive to the bond-dimension limit.We note that the truncation error jumps at the time of the first kink-antikink collision (t ≈ 250), consistent with the entanglement jump observed in Fig. 11 of the main text.Fig. 21 provides a more detailed picture of the entanglement structure, showing the full entanglement (Schmidt) spectrum (up to truncation) before and after the first collision at the cut with the largest entanglement entropy.
The computational cost of simulating up to some fixed time t scales as O D 3 δt .We are therefore eventually forced to trade accuracy for computational cost.For this particular simulation, we judge D ≤ 128 and δt = 0.05 to provide sufficient accuracy, while keeping the computational requirements manageable, to enable us to study the outcomes of at least the first kink-antikink collision event in detail.

Appendix F: Comparison with quench approaches
We construct our initial false-vacuum bubble from individual kink and antikink quasiparticles, separated by a region of metastable false vacuum.Similar states can be constructed via a simpler approach: act on the uniform true vacuum |Ω⟩ with a suitable string operator that, in the case of the Ising-type model we study, flips all the spins over a range of sites: For small longitudinal field h, we know that flipping all the spins gets us from the vacuum to a state close to the false vacuum (see App. B 2), so if d := k −j is sufficiently large, the reduced state will be close to that of the false vacuum in the middle of the flipped region and topological excitations will be created at the edges of the string.In general, these excitations will be a combination of many topological quasiparticles of varying energy, hence the walls of a bubble created in this way will be unstable to interactions between these quasiparticles.The walls are also highly localized in position and hence have large momentum uncertainty.By smearing the edges of the string in space, their momenta can be focused.This results in states of the where f j and f k are wavepacket functions for the left and right edge, respectively.If we choose these to be Gaussian, as for the wavepackets in the main text, we get a bubble state similar to those used in the main text, except that the walls have undetermined quasiparticle content.
In Fig. 22, we compare the dynamics of three different initial states for the Hamiltonian parameters (ii) of the main text: (a) the fully localized string of (F1), (b) the smeared string of (F2), and (c) the tuned quasiparticle kink-antikink wavepackets used in the main text.We choose the same kink-antikink separation in all cases, and the same wavepacket widths in the latter two.The TDVP step size and the maximum bond dimension were also the same in all three cases.Simulations (a) and (b) both exhibit clear ballistic spread of energy from the initial bubble edges, indicating their instability, whereas simulation (c) only shows ballistic spread after the initial bubble walls have collided, consistent with the walls consisting of individual quasiparticles.
It is noteworthy that the evolution of state (a) encounters catastrophic numerical errors at around t = 240, unlike simulations (b) and (c), suggesting that scenario (a) is much harder to simulate accurately 15 .In Fig. 23 we show the energy drift and estimated truncation error (smallest retained Schmidt coefficient for the most entangled cut) for the same three simulations.This data clearly indicates that simulation (c) is easiest to simulate, since both the energy drift and the truncation error remain smaller with the same evolution parameters.The large truncation error at early times in cases (a) and (b) is consistent with a large amount of entanglement being generated early on, coming from the multiple excitations created at the ends of the string operator.

Appendix G: Velocity and Bloch oscillations 1. Single-kink evolution
In case of explicit symmetry breaking (h > 0), an isolated kink wavepacket, with true vacuum on one side and false vacuum on the other, will accelerate toward the false vacuum, absorbing the excess energy density of the latter.On the lattice, however, the kinetic energy cannot increase indefinitely.Instead, as discussed in the main text, the kink begins to undergo Bloch oscillations, eventually decelerating and reversing its direction of travel, as shown in velocity derived from the momentum (via the dispersion relation).

Bubble evolution
In the case of a bubble state, the initial kink and antikink behave as their isolated counterparts, accelerating into the false vacuum until they near each other and interact.Their momenta increase linearly until the collision, as shown in lattice units (−π < p ≤ π) for the Ising model (parameter set (i) of the main text) in Fig. 24.It is interesting to note that the momentum variance is significantly larger after the collision than it is before, indicative of an (in this case elastic) interaction.
The velocity of the kink and antikink, defined here as the velocity of the point in space at which the (inter-  polated) spin expectation value crosses zero, evolve as shown in Fig. 25, in accordance with the dispersion relation of the kink quasiparticle excitation.In this simulation, the kink achieves its maximum velocity well before the collision, and begins to decelerate as part of a Bloch oscillation.The pre-collision velocity can be kept from reaching its maximum by reducing the kink-antikink separation in the initial bubble state, as we have confirmed with other simulations.In Fig. 26, we show the kink and antikink velocity evolution for simulation (ii) of the main text.In this case, we have set the initial kink-antikink separation so that the velocity does not reach a maximum prior to collision (indicating that Bloch oscillations have not yet begun).be either above or below the threshold for quasiparticle pair production, which we numerically estimate to be 2m µ = 1.88 (relative to the vacuum energy).In Fig. 27, we show the cut entropy as a function of space and time and, separately for clarity, the timedependence of the cut entropy at the midpoint between the quasiparticle wavepackets.We choose three different initial momenta: p = 4π 32 , p = 5π 32 , and p = 6π 32 , corresponding to bubble energies of E = 1.65,E = 1.92, and E = 2.20, respectively.We observe that the postcollision mid-chain entropy returns to its vacuum value for p = 4π 32 , suggesting a trivial scattering event (see App. I).For p = 5π 32 and p = 6π 32 , we observe a residual entropy surplus after the collision, suggesting nontrivial scattering.Since the onset of this extra entropy contri- The interpretation of the 0-intercept as the position breaks down both during and, to some extent, after the collision: During the collision, the zero intercept disappears altogether as the kink and antikink merge and all spin expectation values are > 0. After the collision, there are in this case (see Fig. 10) at least two different bubble "branches" of the wavefunction, both contributing to the spin expectation values.The onset of the first collision is indicated by the dotted line, which is the time at which the maximum cut entropy begins to grow rapidly.
bution coincides with the energy crossing the two-meson threshold 2m µ , it is likely due to an increasing probability of meson pair production.It is interesting to note that, if we set h > 0 while keeping the other Hamiltonian parameters the same, we observe nontrivial, albeit elastic, scattering of kinkantikink pairs when the energy is below the two-meson threshold.Turning off the longitudinal field appears to turn off this nontrivial elastic contribution, so that kinks and antikinks scatter trivially.We assume the initial state factorizes as the product ψ(p, q) = ψ A (p)ψ B (q) of two widely separated wavepackets, but after elastic scattering, the wave packets becomes correlated due to the momentum dependent phase shift ϕ(p, q): ψ(p, q) = ψ A (p)ψ B (q)e iϕ (p,q) .(I3) Tracing out particle B, we obtain the reduced density matrix for particle A ρ A = ˆdp 1 dp 2 |p 1 ⟩ ρ A (p 1 , p 2 ) ⟨p 2 | , ρ A (p 1 , p 2 ) = ˆdq ψ (p 1 , q) ψ * (p 2 , q) , (I4) where * denotes complex conjugation.
To quantify the entanglement of particles A and B, we compute the Rényi entropies of ρ A , Now suppose that the wave packets for particles A and B are Gaussian: If the phase shift were slowly varying over the range in p and q where the wave packets have significant support, we could approximate trρ n A by expanding ϕ(p, q) to quadratic order about (p, q).But in that case the scattered wave packets are only slightly entangled.In order to do an analytic computation, we will assume that ϕ(p, q) is exactly quadratic even if the phase shift Figure 1.Cartoon illustrating the relaxation of a falsevacuum bubble in a spin chain.The magnetization ⟨Z⟩ is positive in the true vacuum, but negative in the false vacuum.A bubble-wall collision is a scattering process, which may be (a) free (no interaction), (b) elastic (no particle production), or (c) inelastic (particle production).Note that in free and elastic scattering of topological particles the left (right) particle always remains a kink (antikink).

Figure 2 .
Figure 2. Evolution of the excess energy density e (relative to vacuum), as a fraction of total excess energy E, in a spin chain for two initial states: (a) created by applying a spatially smeared string operator to the vacuum and (b) constructed from MPS tensors to contain kink and antikink quasiparticle wavepackets.In (a) meson pairs are produced immediately at the string edges, whereas in (b) there is no particle production until the initial kink and antikink collide.The dynamics are restricted to a window of ∼ 1000 sites, leading to boundary effects in (a).For more details, see App.F.

Figure 3 .
Figure 3. Partial phase diagrams of the extended Ising chain[51,52], both without and with a small longitudinal field h that breaks the Z2 symmetry.At h = 0, g = 1 there is a continuous (symmetry-breaking) phase transition described by the Ising CFT for λ ≲ 0.428, and by the Tri-Critical Ising (TCI) CFT at λ ≈ 0.428.By studying the spin chain near to these transitions (g → 1, h → 0), we can access emergent, relativistic quantum field theories with confined kinks[53,54].Points (i), (ii), and (iii) correspond to the data shown in the figures below.

Figure 4 .
Figure 4. Diagram illustrating the various types of states, and their spin profiles, relevant for simulations.For example, our initial states are wavepackets constructed from kink-antikink pairs.The product states listed are eigenstates of H when g = 0, λ = 0. Away from this regime, we use MPS to accurately capture fluctuations in the vacua, as well as the quasilocal nature of the excited quasiparticle states.

Figure 8 .
Figure 8. Portion of state (by probability) outside of the MPS kink-antikink subspace κκ (0,0) for simulations (i) and (ii) of Fig.7.Here we fully account for momentum dependence of the basis states |κκ j,k ⟩ via a Fourier analysis and count only contributions with k − j ≥ 60 (see App. D).For Ising (i), the small probability after the first collision of t ≈ 140 indicates elastic scattering of kinks, in stark contrast with the TCI case (ii), where the probability remains high after the first collision at t ≈ 250.In (i), the growth of the post-collision probability with subsequent collisions is consistent with delocalization of the wavepackets, since contributions from kink-antikink pairs with k − j < 60 are not counted.Hence the larger values for (i) after the second collision at t ≈ 430 should not be taken as evidence of inelastic scattering.

Figure 9 .
Figure 9. Dispersion relations (numerical, using MPS) of kinks κ and mesons µ for λ = 0.41, g = 0.98 in lattice units for the Hamiltonian(1).Momentum ranges from −π to π.For mesons, energies are shown with and without a weak longitudinal field.Individual kinks do not have a finite energy for h > 0. Threshold energies for pair production are shown (computed assuming h = 0 for kinks and h = 0.001 for mesons), as is the energy (labelled Ψ) of the simulation shown in Fig.7for parameter-set (ii).The dispersion relations reach their maximum gradient at momentum p(vmax) ≈ 0.6.

Figure 11 .
Figure 11.Entanglement entropy (base 2) for cuts (left-right bipartitions) of the spin chain as a function of time for the simulations (i) and (ii) of Fig. 7. Convergence with the bond dimension D slows as time goes on.For example, in (i) the max.cut entropy at D ≤ 128 is very likely not converged after t ≈ 700.

Figure 12 .
Figure 12.Peak maximum cut entropy during the first collision as a function of energy (ii) close to (λ = 0.41 g = 0.98, h = 0.001) and (iii) further from the TCI point (λ = 0.3, g = 0.9, h = 0.0069, mµ ≈ 0.97, vmax ≈ 1.43), with the energy controlled by the initial kink-antikink separation.The vacuum correlation length ξ is 3.6 sites for (ii) and 1.8 sites for (iii), indicating that (ii) is closer to criticality.The initial wavepacket width is σ = 40 for (ii) and σ = 19 for (iii).The kink velocity at the start of the first entropy jump (see Fig.11), normalized so that the maximum is 1, is also shown.Decreasing velocity with energy indicates the onset of Bloch oscillations.Simulation (ii) of Fig.7corresponds to the leftmost point of the black curves.

Figure 13 .
Figure 13.Error made (1−⟨ϕj(B)|ϕj(B(p))⟩) in ignoring the momentum dependence of the tensor B used to construct MPS quasiparticle states, for both kink and meson excitafor the Hamiltonian parameters used in simulations (i) and (ii) of the main text.The momentum-eigenstate freedom on B(p) is fixed so that ⟨ϕj(B)|ϕ k (B ′ )⟩ = δ jk .

Figure 15 .
Figure 15.Spin expectation values of a kink position MPS |κj⟩ for (i) the Ising model and (ii) close to the Tri-Critical-Ising (TCI) point.The bond dimension is D = 8 for the Ising data, and D = 18 for the TCI data.We plot the spins for various ways of fixing the momentum-eigenstate freedom: the left and right orthogonal conditions, (C4) or (C6), and the reflection-symmetric conditions (C17), beginning from a B tensor optimized using either the left or right conditions (since this makes a small physical difference to the resultsee App.C 1).
(a,b) jk |ψ(t)⟩ is sensitive to the presence of a pair consisting of a type-a kink quasiparticle at position j and a type-b antikink quasiparticle at position k.
jk ⟩ of (C19) (where ϕϕ is generic notation for either a kink-antikink pair or a meson pair), we can achieve orthogonality by enforcing the left orthogonality conditions on the left excitation tensor, B L , and the right orthogonality conditions on right excitation tensor, B R .With these conditions, the twoparticle states are orthogonal for one pair of species a, b ⟨ϕϕ (a,b) jk |ϕϕ (a,b) lm ⟩ = δ jl δ km .(D1) Furthermore, by limiting the subspace to states with separation d := k − j large enough so that interaction effects are negligible (see App. C 4), we can achieve an approximately orthonormal basis across species as well as positions: ⟨ϕϕ (a,b) jk |ϕϕ (c,d) lm ⟩ ≈ δ jl δ km δ ac δ bd ∀ k ≫ j, m ≫ l. (D2) If the left and right bulk vacuum tensors, A L and A R of the evolved state |ψ(t)⟩, which in our simulations have the form (A3), match the left and right bulk tensors

Figure 18 .
Figure 18.Portion of the evolved bubble wavefunction outside of the kink-antikink "sector" for simulation (i) of Fig.7(λ = 0, g = 0.8, h = 0.007) at time t = 240, as a function of the minimum separation dmin := min(k − j) permitted in the two-quasiparticles basis states |κκ (0,0) j,k ⟩.Probabilities are computed via a Fourier analysis, taking into account the momentum-dependence of the basis states.The basis error due to interaction effects is estimated using (C22).

Figure 19 .
Figure 19.Quasiparticle "sector" projection probabilities for simulation (ii) of Fig.7(λ = 0.41, g = 0.98, h = 0.001) at time t ≈ 426, as a function of the minimum separation dmin := min(k − j) permitted in the two-quasiparticles basis states | • • j,k ⟩.Probabilities are computed via a Fourier analysis, taking into account the momentum-dependence of the basis states.The basis error due to interaction effects is estimated using (C22).

Figure 21 .
Figure 21.Schmidt spectra for the maximum-entropy cut before and after the first collision in simulation (ii) of the main text.The bond dimension is 128. form

Fig 5 forFigure 22 .
Figure 22.Evolution of spin and energy density expectation values and the cut entropy for parameter set (ii) of the main text (λ = 0.41, g = 0.98, h = 0.001), with three different initial states.State (a) is prepared by acting on the vacuum with a string operator x R −1 j=x L Xj, which flips the spins to form a bubble-like state with energy E/mµ = 8.69.State (b) is similar to (a), but with the ends of the string smeared out using Gaussian packets of width σ = 40, reducing the energy to E/mµ = 4.02.State (c) is the initial state discussed in the main text with E/mµ = 2.62, using quasiparticle wavepackets for the kinks and the false vacuum for the middle region.The evolution parameters are the same in all cases: The maximum bond dimension is 64 and the RK4 step size is ≈ 0.08 (0.05 in unscaled lattice Hamiltonian units).In (a), dramatic errors in the simulation occur at t ≈ 240, indicating the difficulty of simulating these dynamics versus (b) and (c).In both (a) and (b), ballistic energy-spread emanating from the initial kinks indicates that they have complex quasiparticle content, resulting in immediate inelastic scattering.In contrast, the tuned quasiparticle kinks of (c) do not produce appreciable ballistic spread until the bubble walls have collided.

Figure 23 .
Figure 23.Evolution of the energy expectation value and MPS truncation error for the simulations of Fig. 22.Energy drift (|1−E(t)/E(0)|) indicates deviation from unitary evolution and results from restriction to a maximum bond dimension of 64 as well as from numerical integration errors (RK4 step size ≈ 0.08).Truncation error (estimated as the maximum over cuts of the smallest Schmidt coefficient) results from the limited bond dimension and increases as entanglement is produced.

Figure 24 .Figure 25 .
Figure24.Momentum of the kink in simulation (i) of the main text (Ising), computed from the projected wavefunction in the κκ (0,0) basis.The bond dimension is 128.Since the κκ (0,0) basis does not accurately capture the state when the kink is very close to the antikink, this data is not complete during the collision (t ≈ 140).The error bars indicate the standard deviation of the momentum distribution.

Figure 26 .
Figure 26.Kink and antikink wavepacket velocity for simulation (ii), computed as the finite difference of the interpolated position of the 0-intercept of ⟨Zj⟩.The bond dimension is 128.The interpretation of the 0-intercept as the position breaks down both during and, to some extent, after the collision: During the collision, the zero intercept disappears altogether as the kink and antikink merge and all spin expectation values are > 0. After the collision, there are in this case (see Fig.10) at least two different bubble "branches" of the wavefunction, both contributing to the spin expectation values.The onset of the first collision is indicated by the dotted line, which is the time at which the maximum cut entropy begins to grow rapidly.

Figure 27 .
Figure 27.Cut entropy for kink-antikink collisions, in the absence of a longitudinal field, with initial kink momentum p and antikink momentum −p, for three different values of p.The Hamiltonian parameters are g = 0.9, λ = 0.3, h = 0, and the wavepacket width is σ = 19.0.The vacuum bond dimension is D = 14, with a limit D ≤ 64 imposed during evolution.The integration time-step size was δt = 0.05 in unscaled lattice Hamiltonian units.

log 2
varies rapidly.Then the only term that matters is the cross termϕ(p, q) = ϕ 2 pq + . . ., (I8)because the exponential of the other terms factorizes into a function of p times a function of q, which does not contribute to the entanglement of particles A and B. By evaluating a Gaussian integral, we find trρ n A = 1 + 4α 2 sin 2 kπ n , (I9) Figure 20.Evolution of the energy expectation value and MPS truncation error for simulation (ii) of the main text.We compare different maximum bond dimensions D as well as different RK4 time-step sizes δt (shown here in unscaled lattice Hamiltonian units, for D ≤ 64), observing that the effects of the time step are most noticeable in the energy drift, whereas the bond-dimension most obviously affects the truncation error at around the time of the first collision (t ≈ 240).The D ≤ 128 simulation is initialized from D ≤ 64, δt = 0.05 at t ≈ 142.The uptick in the energy drift at t ≈ 750 is due to unconfined wavepackets hitting the boundaries of the simulation window.