Entanglement generation in $(1+1)D$ QED scattering processes

We study real-time meson-meson scattering processes in $(1+1)$-dimensional QED by means of Tensor Networks. We prepare initial meson wave packets with given momentum and position introducing an approximation based on the free fermions model. Then, we compute the dynamics of two initially separated colliding mesons, observing a rich phenomenology as the interaction strength and the initial states are varied in the weak and intermediate coupling regimes. Finally, we consider elastic collisions and measure some scattering amplitudes as well as the entanglement generated by the process. Remarkably, we identify two different regimes for the asymptotic entanglement between the outgoing mesons: it is perturbatively small below a threshold coupling, past which its growth as a function of the coupling abruptly accelerates.

We study real-time meson-meson scattering processes in (1 + 1)-dimensional QED by means of Tensor Networks.We prepare initial meson wave packets with given momentum and position introducing an approximation based on the free fermions model.Then, we compute the dynamics of two initially separated colliding mesons, observing a rich phenomenology as the interaction strength and the initial states are varied in the weak and intermediate coupling regimes.Finally, we consider elastic collisions and measure some scattering amplitudes as well as the entanglement generated by the process.Remarkably, we identify two different regimes for the asymptotic entanglement between the outgoing mesons: it is perturbatively small below a threshold coupling, past which its growth as a function of the coupling abruptly accelerates.
The investigation of fundamental interactions is one of the more challenging research fields in physics from a theoretical, experimental and computational point of view.The Standard Model (SM) of particles physics describes the fundamental components of matter as quantum fields, whose interactions are set by the system invariance under specific gauge transformation groups [1].Many unveiled phenomena, such as the matter-antimatter asymmetry, the origin of dark matter and dark energy, as well as a complete understanding of the Quantum Chromodynamics (QCD) phase diagram, motivate intense research for physics in and beyond the SM [2].In particle accelerators, scattering events are extensively exploited to investigate the properties of fundamental particles and their interactions.In this scenario, numerical simulations play a fundamental role to compare experimental results and theoretical predictions.The most common approach for the numerical simulation of quantum gauge theories relies on Lattice Gauge Theory (LGT) models [3,4], in which quantum fields are defined on a discrete space-time lattice.Within this framework, Monte Carlo techniques are the established tool to predict scattering amplitudes [5].Even though Monte Carlo simulations of LGTs have been able to correctly interpret an impressive number of experimental observations, they suffer from the notorious sign problem, that severely limits the study of the finite-density region of the QCD phase diagram and the real-time dynamics [6][7][8][9].In the last decades Tensor Network (TN) methods [10][11][12][13][14][15], thanks to their capability to efficiently compress the information contained in quantum many-body states [16][17][18][19][20][21][22][23][24][25][26], have emerged as a sign-problem free tool to perform numerical simulations of fermionic LGT models [27][28][29].TNs have been successfully applied to study (1 + 1)-dimensional LGTs, characterising the phase diagram and the dynamics in Abelian [30][31][32][33][34][35][36][37] and non-Abelian [38][39][40]  The ground state of the Schwinger model is computed and expressed as an MPS (blue bottom squares).Then a set of MPOs (red squares) is applied to prepare a gauge invariant initial state made up of two spatially separated meson wave packets with opposite electric field orientations.Indexes i and j indicate lattice sites on which matter and antimatter field components are defined, while the electric field is defined on the links between sites.The two mesons evolve under the Hamiltonian of Eq. ( 4) with TEBD (pictorially represented by yellow boxes) and their scattering is studied.Note the persistent entanglement between them after the scattering process.Numerical parameters: a = 1, m = 1.0, g = 0.14.
Here, we perform TN simulations by using Matrix Product States (MPS) to investigate the real-time dynamics of scattering events in a (1 + 1)-dimensional Abelian LGT in the Hamiltonian formulation [49,50].In particular, we set up a protocol for the initialization of wave packets of asymptotic isolated states and observe their scattering dynamics.We focus on the Schwinger model, i.e.QED 2 , which describes electrons and positrons interacting on a 1D lattice via a scalar electric field [49].Since the Schwinger model does not admit stable isolated-charge excitations [51], we compute the scattering between two meson wave packets, each of them composed by an electron and a positron, sufficiently separated in space to be initially uncorrelated, as pictorially shown in Fig. 1.In order to prepare the initial state for the dynamics, we compute the ground state of the model Hamiltonian by using the Density Matrix Renormalization Group (DMRG) algorithm with MPS representation [10,[52][53][54].Then, we apply a set of Matrix Product Operators (MPO) which create the mesons [55].This approach allows us to fix the space and momentum centerof-mass position for each meson independently, as well as their widths.After the preparation of the initial state, the dynamics is computed via the Time Evolving Block Decimation (TEBD) algorithm [56,57].We analyze the scattering products in terms of the final momenta, by computing the structure factor of the final state and some scattering amplitudes, in order to better characterize the final state and compare it to the initial one.Finally, we focus on the entanglement production during the process: the initial entanglement between the two mesons is zero before the scattering event, as expected.After the scattering, we observe two very well-separated behaviors as the Hamiltonian parameters are varied: one in which the entanglement at long times is perturbatively small with respect to the electric field coupling, and another in which a significant amount of entanglement persists asymptotically in time.
Our work paves the way to future insights into the role played by entanglement in scattering processes in LGTs, extending previous analytical studies and numerical results obtained for spin models [58][59][60][61][62][63].Moreover, TNs provide the ideal framework and language to make a link between our protocol and experimental quantum simulations and computations of LGTs.Indeed, motivated by the improvements in experimental techniques for the manipulation of isolated quantum many-body systems [64][65][66][67][68][69][70][71], a plethora of analog and digital quantum simulators proposals has come out to implement Abelian and non-Abelian LGTs into experimental setups [72][73][74][75][76].In this framework, our TNs protocol could also serve as a benchmarking toolbox for realistic scattering implementations on quantum hardware.
The manuscript is organized as follows.In Section I we illustrate the models involved in our simulationsnamely, free lattice fermions and QED 2 .Section II contains an overview of the generic tools employed in our investigations.In Section III we present the protocol that we use to model the initial state of a scattering experiment and we introduce operators that prepare wave packets of fermion and antifermions (free theory) and mesons (QED 2 ).In Sections IV and V we analyze the simulations of some meson-meson collisions.Section IV focuses on the collision phenomenology and includes the computation of some scattering amplitudes.Section V focuses on the entanglement content of the system during the scattering process and includes an estimation of the entanglement between the final products that is generated by the interactions.

I. MODELS
We adopt units where ℏ = c = 1 and consider a discretization of one-dimensional space as a uniform chain Λ of an even number L ∈ 2Z of sites separated by lattice spacing a, Λ = a {0, . . ., L − 1}.Lattice positions are denoted by x, y, z ∈ Λ; lattice momenta are k, p, q ∈ Λ * , We now give a possible lattice definition of the (1 + 1)dimensional theories of free relativistic Dirac fermions and of quantum electrodynamics; the latter being the main subject of this work.

A. Free fermions
The Kogut-Susskind discretization of a mass m free relativistic Dirac field in (1 + 1)-spacetime dimensions is described by the Hamiltonian [49,77,78] The staggered fermion ξ x degrees of freedom on even and odd lattice sites correspond respectively to the upper and lower components of the (1 + 1)-dimensional Dirac spinor field; hence the alternating sign in front of the mass term in Eq. ( 2) [78].The fields ξ x and ξ † x satisfy canonical anticommutation relations [77]: while other fundamental anticommutators vanish.The Jordan-Wigner transformation [78,79] provides an irreducible matrix representation of this algebra in the local occupation number As its continuum counterpart, the Hamiltonian in Eq. ( 2) has a global U(1) symmetry, generated by the particle number conserved charge Q = a x ξ † x ξ x [78].In the thermodynamic limit or for periodic boundaries the theory is also translation invariant.Valid translations, however, are generated by shifts by two lattice steps [78].As a consequence, the effective momentum space of staggered fermions is Λ * ′ = Λ * ∩ [−π/2a, π/2a[; we denote momentum sums restricted to this sublattice by ′ .Moreover, it is sometimes convenient to isolate the even and odd sublattices, namely E = Λ ∩ a(2Z) and O = Λ ∩ a(2Z + 1).

B. QED
A discretization of (1 + 1)-dimensional Quantum Electrodynamics (QED 2 , also known as massive Schwinger model [80,81]), is obtained promoting the global U(1) symmetry of Eq. ( 2) to a gauge symmetry.In the Kogut-Susskind formalism, the theory is defined by the Hamiltonian [49,77,78] Equation ( 4) involves matter degrees of freedom, ξ x and ξ † x , living on the lattice sites; as well as unitary gauge parallel transporters U x, x+a and Hermitian electric field gE x,x+a operators, acting on local Hilbert spaces associated to the links between neighboring sites [77].Both the bare mass m and coupling g parameters have the dimension of a mass, thus the ratio g/m is adimensional and can be used to quantify the strength of the interaction and interpolate between the weak (g ≪ m) and strong (g ≫ m) coupling limits [82].In addition to the matter anticommutator from Eq. (3), we have the non-vanishing fundamental commutator [77] [U y,y+a , E z, z+a ] = δ yz U y, y+a , showing that parallel transporters act as lowering operators for the electric field on the same lattice link.An irreducible representation of Eq. ( 5) on a given link (left implicit) reads [77,83] where σ(E) denotes the spectrum of the operator E.
There exist infinite other unitarily inequivalent representations related to Eq. ( 6) by a shift of the electric field: Here we assume 0 ∈ σ(E) on all links.The representation in Eq. ( 6) extends trivially to the whole chain.Time independent U(1) gauge transformations are implemented via exp(−ia x θ x G x ), with θ x ∈ R L parametrizing the local transformation generated by the Gauss operator G x [84,85], which commutes with the Hamiltonian.As a boundary condition, we identify the electric field at the left and at the right of the chain with the null operator, thus restricting our analysis to the sector with total charge zero.Physical states |Ψ phys ⟩ ∈ H phys ⊂ H are required to be gauge invariant [49,85], namely to satisfy The physical state condition is equivalent to Gauss law, i.e., aG x |Ψ phys ⟩ = 0 for all x.
In 1 + 1 dimensions and with E −a,0 ≡ 0, given an arbitrary configuration N x of the matter fields there is one and only one configuration E x, x+a of the link degrees of freedom complying with Gauss law; precisely where H free is the free staggered fermion Hamiltonian from Eq. (2) [86].With this procedure the link degrees of freedom have been removed -namely, integrated outand the gauge redundancy of the model has been eliminated.
The reformulation of the theory in Eq. ( 10) is convenient for some derivations carried out in this work.In numerical simulations, instead, we keep the (redundant) link degrees of freedom to avoid long-range interactions.Their infinite dimensional local Hilbert spaces are truncated introducing a cutoff E max ∈ N in the |E| spectrum and imposing U |−E max ⟩ = |+E max ⟩ on each link; in this way σ(E) is identified with Z n , n = 2E max + 1.This truncation spoils the commutator in Eq. ( 5), explicitly breaking the U(1) gauge invariance of the model down to the Z n residual symmetry group of transformations with parameters θ x ∈ (2π/n)Z [84,[87][88][89][90].The physical state condition in Eq. ( 8) is accordingly weakened and the Gauss law only holds modulo n.Namely, the physical Hilbert subspace H phys is spanned by the occupation number and electric field eigenstates satisfying aG y |N x , E x, x+a ⟩ ∈ (nZ) |N x , E x, x+a ⟩ for all y.Our simulations rely on an n = 7 truncation.It has been shown that n = 3 gives already an excellent approximation of the exact ground state of QED 2 [84,91,92]; even though the Z 3 model is also capable of reproducing accurately some dynamical processes of the untruncated theory, the quality of the approximation depends on the value of the model parameters and the specific process under consideration [93,94].Monitoring the system state |Ψ⟩ during Z n (n = 3, 5, 7) scattering simulations, we find that the fraction χ trunc of configurations affected by truncations of the electric field is below the numerical precision of the simulation for n = 7. Precisely, χ trunc = 1 − ⟨Ψ|Π U(1) |Ψ⟩ ≤ O(10 −10 ), where Π U(1) is the projector on the Hilbert subspace of states complying with the U(1) Gauss law.

II. METHODS
In this work, Tensor Network (TN) methods [12][13][14][15]95] are employed to tackle the exponential growth of the many body Hilbert space of the lattice theories introduced in Section I. Specifically, we use the Matrix Product State (MPS) [96] and Matrix Product Operator (MPO) [55] ansätze to represent states and operators.In conjunction, we use the MPS implementation of (i) the (two-site) Density Matrix Renormalization Group (DMRG) algorithm to perform variational optimizations; and (ii) the (fourth-order) Time Evolving Block Decimation (TEBD) algorithm to compute Trotterized time evolutions [10, 52-54, 57, 97].In order to recast the QED 2 Hamiltonian of Eq. ( 4) in a nearest-neighbor Hamiltonian, as required by TEBD, we fuse the Hilbert spaces associated to a site H x and to the subsequent link H x,x+a in a single (larger) local computational Hilbert space Finally, we exploit some of the symmetries of the problem to improve the efficiency and the numerical precision of our simulations.Using a symmetric MPS ansatz [14,15,28], we constrain the ground state search and the dynamics in the desired charge sector.Specifically, we impose -exactly -the conservation of the U(1) global charge Q and, for the QED 2 case, also the vanishing of the total Z n Gauss charge on even and odd sublattices , i.e., a x∈E G x and a x∈O G x [98,99].
Schematically, the presented results are obtained as follows (as depicted in Fig. 1 For QED 2 we also measure the electric field g ⟨E x,x+a ⟩ t and its energy density h elec , Finally, we characterize the entanglement content of the system by computing the Von Neumann entanglement entropy S associated to every bipartition of the chain in two subsystems, L = {y < x} and R = {y ≥ x}; namely, where ρ(t) is the reduced density matrix of one subsystem, e.g., the partial trace of |Ψ(t)⟩⟨Ψ(t)| over y∈L H y .
Unless otherwise stated, we report the deviation from the vacuum (ground state) expectation value of all the above quantities.Moreover, we use the lattice spacing to fix the overall length scale, i.e., we set a = 1, and all simulations are carried out with open boundary conditions.
Numerical errors have been kept under control.The MPS compression is achieved truncating singular values below 10 −6 , resulting in an overall maximum MPS bond dimension of 991; the Trotterization of the evolution (we exploit a fourth-order Suzuki-Trotter decomposition [100] with time step δt = 0.05) is responsible for an overall O(10 −3 ) error on the final state; while the truncation of the electric field spectrum (|gE| ≤ 7g) is discussed in Section I.The convergence of some relevant quantitiese.g., of the midchain entanglement entropy -with the MPS bond dimension has been verified.

III. INITIAL STATE PREPARATION
The state of a system undergoing a scattering process long before the collision is populated with wave packets of stable particle excitations localized in far distant space regions [101,102].A prerequisite for the description of a scattering experiment is thus the identification of the model's particles and their asymptotic-times dynamics [103,104].In the absence of bound states and long range interactions, the aforementioned task consists in the computation of the plane wave solutions of the free theory associated to the degrees of freedom of the model [104].The previous assumptions are trivially verified by a free theory, e.g.free fermions, but do not hold for QED 2 [51,81,105,106].
In this Section we present and motivate our protocol for preparing the initial state of a QED 2 scattering simulation, providing an expression for the wave packet amplitudes and operators involved.An MPO representation of these operators is given in Appendix B. We consider initial states of uncorrelated particles -i.e., particle wave packets described by independent amplitudes.We start by illustrating how localized particle excitations are prepared in the lattice theory of free staggered fermions, as it provides useful insights into the QED 2 case.

A. Free fermions
In Appendix A we solve exactly the theory of free staggered fermions, presenting the Fock space structure of the Hilbert space.Periodic lattice boundary conditions are assumed in the derivation in order to neglect open boundary effects.Energy-momentum eigenstates with definite particle and antiparticle number are obtained acting with c † k and d † k operators on the vacuum |Ω⟩, i.e., the ground state of the model.In Eq. (A6) we express the anticommuting creation operators c † k and d † k in terms of ξ † x and ξ x .However, these states are completely delocalized in space (and time) while realistic particle states are characterized by some space localization and are prepared acting with wave packets of c † k and d † k , defined as for a fermion and an antifermion respectively, with momentum space probability density Combined with the orthonormality condition in Eq. (A14), Eq. ( 15) enforces the normalization of the prepared state.The functions φC,D x in Eq. ( 14) are specified in terms of ϕ k through Eq. (A6).The typical momentum space amplitude ϕ k is that of a Gaussian wave packet, namely where the phase centers the wave packet in µ x in position space and the normalization constant N ϕ is fixed by Eq. ( 15).In the continuum (σ k ≫ a) and thermodynamic limits, the |ϕ k | 2 given in Eq. ( 16) approaches the probability density of a Gaussian momentum space distribution with mean µ k and standard deviation σ k [107].
The time evolution of a state consisting of two Gaussian fermions wave packets with different µ x , µ k and σ k is depicted in Fig. 2. The propagation speed of the two wave packets matches the expected result from the lattice group velocity, ω ′ k , given by Eq. (A5c).As Fig. 2 shows, the free theory phenomenology consists only of pure kinematics; in order to observe a nontrivial dynamics an interacting theory, such as QED 2 , has to be investigated.

B. QED
Continuum QED 2 has no free asymptotic charged states [81,105,106]: due to Gauss law and the linear rise of the Coulomb potential, typical of 1 + 1 dimensions, the model exhibits a confining force and isolated charges correspond to states of infinite energy [106].The stable QED 2 particle states are thus neutral mesons, i.e., fermion-antifermion bound states.In general, a creation operator for a lattice meson of momentum for some function ηpq .Here c † k and d † k are gauge invariant fermion and antifermion creation operators.They are obtained (i) writing c † k and d † k in terms of the position space fields ξ † x and ξ x and (ii) dressing the latter with strings of unitary electric field rising or lowering operators acting on their right: to comply with Gauss law.Unless appropriate conditions are imposed on ηpq , the operator b † k in Eq. ( 17) creates excitations of definite momentum k but not of definite energy.Rather than seeking an approximate solution for QED 2 of the (notoriously difficult [102,108]) bound state problem, in this work we fix the functional form of ηpq with an ansatz, namely ηpq = η p−q = N η e −(q−p) 2 /4σ 2 ∆k e −i(q−p)µ∆x/2 , ( where In Eq. ( 19) we switched to center of mass and relative coordinates of the fermion and antifermion constituents of the meson and assume ηpq depends only on the latter.Moreover, we require: (i) that the fermion and antifermion are located in close real space positions, y and z, with average separation ⟨z − y⟩ = µ ∆x ; (ii) that (q − p) follows a Gaussian probability distribution centered in µ ∆k = 0 and with standard deviation σ ∆k .Using Eq. ( 19) in place  of the exact ηpq is equivalent to introducing some excitation in the bound state created by b † k .As a consequence, some internal dynamics is to be expected.We monitor this approximation a posteriori and set the simulation timescale shorter than the lifetime of the mesons.
Meson wave packets are prepared starting from the QED 2 ground state |Ω⟩ and acting with operators with ψ pq = ϕ p+q η q−p .As for the free theory, ψyz is implicitly defined in terms of ψ pq via Eq.(A6).The functions ψ pq and ψyz for a pair of Gaussian meson wave packets, with ϕ k given by Eq. ( 16), are plotted in Fig. 3.The states prepared with B † ϕ operators are normalized  a posteriori because, in our approximation, no exact orthonormality condition analogous to Eq. (A14) holds a priori for the states created by b † k , when g > 0. Before focusing on meson-meson scatterings, we simulate the free propagation of one meson and test the approximation in Eq. ( 19) in the parameters region explored hereafter.The time evolution the meson wave packet in Figs.3a and 3b is shown in Fig. 4 for mass m = 0.8 and couplings up to g = 0.15.Similar results are observed for masses m = 0.6 and m = 1.0.As antici- pated, the meson does not decay during the simulations and, crucially, its entanglement track remains confined in the region where the wave packet is localized.On the other hand, we do observe some internal dynamics in the wave packets, taking the form of periodic inversions of the meson polarization (sign of its electric string), whose frequency increases with the coupling g [31,109].This behavior is also captured by Fig. 5, where the voltage drop across the chain, ∆V = a x gE x, x+a , is plotted as a function of time.Looking at the relevant rows of Fig. 4 we see that the inversions are accompanied by a drop in the electric field energy h elec and a concentration of the mass energy density h mass .We infer that the inversions correspond to damped oscillations of the fermion and antifermion constituents of the meson around their center of mass.
As far as the interpretation of the results is concerned, hereafter we consider as exact the ansatz for the meson creation operator in Eq. (19).For instance, in the analysis of Section V we assume b † k creates monochromatic mesons, neglecting the consequences of the internal dynamics on the entanglement.

IV. DYNAMICS
In this Section we present some tensor network simulations of lattice QED 2 meson-meson scatterings.We explore a region of model parameters ranging from weak (g ≪ m) to intermediate (g/m ≈ 1/4) coupling.All the reported simulations follow the scheme outlined in Section II.They are carried out in the center of mass frame of reference and for parity-symmetric initial configurations.Accordingly, the initial scattering state is prepared by using a pair of uncorrelated meson wave packet creation operators form Eqs. ( 19) and ( 20) with identical σ k and σ ∆k , opposite mean momenta µ k , and symmetric µ x and µ ∆x .

A. Scattering phenomenology
An overview of the dynamics is shown in Fig. 6.Figures 6a to 6c show the scattering of two initial mesons for increasing bare mass m = 0.6, 0.8, 1.0 and fixed coupling strength g = 0.08.Increasing m clearly affects the kinematics by slowing down the propagation of the mesons before and after the collision but it also results in the generation of a larger entanglement between the scattering products.To quantify the first effect we focus on the initial stage of the evolution (t < 30) and linearly interpolate the trajectories µ x (t) of the mass energy density peaks: The resulting group velocities ω ′ are plotted in Fig. 7.
The second effect, namely the increased entanglement after the scattering, may be interpreted as an indirect consequence of the slow down of the colliding particles, as they effectively interact for a longer time.Figures 6d to 6f show again the scattering of the initial mesons for fixed bare mass m = 1.0 and coupling g = 0.02, 0.08, 0.14.We observe a drastic increase of the post-collision entanglement with the strength g of the interactions.A detailed discussion of this phenomenon is given in Section V. Furthermore, in Figs.6d and 6e the polarizations of the outgoing mesons are inverted, as indicated by the sign of the electric field.The string inversion, however, is not a consequence of the collision.As discussed in Section III, inversions are an (accidental) internal feature of our meson states and are observed also in Fig. 4, during the free propagation of a single meson.

Energy balance
The total energy ⟨H⟩ of a QED 2 state can be partitioned in mass, electric and kinetic energy fractions χ γ , γ ∈ {mass, kin, elec}.In the simulation frame of reference these are defined as We monitor χ γ (t) during various meson-meson scatterings.The time averages over the whole simulation time span, denoted χγ , are reported in the stacked area plots of Fig. 8, together with the standard deviations (time fluctuations) ϵ γ , (ϵ γ ) 2 = (χ γ (t) − χγ ) 2 , of the mass and kinetic fractions.The mass and kinetic energies are the dominant contributions in the parameter region we explored, while the electric energy fraction never exceeds 5%, decreasing mildly with the mass and growing as g 2 for large enough couplings (g ≥ 0.05).The relative weight of the mass and kinetic energies is mostly controlled by the mass parameter.On the other hand, energy transfers (namely, amount of fluctuations ϵ γ ) are boosted for larger couplings, i.e., stronger interactions.We stress that, even for the largest simulated couplings, the mass energy is approximately constant during the evolution, with time fluctuations amounting at most to 0.5% of its value.This behavior strongly hints that the simulated processes are elastic collisions.In (1 + 1)dimensions, energy-momentum conservation implies that the products of an elastic collision of two particles (of equal mass) have the same momenta of the incoming particles.In the next Section we verify this kinematical constraint by analyzing the momentum content of the initial and final states.

Momentum space analysis
Among the most important observables for a scattering process are the species and the momenta of the incoming and outgoing particles involved in the collision.In experiments, the momenta of the incoming particles are tuned by collimating the colliding particle beams.The species and momenta of the scattering products, instead, are inferred from detector data.Our simulations follow a similar procedure: the momenta of the incoming particles are set by the initial state, while those of the scattering products are identified analyzing the final state.The momenta of the excitations contained in a state are reflected by its spatial periodicities.To detect the correlations between pairs of meson excitations in different positions we evaluate the connected 2-point correlation function for the meson composite operator η Then, we compute the translation invariant momentum space connected 2-point function Before evaluating Eq. ( 23) on the initial and final scattering states of our simulations, let us discuss the case of the ground state.Due to the nonzero correlation length λ of the full QED 2 vacuum, in the thermodynamic limit [110], yz ∝ e −|y−z|/λ , G (Ω) k ∝ sinh(2a/λ) cosh(2a/λ) − cos(2ak) .
(24) Fitting the numerical G k via Eq.( 24) we find that the thermodynamic limit provides an excellent approximation of the numerical results (up to a constant shift) see Fig. 9a.Moreover, for all the simulated m and g values we find λ/(aL) ≈ O(10 −3 ) ≪ 1; signaling that we are indeed at the thermodynamic limit and justifying a posteriori the usage of the expressions in Eqs. ( 23) and ( 24).
The 2-point function G k of initial and final scattering states also presents a background of the type in Eq. ( 24) (up to a shift).On top of this background, we observe peaks detecting the momenta of the incoming and outgoing mesons.The momentum space correlators of some initial and final scattering states are plotted in Figs.9b  and 9c, with the background removed.All initial and final correlators are peaked around the mean momentum of the initial wave packet, as expected for a 2 → 2 elastic scattering in 1+1 dimensions.The distortions appearing in final state correlators are likely caused by the inexact modelling of the meson excitations.Time-independent momentum space connected meson-meson correlation functions, G k , at m = 0.6 and for positive momenta k (the k < 0 branch is symmetric).Evaluated: at g = 0.08 on the full QED2 vacuum (a) and the initial meson-meson scattering state (b); as well as on final scattering states for various couplings g (see color bar) (c).The vacuum and initial G k are almost coupling g independent (up to percent order deviations).In (a) the fitted G (Ω) k from Eq. ( 24) is also reported, while in (b) and (c) it has been subtracted.The dashed vertical line corresponds to the mean momentum of the (left) initial wave packet.

B. Towards S-matrix elements
The central quantity in scattering theory is the S operator or S-matrix [102].In this Section, after having briefly introduced the problem in the continuum, we lay out a prescription for extracting S-matrix elements from dynamical lattice simulations.
S-matrix elements are essential to verify the predictions of a theoretical model, as they relate the particle content of the initial (infinite past) and final (infinite future) states of a scattering experiment [102].Here we set coordinates in which the collision takes place at x = t = 0. Ideally, the S-matrix reads S = lim t→∞ e −2itH , where H is the Hamiltonian of the model.However, this limit involves infinitely oscillating phases and thus does not exist; nor does lim t→±∞ |Ψ(t)⟩, where |Ψ(t)⟩ = e −itH |Ψ⟩ describes the state of a system undergoing a scattering process.To overcome this problem, the asymptotic evolution has to be factored out in the definition of the S-matrix [103,104].Explicitly, where H 0 is the Hamiltonian describing the free kinematics of the stable particle states of the theory de-fined by H.The definition in Eq. ( 25) is motivated by the assumption [103,104,111,112] that the interaction decays rapidly enough so that, at asymptotic times t → ±∞, when particles are far apart, the evolution specified by H coincides with that of H 0 and the scattering solution |Ψ(t)⟩ approaches the trajectories of some freely evolving particle states |Φ ± ⟩; namely e −itH |Ψ⟩ ∼ e −itH0 |Φ ± ⟩. Parametrizing a complete set of asymptotic configurations as {|Φ α ⟩}, the S-matrix elements read S α ′ α = ⟨Φ α ′ |S|Φ α ⟩.The index α typically runs over momenta, spin projections and possibly other discrete labels [102].
Two observations are in order.As just mentioned, S-matrix elements are usually specified for momenta (and energy) eigenstates.However, these states are completely delocalized in space (and time) and thus cannot describe noninteracting particles localized in far-distant space regions.Indeed, a limit of narrow momentum space wave packets is implied in the previous construction [102,104].Moreover, the infinite time limits are a useful idealization [102,111]: in real world experiments, measurements are carried out at macroscopic times that precede and follow the collision by a time lapse T , much larger than the microscopic timescale of the collision itself but still finite.These measurements should reveal a state that approximates a free multiparticle wave packet with a degree of precision related to that of the time limit.It follows that the transition amplitudes measured experimentally are with |Φ(−T )⟩ and |Φ ′ (+T )⟩ wave packets of approximate free multiparticle energy-momentum eigenstates.
With tensor network scattering simulations, the evaluation of the transition amplitudes in Eq. ( 26) is straightforward.To this aim, the initial state |Φ(−T )⟩ is evolved until some scattering products emerge, after the collision, in the form of well separated particle wave packets.An estimate of how well the evolving state resembles a state of isolated particles is obtained, e.g., inspecting its mass energy density profile, in particular the location and width of its peaks.Once the desired separation is reached, the state |Φ(T )⟩ = e −2iT H |Φ(−T )⟩ is stored and the simulation is terminated.At this point, Eq. ( 26) provides the amplitude of the transition -due to the collision -from the initial state, |Φ(−T )⟩, to any |Φ ′ (+T )⟩ we are interested in.
Here we identify T with the time at which the mass energy density peaks of the outgoing wave packets are separated by ∆x > ∼ 100 and focus on transitions amplitudes to final states |Φ ′ (+T )⟩ of two mesons.In order to study the distribution of their momenta, we consider meson wave packets with amplitudes θ (p,L) k and θ (q,R) k peaked at momenta p and q and completely delocalized in the left and right half of the chain respectively.That is, we compute the overlap of the final state with states of the form with N pq enforcing the normalization and Θ(x) being the step function.Similarly for θ (p,L) .Exemplary transition probabilities are plotted in Fig. 10.As for Fig. 9 and as expected for kinematical reasons, the momentum distributions of the final meson are concentrated exactly on the momentum space support of the initial wave packets.As a final remark, let us stress that it is in principle possible to extract the S-matrix defined in Eq. ( 25) from lattice tensor network simulations.The starting point are the finite time amplitudes in Eq. ( 26): in order to correctly identify which S-matrix element is computed, the free H 0 evolution of the wave packets |Φ ′ (+T )⟩ and |Φ(−T )⟩ from time t = 0 to t = ±T has to be compensated in Eq. (26).In this way approximate (due to the finite T ) and smeared (due to the finite momentum spread of the wave packets) S-matrix elements are obtained.Finally, the continuum and thermodynamic limits have to be performed.The thermodynamic limit reads L, T → ∞, that is, the limit is approached not only increasing the lattice size (aL), but also, simultaneously, the duration of the simulation (T ), while keeping the reference wave packet, |Φ ′ (0)⟩ and |Φ(0)⟩, unchanged.In the T → ∞ limit, |Φ ′ (+T )⟩ and |Φ(−T )⟩ approach states of truly noninteracting particles.S-matrix elements between initial and final states of definite momentum are then extrapolated taking a limit of vanishing momentum spread of the reference wave packets.In Appendix C we show that the wave packets θ (p,L) and θ (q,R) behave as momentum projectors in the continuum and thermodynamic limit.

V. ENTANGLEMENT GENERATION
The Von Neumann entanglement entropy S observed in the scattering simulations of Section IV can be traced back to three major sources: the ground state (S grn , background), the particle -either fermion, antifermion or meson -wave packets (S wps , intraparticle entanglement) and their interactions (S int , interparticle entanglement).As a first approximation, we consider these contributions as additive.
The models in Section I are translationally invariant in the thermodynamic limit or for periodic boundaries.Hence, the ground state |Ω⟩ is responsible for a uniform (up to boundary effects) background entanglement S grn .In what follows we estimate a priori the wave packet contribution to infer the interaction entanglement generated by the dynamics.

A. Wave packet entanglement
Acting on |Ω⟩ with the wave packet creation operators from Eqs. ( 14) and ( 20), additional entanglement, S wps , appears on top of S grn in the space region where the particles are localized.This entanglement contribution can be characterized for freely propagating wave packets in the m → ∞ limit, assuming that the contributions of different wave packets are additive, as expected for uncorrelated wave packets.We checked this assumption for a few overlapping wave packets of either fermions or antifermions, and found additivity violations (due to the exclusion principle) of a few percent order.

Fermion and antifermion entanglement
We consider a single fermion wave packet |Ψ⟩ = C † ϕ |Ω⟩, prepared via Eq.(14a); analogous results hold for an antifermion.In the infinite mass limit, the Hamiltonian in Eq. ( 2) is dominated by the mass term and the ground state approaches the bare vacuum -i. with where S ∈ {L, R}.Furthermore, exploiting the expression of φC x in Eq. (A17), Φx is identified with the cumulative distribution function ( Equation ( 29) represents the coherent superposition of the overall fermion excitations in each side of the chain.In a basis obtained extending {|Ω⟩ S , |Ψ⟩ S } to a complete orthonormal system, the reduced density matrix of the S subsystem reads ρ = diag( Φx , 1 − Φx , 0, 0, . ..), yielding the entanglement entropy profile Note how, at the median x of the wave packet distribution Φx = 1/2 and |Ψ⟩ in Eq. ( 29) becomes a Bell state, thus S(x) = 1.Specializing to a Gaussian (fermion or antifermion) wave packet, it follows from Eq. (A19) that in the thermodynamic and continuum limits Φx becomes the CDF of a Gaussian distribution with mean µ x and standard deviation 1/2σ k : An analytic expression for S(x), denoted Ŝ(x; µ x , σ x ), is obtained substituting Φx from Eq. ( 32) in Eq. ( 31).

Meson entanglement
The previous results are straightforwardly extended to meson wave packets exploiting the formulation of QED 2 with the link degrees of freedom integrated out from Eq. ( 10).In the infinite mass limit the ground state is again the Néel state.Suppose that the meson wave packet creation operator in Eq. ( 20) factorizes as a product of fermion and antifermion creation operators, for some amplitudes ϕ ± .Then, by the additivity assumption, the entanglement profile of the meson is the sum of a pair of Eq. ( 31) contributions, S = S ϕ + + S ϕ − .Figure 3 shows that this is a good approximation for the Gaussian mesons of our simulations, which are characterized by σ k ≈ σ ∆k .Indeed, in this case, Eq. ( 33) holds up to O(σ 2 k − σ 2 ∆k ) terms, with ϕ ± k Gaussian amplitudes whose parameters are related to those of the meson wave packet by the substitutions Moreover, the µ ∆x /2 phase shift can be neglected if the fermion and antifermion are almost overlapped, i.e., if µ ∆x ≪ 1/σ k , as is the case for the mesons in Fig. 3.We conclude that the entanglement due to a Gaussian meson wave packets is approximated by 2 Ŝ(x; µ x , σ x ), where Ŝ is the result derived for fermions and antifermions and σ x = 1/ σ 2 k + σ 2 ∆k .

Entanglement propagation
The evolution of a particle wave packet k ϕ k |k⟩, i.e., a wave packet of eigenstates |k⟩ of momentum k and energy ω k , is given by ϕ k (t) = e −iω k t ϕ k .This is true for fermions and antifermions in the free theory, as well as for (exact) mesons in QED 2 , even though the dispersion relation ω k of mesons is not known analytically if g > 0.
Let us focus again on a Gaussian wave packet in the continuum with ϕ k (t = t 0 ) given by Eq. ( 16).A stationary-phase approximation shows that the absolute square of the inverse Fourier transform of ϕ k (t) is a Gaussian PDF also for t ̸ = t 0 , but translated and widened according to Combined with Eq. (A19) and the results of the previous paragraphs, Eq. ( 35) provides the entanglement of freely propagating fermions, antifermions and (factorizable) mesons.We model the wave packet entanglement due to a pair of parity-symmetric scattered mesons introducing where w = (N = 2, t 0 , µ x , σ x , ω ′ µ k , ω ′′ µ k ) collects all the parameters entering Eqs. ( 35) and (36).Having been derived for freely propagating particles, Eq. ( 36) applies only for g = 0 or when the mesons are distant enough so that their interaction can be neglected.In Fig. 11 we compare the numerical wave packet entanglement S wps (t, x) from a g = 0 simulation with the prediction from S th wps (t, x; w), for both the expected parameters w = w th , and the values w fit that best interpolate the numerical data (setting t 0 = 0).The function S th wps models accurately the numerical entanglement; comparing w th and w fit , the prominent distortion is a slight squashing of the expected entanglement profiles due to the lack of exact additivity between the various entanglement contributions.Comparison of the numerical results with the S th wps (t, x; w) prediction from Eq. (36).The contour plot reproduces the numerical data Swps.The bottom and top panels show the initial and final time spatial entanglement profiles; the left panel shows the midchain time profile.In each side panel, along with Swps, we plot the relevant section of S th wps (t, x; w) for (blue lines) w = w th and (red lines) w = w fit .

B. Interaction entanglement
The ground state and the wave packets completely explain the entanglement observed in our simulations of the free kinematics.In particular, if the system is cut outside the support of the outgoing wave packets, the entanglement entropy of the bipartition comes from the ground state only.When the interaction is switched on (g > 0) this is no longer true, because additional entanglement is generated by the dynamics.In the elastic scattering regime explored with our simulations, the final time entanglement entropy profiles are characterized by a uniform plateau in the region enclosed between the two outgoing wave packets; we thus interpret the dynamically generated entanglement as entanglement between the scattering products.The entanglement plateau is clearly visible in Fig. 12a and especially Fig. 12b, where we subtract the contribution from the g = 0 simulation.For g > 0, entanglement is present also in the middle of the chain, where the mass energy density, shown in Fig. 12c, vanishes.
The time profile of the midchain entanglement in Fig. 12d shows that the correlation between the outgoing mesons is produced by the interactions, during the collision.Indeed, in the initial stage of the evolution, when  the incoming mesons approach one another, no entanglement is present between them.However, once their wave packets overlap, the associated entanglement is detected at the middle of the chain.This process gets delayed as the coupling is increased (see Fig. 13a), likely due to a combination of factors, such as a deflection of the meson trajectories and the takeover of the interaction-generated entanglement.In Fig. 13b we plot the peak midchain entanglement S(t ⋆ , L/2) as a function of u = m α k β g, where k = 0.81 is the absolute value of the mean momentum of the initial mesons and the exponents α, β are reported in Fig. 14.S(t ⋆ , L/2) is decreasing at small u, but above the threshold value u ⋆ ≈ 0.08 it rapidly increases again.
Nevertheless, it is only after the collision that the fundamental distinction between the free and interacting cases -namely, the entanglement of the scattering products -emerges.Figure 14 shows the the coupling dependence of the entanglement between the outgoing mesons generated by the interaction, S int , for different values of the bare fermion mass m and of the mean momenta ±k of the initial mesons.To quantify S int we assume that all the entanglement due to the interaction is produced during the collision and remains approximately constant afterwards.Thus, in the final stage of the evolution, the midchain entanglement can be decomposed as a constant contribution S grn + S int , due to the ground state and the interaction, plus a decaying component S wps (t), attributed to the tails of the outgoing meson wave packets.We thus evolve the system until the mesons are spatially separated by ∆x > ∼ 100 and fit the midchain entanglement values sampled in the final ∆t = 40 time interval using the expression provided by Eq. ( 36) for x = L/2, plus a constant background.Since we always remove the ground state contribution, the value of the background is precisely the extrapolated S int .The fits are reproduced in Fig. 12e.According to Fig. 12c, for the highest g values the outgoing wave packets are not Gaussian, as assumed in Eq. ( 36), thus the fits are only partially justified.Yet, especially for stronger couplings, towards the end of the evolution the wave packets give a small relative contribution to the midchain entanglement and its final value is already an accurate estimate of S int .
We find that the entanglement produced by the interaction -i.e. the amount of of quantum correlations between the outgoing mesons -increases with the coupling, as expected.Moreover, we find that the collapse of the interaction entanglement curves on a unique curve F (u) shown in Fig. 14c is described by The optimal exponents α, β, γ, δ in the least-squares sense are reported in Fig. 14c.They are obtained minimizing the residual sum of squares of the rescaled data points from a common interpolating polynomial of degree 10.Equation (37) allows us to express our findings in terms of u only: after an initial gentle growth for small u, at u ⋆ we observe an abrupt increase in the slope of S int , which then stays constant up to u ≈ 0.12, at which point the slope diminishes again.

VI. CONCLUSIONS AND OUTLOOK
We have shown that Tensor Network (TN) methods are viable numerical techniques for the simulation of the dynamics of (1 + 1)-dimensional lattice Quantum Electrodynamics (QED 2 ).We identified a Matrix Product State (MPS) representation of the stable particle states of the theory, namely meson bound states.We performed real-time TN simulations of the elastic collision of a pair of mesons for moderately week couplings (g/m < ∼ 1/4) and large lattice sizes (aL ≫ λ, λ being the correlation length of the model).In the simulated processes the entanglement growth is perfectly sustainable, attesting that classical TN algorithms are capable of attacking this problem efficiently and accurately.
In Section III and Appendix B we elaborated a protocol for the preparation of an MPS consisting of some approximate meson wave packets.The protocol is based on the exact solution of the theory of free staggered fermions, given in Appendix A. With this regard, let us stress that rather than the breakdown of the numerical tools employed, it is only our approximation of the meson states that prevented us from studying the QED 2 dynamics at stronger couplings.Recently developed numerical techniques, such as tangent space MPS methods [60,[113][114][115][116][117] and in particular [118], should furnish a valid replacement for our (perturbative) initial state preparation protocol.We expect the identification of the exact particle excitations of the model to provide access to stronger coupling regimes and thus, likely, more diverse dynamical processes.Given the resources required by our simulations, we expect the dynamics of more complex gauge theories to be accessible via tensor network simulations as well.Examples that are worth looking into are theories involving multiple matter field flavors and non-Abelian gauge groups [27,28,38,72].
In Section IV we put forward two strategies for analyzing the momentum content of the system and applied them to identify the momenta of the collision products in our scattering simulations.The first strategy involves the computation of a specific 2-point connected correlation function.The second relies on the evaluation of a set of finite-time scattering amplitudes, obtained projecting the system state on a family of carefully constructed wave packets whose properties are discussed in Appendix C. Building upon these scattering amplitudes, we gave a prescription for extracting continuum S-matrix elements from tensor network lattice simulations.We believe that its implementation is within the reach of the currently available numerical tools, with reasonable computational resources.For an alternative approach, based on the Lippmann-Schwinger formalism, see [119].
In Section V we analyzed the evolution of the entanglement content of the system during our simulations of elastic meson-meson collisions.We found that the entanglement entropy observed in these scattering processes can be decomposed in three, approximately additive, distinct contributions, namely the following: 1. Vacuum entanglement: a uniform layer of entanglement produced by the correlations in the ground state of the theory.
2. Intraparticle entanglement: entanglement bumps in the regions where particle wave packets are localized.
3. Interparticle entanglement: a dynamically generated entanglement string correlating the products of the collision.
After having derived approximate analytical expressions that model a particle internal entanglement, we focused on the inter particle entanglement S int .We studied its growth with the interaction strength g for different values of the bare mass parameter m and mean momentum k of the initial meson wave packets.In the explored parameters region, we found a phenomenological relation describing the interplay of the mass, the coupling, the meson momenta, and the interaction entanglement S int .Precisely, by Eq. ( 37), where F is the function plotted in Fig. 14c.At fixed m and k, an initial slow growth of S int for small values of the coupling g is followed by a rapid and steady rise, for u > u ⋆ ≈ 0.08.The above relation offers a testbed for the investigation of the relation between entanglement and scattering amplitudes, motivating additional real-time dynamical investigations of scattering events.
the Hamiltonian becomes In the last step we recognized that the quadratic forms of the k and k + π/a modes in Eqs.(A5a) and (A5e) are degenerate and restricted the sum to Λ * ′ .Explicit expressions for c † k and d † k in terms of the position space fields, ξ † x and ξ x , are obtained combining the previous results.We report them below for reference: On Λ * ′ only the first Kronecker delta in Eq. (A9) contributes and Eq.(A8) follows from the unitarity of V k .Recalling the d k , d † k anticommutator in Eq. (A8), Eq. (A5e) can be rewritten as a positive semidefinite operator minus a constant.The constant term, namely ′ k ω k , incorporates both ultraviolet and infrared divergencies when the respective regulators are released, making the continuum limit ill defined.Dropping it by a normal ordering prescription, the free staggered fermions energy-momentum operator P µ = (H, P ) reads Equation (A10) shows that c † k and d † k create excitations of energy-momentum k µ : It also identifies ω k in Eq. (A5c) as the staggered fermion dispersion relation.Its derivative ω ′ k gives the group velocity of fermion wave packets on the lattice.

Fock space
The ground state of the theory, |Ω⟩, is the vacuum P µ |Ω⟩ = 0, i.e., the state with no excitations to destroy: The Fock space of all the (normalized and antisymmetrized) single and multiparticle states of the theory is generated by acting on |Ω⟩ with products of creation operators c † k and d † k : |q N . . .q 1 ; p M . . .
A13) with p i , q j ∈ Λ * ′ .For instance, by Eqs.(A8) and (A12), In terms of c k and d k , the global U(1) charge Q reads (A15) where in the last step we have normal ordered, discarding the half-filling constant L/2.Equation (A15) shows that d-type excitations (antifermions) are the antiparticles of c-type excitations (fermions).

Infinite mass limit
The theory of free staggered fermions simplifies greatly when m ≫ 1/a.For instance, recalling the derivation in Eq. (A5), in the m → ∞ limit Eq. (A6) reduces to We conclude that fermions (antifermions) excitations are supported on the even (E) and odd (O) sublattices.Moreover, inserting Eq. (A16) in the expressions of C † (A18) has been introduced for formal convergence.This θ q (k) should be interpreted in the sense of distributions, it is unnormalizable and the variance of |θ q (k)| 2 is undefined.Nonetheless, as we now show, ⟨θ q | projects (sufficiently well behaved) wave packets peaked at x > 0 in position space, on their component, i.e., wave packet amplitude, of momentum q.
Let |ϕ y ⟩ = dk e −iky ϕ(k) |k⟩ be a wave packet peaked at x = y such that, as a complex function, ϕ(k) has no singularities and |ϕ(k)| → 0 for |k| → ∞ (these criteria are satisfied, e.g., by a complex Gaussian e −|k| 2 ).We want to compute and the residue theorem yields ⟨θ q |ϕ y ⟩ = Θ(y)e −iqy ϕ(q) , (C5) which is exactly the anticipated claim.The analogous result for wave packets delocalized in the x < 0 region follows by parity symmetry.
FIG.1.The ground state of the Schwinger model is computed and expressed as an MPS (blue bottom squares).Then a set of MPOs (red squares) is applied to prepare a gauge invariant initial state made up of two spatially separated meson wave packets with opposite electric field orientations.Indexes i and j indicate lattice sites on which matter and antimatter field components are defined, while the electric field is defined on the links between sites.The two mesons evolve under the Hamiltonian of Eq. (4) with TEBD (pictorially represented by yellow boxes) and their scattering is studied.Note the persistent entanglement between them after the scattering process.Numerical parameters: a = 1, m = 1.0, g = 0.14.

FIG. 2 .
FIG. 2. Charge density during the free propagation of two different Gaussian fermion wave packets.The black lines in overlay are the trajectories of the wave packet peak predicted by the group velocity ω ′ k .

FIG. 4 .
FIG. 4. Test of the stability of QED2 mesons prepared using Eqs.(19) and (20).The electric field ⟨E x, x+a ⟩, its energy density h elec x , the mass energy density h mass x and the entanglement entropy Sx (rows) are measured during the propagation of a single meson, for mass m = 0.8 and increasing values of the coupling g (columns).

15 FIG. 5 .
FIG.5.Damped oscillation of the voltage drop across the chain ∆V during the propagation of the a single meson (simulations in Fig.4) for different values of the coupling g.Notice the increase of the oscillation frequency with the coupling g.

FIG. 6 .
FIG. 6. Simulations of meson-meson collisions in QED2.Measurements of the electric field ⟨E x, x+a ⟩ t , its energy density h elec x (t), the mass energy density h mass x (t) and the entanglement entropy S(t, x) are reported (rows).Various values of the mass m (left block) and coupling g (right block) parameters are considered.The initial meson wave packet distributions are those of Fig. 3.

FIG. 7 .
FIG.7.Mass m and coupling g dependence of the meson group velocity ω ′ , as defined in Eq. (21) and measured during the first (approximately free) stage of the propagation.Error bars represent three standard deviations.The uncertainty increases with the coupling g because stronger interactions anticipate the deflection of the meson trajectories.

FIG. 8 .
FIG.8.Total energy partition in the mass, kinetic and electric energy contributions from Eq. (22): time averages (hatched regions) and fluctuations (filled regions) during various meson-meson scatterings.Masses m = 0.6 (a), 0.8 (b) and 1.0 (c) and couplings ranging from g = 0 to g = 0.14 are considered.The thickness of the filled regions reproduce a ±2ϵ γ confidence belt around the mean value χ γ of the kinetic and mass energy fractions.
FIG. 9.Time-independent momentum space connected meson-meson correlation functions, G k , at m = 0.6 and for positive momenta k (the k < 0 branch is symmetric).Evaluated: at g = 0.08 on the full QED2 vacuum (a) and the initial meson-meson scattering state (b); as well as on final scattering states for various couplings g (see color bar) (c).The vacuum and initial G k are almost coupling g independent (up to percent order deviations).In (a) the fitted G

FIG. 10 .
FIG.10.Transition probabilities |Apq| 2 to the states defined in Eq. (27), for three meson-meson scattering simulations with bare mass m = 0.8 and coupling g = 0.02 (a), 0.08 (b), and 0.14 (c).Namely, probability of the two mesons in Fig.3to evolve, after a collision, into pairs of meson wave packets peaked around momenta p and q and delocalized in the left and right side of the chain, respectively.The dashed lines correspond to the initial meson momenta µ k ≈ ±π/4.The resolution of the above images is related to the momentum space standard deviation ε ≈ 0.14 of these wave packets and thus, indirectly, to the lattice size.The transition probabilities in Figure are computed projecting on a family of final states with p and q values spaced by ε/2.
e., the Néel product state |Ω⟩ = |0101 . ..⟩.Consider a bipartition of the chain Λ in the subsystems L = {y < x} and R = {y ≥ x} with x ∈ E (even sublattice), so that |Ω⟩ = |Ω⟩ L ⊗ |Ω⟩ R and |Ψ⟩ decomposes as FIG.11.Wave packet entanglement entropy during the free propagation of the mesons in Fig.3for m = 1 and g = 0.Comparison of the numerical results with the S th wps (t, x; w) prediction from Eq.(36).The contour plot reproduces the numerical data Swps.The bottom and top panels show the initial and final time spatial entanglement profiles; the left panel shows the midchain time profile.In each side panel, along with Swps, we plot the relevant section of S th wps (t, x; w) for (blue lines) w = w th and (red lines) w = w fit .

FIG. 12 .
FIG. 12. Entanglement entropy S and mass energy density h mass space and time sections (with ground contribution removed) for elastic meson-meson scatterings with m = 0.6 and various couplings g (see the color bar in the last row).Panels (a)-(c) plots are obtained at a fixed (g dependent) time t = Tg, corresponding to outgoing wave packet peaks separated by ∆x ≈ 100.Panels (a), (b) reproduce the final time Tg entanglement profile, with the free contribution (g = 0) subtracted in (b); while (c) shows the final time mass energy density.Panels (d), (e) show the midchain entanglement as a function of time; (e) in an enlargement of its decaying tail, approaching the extrapolated asymptotic value, S(∞, L/2) = Sint + Swps(∞) ≈ Sint.

FIG. 13 .FIG. 14 .
FIG.13.Time t ⋆ (a) and entanglement S(t ⋆ , L/2) (b) corresponding to the peak of the midchain entanglement profiles, for different masses m and couplings g.The ordinate represents the relative deviation from the values obtained for the free (g = 0) simulation, namely t ⋆ 0 and S0(t ⋆ 0 , L/2).

FIG. 15 .
FIG. 15.Wx matrices of the meson wave packet creation MPO.Empty entries represent null operators.

⟨θ q |ϕ y ⟩ = − 1
|k|, ϕ(k) k − q + iϵ ≤ |ϕ(k)| |k| − |q − iϵ| < |ϕ(k)| ; (C3)therefore, we can invoke Jordan's lemma closing the contour of integration in the lower half (upper half) of the complex plane for y > 0 (y < 0).The integrand has a single pole at k = q − iϵ with limϵ→0 + Res k=q−iϵ e −iky ϕ(k) k − q + iϵ = e −iqy ϕ(q) (C4) It follows that |N x ⟩ provides a basis for H phys , which is thus unitarily equivalent to the Hilbert space of free staggered fermions H free .Consequently, a free theory operator O also defines an operator on H phys that can be extended to a (gauge invariant) dressed operator Ō on the whole H . Conversely, as far as gauge independent properties of the model are concerned, gauge invariant QED 2 operators can be expressed in the |N x ⟩ basis.For instance, in this basis the Hamiltonian becomes