Non-Abelian string breaking phenomena with Matrix Product States

Using matrix product states, we explore numerically the phenomenology of string breaking in a non-Abelian lattice gauge theory, namely 1+1 dimensional SU(2). The technique allows us to study the static potential between external heavy charges, as traditionally explored by Monte Carlo simulations, but also to simulate the real-time dynamics of both static and dynamical fermions, as the latter are fully included in the formalism. We propose a number of observables that are sensitive to the presence or breaking of the flux string, and use them to detect and characterize the phenomenon in each of these setups.


Introduction
The phenomenon of string breaking, as a consequence of confinement, is one of the most fundamental aspects of the gauge theories that lie at the basis of our understanding of high energy physics. While a comprehensive understanding in many cases, such as QCD, is still lacking, a lot of insight has been gained thanks to lattice gauge theory (LGT). In LGT, originally pioneered by Wilson [1], the theory is formulated on a discretized space-time lattice in such a way that the gauge symmetry is preserved and the continuum model is recovered in the limit of vanishing lattice spacing. This allows powerful Monte Carlo simulations which besides mass spectra [2] and phase diagrams [3] have been also successfully used to study static properties of confinement and of string breaking [4][5][6][7][8][9]. Despite the great success of this method, it suffers from the sign problem [10], which limits the kind of feasible simulations, and in particular real-time dynamics (although recently introduced techniques have allowed some real-time calculations in certain regimes [11,12]). It is then of undeniable interest to explore other techniques that can overcome the limitations of standard LGT techniques.
Another direction explored in recent years is quantum simulation of gauge theories [17,[28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47]. The proposals for quantum simulators typically work as well with the Hamiltonian formulation of the LGT and map the gauge theory to the physical degrees of freedom of another quantum system which can be controlled experimentally thus being free from purely numerical limitations. Although a fully fledged quantum simulator is not available yet, this route seems promising for the future as the first proposals might come into reach with nowadays high level of experimental control.
In this work, we consider a non-Abelian lattice gauge theory with dynamical charges in one spatial dimension. The Hamiltonian we simulate realizes an exact SU(2) gauge symmetry using a finite-dimensional representation for the bosonic degrees of freedom [47,48]. This model is hence a suitable candidate for the design of an atomic quantum simulator of SU (2) LGT. The finite dimension of the link Hilbert spaces allows the simulation of the model via direct application of MPS techniques, which we employ here in order to numerically explore the statical and dynamical aspects of the string breaking phenomenon.
In order to detect and characterize the breaking of the string, we propose three different observables. They can be monitored during the MPS simulation, and should in principle also be accessible in an experiment, so that our results will be measurable in a potential future quantum simulator.
More specifically we study three different scenarios. First, we consider the static aspects of string breaking by determining the ground state of the system, which includes fully dynamical fermions, in the presence of two additional external static charges. This calculation relates closely to the ones traditionally accessible by lattice Monte Carlo methods. We show that we can reliably determine if the string is present, and therefore identify the regions where we expect string breaking to occur. These calculations additionally demonstrate the suitability of the proposed observables for the detection of string breaking in dynamical scenarios. In the second place, we investigate the dynamics of string breaking by introducing the external static charges on top of the interacting vacuum and evolving the state in real time. When the string breaks, we can explicitly observe the screening of the charges via the creation of dynamical particles. Finally, we analyze how the picture changes when the charges added to the vacuum are themselves dynamical, a scenario which is closer to more realistic out-of-equilibrium situations. Also in this case we can recognize the string breaking if the fermion mass is small enough.
The rest of the paper is organized as follows. In section 2 we review the model and describe the numerical methods we are applying, with a special focus on the MPS techniques we are using. In section 3 we present our numerical results for the static calculations. Section 4 contains the results for the real-time dynamics of string breaking with two static external charges. In section 5 we present the results for the real-time dynamics of string breaking when a pair of dynamical charges is added to the interacting vacuum. To conclude, we summarize our findings in section 6.

Truncated theory
The model we consider is the Hamiltonian formulation of a 1+1 dimensional SU (2) LGT with dynamical fermions, in which the gauge symmetry is exactly realized with finitedimensional link variables [47]. The model can be understood as a truncation of the complete SU(2) gauge theory, formulated as the Kogut-Susskind Hamiltonian with staggered fermions [49] H =ε where ψ n = ψ r,n ψ g,n is a two component spinor taking into account the two "colors" of fermions ("red" and "green") on site n. U n is a matrix on color space in the fundamental representation. Its entries, U ij n , act on the link between the fermionic sites n and n + 1 and change the corresponding electric flux. The operator J 2 n is a group scalar, as it does not carry any color index, and gives the flux energy on link n. Physical states, |φ , have to be eigenstates of Gauss Law, i.e. G α n |φ = q α n |φ for all sites, n, where G α n = L α n − R α n−1 − Q α n , α ∈ {x, y, z}.
Q α n = 1 2 ψ † n σ α ψ n is the non-Abelian dynamical charge at site n, where σ α are the usual Pauli matrices. A nonzero value for q α n indicates a non-vanishing static external charge. 1 L α n , R α n are the generators of left and right gauge transformations acting on link n, which correspond to the left and right electric field on the link and fulfill the commutation relations [L α n , L β m ] = −iδ nm γ ε αβγ L γ n , [R α n , R β m ] = iδ nm γ ε αβγ R γ n and [L α n , R β m ] = 0. They are related to the operator for the electric flux energy as J 2 n = α L α n L α n = α R α n R α n . The Gauss Law components do not commute among themselves, as [G α n , G β m ] = −iδ nm γ ε αβγ G γ n , and thus cannot be diagonalized simultaneously. However, they all commute with the Hamiltonian, [H, G α n ] = 0, so that the Hilbert space is a direct sum of sectors characterized by their configuration of external charges, {q α n } [47]. In the strong coupling limit, corresponding to ε = 0, the ground state is given by links carrying no flux, odd sites occupied with two fermions, and empty even sites. This state can be written in a product basis as In the previous expression bold numbers represent the occupation number of a fermionic site, while |0 represents the state of a link with no flux. The fermionic part of this state corresponds to the Dirac sea, and it is easy to check that it fulfills Gauss Law with G α n |φ SC = 0 ∀n, α. Applying on this state the gauge invariant string operator or its adjoint, S † nl , for odd l, 2 generates a particle-antiparticle pair at sites n and n + l, connected by a flux tube of length l. Such configurations will have an excess of energy of 2m due to the particle-antiparticle pair, plus a flux energy proportional to the string length, l. Consequently, from a certain length on, it will be energetically favorable to reduce the flux energy by creating extra particle-antiparticle pairs, leading to configurations with a broken string.
In the complete SU(2) model, the flux on a link is not bounded, and the dimension of the Hilbert space for each link is infinite. It is nevertheless possible to consider a theory truncated in a gauge invariant manner [48], where the maximum flux a link can carry is limited and the Hilbert spaces of the links are finite dimensional. Following the method in Ref. [48], each matrix element U ij n can be decomposed as a sum over all irreducible representations, which may be separated to summands that are gauge invariant themselves.
Here we consider the model corresponding to the simplest non-trivial truncation of the full theory, meaning that only the trivial and the fundamental representation are kept, resulting in dimension 5 for the links.

Numerical Methods
We use the MPS ansatz to describe the system [50]. A MPS with open boundary conditions for N sites is given by where the A i k k are complex matrices in C D×D for 1 < k < N and The string breaking phenomenon can be studied statically or dynamically. Different well known numerical algorithms exist to find MPS descriptions of stationary or timedependent states [26,51]. In this work, we employ simulations of time evolution [25,26,52] where the evolution operator is approximated via a first order Taylor expansion. This allows us to study both static and dynamic scenarios by respectively using imaginary or real time. Additionally, the Taylor expansion preserves the symmetries of the Hamiltonian, such as gauge invariance. Then it is possible to explore a specific sector of the external charge distribution without explicitly implementing symmetries in the tensors.
The numerical simulations have three main sources of error, each of them controllable by a suitable choice of parameters. The first one is due to the approximation of the evolution operator via a Taylor series. This error can be controlled by choosing the time step suitably small. Another source of error is due to the limited bond dimension. By obtaining results with different bond dimensions, the size of this truncation error can also be estimated and controlled. As we are working with finite systems, a third source of error arises from finite size effects which can be avoided by using sufficiently large systems.
Although TN and in particular MPS can be formulated in terms of fermionic degrees of freedom, we choose to translate the fermionic degrees of freedom in Hamiltonian (2.1) to spins via a Jordan-Wigner transformation for our numerical simulations (details about the procedure and the relevant operators in spin formulation can be found in the appendix A) and in the following we work with the spin formulation.

Detection of string breaking
Detecting the presence or absence of the flux string during the evolution requires observables that are suitable for the dynamical case. The Wilson loop and more general correlation functions, widely used in lattice calculations to determine the static potential, are typically evaluated in the limit of large Euclidean time and therefore not suitable for our setup [5,8,53]. Instead, to detect strings and string breaking in the dynamical setup, we propose three different observables. We monitor their values throughout the evolution, taking advantage of the fact that we have access to the MPS wave function at all times.
In the first place, the spatially resolved spin and flux configurations in the system allow us to visualize the change with respect to the initial configuration (see e.g. figure 3).
A second observable is the local imbalance between red and green species, also spatially resolved. A string operator (2.2) changes the fermionic content of the sites at its end-points, such that, on the strong coupling vacuum, it will produce a superposition of states having a single red or a single green fermion at the beginning or at the end of the string. This can be detected by the operator Q α n 2 = 1 4 (n r,n − n g,n ) 2 , where n r,n and n g,n are the occupation numbers for the two fermionic species on site n. 3 A third observable can be proposed that looks only at the flux content of the links. Applying a string of length l, starting at site n, on the strong coupling vacuum produces a state in which links between sites n and n + l carry non-vanishing flux, whereas outside the region there is no flux. We can construct projectors P nl on this kind of configurations. However, we are not interested in a single string but rather in the statistics of string lengths 3 One should note that the index α in Q α n 2 is not summed. This might look puzzling at first, as it seems that we are using a non color-neutral object as observable. However, a little calculation shows (see appendix A) that Q α n 2 is identical for all α. Therefore summing α would only yield an additional factor of 3 which we are not taking into account. in a (time dependent) state |φ(t) . Thus, we bin all strings of a certain length together and normalize by the number of possible strings of length l: In this way we can obtain histograms for the distribution of the string lengths in the state at a given time t.
As we show in the next section, these three observables allow us to reliably determine whether our system still contains the initially imposed string or the string is broken.

Ground state with static external charges
The usual way of probing for string breaking in lattice Monte Carlo studies is the analysis of the static quark-antiquark potential. Using the MPS method, we can also determine the energy of an extra pair of external charges as a function of their distance by simply computing the ground state in a sector, where two external static charges are placed at the desired separation.
In particular, we choose two static external charges q y = ±1/2, located in the central region of the system (to minimize finite-size effects) at a distance l. Gauge invariance requires that there is a flux tube connecting them. In the strong coupling vacuum the charge at each site is zero, as well as the left and right electric field on each link. Consequently, we can prepare a state with the desired external charge distribution by applying the (noninvariant) operator U n . . . U n+l−1 on the strong coupling ground state. This operator only acts on links from n to n + l − 1, on which it creates a finite electric field. Hence the Gauss Law at the beginning and the end of the string yields a non-zero value corresponding to a state with q y n = ±1/2 and q y n+l = ∓1/2. Subsequently we evolve this state in imaginary time to determine the ground state of (2.1) for a chosen set of parameters, (ε, m, g). For the results presented in this section we use a time step ∆t = 1.0 × 10 −3 and a bond dimension D = 100, parameters which turn out to be sufficient to avoid noticeable numerical errors (see appendix B for a more detailed error analysis).
In figure 1 we compare the ground state energy E for ε = 3.0, g = 1.0 and various masses, where we subtracted the energy of the interacting vacuum without external charges E vac . For m = 3.0 and m = 5.0 we observe three regions as functions of the string length. For short strings the energy grows linearly with the string length, indicating the stretching flux tube between the charges, namely the presence of the string in the ground state. From a certain value l c on, the ground state energy does not depend on l. This is the signature for string breaking, as after reaching the threshold for creating particle-antiparticle pairs, reducing the flux is energetically favorable, the string breaks and the energy is independent of the initial string length. Finally we observe a third region, when the string length is already close to the system size, and finite size effects become noticeable. In our plots we clearly see that the values of l c , where one transitions from the string region to the breaking region, are independent of the system size. For m = 10.0 the mass is large enough that one does not leave the linear scaling region, even if we create the longest string that fits in our system. Hence, we do not expect string breaking to occur in this system. The difficulty of numerically detecting string breaking with Monte Carlo techniques has been traced back to the mixing of two-meson states being hard to capture by Wilson loops [5,6,54]. With our method we can explicitly see this mixing happening, as illustrated by figure 2. The figure shows the energy as a function of imaginary time for several cases in the breaking region. At the beginning of the imaginary time evolution, we observe the energy going down until it reaches a metastable plateau, corresponding to the interacting string state, evidenced by the spin and flux configurations at that point, shown in the left inset panels. For later times there is then again a significant decrease in energy corresponding to the breaking of the string. This can also be seen in the spin and flux configuration (right inset panels of figure 2), where the region of high flux between the static charges is going away after the decrease in energy, and only a peak in the flux around the external charges remains. To better compare the string breaking scenario with that of a string ground state, the whole imaginary time evolution of spin, flux and charge square configurations for masses m = 3.0 and 10.0 and system sizes N = 22 and 30 is shown in figure 3. The flux configuration is giving an indication that the particle-antiparticle pairs created during the breaking process are clustering around the heavy external charges and screen the electric flux. Furthermore figure 3 reveals that there is essentially no difference between both system sizes.
In figure 4 we plot the charge squared Q α n 2 and the histograms P l for the initial configuration and for the final ground state respectively. We subtract the configuration of the interacting vacuum in the sector without external charges, in order to better visualize the difference. As one can see, in the breaking case two peaks around the external charges form in the charge square configuration, thus verifying that the particles created during the breaking process indeed cluster around the external charges and screen the flux. By contrast, in the nonbreaking case the charge square configuration only changes slightly. The histograms for the string lengths show a similar picture. In the nonbreaking case the clear initial peak at l = 11 is preserved, whereas for m = 3.0 it vanishes and peaks emerge around smaller string length. From figure 1 we can identify the parameter regions in which we expect string breaking to occur. To study the real-time dynamics of the string breaking we select two distinct situations. We choose l = 11, which is deeply in the breaking region for m = 3.0 but still far enough from the point where finite size effects are noticeable, and for which we do not expect any string breaking with m = 10.0. For all the following real-time cases we set ε = 3.0 and g = 1.0, as in the imaginary-time setup.

Real-time evolution with static external charges
MPS methods give us access not only to the static properties, but also to the real-time dynamics of the system. We may then investigate how the string breaking process manifests dynamically when the static charges are introduced in the interacting vacuum. To this end, we first compute a MPS approximation to the interacting vacuum of the theory using variational energy minimization [51,55,56]. Starting from this state, we apply again the non gauge invariant operator U n . . . U n+l−1 , which effectively creates two static external charges q y = ±1/2 separated by a distance l = 11, in a gauge invariant manner, thus starting the connecting flux tube. Subsequently we evolve this state in real time with ∆t = 1.0 × 10 −4 , D = 100 and determine the time dependent spin, flux and charge square configurations along the chain, as shown in figure 5. As in the imaginary time case, one can see that a system size N = 22 is sufficient to avoid finite size effects. In order to better appreciate the dynamics, we show the details of the three proposed observables at fixed   antiparticle pairs are created inside the string region, leading to a significant increase in Q α n 2 there. This is accompanied with a considerable decrease of the initial peak in the string length histogram around l = 11 and a change in the spin and flux configurations. This growth is continuing up to t = 0.5 where there has been a large amount of particles created in the string region and the peak in the histogram at l = 11 is already gone. At later times, t = 2.0, these particles are clustering around the region of the external charges and, as the flux configuration reveals, the external charges are screened, leading to a reduction of flux in the center region of the original string. For m = 10.0 the picture is significantly different, as there is essentially no change during the evolution in either the charge, spin and flux configurations or the string length histograms, which show a single dominant peak at l = 11 at all times, therefore indicating that the initial string is preserved during the evolution. The minor changes over time present in figure 7 result from the fact that the starting state of the evolution is not an eigenstate of the Hamiltonian, and consequently it is not perfectly steady.

Real-time evolution with dynamical charges
The scenarios with external static charges allow us to isolate and study the string breaking phenomenon. However, it is also worth investigating the alternative scenario in which the charges added to the vacuum are themselves dynamical, as this is arguably a more realistic situation. This enables additional configurations, in which the charges can move, to also play a role in the evolution, so that the string breaking phenomenon may be displayed differently.
Using the same MPS techniques, we can also explore this setup. Thus we repeat the simulations described in section 4 but applying the gauge invariant string operator from eq.  (2.2) on the interacting vacuum to construct our initial state. This again results in a state with a string between two charges. Different to the systems studied earlier, the charges are not external, but they are now created on a site and are fully dynamical. Again, we study the time evolution of the spin, flux and charge square configurations for different fermion masses and system sizes, as shown in figure 8, where we use ∆t = 1.0 × 10 −4 and D = 100 as in the previous section. We compare the case m = 10.0, in which the string does not break (see also figure 11), with m = 3.0. In the latter case, which exhibited clear breaking for static external charges, the fully dynamical situation shows differences (compare figures 6 and 10). To better identify the features of string breaking, we look at an even smaller fermion mass, m = 1.0, shown in figure 9.
In all cases we see a clear initial peak in the charge square configuration at the beginning and at the end of the string. For m = 1.0 and 3.0 new charges emerge especially in the string region and these peaks are quickly decaying. Also the flux in the string region is decaying while two small peaks are preserved roughly around the start and end point of the original string. The histograms for the distribution of string lengths show a similar  picture: at the beginning of the evolution there is a clear peak at around l = 11 which is gone at around t = 0.25. In the m = 1.0 case the initial peak at l = 11 is less dominant as in the m = 3.0 case, due to the fact that in this case the interaction strength ε is more important, leading to a state less close to a strong coupling string. For later times one sees that in both cases smaller string lengths are dominating in the system, therefore indicating that the string is broken. By contrast, in the case of m = 10.0 the dominant peak at l = 11 is preserved during the entire evolution. Although the magnitude decreases, one can see from the flux and charge square configurations that there is not a lot of change, which indicates that the original string is still present. The slight changes have the same origin as in the previous cases. As we are not starting with an eigenstate of the Hamiltonian but rather do a local quench, the state is not perfectly steady. Furthermore in this case the charges are fully dynamical, what leads to richer dynamics.

Conclusion
We have studied the real-time dynamics of string breaking with MPS in a truncated 1+1 dimensional SU(2) lattice gauge theory. By looking at the charge square, the flux config- uration and the statistics of string lengths, we are able to detect and characterize string breaking in statical and dynamical scenarios.
We calculate the ground state in the sector with two static external charges. This allows us to clearly pick up the static signature of string breaking, and to identify parameter regions where we expect string breaking to occur. Furthermore our calculation shows that the string state is metastable during the imaginary time evolution, until configurations with a broken string mix in and reduce the energy, as expected. We have also shown that the particle-antiparticle pairs resulting from the breaking string cluster around the external charges and subsequently screen the electric field. This can be detected with the spin and flux configuration as well as in the charge square and the distribution of string lengths in the system.
Using real-time simulations with MPS, we have studied the dynamics of string breaking in a setup with two static external charges. With the observables proposed before, the breaking and non-breaking cases are clearly distinguishable. We can explicitly observe that, in case the string is breaking, dynamical fermions are created that cluster around the external charges and screen them, therefore reducing the flux in the system.
Finally, we have also simulated the time evolution of a string between fully dynamical fermions. We also identify situations in which the string breaks in this case. In particular, the decay of charges and the reduction of the flux in the middle of the string region indicate string breaking. Due to the fact that the initially created particle-antiparticle pair is now dynamical, we do not observe the clustering of the dynamically created charges, but rather a distribution along the system.
In our study we have used a 1+1 dimensional SU(2) gauge model. This model could be a suitable candidate for the implementation of a quantum simulator for a SU(2) lattice gauge theory. The observables and setups that we describe would also be viable for such a potential experimental realization. Furthermore the applied techniques are not limited to the specific model studied here. The truncation method used from Ref. [48] works for arbitrary compact (or finite) Lie groups and is not limited to the case of one spatial dimension. Therefore also other gauge groups could be studied and with more general TN such as PEPS [57] this study could also be extended to higher dimensions. While completing the manuscript we became aware of a related work on TN and string breaking [58].

B Analysis of the numerical errors
In this appendix we analyze the numerical errors of the results presented in sections 3 -5. In the main text we have already presented data for different system sizes thus showing that N = 22 is enough to avoid noticeable finite size effects for the considered range of parameters. Here we focus on the influence of the time-step size and the bond dimension on our results.

B.1 Ground state with static external charges
In figure 12 we show how the spin and flux configuration changes for a system of size N = 22, if a larger bond dimension or a smaller time step is used. As the figure shows, there is no significant improvement with larger bond dimension. Reducing the time step by a factor of two to ∆t = 0.5 × 10 −3 leads to a slightly changed spin and flux configuration on the order of 10 −2 for m = 3.0. However, these changes are predominantly at an early stage of the evolution and difference in the final configuration is much smaller. For m = 10.0 the picture is qualitatively similar, at a very early stage there are small differences which are roughly one order of magnitude less than in the m = 3.0 case. All in all the changes for reduced time step and larger bond dimension are rather small compared to the data for ∆t = 1.0 × 10 −3 , D = 100 presented in figure 3 hence justifying our choice of time step and bond dimension.

B.2 Real-time evolution with static external charges
Also for the real-time evolution case with static external charges, we analyze our error analogous to the imaginary time case. The difference in results for a smaller time step ∆t = 0.5 × 10 −4 and larger bond dimension D = 130 is presented in figure 13. A reduction of the time step shows a rather similar effect to the imaginary time case and leads to differences on the order of 10 −3 for m = 3.0 and 10 −4 for m = 10.0. For m = 3.0, the differences comparing to results with enlarged bond dimension are bigger than in the imaginary time case, whereas for m = 10.0 there is essentially no change up to machine accuracy. Overall, we see the same picture as in the imaginary time cases that our choice of time step and bond dimension controls the error well enough to avoid considerable influences on the effects observed in figure 5.

B.3 Real-time evolution with dynamical charges
For the real-time evolution with dynamical charges the same error estimation as in the cases with static charges yields the results in figure 14. A reduction of the time step from ∆t = 1.0 × 10 −4 to 0.5 × 10 −4 yields a change in the spin and flux configuration on the order of 10 −2 . Contrary to the cases with static charges, the change is more pronounced for a large fermion mass m = 10.0, whereas for smaller fermion masses there is less change. For the bond dimension however, we see again the same effect that the changes are more pronounced if the mass is smaller, whereas there is almost no difference for larger mass. Nevertheless, also in this case the absolute changes in the spin and flux configuration are