Real-time Dynamics in U(1) Lattice Gauge Theories with Tensor Networks

Tensor network algorithms provide a suitable route for tackling real-time dependent problems in lattice gauge theories, enabling the investigation of out-of-equilibrium dynamics. We analyze a U(1) lattice gauge theory in (1+1) dimensions in the presence of dynamical matter for different mass and electric field couplings, a theory akin to quantum-electrodynamics in one-dimension, which displays string-breaking: the confining string between charges can spontaneously break during quench experiments, giving rise to charge-anticharge pairs according to the Schwinger mechanism. We study the real-time spreading of excitations in the system by means of electric field and particle fluctuations: we determine a dynamical state diagram for string breaking and quantitatively evaluate the time-scales for mass production. We also show that the time evolution of the quantum correlations can be detected via bipartite von Neumann entropies, thus demonstrating that the Schwinger mechanism is tightly linked to entanglement spreading. To present the variety of possible applications of this simulation platform, we show how one could follow the real-time scattering processes between mesons and the creation of entanglement during scattering processes. Finally, we test the quality of quantum simulations of these dynamics, quantifying the role of possible imperfections in cold atoms, trapped ions, and superconducting circuit systems. Our results demonstrate how entanglement properties can be used to deepen our understanding of basic phenomena in the real-time dynamics of gauge theories such as string breaking and collisions.


I. INTRODUCTION
The mechanism of quark confinement stands as a key concept in our understanding of the fundamental interactions in high energy physics [1][2][3][4]. As a quark and an anti-quark are pulled apart, the energy stored in the gluon string connecting them grows linearly with distance, due to the confining nature of strong nuclear forces described by quantum-chromodynamics (QCD). In gauge theories hosting dynamical charges, there exists a critical length scale at which the confining string breaks, creating particle-antiparticle pairs which reduce the energy density in the string and give rise to the hadrons at the string edges [5].
The static properties of string breaking have been widely explored using a variety of lattice methods, wherein the effective string potential separating static charges can be extracted by the Polyakov or Wilson loops [6][7][8]. However, the real-time dynamics of gauge theories are usually biased by a severe sign problem, and as such cannot be accessed using lattice Montecarlo simulations [9][10][11]. In this paper, we apply tensor network (TN) methods in order to study the real-time dynamics of a lattice gauge theory (LGT) with dynamical charges and quantum gauge degrees of freedom in one dimensional systems. In particular, we investigate the realtime string-breaking dynamics in Abelian U (1) LGTs in (1+1)d, which share with QCD the basic feature of con-finement.
In recent years, efficient numerical methods based on TNs have found widespread application to the realtime dynamics of strongly correlated low-dimensional systems [12]. They are nowadays routinely used to tackle a variety of condensed matter and atomic physics problems, such as the evaluation of spectral functions of low-dimensional magnets and the quench or controlled dynamics of ultra-cold quantum gases in optical lattices [13][14][15][16][17][18][19][20][21][22][23][24][25][26]. While TN methods have been extensively applied to spin and Hubbard-type models, only recently it has been shown how TNs can provide an ideal platform for the investigation of gauge theories, for example in the study of the static properties of the Schwinger model, the low-energy mass excitation spectrum, and the dynamics of deconfinement in 2D Z 2 LGTs [27][28][29][30][31][32][33][34][35][36][37]. In particular, the quantum link model (QLM) formulation of LGT [38][39][40] -gauge theories whose link Hilbert space is finite-dimensional -has been used to develop efficient general-purpose TN algorithms to describe static and real-time dynamical properties of Abelian and non-Abelian LGTs including generic forms of matter fields, and to present different possible quantum simulator implementations on different platforms: trapped ions, cold atoms in optical lattices and circuit quantum electrodynamics (QED) [41][42][43][44][45][46][47][48][49][50][51][52][53][54][55] (see also [56,57] for recent reviews and references therein). Notice that the Gauss law depends on the lattice site due to the staggered fermions. Middle panel: cartoon states for the different stages of the string breaking dynamics (see text). Here the leftmost site is an odd one. Right panel: sample simulation for the electric field dynamics when quenching an initial string state (B region) connecting two charges, and surrounded by the vacuum (A regions) for m = 0 = g. Primary string breaking takes place in four stages (C-F), until an anti-string is created in place of the original string. The latter also decays during the secondary string breaking. The shaded areas represent the wave-fronts estimated from entanglement entropies (see Sec. IV), which are directly related to the electric field evolution.
ansatze for the many-body wave function of the quantum system of interest: the tensor structure is chosen to best accommodate some general system properties, e.g., dimensionality, boundary conditions and symmetries, while a controlled approximation is introduced in such a way that one can interpolate between a mean field and an exact representation of the system. Being a wave function based method, one has direct access to all relevant information of the system itself, including quantum correlations, i.e., entanglement. In one-dimensional systems, an efficient tensor structure is given by the Matrix Product State (MPS) ansatz [12,14], defined as, where the tensor A contains the variational parameters needed to describe the system wave-function, α i = 1, . . . , d characterize the local Hilbert space, and β i = 1, . . . , m account for quantum correlations or entanglement (Schmidt rank) between different bipartitions of the lattice. Indeed, setting m = 1 results in a mean field description, while any m > 1 allows for the description of correlated many-body states. Given the tensor structure, the tensor dimensions and coefficients are then optimized to efficiently and accurately describe the system properties by means of algorithms that scale polynomially in the system size and m. Usually, these algorithms exploit the system Hamiltonian tensor structure, naturally arising from the few-body and local nature of the interactions, to efficiently describe the system ground state or low-lying eigenstates, or to follow the real-or imaginary-time evolution of the system itself. Indeed, in the TN approach, real-and imaginary-time evolution have no fundamental differences at the computational level as there is no sign problem, and limitations arise -only in some scenarios depending on the specific dynamics of interest as witnessed by the fast increasing literature appearing based on this approach [12] -from the amount of quantum correlations present in the system wave function.
Here, we show how TN algorithms allow for the study of the real-time dynamics of LGTs, focusing on the string breaking in a paradigmatic confining theory -the Schwinger model [59][60][61] in a quantum link formulation. We characterize the real-time dynamics of the primary and secondary string breaking and show that string breaking is intimately related to entanglement production in the system. A qualitative picture for string breaking in our models, together with a typical result for our time-dependent simulations on a system of L = 100 lattice sites, is illustrated in Fig. 1. Even more importantly, our simulations allow us to track the entanglement evolution during string breaking: as we will show below, string breaking and the so-called Schwinger mechanism are intimately connected to entanglement propagation, which we address by evaluating the so called von Neumann entanglement entropy. Finally, we show that TN methods can be used to study scattering processes between bound states of LGTs: we develop a scheme to engineer meson collisions [62], and we show how, very surprisingly, the scattering not only reflects into an enhanced rate of particle-antiparticle creation, but it does affect drastically the entanglement properties of the system, which stays significantly correlated well beyond the scattering time-window.
The paper is structured as follows: in Sec. II we present the system Hamiltonian and recall the TN algorithm we are using throughout this work. In Sec. III we present the results on string breaking and mass produc-tion dynamics, including a discussion on how this phenomenon can be observed using a quantum simulation platform. In Sec. IV we show how entanglement follows the string breaking dynamics, providing a quantitative picture which underlines how entanglement entropies are directly tied to string-breaking. Finally, we present our result on scattering in Sec. V, and draw a summary of our results in Sec. VI.

II. MODEL AND METHODS
A. Model Hamiltonian: QED in (1+1)d QED in (1+1)d, also known as the Schwinger model [59], represents an ideal testing-ground for the benchmarking and development of new computational methods. Despite its relative simplicity, this model captures fundamental aspects of gauge theories such as, e.g., the presence of a chiral symmetry undergoing spontaneous symmetry breaking [59][60][61][63][64][65][66][67][68][69][70]. Even more importantly, this theory, like QCD, displays confinement: differently from (3+1)d QED, in (1+1)d electrons and positrons are confined, and interact via a long-range potential which increases linearly with distance. Due to the large energy cost associated with the electric flux between charges at large inter-charge distances, the electric flux string is unstable to particle-antiparticle creation, as in QCD, and string breaking takes place [58]. While this phenomenon, directly connected to the Schwinger mechanism of mass production out of a vacuum, has long been debated, and notable insights have been provided using a variety of approximate methods, a full quantum mechanical understanding of the complex real-time dynamics taking place during string breaking is lacking due to the computationally complexity of the many-body problem [58,[71][72][73].
In the Hamiltonian formulation, its dynamics are defined by the following form: where ψ † x , ψ x are fermionic creation/annihilation operators describing Kogut-Susskind (staggered) fermions (see Fig. 1), U x,x+1 are the gauge fields residing on the (x, x + 1) link, and we denote the strength of fermionhopping (the kinetic energy of electrons and positrons) with t, the staggered mass of the fermions with m, and the electric coupling strength with g, where E x,x+1 is the electric-field operator. The gauge generator is given bỹ and, as such, act on a finite-dimensional link Hilbert space [45]. In particular, the electric field operator allows three possible states for the electric flux, constraining the physical states per site as described in Fig. 1. A detailed discussion of the quantum link formulation can be found in Ref. [38][39][40], while in Ref. [31] it was shown how such quantum link formulation reproduces the phase diagram and quantum criticality of the continuum theory.

B. String breaking and classical cartoon states
String breaking is the process of cutting and shortening the electric flux string that connects a particleantiparticle pair by creating a new charge-anticharge pair [45]. Within our framework, a string consists of two charges creating non-zero electric flux between them. The charges are represented by appropriate boundary conditions (static charges) or by excitations of the mass field at the site of the fermion (dynamical charges). This is realized by an effective jump of a fermion from the site of one charge to the site of the second charge satisfying the Gauss law. The string of electric flux then follows from Gauss' law. The charges force the links into a non-zero flux state, according to the configuration of the charges in either one direction or the other. Before embarking in a full quantum mechanical investigation of string breaking, we now discuss the classical (t = 0) static picture which provides a simple yet informative illustration of the different stages of the string breaking mechanism. A set of cartoons of the classical states is provided in Fig. 1: Vacuum.-In the vacuum (A), neither mass nor electric field excitations are present. The vacuum energy is thus E 0 = − L 2 m. String.-In the string state (B), two mass excitations are present at the boundaries, and all electric fields connecting the two are also in the | + 1 state. The resulting string energy then takes the form Pairs.-In the pairs state (C) all the masses are excited forming charge-anticharge pairs with an energy E pairs = g 2 L 4 + mL.
Mesons.-In a confined phase, particle-antiparticle pair production can favor the establishment of a vacuum state between two static charges, with a pair of mesons at the boundary of the string (see (D)). The resulting energy is: At the static level, string breaking takes place at a critical distance L c , above which the meson state is energetically favored over the string state (E string (L c ) = E mesons ): Antipairs and Antistring.-The antipair-state (E) and the antistring (F) denote the pair-state and the string with the electric flux having the opposite sign.
At a dynamical level, string breaking takes place as a consequence of the Schwinger mechanism: in (1+1)d, the vacuum between two charges of opposite sign is unstable against particle-antiparticle creation [58,[71][72][73], eventually leading to the electric field in the system flipping sign and to mass production. In real-time, this process develops following intermediate consecutive steps, and is schematically illustrated in Fig. 1: at very short timescales, particle-antiparticle creation takes place in the middle of the string, creating the pair state depicted in (C). Subsequently, the electric field in the center of the string relaxes to 0, and the external charges get screened, effectively forming mesons (D). At this point the process reverses, first establishing a state with anti-pairs (E), which finally decays into a string of oppositely signed electric field with respect to the initial state, an antistring (F).
String-breaking is a direct consequence of confinement: while in a deconfined phase such as QED in (3+1)d, it is possible to separate opposite charges at large distances due to Coulomb's law, in a confined phase the corresponding electric field string breaks due to the effective potential between charges increasing as a function of distance. However, the exact real-time dynamics of string breaking are inaccessible to classical simulations based on Montecarlo sampling due to a severe sign-problem. While classical-statistical approaches can provide remarkable insights in some parameter regimes [58,[71][72][73] such as small masses, unbiased numerical simulations for arbitrary parameter regimes have been lacking. In the following, we present a systematic study of the string breaking dynamics using MPS techniques.

C. Tensor networks for lattice gauge theories
Tensor network algorithms are one of the paradigms for simulating quantum many-body systems in lowdimensions, both in and out of equilibrium, via a representation of the quantum state with a variational ansatz for wave-functions and/or density matrices [12].
For one-dimensional pure states, on which we focus here, the starting point is to consider a class of states of the form given in Eq. (1) and depicted in Fig. 2 with some fixed auxiliary dimension m and physical dimension d: for example, for spin one-half systems one has d = 2, while the dimension m depends on the states studied and on the desired accuracy. For ground states of onedimensional gapped Hamiltonians, m is in general independent of the system size N , while for critical systems, due to area-law logarithmic violations, the auxiliary bond dimension scales as m ∝ c log N , where c is the central charge of the system (see Ref. 75 for a review). The case of out-of-equilibrium dynamics as considered here is more challenging, and no general picture is known, thus the bond dimension m has to be adapted to each specific case and convergence has to be checked comparing the results at increasing bond dimension [76,77].
Global symmetries of the system of interest, such as particle number or global magnetization and even more complex non-Abelian global symmetries, can be embedded in the wave function ansatz given in Eq. (1) in an elegant way by promoting each tensor A βj ,βj+1 αj to a symmetry-sector-preserving tensor to a tensor that preserves every symmetry sector; that is, every index of the tensor is dressed with the corresponding symmetry charge number [78][79][80][81]. This symmetric formulation of tensor networks allows one to describe wave functions exactly and more efficiently with the desired quantum numbers, addressing each symmetry sector separately.
Recently, it has also been shown that local symmetries -namely gauge symmetries -can be embedded in such description by generalizing the wave function ansatz to a gauge invariant one [31,32,[34][35][36]. We briefly recall such a construction in the rest of this Section, and tailor it to the model we will study in the rest of the paper.
The Hilbert space of a gauge invariant system is given by the direct sum of every sector with different values of the "Gauss law",G x |Ψ QLM = g x |Ψ QLM and H = ⊕ gx H gx . Due to the gauge symmetry, there is no physical (e.g., gauge invariant) operator, including the Hamiltonian, that connects any two different sectors. Hence, starting with a quantum state partially defined by the values of the Gauss law, or a set of local (gauge) constants of motion, the time-evolved wave-function will remain in this sector under the action of the Hamiltonian. In the QLM formulation of LGTs, a gauge invariant tensor description is immediately obtained if we use a "rishon" [40] or Schwinger representation of the gauge degrees of freedom (independent of the nature of the local continuous symmetry -Abelian or non-Abelian -and the dimensionality of the lattice [31]).
In our case, the U (1) QLM in (1+1)-d with a spin-1 per link, the spin operator or electric field E is represented by a pair of Schwinger bosons E x,x+1 = 1 2 (n R,x+1 − n L,x ) with a total occupation of two bosons per link, i.e. n R,x+1 + n L,x = 2. The first step to build an efficient tensor network representation of such states is to identify a local Hilbert space spanned by the states |α defined on the tensor product of the fermionic matter field on a lattice site and of the rishon states on its left and right, that is |α = |n R,x , n Ψ,x , n L,x . This allows for the projection of the state into the gauge invariant subspace, restricting the local bases only to the "physical" states: in our model, this process results in the five gauge invariant states (uniquely identified in terms of the Schwinger rishons for even and odd lattice sites respectively and depicted in Fig.1 left) given by [45]: These are the only states allowed locally by gauge invariance (notice that the local Hilbert space is different on even and odd sites as a consequence of the staggered nature of the fermions). A gauge invariant many-body state clearly exists in the tensor product of such local basis states. However, not all tensor product combinations are compatible with the spin representation S: with our choices, only the states with n R,x+1 + n L,x = 2 are allowed. This additional constraint can be satisfied by modifying the tensor structure of the many-body state ansatz, that is by introducing a projector which acts on nearest-neighbor lattice sites and enforces the correct number of rishons on the link [31,34].
Hence, the gauge structure of systems with a local continuous symmetry can be encoded in a Matrix Product Operator (MPO) as depicted in Fig. 2, such that: For the model we study, the MPO has bond dimension three (γ = 1, 2, 3) and the non-zero elements of the tensors B for the odd sites can be expressed as: The tensors for the even sites can be computed in a similar way. Given the gauge invariant tensor structure introduced above, one can reformulate the standard algorithms for tensor networks to compute the ground state of the model, or, as we will do in the next Sections, compute the real time evolution of some initial state, e.g., by means of a Suzuki-Trotter decomposition of the time evolution operator acting on a pair of neighboring sites [76,77]. The parameters used in our calculation, together with a discussion of the errors involved in the approximations, are presented in Appendix A. Finally, we mention that one can further simplify the tensor structure in Eq. (9) to increase the algorithmic efficiency and that this construction is completely transparent with respect to the Abelian or non-Abelian nature of the gauge symmetry, resulting in a drastic simplification of numerical analysis of non-Abelian lattice gauge theories [31,34].

III. STRING BREAKING
In this section we present our results of the lattice gauge TN numerical simulations for the real-time dynamics of string breaking: we focus on the time evolution of the electric and matter fields, quantitatively studying the properties of string breaking and of the Schwinger mechanism. In the following section we analyze the time evolution of quantum correlations during this processan analysis enabled by the TN approach -and show that the two figures of merit are intimately related.
The setup for our simulations is a dynamical string surrounded by the vacuum. The total lattice length N = 100 is chosen such as boundary effects are negligible for the timescales we investigate. Some typical results of the electric field time evolution are shown in the left column of Fig. 3 (Column A) for different values of the fermion mass m and electric field coupling strength g (hereafter  Fig.1). The electric flux real-time evolution (Column A) is shown with the evolution of the mass excitations (Column B) and the evolution of the bipartite von Neumann entropy (Column C) for m = 0, g = 0 (Line 1), m = 0.25, g = 1.25 (Line 2), and m = 3, g = 3.5 (Line 3). The von Neumann entropy is calculated using a bipartition of the system defined via a cut between lattice sites x and x + 1.
we set t = = 1 and times are given in units of /t). In the top row, the string freely breaks as for g = m = 0 no energy-cost is needed for such process to occur and the system evolves according to the free hopping Hamiltonian: the electric field starts to oscillate between string and anti-string displaying primary, secondary (marked by vertical lines) and subsequent string breakings. Moreover, the string propagates into the vacuum creating two diverging electric field excitation wavefronts. In the middle row, representing the same scenario for non-zero mass and electric field coupling strengths, string breaking still occurs; however, the string evolution is damped and stabilizes the mean electric flux around zero. Finally, for large mass m and g, the large mass suppresses vacuum fluctuations while the large electric coupling prevents the occurrence of string breaking and thus the string does not decay. However, some dynamics still occur as we will see in more details in subsection III B, as fermion-antifermion pairs are created using the energy of the electric field and then annihilated, thus restoring the electric field excitations.
To perform a quantitative analysis of string breaking, we repeat the same simulation for different values of m and g and analyze the mean electric field of the central six lattice sites of the string as a function of time. The results are reported in Fig. 4 where one can clearly see that different scenarios might occur: either the string breaks when the mean electric field drops below zero, or the electric flux remains positive throughout the whole evolution and no string breaking occurs. In the limit of m = g = 0 the oscillations show the highest amplitude, which is reduced by changing at least one of the two parameters; as previously noticed, for high values of either system parameter the string never breaks. In the regime between these two extreme cases, we observe a third type of behavior: the electric flux tends to zero, however no anti-string is formed: the oscillation is strongly damped and the system remains in the broken string state (compare with the middle row in Fig. 3). These findings are summarized in Fig. 5: For g, m 1 we observe the full string breaking dynamics with at least the partial formation of a string with opposite electric field flux after reaching the broken string state (red area). For g, m 1 string breaking was not observed and the dynamics are dominated by the interplay of the state of maximum pair creation and the original string (green area). Finally, the white area in between represents the region of parameters where we observe the string breaking with over-damped oscillations.

A. String wavefront spreading
During the string breaking process, a wavefront of electric flux spreads outwards, as can be clearly seen in Fig.  3 (Panel A1 and A2). In this Section, we quantitatively characterize the wavefront spreading by analyzing its spreading velocity and the oscillation intensity.
In Fig. 6, we show the wavefront propagation as a function of time for different electric coupling g for the zero mass case. The lower inset illustrates how we calculated such propagation: we follow the electric field excitation on one side of string by means of tracking the difference ∆E between the gauge field at some position x and the next nearest neighbor site x + 2. Further, we define the  −0.15Emax). The white area represents the parameters where the mean electric field approaches zero and stays around that value (|Emean| < 0.15Emax). And finally, the green area represents the parameters without string breaking (Emean > 0.15Emax). time when this difference displays a maximum s arrival of the wave front [90]. As can be seen from the lower inset of Fig. 6, where different colors represent different coordinates x, a wave-like propagation can be clearly identified.
Following this scheme, we plot the position of the wavefront as a function of time in the main panel of Fig. 6. The result shows an approximatively linear spreading after an initial transient time of about τ ≈ 2, with a velocity almost independent of the values of the electric coupling for g < 1. Increasing g starting from g = 1, the velocity increases as well, until for g > 1.5, where the results are inconclusive as increasing g leads to smaller wavefront amplitudes and consequently the errors bars prevent an accurate analysis to be carried out. However, for small g, the spreading velocity can be extracted directly from Fig. 6. By fitting the values for τ > 2/t and m = 0 we obtain a value of v E = 1.96 ± 0.02.
Finally, in the upper inset of Fig. 6, we repeated this analysis for different masses and g = 0. The results clearly show that, for sufficiently large m/t 4, the wave front spread velocity has an inverse linear dependence on the mass. All these results are in agreement with a theoretical estimate obtained assuming the ends of the string as sources of excitations: in a quasi-free or weak coupled model, the speed is related with the band-width of the kinetic term resulting on an excitation spreading velocity proportional to v th = 2/m.

B. Schwinger mechanism
During string breaking, the Schwinger model dynamics exhibit particle-antiparticle pair production as a consequence of the energy released from the external electric field string. This phenomenon is usually referred to as the Schwinger mechanism, and has been studied extensively  since its first presentation in 1951 [59]. In the following, we provide a systematic investigation of the Schwinger mechanism in the context of the U(1) QLM. In our analysis, in the initial state defining the string, the only two mass excitations present are the two dynamical charges which create the string itself. However, during the dynamics, the energy of the string is transformed into mass excitations. When the maximum mass production is reached, the particles start to annihilate, either to break the string or to restore it. The time needed to reach the maximum mass production τ max depends on the mass and electric coupling as shown in Fig. 7. It clearly displays two different behaviors, depending whether the electric coupling g is greater or less than one. For small-g, τ max is maximal for m = 0 and decreases monotonically with m. This occurs, recalling the results of Fig. 5, in the regime where string breaking is observable. Indeed, these are the cases comparable to the top row in Fig. 3. On the contrary, for larger values of g, we observe the maximum to occur for m > 0. In particular, the maximum is obtained at a point where the energy of the electric field matches approximately the energy needed to fully convert the string into particle pairs, i.e., m = g 2 /4. This corresponds to the dynamics as displayed in the bottom row in Fig. 3. In the regime of high masses we can use second-order perturbation theory to estimate the mass-dependence of the mass-production time-scale analytically, as the model can be approximated as decoupled double-wells potentials, resulting in τ max = π/2 m ≈ 1.57 m . We checked this approximation comparing the analytical estimate with the numerical results in Fig. 8: The best fit results in τ max = (1.54 ± 0.02)/m, in good agreement with the theoretical prediction.

C. Observability of string breaking in synthetic platforms
Recently, the implementation of Abelian quantum link models has been envisaged on different platforms, such as ultra cold atom gases in optical lattices [45,46,52], trapped ions [48], and circuit QED architectures [47,51]. In this context, the possibility of investigating the realtime dynamics using TN methods provides an invaluable tool to benchmark experiment against theory in (1+1)d, and to address the role of possible imperfections in quan-tum simulators.
Unavoidable imperfections which will be present in any implementation can be detrimental to the observation of both ground state physics properties and real-time dynamics. The former case has recently been investigated in the context of adiabatic state preparation, where it was shown how gauge-variant perturbations weakly affect the fidelity of the loading process [33]. Here, we focus instead on the effect of gauge-invariant imperfections on the string-breaking dynamics. Different from gaugevariant terms, the role of gauge-invariant imperfections cannot be systematically addressed in an experiment using, e.g., post-selection over the experimental data.
Following the implementation schemes in Ref. [45,47,48], one of the most common form of gauge-invariant imperfections are nearest-neighbor interactions between matter and gauge fields: with n x = ψ † x ψ x . This form of the imperfection is usually generated as a resonant term in perturbation theory to next-to-leading order with respect to t. While this implies t ξ, for realistic implementations the difference in magnitudes between these two terms cannot be made arbitrary large: this will require small absolute energy scales, thus making other sources of more detrimental imperfections such as temperature (for the cold atom implementation) and disorder (in the circuit QED implementation) dominant. At a qualitative level, this interaction term can freeze the system into a configuration where the matter fields remain pinned. This is due to the effective attraction generated by the nearby electric field configuration.
To estimate the effects of these imperfections on a quantum simulation of the dynamics considered in this work, we repeat the numerical simulations including realistic imperfections expected in a first generation of experiments. In the top row of Fig. 9 we show the string breaking evolution of the electric field, in the presence of H I , with the same system parameters as used for Fig. 3, with imperfections of the order of 10% ξ = 0.1. The results including the imperfections still clearly exhibit the physics observed in the imperfection-free results. In general, even the quantitative dynamics is very well captured up to long-timescales, as illustrated in the bottom row of Fig. 9. The only exception are intermediate g and m values, where significant discrepancies (up to 50%) are observed for intermediate timescales (middle panel of the last row). In all other regimes, we could observe deviations up to a maximum of 15% caused by typical imperfections ξ = 0.1.

IV. ENTANGLEMENT DYNAMICS
One of the key aspects of MPS-based methods is that they give full access to the wave-function during the real- In the last decade it has been shown that entanglement play a fundamental role in many-body quantum processes, from quantum critical phenomena, to quantum information theory and other fundamental aspects of quantum physics. Moreover, several aspects of quantum field theories, such as the static properties of conformal field theories, have also been extensively studied using entanglement measures [74,75]. Additionally, entanglement was shown to play a crucial role in the limits of classical simulations of quantum systems calling for the need of the development of quantum simulators to overcome such limitations [12].
A common way to quantify entanglement for pure quantum states is by using the so-called Renyi entanglement entropies, and, in particular, the von Neumann entropy. Given a pure state with the density matrix ρ = |ψ ψ|, the entanglement entropy is given by the von Neumann entropy of the reduced density matrix ρ(x) = Tr L−x ρ [75]: S(x) is a measure of the entanglement of a bipartition at the lattice site x. The entanglement entropy takes values between S(x) = 0 for a separable state (product state) and S(x) = log 2 d, with d being the size of the Hilbert space, for a maximally entangled state.

A. Von Neumann entropy after string breaking
In Fig. 3 (Column C) we plot the time evolution of the entanglement entropy for the three different cases considered before (Panels 1-3): as it can be clearly seen, the entanglement evolution resembles the mass and electric field excitations, indicating how the two phenomena are related. Firstly, the vacuum fluctuations for small g and m generate not only mass and electric field fluctuations, but also a high amount of correlations. Moreover, the electric field wavefront is mimicked by the entanglement behavior, once again showing that not only do excitations propagate from the string, but also correlations do as well. Secondly, the string breaking process is clearly a two-steps process: first a correlated pair is created in between odd-even sites, and later in between the evenodd sites (blue-yellow checkerboard pattern inside the string in Panels C1 and 2). Finally, for large g and m, when the string does not break, the correlations behavior drastically changes as well. There, within the string, the correlations are built periodically only between odd-even sites, while in the vacuum the correlations are drastically suppressed (Panel C3).
A transparent picture on how entanglement is generated during the quench dynamics can be gathered by monitoring the entanglement growth at different points in space. In Fig. 10 we show the entanglement time evolution for a set of system bipartitions and parameters, e.g., cutting the system in the middle of the string or in the vacuum.
The solid lines are results obtained for m = 0 and g = 0, while the dashed lines correspond to the case m = 3 and g = 3.5. The orange line represents the entanglement growth for a partition at lattice site x = 20, therefore in the center of the lattice region starting in the vacuum. In the zero-mass case, one can see that the entanglement entropy grows almost linearly for most part of the evolution -signaling the vacuum instability against resonant particle-antiparticle pair production. Towards the end of the evolution, the linear growth breaks down as boundary effects start to play a role. The remaining solid lines in Fig. 10 represent the behavior of a partition between an even-odd lattice site (violet) and between an odd-even lattice site (blue) and display a more complex behavior. These results were obtained in the center of the string while it breaks up: the counter-phase oscillations indicate the competition of two states, together with an overall growth of the entanglement entropy, though not as fast as in the vacuum. As we have seen, in this regime the string breaking is a result of consecutive hopping processes: fermions hop on the lattice creating mass excitations followed by annihilations, which result in the string breaking with two remaining dynamical mesons. This dynamical process goes on and after two hopping processes the anti-string is created. Fig. 10 shows that these dynamics are well captured by the oscillations of the entanglement entropy. Each maximum represents one hopping process: at the first maximum of the blue line the fermions are about to hop for the first time, while at the first maximum of the violet line the second hopping process is at its peak resulting in the string broken state. At around τ ≈ 4, the blue line does not display a maximum: this 'depleted region' signals the anti-string state, where the last hopping event creating the anti-string is again the first hopping to break the negative string (and thus, it takes twice as long for the entropy to reach the next maximum).
Differently from before, in the massive scenario (dashed lines), we see that the entanglement entropy for the vacuum stays close to zero as the large mass and electric coupling strongly suppress the particle-pair creation that triggered the strong growth of the entropy in the previous case. Moreover, in the middle of the string the entanglement entropy is drastically affected: the blue dashed line initially behaves as the solid line in the massless case, reflecting the same mass excitation by pair creation. However, the violet dashed line always remains close to zero as further evolution into the string broken state is energetically forbidden: the state evolves back into the string and the correlations between the even-odd sites cannot be created. The system is then oscillating between two almost degenerate states: the initial string state and the state made out of pairs. This results in the oscillating behavior of the entanglement entropy between zero and one. Finally, the third case with m = 0.25 and g = 1.25 (dot-dashed lines) lies between the two previous limiting cases: here the string breaks, but does not evolve into an anti-string. In the vacuum, the entanglement evolution is very similar to the first case (m = 0 = g) as the entropy grows almost linearly after a transient regime.
However, the slope is reduced by the nonzero mass. The entanglement in the center of the string initially evolves as for the massless case, but after the first two hopping processes for τ 2, the oscillation turns into vacuumlike growth. This is a strong indication for non-periodic string breaking: the dynamics, although being unitary, resemble a dissipative process where the electric field energy irreversibly disperses into the vacuum. This irreversible behavior directly resembles what we observe in the electric field dynamics, which does not display any clear periodic signature.

B. Entanglement propagation and wavefront
Even more remarkably, the real-space particle creation and the entanglement dynamics are quantitatively tied. We concentrate on the signatures of the wavefront of the string imprinted on the evolution of the entanglement entropy. We consider the case m = g = 0 as it is characterized by the most pronounced wavefront, where the string with its slow entanglement growth is embedded in the fast growing vacuum (see Fig. 3, panel C1).
To characterize the entanglement spreading due to the wavefront, we exploit the fact that the entanglement entropy in the vacuum is constant in space even though it evolves in time. Therefore, far enough from both sides of the string there is a plateau of constant entropy much higher than the entropy in the middle of the string. Thus, to define the wavefront of entanglement spreading due to the string, one can look for the lattice site at which the entropy plateau starts to decrease. We identify this point computing the difference of entropy between nearest neighbor bipartitions: tracking when this quantity become bigger than a given threshold allows us to characterize the entanglement wavefront spreading.
In Fig. 11 we show the estimated spreading velocity for different values of the threshold: the limit for the threshold value going to zero gives an estimate of the spreading velocity. A power law fit results in a spreading velocity of v S = 2.0 ± 0.2 in very good agreement with the analytic estimate of v T 2 and the result from the electric field of v E = 1.96 ± 0.02 demonstrating the intimate connection between entanglement and electric field spreading.

V. MESON SCATTERING AND ENTANGLEMENT GENERATION
Bound states are a fundamental component of gauge theories, and understanding their complex internal structure is one of the most ambitious goals of computational physics [83]. In experiments, such internal structure is usually explored by colliding heavy ions, so that the energy released during the process can be released via particle-antiparticle creation. This makes the ab initio numerical simulations of such scattering processes chal- lenging, as MC simulations suffer from a severe sign problem when tackling real-time dynamics.
Here, we show how TN simulations allow to investigate meson scattering for the (1+1)d QLM using TNs. First, we present a general procedure to implement scattering processes between composite particles, and discuss the electric field dynamics after the collision. Then, we present results for the entanglement dynamics during and after the collisions showing that the meson collision is accompanied by the creation of entanglement between the mesons themselves. Indeed, as we will show, the entanglement is bounded by the propagation wavefronts of the particles after collision, and is characterized by a constant plateau of the entanglement entropy within the region.

A. Electric field patterns during meson collisions
In order to produce the scattering process, we shall start with two composite particles. Each particle is charge/anti-charge pair, and divided only by one link, namely a meson, with opposite momentum. For the twomeson problem, there is a simple picture in the strong coupling limit: the massless theory is a free massive boson (meson) theory that is expected to become weakly interacting once a small mass term is included. Hence, in the strong coupling region, a possible two-meson bound state is loosely bound, while in the weak coupling region it is tightly bound.
We start the numerical simulation with the state represented in the cartoon (D) in Fig. 1: two mesons separated by a vacuum state of ten sites, which can be straightforwardly be written in a simple, separable matrix product state with t = 0. We provide momentum to the mesons by adiabatically moving them from the boundaries toward the center of the system: this is done by introducing a deep box-shaped potential that decouples the mesons from the rest of the system, with the only dynamics allowed being the oscillate between its position and a neighboring site. The box-potential is removed at time τ i = 17.4 when the meson is exactly at half oscillation: from that point on the mesons evolve freely with an effective momentum mostly in one direction, one towards the other and eventually colliding [91]. In order to avoid vacuum fluctuations during the process, we choose a large value of g = 8. Fig. 12 shows an example of such a scattering process. In particular, it shows the absolute value of the electric field of two mesons approaching each other, colliding in the center and there parting again. While before the collision the meson are tightly bound, after the scattering process the electric field diffuses, and the corresponding wavefront has a significantly attenuated signal. In the lower panel of Fig. 12, we monitor the time-evolution of the total particle number (blue), clearly indicating that this quantity is approximately conserved over the entire time-evolution, due to the large electric field strength, which suppresses particle-antiparticle creation.

B. Post-collision entanglement generation
A classical-like picture of the scattering process presented above, would have that two particles move against each other and then bounce back as there is not enough energy available to generate a more complex inelastic scattering. However, this picture is oversimplified, as this is a fully quantum process. Indeed one can, once more, monitor the quantum correlations generated during the scattering process. This is done in Fig. 13, where we show the evolution of the bipartite entanglement entropy: one sees that entanglement is created and that it is mostly carried by the two mesons. In this parameter regime, the vacuum does not generate entanglement due to the very large value of g 2 . Studying the bipartite entanglement entropy for different bipartitions and times, one clearly sees that there are two regimes: before the scattering occurs, the entanglement is present only in the bipartition that cuts the mesons wave packets, indicating two electron-positron wave packets internally correlated but not sharing any quantum correlations among them. To the contrary, after the scattering, the two wave packets become highly correlated even when their two centers of mass are clearly separated (see Fig. 12 for times τ > 100).
The values of the entanglement entropy indicate that one ebit of quantum information has been created dur-ing the scattering process. In the right panel of Fig. 13, we present various cuts of the entanglement entropy profile taken at different times, together with a comparison with the entanglement generated by a single meson moving through the lattice (dashed line). The difference of ∆S ≈ 1 between the two cases (highlighted in Fig.  13, green bar) clearly shows that one additional ebit of entanglement has been generated during the scattering process: that is, that a singlet state has been created between the two indistinguishable mesons. The entropy has increased as the information on which process occurs (either mesons bouncing and moving back, or each one moving freely through the other one and keeping moving with its own momentum) is completely lost. Notice that this is a direct consequence of the wave nature of the mesons: the classical case of two scattering particles in one dimension -with the constraint that no double occupancy might occur in a single matter site -would necessary results only in backscattering, since the two hard classical particles cannot pass through each other.

VI. CONCLUSIONS
We presented a detailed tensor network study on the real-time dynamics of a lattice gauge theory in the presence of dynamical charges and quantum gauge fields. Within this approach, we have shown that one has direct access to all local quantities of interests -the time evolution of the mass, charge and gauge fields -and to the quantum correlation between bipartitions of the system by means of the von Neumann entropy. We investigated the primary and secondary string breaking in QED in (1+1)d represented by an S = 1 quantum link model with staggered fermions. In this context, we studied the real-time evolution of the Schwinger mechanism leading to mass creation and annihilation by means of the interplay with the electric energy released by the string. We quantified key properties of these effects such as the mass production rate of the Schwinger mechanism and the velocity of the electric field spreading. Moreover, we unveiled the relation between string breaking dynamics and the entanglement spreading in the systems and we showed that it is possible to study scattering dynamics, characterizing not only mass and charge real-time evolution but also the creation of quantum correlations among scattered particles. Finally, we showed that the presented results can be in principle verified experimentally in possible future quantum simulations, as they appear to be robust with respect to the most common sources of gauge-invariant imperfections appearing in most of the proposed implementations.
This work paves the way to systematic studies of real time dynamical phenomena in Abelian and non-Abelian LGTs in low dimensional systems. Indeed, the present approach can be straightforwardly generalized to more complex LGT and geometries, e.g., ladders or cylinders, and also can be studied in presence of an external en-vironment by means of, for example, the tensor network approach presented in [84]. Moreover, one can also study the continuum limit of the LGT (as already discussed for Wilson theories [30]): for QLMs, this can either be done using dimensional reduction [56], or by increasing the quantum link representation, similarly to what has been done in [29]. Notice that the gauge invariant tensor network formulation behaves favorably as the speedup it grants scales as the number of rishons per link squared [34]. The unprecedented access to the entanglement dynamics in LGTs will allow investigations on the role of quantum correlations in different contexts enabling a deeper understanding of the quantum real-time dynamics of lattice gauge theories. Finally, exploiting the capability to prepare a wide class of complex states granted from recent developments in quantum optimal control of many-body quantum systems [85,86], more complex dynamics could be investigated. One example would be to perform extensive studies of the scattering at different energies.
Alongside, these methods can also be applied to study condensed matter systems to compute, e.g., response functions of antiferromagnets described by LGT (spin ices, Resonating Valence Bond models, etc.) which are very hard to evaluate using MC due to analytic continuation [87]. Moreover, real time dynamics of gauge theories is fundamentally interesting to ab initio investigations of scattering equilibration and pre-thermalization [56]. Finally to benchmark and verify small quantum simulations whose proposal have recently appeared for different platforms (ranging from cold atoms to superconducting circuits and trapped ions) and that can be foreseen to be experimentally implemented in the next years [41][42][43][44][45][46][47][48][49][50][51][52][53][54][55].
While writing this manuscript, we became aware of a related work on string breaking using TN methods [88]. We investigated String-breaking in a one-dimensional system with up to 100 lattice sites. The simulation was performed using a matrix-product-state (MPS) algorithm using a second order Suzuki-Trotter decomposition for the time-evolution. In this appendix we present the parameters we used for our simulations and provide a discussion on the relative errors.
The main sources of numerical error in our calculations are the finite bond dimension of the MPS-state representation and the Suzuki-Trotter time-step. As a bond dimension we used a value up to χ = 200, which ensures a truncation error on the corresponding wave function of maximum order 10 −8 and 10 −3 for τ = 5 and 9 respectively for the string breaking calculations, and 10 −5 for τ = 250 for the scattering processes. Examples of the change of the truncation error during the evolution can be seen in the inset of the left panel in Fig. 14.
In order to make sure that these errors lead to small changes in our main observables, we tested the convergence of the mean electric field in the center of the chain (cf. Fig. 4) for different values of χ and typical system parameters. The results can be seen in the left panel of Fig. 14 where we plot the mean electric field at the end of the evolution (τ = 8) as a function of the inverse bond dimension 1/χ. The mean electric field was subtracted by a fitted offset E lim to allow a better comparison between the two sets of system parameters. As we see from the plot, even for rather small bond dimensions χ < 100 the change to the largest bond dimension used (χ max = 220) is in the order of 0.1 − 0.01%. For the bond dimension used in our simulations (χ = 200) Two numerical lattice sites are combined to form one physical lattice site due to our staggered fermions formulation. The dynamical charge density decays exponentially with the distance from the charges, the black line represents the theoretical model with q(x) = a exp(−bx) − a exp(bx) where the decay parameter is b = g/ √ π and the scale parameter is a = exp(b) − 1.
the difference to the extrapolated correct value is in the order of E − E lim ∼ 10 −5 . As said before, we used a second order Suzuki-Trotter decomposition with a time-step of δτ = 0.01 to simulate the time-evolution. In the right panel of Fig 14 we report the convergence test obtained repeating the same simulation with different time-steps. As expected, we find a clear E − E lim ∼ δτ 2 dependence: This clean powerlaw behavior allows us to give a very good estimate of the correct value of E lim . In summary, the error from the Suzuki-Trotter decomposition using a time-step of δτ = 0.01 is of the order of E − E lim ∼ 10 −5 , that is, of the same order as the error introduced by the truncated bond dimension. we limit our system size to 20 lattice sites and encode the two static charges which build up the initial string in the boundary conditions on both sides of the lattice. In such a way, we avoid to simulate the vacuum regions where quantum correlations build up very quickly (see Fig. 3) and are less interesting for the phenomena under study. Thus, we simulate the time evolution of the string of electric field among the charges with a small, but finite mass: after a long time compared to other timescales of the system (τ = 400) when the largest fluctuations in the lattice damped out, and report the final charge distribution q x (averaged over the last τ = 10 to remove remaining oscillatory effects) in Fig. 15. The net results is an exponential decay of the charge density with the distance from the charges, as discussed in Ref. [89] where the authors showed that the screening also at m = 0 is equal to that experienced by the massless Schwinger model, as q x = a exp(−bx) with a decay rate b = g/ √ π ≈ 0.5642 [61]. The parameter a = exp(b) − 1 is derived from the normalization condition ∞ x=1 q x = 1.
It can be clearly seen that our findings are compatible with the theoretical results from Ref. [89] (black solid line in Fig. 15).
We stress that the goal of the presented analysis is to show that it is possible to perform a study of the screening of charges using tensor networks, not to perform a thorough discussion which will be presented elsewhere. Indeed, an extensive numerical analysis is needed to individuate or rule out possible deviations of the massive from the massless case in a full quantum mechanical model, the influence of the spin representation in the quantum link formulation, possible finite size effect, and also more challenging regimes of parameters, e.g. m ∼ 1 as shown in [89]. However, the agreement between the presented results and the semiclassical approach, indicates that both analysis are capable of describing the physics of the system and might be alternative approaches which can complement and benchmark each other.