Real Time Dynamics and Confinement in the $\mathbb{Z}_{n}$ Schwinger-Weyl lattice model for 1+1 QED

We study the out-of-equilibrium properties of $1+1$ dimensional Quantum Electrodynamics (QED), discretized via the staggered-fermion Schwinger model with an Abelian $\mathbb{Z}_{n}$ gauge group. We look at two relevant phenomena: first, we analyze the stability of the Dirac vacuum with respect to particle/antiparticle pair production, both spontaneous and induced by an external electric field; then, we examine the string breaking mechanism. We observe a strong effect of confinement, which acts by suppressing both spontaneous pair production and string breaking into quark/antiquark pairs, indicating that the system dynamics deviates from the expected behaviur toward thermalization.


I. INTRODUCTION
Gauge fields coupled to matter are at the heart of the microscopical description of nature, representing the key ingredient to describe the physics of fundamental particles and interactions and appearing in many problems of condensed matter physics, such as superconductivity or the Quantum Hall Effect [1,2]. Quantum field theory has developed by tackling more and more complicate issues, and by now we have not only a comprehensive description of the low energy sector of the model, but also techniques to study strong coupling, non-perturbative and topological effects. Numerical simulations, especially the ones based on Monte Carlo techniques [3], have been able to confirm many phenomena predicted with analytical methods and shed light on several non-trivial aspects. The possibility to describe a given field theory through a lattice theory is an important problem [4,5], which has been investigated by outstanding physicists since the 70' [6][7][8][9].
It is clear that, from a theoretical point of view, it is important to determine how effectively a lattice model is able to reproduce the phenomenology of a given quantum gauge theory. Very often this offers a formidable computational problem, but recently techniques developed from quantum information insights, such as tensor networks [10,11], have been used to tackle these questions, thanks to their intrinsic ability to restrict the dynamics to relevant subspaces of the Hilbert space, exploiting the entanglement of the states that contribute to the dynamics. From an experimental point of view, recent advances in the control of ultra-cold atomic systems have paved the way to a whole new way of studying fundamental theories in the spirit of quantum simulations, as suggested by Feynman [12]. At present, there are many proposals in the literature to use ultra-cold atoms or trapped ions in optical lattices to simulate many body Hamiltonians [13][14][15][16][17], as well as Abelian and non-Abelian lattice gauge theories [18][19][20][21][22][23][24][25][26]. These ideas led to the first experimental realization of a quantum simulation of 1+1 QED with four trapped ions [27]. On the other hand, the implementation of classical "synthetic" fields has been proposed and realized in a number of experiments [28][29][30][31][32].
In general, an approach based on quantum simulators opens the possibility to study in a more efficient way problems that have always been difficult to tackle, including non-perturbative effects as well as genuine dynamical behavior. These phenomena are necessary to understand the essence of Abelian and non-Abelian quantum field theories in applications to particle physics. Two examples on which we will elaborate in this article are vacuum stability in QED with respect to particle/antiparticle pair production [33] and confinement in QCD [34], together with the associated string breaking mechanism. In recent years, much attention has also been devoted to the interplay between correlations and dynamics in statistical mechanics and many-body systems in the quantum realm, which may exhibit unusual thermalization and equilibration properties. It has been stressed that, in a closed quantum system, non-equilibrium dynamics might display behaviors that cannot be understood by means of the standard theory of ergodic or integrable models and the corresponding (Gibbs ensemble) thermalization hypothesis. This is the case of the so-called many-body localization phenomenon (for reviews of this topic see, for example, [35,36]) or of the recently discussed quantum scars [37,38]. The interest in these systems and phenomena has fostered the development of novel paradigms.
In the spirit of such research developments, in this article we tackle the problem of the out-of-equilibrium real-time dynamics of (1 + 1)-dimensional QED, which represents the simplest model of fermionic matter coupled to (Abelian) gauge fields. In spite of its simplicity, the model shares many phenomenological features with arXiv:1909.04821v2 [quant-ph] 13 Sep 2019 more complicated theories, such as QCD, thus representing an interesting instance to test quantum field theory and quantum simulations. The model has been widely investigated since the 70's by a number of outstanding physicists [6][7][8][9] and has been recently proposed as a possible viable experimental option in the not-too-distant future [39][40][41][42][43]. Indeed, a first experiment reproducing such a model with few qubits in a ion trap has already been reported [27] and new ideas have been proposed for possible experimental realizations with Rydberg atoms [39,44].
The first problem to face when trying to encode a lattice-gauge theory in a quantum simulation with cold atoms or a classical numerical simulation is the fact that the number of states of the gauge field, associated to links between lattice sites, must be finite. This problem has been tackled in different ways in literature. Sticking to U (1) theories, it is possible to represent gauge fields by means of spin variables [45,46], considering the so called quantum link model [47][48][49][50]. Otherwise, one can try to discretize the gauge group by replacing U (1) with Z, which however admits only infinitedimensional representations. In such a case, the Hilbert space describing the gauge degrees of freedom must be truncated [42,51]. These approaches allows for the investigation of questions that are traditionally difficult to analyze, such as the emergence of quantum phase transitions, non-perturbative phenomena and dynamical aspects [22,45,46,[51][52][53][54][55].
We adopt here the scheme presented in [56], where we considered a lattice version of the staggered-fermion Schwinger model [57], in which the U (1) gauge degrees of freedom are given in terms of the possible finite dimensional representations of the corresponding Weyl group [40], thus yielding a discrete and finite implementation through the group Z n . In this approach, a systematic control of finite-n effects is possible and a rigorous approximation of lattice U (1) quantum electrodynamics in the large-n limit becomes available to investigation. The static phase diagram of this model has been studied in [56], while in [58,59] a variant of it has been shown to admit topological phases. In this paper, our goal is to study out-of-equilibrum properties of this model, by analysing two major phenomena: i) the stability of the Dirac vacuum with respect to the production of virtual particle/anti-particle pairs, induced by quantum fluctuations and relevant to understand striking behaviors of the quantum vacuum such as the Casimir effect [60]; ii) dynamical effects of confinement, such as the string breaking mechanism, which is strictly connected to the question of the asymptotic freedom of quarks [34].
We will see that both these phenomena depend on the values of the two parameters that enter the Hamiltonian, namely the fermionic mass m and the gauge coupling g. In particular, we will show that, in the strong coupling regime, the dynamical behaviour of the model strongly deviates from the usual thermalization and relaxation properties which are expected to be found in a many-body non-integrable system, resulting in stable and/or recurrent evolution of interesting physical quantities. We remark that similar results have been obtained in different many-body systems, such as non-integrable spin-chain models [61] or constrained Hamiltonians [38] .
This article is organized as follows. In Sect. II, we review the discretized version of the Schwinger model for 1+1 dimensional QED, that is accommodated on a onedimensional lattice and endowed with a Z n symmetry. This model, which depends on two physical parameters, namely the fermionic mass m and the gauge coupling g, exhibits a quantum phase transition at m c = −0.33, belonging the Ising universality class. For m > m c the ground state of the model is in a confined phase, in which elementary excitations above the Dirac sea vacuum are of mesonic type. In Sect. III, we set up the quench protocol to simulate spontaneous pair production, occurring in absence of an external electric field. We analyze this phenomenon by looking at the dynamical evolution of several physical quantities of interest, such as particle density, entanglement entropy, density correlation functions. We will see that, contrary to what is found in many other integrable and non-integrable models, the production of correlated particle/antiparticle pairs is strongly suppressed when we consider system parameters deep in the confined phase. In Sect. IV, we examine the phenomenon of pair production induced by the presence of an external electric field, finding an agreement between our simulations and old predictions by Schwinger [33] for the rate of pair production. Also, by examining the time evolution of entanglement, we conclude that the formation of mesonic excitations is stimulated by a dynamical effect due to the presence of the external field, of a different nature from the one emerging in the spontaneous case, examined in the Sect. III. Finally, in Sect. V, we will consider the real time evolution of a string excitation. We observe that the string breaks into mesons, thus giving rise to the so-called string-breaking mechanism, only when interactions are sufficiently weak. Vice versa, deep in the confined regime, the strings remain localized and are apparently stable.

II. THE MODEL
We start from the Hamiltonian of quantum electrodynamics (QED) in one spatial dimension: which describes a U (1) gauge theory, where a charged particle (electron), represented by the spinor field ψ(t, x), interacts with the electromagnetic field F µν = ∂ µ A ν − ∂ ν A µ (µ, ν = 0, 1), associated with the potential A µ .
Here m and g represent the electron mass and charge, respectively, while the Dirac matrices γ satisfy {γ µ , γ ν } = 2η µν with η = diag(1, −1). We work in the canonical gauge, in which the temporal component A 0 is set to zero, while the spatial component A := A 1 is taken as the conjugate variable to E: The Hamiltonian (1) must be complemented with the enforcement of Gauss law in a weak sense, meaning that it is possible to select a subspace, called the physical subspace, of states |ψ for which G(x)|ψ = 0.
As described in Ref. [40], to which we refer for further details, it is possible to discretize this model following two steps: i) first, we perform a spatial discretization, in which the space continuum is replaced by a linear lattice of points with spacing a; in order to avoid the fermion doubling problem we adopt Schwinger's approach [4][5][6][7] of staggered fermions; ii) second, by following the Weyl-Wigner quantization scheme, we approximate the gauge group U (1) with the finite group Z n ; this step is essential in order to work with a finite number of local degrees of freedom also for the gauge variables.
The discretized Hamiltonian reads: with x (1 ≤ x ≤ N ) labelling the sites of a onedimensional lattice of spacing a = 1. Here, the onecomponent spinor is represented by the creation and annihilation operators ψ † x and ψ x , defined on each site x and characterized by a staggered mass (−1) x m , while the gauge fields are defined on the lattice links (x, x + 1) through the pair of operators E x,x+1 (electric field) and U x,x+1 (η) (unitary comparator). The gauge field operators act on the n-dimensional Hilbert spaces H (x,x+1) , attached to each link and spanned by the orthonormal bases {|v k x,x+1 } 0≤k≤n−1 , that diagonalize the local electric field: Here, a non-zero value of the angle φ entails the presence of a constant background field, which in turn corresponds x,x+1 living on links between adjacent sites, described by a periodic discrete variable with a step 2π/n. to charges at the boundary of the chain. The unitary comparator, instead, acts as a cyclic ladder operator: Gauss law is implemented by requiring that the physical states belong to the null space of the operator Let us remark that, in one spatial dimension, the fermionic density and the local electric field completely determine each other up to a constant, which corresponds to the value of the electric field on one boundary. Thus, one can integrate out the gauge field in order to obtain an effective Hamiltonian in which only the matter fields appear, corresponding to a spin 1/2 model with longrange interactions [44]. On the other hand, one can eliminate the fermionic field to get an effective Hamiltonian that contains only the gauge variables. Within our approach, one thus obtains a local Hamiltonian with a Z nsymmetry, acting as (4) when restricted to the physical subspace. The study of such Hamiltonian will be presented in a future work. Here, it is worth noticing that these ideas can be generalized also to higher dimensions [62,63].
In Ref. [56], we analyzed the phase diagram of the discretized Hamiltonian (4) showing that, for any n, it exhibits a phase transition at a critical value of the mass m c between a confined phase for m > m c , in which the ground state is given a dressed Dirac sea vacuum, and a deconfined phase for m < m c , in which the ground state is characterized by a nonvanising average electric field, which entails a symmetry breaking.
In our representation, the Dirac sea vacuum is obtained by filling up all odd sites (negative mass fermions) and leaving the even ones (positive mass fermions) empty, as shown in the first line of this case, Gauss law requires that the electric field on the connecting link is equal to + 2π/n (− 2π/n), represented in the picture by an oriented link pointing to the right (left). In the fourth line of Fig. 2 we also show another gauge invariant configuration, representing a string, in which the particle/antiparticle excitations that constitute a meson are moved farther away from each other, with the electric field between constant and nonvanishing on all the intermediate links. The real-time dynamics of such string configurations will be studied in the second part of this article.
The phase transition belongs to the Ising universality class for all n, with the confined/deconfined cases corresponding respectively to the paramagnetic/ferromagnetic phases of the Ising model in a transverse field. Here, the quantum phase transition is driven by the value of the fermion mass (playing the role of the external transverse magnetic field) and the order parameter (the magnetization in the Ising case) is represented by the mean value of the electric field operator: or, equivalently, by the mean fermion density which, in the thermodynamic limit, are nonvanishing only in the deconfined phase. The value of m c depends on the dimension n used to discretized the U (1) gauge group, but it approaches m c −0.33 in the large-n limit.
In addition, as expected in the Ising model, close to the transition point the first excitation has conformal dimension d = 2, corresponding to a kink-like (or domain wall) solution.
Let us remark that while the Ising model in a transverse field is integrable for any value of the magnetic field, our model is more complicated and never integrable (except for the trivial case g = 0), because of the gauge coupling between fermionic matter and electric field. In the Ising case, integrability is lost and effective interactions between domain-wall excitations are present only if one adds the coupling with an external uniform longitudinal field [64]. This effect can be mimicked in our model by introducing a background constant electric field.
The effects of such gauge-mediated interactions might be quite strong and will be studied by looking at realtime dynamical properties of our model in the next sections. The analysis will be performed by numerically studying the Hamiltonian (4) with a t-DMRG algorithm, whose dynamical evolution is implemented through the Runge-Kutta method. Details about the precision of our algorithm are given in the Appendix.

III. SPONTANEOUS PAIR PRODUCTION
In this section, we will examine the phenomenon of spontaneous pair production in the Schwinger model for (1 + 1)-dimensional QED, by simulating with our model the real-time dynamics of the Dirac sea vacuum, in absence of an external electric field. This effect has also been considered in other approaches [51,[53][54][55] and experimentally analyzed in a small system (4 qubits) of trapped ions [27].
In order to test the stability of the Dirac vacuum with respect to spontaneous pair production, we prepare the system in the m → ∞ ground state corresponding to the the Dirac sea vacuum, and study its evolution under the action of the Hamiltonian with different values of the coupling constants m, g, with either m > m c or m < m c . For completeness, we will perform our analysis for all values of m, both positive and negative, but we stress that the Dirac vacuum we take as initial state is highly excited for m < m c , while it is very close to the true ground state in the confined phase and in particular when m is large and positive. In the language of spin systems, this would be analogous to preparing an Ising system in the paramagnetic ground state, by setting the external transverse field h = 0, and then suddenly quenching the Hamiltonian to a different value of h, with h < 1 (ferromagnetic phase) or with h > 1 (quenching into the paramagnetic phase) [61]. Let us consider the Z 3 case. To test the evolution of the Dirac sea vacuum after the quench, we numerically calculate three different quantities.
1) The first quantity we consider is the mean density of particles where, since = 1, time is measured in units of [energy] −1 . The quantity (12) vanishes on the Dirac sea vacuum and takes the value 1 on the state with the maximal number of mesons.
To perform a first test on our model, we calculate ρ(t) for a chain of N = 4 sites, after a quench to m = 0.5 and g = 6/2π, which correspond to the parameters used in the eperimental set up of [27]. Figure 3 shows our result: the density starts form zero to go very close to the value 1/2, corresponding to the formation of one meson. After reaching a maximum, recombination effects bring the value down again, in perfect agreement with the observations of [27].
More generally, we set g = 6/2π and quench to different values of m. The temporal evolution of ρ is shown in Fig. 4  We clearly see that for very large values of the mass (both positive and negative), ρ(t) oscillates periodically between zero and a rather small value, due to a small pair production rate and to strong recombination effects: thus, the Dirac sea vacuum is essentially stable. On the contrary, for smaller values of |m|, the density increases rapidly up to values close to 1 (corresponding to the mesonic ground state), to start then oscillating because of recombination effects, now around a large value, showing that the Dirac sea vacuum is unstable. This is particularly evident for quenches to values of the mass close to the critical one, m c ∼ −0.33.
Let us analyze better the curves of Fig. 4 as a function of m. Figure 5(a) shows the maximum value of ρ(t), extrapolated form the first peak: it is clear that the largest pair population occurs close to the phase transition. In Fig. 5(b) we show instead the period T of the oscillations,  for m > m c . We see that the data are very well described by the continuous line, which corresponds to the best fit obtained with the functional form T (m) = 1/(am + b). This form is expected from a first-order approximation, according to which the period might be estimated from the inverse of the energy difference between the Dirac sea energy (E Dirac /N = −m) and the energy of the mesonic state (E meson /N = m + g 2 /2).
In the Appendix, we give details on the finite size scaling analysis we performed.
To summarize our results, we show in Fig. 6 the contour plot of the density ρ ∞ (t) extrapolated for N → ∞, in the whole range of the quenched mass m ∈ [−5, +5]. We can conclude that the Dirac sea vacuum is unstable for values of the mass close to the critical value m c , whereas pair production is 2) The second quantity we calculated is the time evolution of the half-chain entanglement entropy where the chain has been partitioned in the two subsytems A and B, consisting of the first and last N/2 sites of the chain, respectively. The initial Dirac sea vacuum is separable, hence S N/2 (t = 0) = 0. 3/π. We see that the entanglement entropy shows an oscillatory behaviour, reaching a maximum at small values (e.g. S max N/2 ∼ 0.5 for m = 2.0 or m = −3.0). The period of these oscillations increases as m decreases and, at the same time, the entropy reaches higher values, of the order of unity, for small m, as one can observe in Fig. 7(c). For m close to the critical value m c (see for example the data for m = 0.0, −0.5), the entropy increases much faster and monotonically. Thus, we can conclude that when the phenomenon of pair production is dominant, quantum correlations between the two parts of the chain are larger due to the fact that the spontaneously created entangled particle/antiparticle pairs spread along the chain. This effect has been extensively studied in free fermionic [65] and/or integrable models [? ], where entanglement can be shown to increase linearly with time, at least as long as one can assume a maximum speed of propagation (the Lieb-Robinson bound [73]). Not considering disorder and/or long-range interactions, this behavior has also been confirmed in several integrable and non-integrable models [66][67][68][69][70][71][72]. In [61], the time evolution of the half-chain entanglement entropy has been studied for the Ising model with both a transverse and a longitudinal field, which admits a confined phase [64], showing that the growth of the entanglement entropy is strongly reduced for quenches within the confined (ferromagnetic) phase. This corresponds to what we observe in our model [74]. The black dashed line in Fig. 7(c) shows the linear behavior of the entanglement entropy in the free case. For small but non-zero values of m, we can clearly recognize a slow-down of the entanglement entropy growth, which is now well described by a logarithm: S N/2 (t) ∼ log t. When the interaction is strong and we are deep in the confined phase, the system does not seem to evolve toward an equilibrium configuration, at least on the time scale of our simulations. Indeed, entropy is strongly suppressed and shows revivals, similar to what happens in models where the quantum scar phenomenon appears [37,38,75].
3) As a last indicator, we consider the time evolution of connected correlation functions (15) shown in Fig. 8. For small values of m, these functions exhibit a light-cone spreading, typical of conformal or integrable theories, as predicted in [76]. However, as we enter in the intermediate and strong-coupling region, we observe a localization effect and an oscillatory behaviour, indicating that particle/antiparticles pairs do not spread, but are cyclically created and recombined. Similar behaviors have been observed in other models [44,61].

IV. PAIR PRODUCTION IN AN EXTERNAL FIELD
In the context of particle physics, pair production is often studied in presence of an external electric field, from which the pair production rate is expected to depend. For the Schwinger model, the critical value E c of the external electric field at which the e + e − production rate should reach a maximum is given by E c = m 2 /g (which, in SI units, reads E c ≈ 1.32 × 10 18 V/m [33]), but this effect has never been observed experimentally since the aforementioned critical value is still out of the range of even the most powerful lasers. In 1 + 1 dimensions, a formula for the production rate has been proposed in [77,78]: whereρ represents the time derivative of the total density of the chain in the infinite length limit, while = E 0 /E c , E 0 being the value of the external field. This formula predicts that the production rate increases linearly for large values of and it is exponentially suppressed for 1. This formula can be checked in our simulations. We start by considering the Z 3 -model. We choose values of (m, g) for which Dirac sea vacuum results stable according to the analysis described in Sec. III. At t = 0, we apply a constant uniform electric field E 0 along the whole chain and run the simulation to obtain ρ(t) for the corresponding value of = E 0 /E c , with E c = m 2 /g. Fig. 9(a) shows the results of our simulations for the case m = 4.5 and g = 3/π. We observe that the vacuum is stable, with some small oscillations, for very small values of . As we increase , we start observing a linear regime, followed by a saturation effect, as predicted by Eq. (16). Figure 9(b) shows the range in which we performed the linear fit in order to evaluateρ and verify Eq. (16).
Since we expect relevant size effects, we repeated the simulations for different values of N , namely N = 50, 60, 70, 80, 90. Also, to obtain the large-n limit and check if we can reasonably approximate the U (1) limit with our model, we have preformed an analogous analysis for the Z 5 and the Z 7 models. Our results, which are summarized in the Appendix, show that the approximation of the continuum limit improves with n. Finally, in order to get some additional insight on the phenomenon, we have calculated the time evolution of the half-chain entropy, for different values of the external magnetic field. A typical behavior is shown in Fig.  10, which shows the data for m = 4.5 and g = 3/π. We notice that the entropy always presents an oscillatory behaviour with not too large maximum values. More interestingly, we observe that an increase of the external field, which results in a rapid increase in particle pair production as shown in Fig. 9(a), does not contribute at all to the increase and spreading of entanglement. Thus, we can postulate that the mechanism that dictates the particle/antiparticle pair production is different in presence and in absence of the external field: in the former case, pair formation is determined by the system charges rearranging as an effect of the external field, while, in the latter case, entangled quark/antiquark pairs are spontaneously created out of the Dirac sea.

V. THE STRING BREAKING MECHANISM
As recalled in the introduction, the string breaking mechanism is a very interesting phenomenon that is expected to occur in theories, such as (3 + 1)-QCD, admitting confined phases, but it is also very hard to prove either analytically or numerically since its effects are mainly dynamical.
In our approach, we are able to investigate this phenomenon, since we can simulate the real-time dynamics of a generic initial state, which can be let to evolve with the discretized Hamiltonian (4). In particular, we prepare our system in the string-excitation state shown in Fig. 2(c), where a particle and an antiparticle separated by a distance l are created from the Dirac sea vacuum, giving rise to a non-zero string of electric field in between.
In the following, we will consider the Z 3 model and a chain of length N = 80. We initialize our system in a state with a string placed at the center of the chain: the particle and the antiparticle are put at a distance equal to 20 lattice sites, so that the electric field is different from zero (and equal to + 2π/3) only on the 19 median links of chain. The evolution of this state is followed by looking at the value of the electric field E x,x+1 (t) on each link. Clearly, this process is also affected by the instability of the Dirac sea vacuum that we studied in the previous sections. We report the analysis of the real time dynamics of such a string for three different values of (m, g), specifically corresponding to (a) weak interactions (m = 0.1, g = 0.1), (b) intermediate interactions (m = 0.3, g = 0.8), (c) strong interactions (m = 3.0, g = 1.4). Clearly, the total evolution is also affected by the instability of the Dirac sea vacuum that we examined in the previous sections. Thus, to improve the interpretation of our numerical results, we subtract from E x,x+1 (t) the value of the electric field that would be obtained starting from the Dirac sea vacuum for the same values of (m, g). The corrected data are shown in Fig. 11, which clearly show three different situations: a) for weak interactions, the string starts to spread and breaks into particle/antiparticle pairs (mesons) that, after a short time in which one can notice a rich process of pair production and recombination, stabilize in a configuration with two mesons localized at the edges of the string; the two mesons are deconfined, since they move away one from each other at a constant speed; b) for intermediate interactions, the string does not spread; still it breaks into particle/antiparticle pairs (mesons) that, as in the previous case, quickly stabilize in a configuration where we can still distinguish two mesons localized at the edges of the string; however, now the two mesons are confined, since they are kept at a fixed distance from each other; c) for strong interactions, the string is completely stable: it does not spread and it does not break into mesons.
The simulations which are reported in Fig. 11 can be repeated for any couple of values (m, g). Our results are summarized in Fig. 12, where we show the contour plot of the large-t (namely t = 4.0 in our units) total value of the electric field at the center of the chain, defined as the sum of the electric field on the 12 central links, as a function of the coupling constants (m, g). The two (white) level curves correspond to 10% (dotted line) and 50% (solid line) of the initial value, respectively. From  this picture, the three regimes described above are clearly identified: (a) the lighter and central part of the diagram corresponds to the breaking of the string into two deconfined mesons; (b) the reddish part of the diagram to the breaking of the string into two confined mesons; (c) the darker part to a stable string configuration.
It is interesting to examine the string breaking phenomenon by looking also at the time evolution of the half-chain entanglement entropy, for different values of (m, g). Figure 13 shows the behaviour of S N/2 evaluated on the time evolution of the string state by first keeping g = 0.1 fixed and let m vary (a) within positive values, (b) within negative values and then (c) keeping m = 0.1 fixed and letting g vary. These graphs confirm what found by looking at the real time dynamics of the electric field configuration, showing a growth of entanglement entropy that is linear for weak interactions (small m and/or small g), sublinear for intermediate interactions and suppressed (with for small oscillations) in the strong coupling regime.

VI. CONCLUSIONS
We have investigated out-of-equilibrium properties of (1 + 1)-dimensional QED, approximated via a Z n Schwinger model. By means of simulations of the stability of the Dirac vacuum with respect to particle/antiparticle pair production and of the string breaking mechanism, we have studied the effect of confinement on the real time dynamics of the model.
We have found that confinement has a relevant effect on the dynamical properties of the model, resulting in oscillatory behaviours and lack of thermalization for pair production as well as in a total suppression of the breaking and spreading of string excitations, with a perfect localization of the latter.
Let us notice that such a reduction of entanglement and slow-down of dynamics have been observed in other systems. This is the case of not only of model with longrange interactions [79] but also of the Ising and Potts models with both a transverse and a longitudinal magnetic fields [61,80]. A similar behaviour has also be seen in constrained models which exhibit quantum scar states [37,38] or for spin-1/2 chain Hamiltonians that can be derived form abelian gauge models in the quantum link approach [44]. It is interesting to notice that, in the two latter cases, the physical states of the system under consideration are constrained to lie in a restricted subspace of the total Hilbert space, a fact that it is shared by our model where the role of Gauss law constraint is crucial. This is also at the heart of the recent idea to experimentally implement these Hamiltonian with Rydberg atom systems in the Rydberg-blockade regime [39]. We plan to further investigate the role of the gauge constraint in future works.
sion [81]. The time-evolution is based on a Runge-Kutta 4th order scheme, with a time step of δ = 0.01. We implemented the initial state for the time-evolution by calculating the ground-state of two different Hamiltonians: i) for the pair production analysis, we obtained the vacuum state (see Fig. 2) by using the Hamiltonian (4), setting t = 0 and large values of m and g, i.e. m = 5, g = 3; ii) for the string breaking mechanism, we added to the aforementioned Hamiltonian localized electric field terms E x,x+1 on each link in which we want to create the initial string. In our simulations, we used a variable number of DMRG-states, up to 1200, in order to keep the truncation error below 10 −6 at each time step.
Spontaneous pair production. To study finite size effects, we have repeated our simulations for different chain sizes (N = 16,20,24,28,32,36,40), in order to extract the infinite size limit of ρ(t). For example, Fig. 14(a) shows the behaviour of ρ(t) close to its first maximum, for m = 2.0 and different chain sizes. We find that, at all instants of time t 0 , we can estrapolate its infinite size limit ρ ∞ (t 0 ) according to the following fit: An example, for m = 2.0 and t 0 = 0.52 is shown in Fig.  14(b), for which we get ρ ∞ (t 0 ) = 0.2175 ± 0.0001 and β(t 0 ) = 0.1703 ± 0.0001. The time evolution of ρ ∞ (t) is given in Fig. 14(c). Also, we have repeated the same simulations in the case of Z 5 and Z 7 . In Fig. 15(a) and Fig. 15(b) we show the contour of ρ ∞ (t) in the whole range of the quenched mass m ∈ [−5, +5], similarly to what we have done in Fig. 6 for the Z 3 -model. We see a very similar behaviour.
Pair production in an external field. To evaluate finitesize effects in the pair production rate, we have repeated the simulations for N = 50, 60, 70, 80, 90 while, to check the large-n limit to see if we can reasonably approximate the U (1) limit, we have performed simulations also for the Z 5 and the Z 7 models. The results are shown and compared with Schwinger formula in Fig. 16(a) and Fig.  16(b). Even at the largest system sizes that we are able to investigate, it is evident that our numerics still suffer from strong finite-size and finite-n corrections. However, the data follow a pattern that is qualitatively in agreement with the predictions of the continuum U (1)-model and shows the correct scaling. Thus we can conclude that indeed we are describing the physics of the Schwinger model.