{'fraction_non_alphanumeric': 0.057982788710356864, 'fraction_numerical': 0.02046097370163086, 'mean_word_length': 4.180018703241895, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Real-time simulation of the Schwinger effect with Matrix Product States

Matrix Product States (MPS) are used for the simulation of the real-time dynamics induced by an electric quench on the vacuum state of the massive Schwinger model. For small quenches it is found that the obtained oscillatory behavior of local observables can be explained from the single-particle excitations of the quenched Hamiltonian. For large quenches damped oscillations are found and comparison of the late time behavior with the appropriate Gibbs states seems to give some evidence for the onset of thermalization. Finally, the MPS real-time simulations are explicitly compared with the semi-classical approach and, as expected, agreement is found in the limit of large quenches.

Gauge theories lie at the heart of high energy physics and hence play an essential role in our understanding of nature. Morever, gauge theories also emerge as low energy effective theories in several condensed matter systems [1]. Lattice gauge theories provide a nonperturbative regularization of such theories that can often be simulated very efficiently by using Quantum Monte Carlo (QMC) methods. However, several most pressing questions in that regard, e.g., the phase diagram of quantum chromodynamics (QCD) at finite chemical potential or the real-time dynamics of relativistic heavy ion collisions have largely remained out of reach [2].
Over the last decade, the Tensor Network States (TNS) approach has become a powerful alternative method to study strongly correlated quantum systems since it does not suffer from the sign problem [3][4][5]. The most famous example of TNS are the Matrix Product States (MPS) [6] in one spatial dimension. Ever since the formulation of Density Matrix Renormalization Group [7] in terms of MPS, the number of algorithms for quantum many-body systems has increased rapidly. Recently, MPS have also been successfully applied to lattice gauge theories [8][9][10][11][12][13][14][15][16][17][18][19][20].
An intriguing effect in the Schwinger model concerns the non-equilibrium dynamics after a quench that is induced by the application of a uniform electric field E 0 = gα onto the ground state |Ψ 0 at time t = 0. Physically, this process corresponds to the so-called Schwinger pair creation mechanism [45] in which an external electric field separates virtual electron-positron dipoles to become real electrons and positrons. Recently, this process has attracted much interest since high-intensity laser facilities like the Extreme Light Infrastructure (ELI) will for the first time be powerful enough to probe this effect experimentally. So far, theoretical investigations have mainly been restricted to the regime in which the fermions are treated quantum mechanically whereas the gauge fields are described classically (quantum kinetic theory [46,47] or phase-space methods [48]), or classicalstatistically (real-time lattice techniques [29,49]).
In this publication we apply the MPS framework to investigate the non-equilibrium dynamics at the full quantum level. We perform real-time simulations for small, intermediate and large quenches. Furthermore, we use MPS computations of ground states, single-particle excitations and Gibbs states to analyse and interpret our results. Finally, we explicitly compare the MPS simulations with those obtained using real-time lattice techniques.

II. SETUP
A. Kogut-Susskind Hamiltonian The massive Schwinger model describes (1 + 1)dimensional QED with one fermion flavor that is described by the Lagrangian density Here, ψ is a two-component fermion field, A µ denotes the U (1) gauge field and F µν = ∂ µ A ν − ∂ ν A µ is the corresponding field strength tensor.
In the following, we employ a lattice regularizationà la Kogut-Susskind [50]. Therefore the two-component fermions are decomposed into their particle and antiparticle components which reside on a staggered lattice. These staggered fermions are converted to quantum spins 1/2 by a Jordan-Wigner transformation with the local Hilbert space basis {|s n n : s n ∈ {−1, 1}} of σ z (n) at site n. The charge −g 'electrons' reside on the odd lattice sites, where spin down (s = −1) denotes an occupied site whereas spin up (s = +1) corresponds to an unoccupied site. Conversely, the even sites are related to charge +g 'positrons' for which spin down/up corresponds to an unoccupied/occupied sites, respectively.
Moreover, we introduce the compact gauge field θ(n) = agA 1 (n), which lives on the link that connects neighboring lattice sites, and its conjugate momentum E(n), which correspond to the electric field. The commutation relation [θ(n), E(n )] = igδ n,n determines the spectrum of E(n) up to a constant: E(n)/g = L(n)+α. Here, L(n) denotes the angular operator with integer spectrum and α ∈ R corresponds to the background electric field. Accordingly, the Kogut-Susskind Hamiltonian reads [50,51] where σ ± = (1/2)(σ x ± iσ y ) are the ladder operators.
Here, we have introduced the parameter x as the inverse lattice spacing in units of g: x ≡ 1/(g 2 a 2 ). The continuum limit then corresponds to x → ∞. We note that H α is only invariant under T 2 (translations over two sites) due to the staggered mass term in the Hamiltonian. The Hamiltonian is invariant under local gauge transformations that are generated by: If we restrict ourselves to physical (i.e., gauge invariant) operators O for which [O, G(n)] = 0, the Hilbert space decomposes into dynamically disconnected superselection sectors, which are distinguished by the eigenvalues of G(n). The sector with G(n) = 0 at every site n constitutes the physical sector of the Hilbert space. The condition G(n) = 0 is referred to as the Gauss law constraint as it is the discretized version of ∂ z E = j 0 , where j 0 is the charge density of dynamical fermions.
B. MPS for real-time evolution.
Similar as in [31,36] we block site n and link n into one effective site with local Hilbert space spanned by {|κ n = |s n , p n n : s n = −1, 1; p n ∈ Z}. In our approach we approximate the states of the lattice system eq. (2) by Here we have A κ (n) ∈ C D×D and v L , v R ∈ C D×1 . The MPS ansatz associates a matrix A κn (n) = A sn,pn (n) with each site n and every local basis state |κ n n = |s n , p n n . The indices α and β are referred to as virtual indices, and D is called the bond dimension. Note that this ansatz is T 2 invariant. As such we can consider the ansatz directly in the thermodynamic limit (N → +∞), bypassing any possible finite size artifacts. In this limit the expectation values of all local observables are independent of the boundary vectors v L and v R .
In our simulations, we start from the ground state of the Hamiltonian H α=0 without background field, for which we found a faithful gauge invariant MPS approximation |Ψ u [A(1)A(2)] by using the time-dependent variational principle (TDVP) [31,36,[52][53][54] (see Appendix Sec. A 1) for a brief review). At time t = 0, we perform a quench and apply a uniform electric field, which is simulated by evolving the ground state with respect to the Hamiltonian H α =0 with non-vanishing background field: ] . The evolution is performed using the infinite time-evolving block decimation (iTEBD) [55] which adapts the bond dimension of this MPS dynamically according to the Schmidt spectrum, see Appendix Sec. A 3. As explained there, the errors introduced by this method are well-controlled and argued to be only of order 10 −3 or smaller. We refer also to [31,36] for a discussion on the systematics of the iTEBD.

C. Observables and their discretization
We focus on the real-time evolution of the following observables Here, E(t) is the expectation value of the total electric field, j 1 = ψ γ 1 ψ t the current which can also be obtained from the electric field via Ampère's law (Ė = −gj 1 ) and Σ(t) is the discrete version of the chiral condensate ψ ψ t . We will use N (t) = Σ(t) − Σ(0) as a measure for the fermion particle number but notice that only in the non-relativistic limit m g we have a clear notion of electron and positron number. We will also show the real-time evolution of the half-chain von Neumann entropy, S(t) = −T rρ log ρ, where ρ is the density matrix of the half-chain subsystem. Finally, as explained in Appendix Sec. B, notice that E(t), j 1 (t), N (t) and ∆S(t) = S(t) − S(0) are UV finite and already close to the continuum limit x → ∞ for x = 100.

A. Weak-field regime
In [31] we found that the quasi-period of the oscillations of the electric field in the linear response regime (α ≤ 0.01) could be traced back to the first single-particle excitation of H 0 . However, for α 0.1 we observed that the quasi-period grows with α, and, hence, cannot be explained by the mass of the same single-particle excitation for each α. It turns out that for α 0.25 the original vacuum |Ψ 0 is well described as a small density coherent state of single-particle excitations of the quenched Hamiltonian H α . This leads to the oscillatory behavior of Fig. 1. Specifically, as we discuss below, this behavior can be explained quantitatively in terms of the matrix-elements of H 0 , and the considered observables E and N in the truncated Hilbert-space consisting of the  ground state and the two single-particle excitations of the quenched Hamiltonian H α .
As explained in more detail in Appendix Sec. A 4, for a given α, we approximate all observables O in terms of a series of the creation a † m (k) and annihilation operators a m (k) of the single-particle excitations |E m (k) of H α with energy E m and momentum k § ; and this up to first order in a m and a † m [56,57]: Here m = 1, 2 labels the two single-particle excitations of H α and the integral runs over the momenta k ∈ [−π, π]. Using the MPS approximations for the ground state and the two single-particle excitations obtained in [33,34], we can extract the coefficients o 1,m,n and o 2,m . For O = H 0 this leads to the approximation of the α = 0 ground state |Ψ 0 , as a coherent state of H α : a m (k) |Ψ 0 = d m δ(k) |Ψ 0 with d m ∈ C. This corresponds to a state with particle densities ρ m = √ x 2π |d m | 2 of the two zero-momentum single-particle excitations on top of the ground state of H α . In table I we display the obtained densities for different α in units of the correlation length ξ = 1/E 1 (0). One would expect our singleparticle approximation to hold as long as ξρ 1 , ξρ 2 1 which is in the line with our results. The approximation on the evolution for E(t) and N (t) is obtained by extracting the coefficients in eq. (6) for the appropriate operators (Eq. 5), and by considering the proper timeevolution a m (t) = a m e −iEmt . As can be observed in Fig. 1, the approximation works very well for α = 0.1, which lies already well beyond the linear response regime. For α 0.2 our approximation still predicts the right quasi-periods, but overestimates the amplitudes of the minima of E(t) and the amplitudes of the maxima of N (t) by approximately 20%. These discrepancies become § Notice that in this approximation we drop the multi-particle scattering states and can therefore consider the creation/annihilation operators a † m (k), am(k) as corresponding to the asymptotic inor out-states. larger when α increases and, eventually, when α 0.4 this approximation also fails in predicting the right quasiperiods (see Appendix Sec. A 4, in particular Figs. 6 and 7).
Finally, let us mention the holographic approach of [58] and the studies of other models with confinement [59,60] that also obtained an oscillatory behavior of local observables, bearing some resemblance with our results.

B. Strong-field regime
Let us now consider larger quenches: α ≥ 0.75. In Fig. 2 we compare the full quantum simulations (full line) with results obtained from real-time lattice gauge theory simulations [29] (dashed line). The latter should give reliable results as long as the classicality conditions are fulfilled, i.e., anti-commutator expectation values for typical gauge field modes are much larger than the corresponding commutators. This regime is characterized by non-perturbatively large field amplitudes [49].
Focussing first on the electric field E(t) and the current j 1 (t) ‡ one can observe good agreement between the MPS and real-time lattice simulations. The agreement further improves for growing α which is a nice cross-check for these two different techniques. However, for the particle number N (t) we find sizeable deviations. We attribute this discrepancy to differences in the initial states: the MPS simulation starts from the full ground state of the Hamiltonian H α=0 and hence incorporates interactions of the fermions with the fluctuating gauge field. On the other hand, the real-time lattice simulations are initialized in the bare Dirac vacuum that does not account for these interactions.
In a semi-classical picture the behavior of E(t), j 1 (t) and N (t) can be attributed to the nontrivial interplay between fermion and gauge field dynamics (backreaction) [29,47]: the electric field creates electron-positron pairs out of the vacuum and then accelerates them almost to the speed of light. This process costs energy, the electric field therefore decreases due to energy conservation so that particle creation terminates and the current saturates. After this initial creation of electron-positron pairs, which essentially occurs during the first oscillation of the electric field, we enter a regime of plasma oscillations, for which the onset at tg 3 can be observed in Figs. 2a and 2b. Also the behavior of the entanglement entropy ∆S(t) fits nicely with the semi-classical picture [61]: after the local production of entangled electronpositron pairs, the pairs will separate, entangling the system over even larger distances. From Figs. 2c and 2d one can indeed observe that the entropy starts increasing linearly after the initial period of pair production.
Even for large quenches we expect that the classicality conditions that underlie the real-time lattice technique are briefly violated during the times at which E(t) crosses zero. We can indeed observe in Fig. 2a that the full quantum MPS results start deviating from the real-time lattice results after the first transit through zero. In particular the MPS simulations predict a stronger damping. We interpret this damping as the onset of equilibration. It is accepted that a state which is brought out of equilibrium relaxes and equilibrates locally at late times [62]. In fact, it is believed that, under some generic conditions, the state thermalizes to a Gibbs state of the quenched Hamiltonian at a certain temperature [63][64][65][66][67][68][69][70][71][72][73]. There are however some exceptions, such as when the state as a whole is not thermal even if some local quantities already indicate thermalization [74][75][76], when the system is integrable and it converges towards a so-called generalized Gibbs ensemble [77][78][79][80][81][82][83][84], pre-thermalization [85][86][87][88] or many-body localization [80,[89][90][91][92]. ‡ Notice that the time evolution of the electric field and the current are connected through Ampère's law: j 1 (t) = −Ė(t)/g which holds at the operator level. We computed both quantities independently for the MPS simulations and found agreement with Ampère's law up to 10 −3 . For the real-time lattice simulations we derived j 1 (t) from E(t) using Ampère's law.
Under the assumption that the state would thermalize, we can determine its inverse temperature β 0 from energy conservation and by using our results from finite temperature simulations [35] (see Appendix Sec. A 5). In fig. 3 we compare E(t) and N (t) (full line) with its predicted thermal values E β0 and N β0 (dashed line). Note that our finite temperature simulations only enable us to determine β 0 numerically up to ∆β = 0.05, therefore we show the intervals E β0±0.05 and N β0±0.05 . Although the electric field seems to oscillate around E β0 , the amplitudes of the oscillations are still too large for a definite conclusion. On the other hand, one might be more tempted to say that N (t) is close to its thermal value for α = 1.25 and α = 1.5, although one should be cautious here as well.
To reach a definite conclusion, we would have to push the MPS simulations further in time. Unfortunately, the linear growth of entanglement, see Fig. 2d, requires the variational freedom of the MPS representation to grow exponentially in time (see Appendix Sec. A 3, in particular Fig. 4). This precludes computations at large tg and hence constrains the maximum time up to which we can reliably track the state.

IV. CONCLUSION
We demonstrated the potential of MPS to solve the real-time simulation of gauge field theories near the continuum limit, based on the paradigmatic example of an electric quench in the massive Schwinger model. For small quenches the real-time dynamics can be explained by using the single-particle excitations of the quenched Hamiltonian. For large quenches α = O(1), which is related to the phenomenon of Schwinger pair production, we compared the MPS simulations with results from realtime lattice gauge theory simulations and found good agreement between those methods. In this regime, we further investigated whether the state thermalizes at late times by using finite temperature simulations. While we found evidence that supports the onset of thermalization, the increase of entanglement prevented us to reach a decisive conclusion yet.
The MPS method provides a unique means to benchmark quantum simulators of the massive Schwinger model or related models using ultracold ions or atoms in optical lattices [40][41][42][43][44]. On the other hand, it is a major goal to extend this type of real-time simulation technique to more than one spatial dimension using projected entangled pair states (PEPS) [4]. The major progress on PEPS algorithms in the last decade [93][94][95][96][97][98][99][100][101][102] in combination with recent promising PEPS and TNS results for higher-dimensional gauge theories [103][104][105][106][107] makes us confident that this will be realized in the foreseeable future. In this section we explain the Matrix Product States (MPS) methods that are used for the Schwinger model. More specifically we discuss: 1. How the time-dependent variational principle (TDVP) is used to find the optimal translational invariant MPS approximation for the ground state in the thermodynamic limit (see Sec. A 1).
2. How we approximate the single-particle excitations using MPS (see Sec. A 2).
3. How we perform real-time evolution within the manifold of MPS using the infinite time-evolving block decimation algorithm (iTEBD) (see Sec. A 3).
4. How we use the MPS approximations for the ground state and the single-particle excitations to approximate the real-time evolution in the weakfield regime (see Sec. A 4).
5. How we determine the temperature of to the equilibrium state given that the state brought out of equilibrium by the quench thermalizes (see Sec. A 5).
More details and results can also be found in our earlier papers [31][32][33][34][35][36] We block site n and link n into one effective site with local Hilbert space spanned by {|s n , p n n : s n = −1, 1; p n ∈ Z}. Writing κ n = (s n , p n ) and a general state on this system of 2N sites takes the form (2)] assumes now a special form for these coefficients: Here we have A κ (n) ∈ C D(n)×D(n+1) and v L , v R ∈ C D(1)×1 . The MPS ansatz associates with each site n and every local basis state |κ n n = |s n , p n n a matrix A κn (n) = A sn,pn (n). The indices α and β are referred to as virtual indices, and D(n) are called the bond dimensions. Note that here A κ (n) only depends on the parity of n, in accordance with the T 2 symmetry of the Hamiltonian. As such we can consider the ansatz directly in the thermodynamic limit (N → +∞), bypassing any possible finite size artifacts. In this limit the expectation values of all local observables are independent of the boundary vectors v L and v R .
As explained in [31], to parameterize gauge invariant MPS, i.e. states that obey G(n) |Ψ(A) = 0 for every n, it is convenient to give the virtual indices a multiple index structure α → (q, α q ); β → (r, β r ), where q resp. r labels the eigenvalues of L(n − 1) resp. L(n). One can verify that the condition G(n) = 0 then imposes the following form on the matrices: [A s,p (n)] (q,αq),(r,βr) = [a q,s (n)] αq,βr δ q+(s+(−1) n )/2,r δ r,p , (A.2b) where α q = 1 . . . D q (n), β r = 1 . . . D r (n + 1). The formal total bond dimensions of this MPS are D(n) = q D q (n), but notice that, as (A.2b) takes a very specific form, the true variational freedom lies within the matrices a q,s (n) ∈ C Dq(n)×Dr(n+1) . To find the optimal ground state of H α within the class of gauge invariant MPS Eq. (A.2) we apply the timedependent variational principle (TDVP) [52][53][54] to the Schrödinger equation in imaginary time dτ = −idt. When τ → +∞ we indeed find the optimal approximation |Ψ u [A(1)A(2)] for the ground state of H α . As the Schmidt decomposition of it follows that to obtain a faithful approximation for the ground state one has to choose D q such that the discarded Schmidt values for each charge sector are sufficiently small. In particular we could take D q = 0 for |q| > 3 which is explained by the first term in the Hamiltonian Eq. (A.1). A proper justification of truncating the charge sectors is provided in [36]. We refer to [33,36] for the details on the TDVP.

MPS approximation for single-particle excitations
Once we have an MPS approximation |Ψ u [A(1)A(2)] for the ground state of H α , see Sec. A 1, we use the method of [53,108] to approximate the single-particle excitations. The ansatz for the single-particle excitations with momentum k that we will use is: where A(1) and A(2) correspond to the ground state where κ n = (s n , p n ) and b q,s1,s2 ∈ C Dq×Dr . The algorithm to find the optimal approximation |Φ k [B, A(1)A(2)] for the excited states is discussed in [31,36,53]: one has to find b q,s1,s2 such that is minimized with respect to b q,s1,s2 . This boils down in a generalized eigenvalue equation for b q,s1,s2 where the smallest eigenvalues correspond to the energies of the single-particle excitations. Only the ones who are stable against variation of the bond dimensions D q are physical. We refer to [31,34,36] for the details.
In [33,36] we found for m/g = 0.25 and α 0.47 two single-particle excitations with masses E 1 and E 2 . The energies at non-zero momentum are in the continuum limit determined by the Lorentz dispersion relations: E m (k) = k 2 + E 2 m . The corresponding MPS approximations at non-zero lattice spacing a = 1/g √ q,s1,s2 ] αq,βr δ p1,q+(s1−1)/2 δ p2,q+(s1+s2)/2 δ r,p2 (A.5) and are normalized such that [53] The delta-Dirac function originates from the infinite lattice length and has to be read as where 2N (N → +∞) is the number of sites on the lattice. For a local observable O = 2N −1 n=1 T n−1 oT −n+1 , where o is a Hermitian operator which acts only nontrivial on sites 1 and 2, we first subtract the ground state contribution such that With this renormalization we have that (2), B] are finite quantities that can be computed efficiently, see [53]. The delta-Dirac distributions have to be regularized according to Eq. (A.7).

iTEBD for real-time evolution
To evolve a state approximated by a MPS Eq. (A.2) at t = 0, i.e. to find we used the infinite time-evolving block decimation (iTEBD) [55]. At the core of this method lies the Trotter decomposition [109] which decomposes e −idtH into a product of local operators, the so-called Trotter gates. Specifically, we did a fourth order Trotter decomposition of e −iHαdt for small steps dt and projected afterwards |Ψ(t + dt) = e −iHαdt |Ψ u [A(1)A(2)] to a MPS |Ψ u [Ã(1)Ã(2)] with smaller bond dimensions D q . Similar as for the ground state, D q is chosen by discarding the Schmidt values smaller than a preset tolerance 2 in Eq. (A.3). In this way the virtual dimensions are adapted dynamically. We refer to [31] for the details on the implementation of the iTEBD. Taking a non-zero value for yields a truncation in the entanglement spectrum and hence a truncation in turn determines the required bond dimensions D q for every charge sector. For instance, in Fig. 4 we show how the maximum of the bond dimension over the charge sectors D max = max q D q varies with time for a given value of . It is this growth of the required bond dimensions, which can be traced back to the growth of entanglement, that makes the computations more costly at later times. Note that to save computational resources we imposed that D max ≤ 2000. As explained in [31] the simulation should be exact as → 0. Therefore the convergence in can be used to control the truncation error for a certain observable. In order to have a rough idea about the error for taking nonzero we compare the results for the simulation for the two smallest values of . We illustrate this in Fig. 5 for the electric field expectation value and the particle number N (t) where we compare the simulations for = 5 × 10 −5 (full line) with the simulations for = 1 × 10 −4 . As can be observed from the inset, where we plot the differences in magnitude of the electric field, and the particle number, the results are in agreement with each other up to at most 8 × 10 −3 . For other (smaller) values of α we found that this error was even smaller. Therefore we can trust that our results are reliable up to at least 1%.

Weak-field regime approximation
If H α0 is the Hamiltonian in an electric background field α 0 and H α is the Hamiltonian in an electric background field α, then we can write (up to an irrelevant constant) n=1 L(n) and = α − α 0 . Consider now the annihilation and creation operators a m and a † m of the single-particle excitations with energy E m (k) and momentum k of H α . In principle, they can obey either the canonical commutation relations for bosons or fermions, but as we will see later we need to impose boson statistics: [a n (k ), a † m (k)] = δ(k − k)δ m,n , [a m (k ), a n (k )] = 0, [a † n (k ), a † m (k)] = 0. (A.10) Using the TDVP, see Sec. A 1, we have a MPS approximation |Ψ u [A(1)A(2)] for the ground state of H α and by using the method discussed in Sec. A 2, we have a MPS approximation |Φ k [B (m,k) , A] for the m-th singleparticle excitation with momentum k and energy E m (k). They are normalized as The delta-Dirac functions originate from the infinite lattice length and have to be read as, see Eq. (A.7), where 2N (N → +∞) is the number of sites on the lattice. Within this approximation we have that We now want to express the ground state |Ψ(0) of H α0 in terms of the ground state |Ψ u [A(1)A(2)] and the single-particle excitations |Φ k [B, A(1)A(2)] of H α . We will expand H α0 in a series of powers of (a m (k), a † m (k)): where λ 0 , c m (k), µ m,n (k) ∈ C. The integrals over k and k run over all the momenta k, k ∈ [−π, π[. Note that we only displayed the operators that are nontrivial within the single-particle subspace. Indeed, in higher-order terms there appear products of the form a m1 (k 1 ) . . . a mn (k n ) or of the form a † m1 (k 1 ) . . . a † mn (k n ) for n ≥ 2 and, as such, these operators become trivial when projected onto the single-particle subspace. As we have only MPS approximations for the ground state and the single-particle excitations we need to restrict ourselves to the terms that are displayed in Eq. (A.15). Physically this means that we ignore the contributions of multi-particle eigenstates of H α .
Because H α0 is Hermitian, µ m,n should also be a Hermitian operator: µ m,n (k, k ) = µ n,m (k , k).
Using the ground state |Ψ u [A(1)A(2)] and the singleparticle excitations |Φ k [B (m,k) , A(1)A(2)] of H α it follows from Eq. (A.14) that As the energy is only determined up to a constant we can renormalize H α0 such that With this convention, it follows from Eq. (A.14) that we can compute the coefficients µ m,n and c m : and as the states are normalized according to Eq. (A.6), it follows from Eq. (A.8) that: where ] are finite quantities that we can compute efficiently (see [53] for the details).
Using Eq. (A.16) and (A.17) we rewrite H α0 , Eq. (A.15), now as H α0 is now diagonalized by the following transformations: where U (k) is the unitary transformation which diagonalize M (k) and E(k) is the diagonal matrix containing the eigenvalues of M (k), i.e. M (k) = U (k) † E(k)U (k). In vector notation we can write this transformation as One easily verifies now that  ii. In the Hamiltonian H α0 there appear terms of the form c m (k)c n (k) which is ill-defined as c m (k) ∝ δ(k). One can regularize this by replacing the Dirac-functions by δ(k) → δ k,0 2N/(2π) and the dk by dk → 2π/2N (2N the number of sites on the lattice, 2N → +∞).
iii. E r (k) should be positive, otherwise the quadratic expansion of H α0 in the creation and annihilation operators a † n (k) and a n (k) is certainly not a valid approximation anymore.
where O 1 ef f and O 2 ef f are finite quantities which we can compute efficiently.
Hence, we find and To perform real-time evolution with H α we will work in the Heisenberg picture. The creation-and annihilation operator a † m (k) and a m (k) satisfy the following differential equationȧ If we restrict the Hilbert space to the vacuum and the single-particle excitations we find that It follows that within this approximation: a m (k, t) = e −iEm(k)t a m (k) and a † m (k, t) = e iEm(k)t a † m (k).
In the Heisenberg picture Eq. (A. is the expectation value per site and is finite. Within this approximation all coefficients appearing above can be computed from the MPS approximations |Ψ u [A(1)A(2)] and |Φ k [B (m,k) , A(1)A(2)] for the ground state and the single-particle excitations of H α . In our case for m/g = 0.25 and x = 100 we have for the values of α we considered here two single-particle excitations. Hence, the sum over m runs from 1 to 2. We expect the above approximation to be true as long as the contribution of the multi-particle excitations of H α is negligible. Physically this means that the ground state |Ψ(0) of H α0 is a coherent state of the creation and annihilation operators of H α . This can be interpreted as the fact that |Ψ(0) is constructed from the ground state of H α with a small density of single-particles of H α on top of it. We can indeed expect that for small values of α and at early times that this is the case. In Figs. 6 and 7 we compare the real-time simulations with iTEBD (full line) with this approximation Eq. (A.34) and find agreement for α 0.2 while for α = 0.4 the difference between both results is quit large. A discussion is provided in Sec. III A.

Predicting the asymptotic thermal values of real-time evolution
In [35] we succeeded to approximate the Gibbs state ρ(β) at temperature T = 1/β by using Matrix Product Operators (MPO) with where P is the orthogonal projector onto the (G(n) = 0)subspace. If the state |Ψ(t) = e −iHαt |Ψ(0) would eventually equilibrate to a Gibbs state then we can estimate its inverse temperature β 0 from the requirement that Ψ(0)|H α |Ψ(0) = tr H α Pe −β0Hα tr (Pe −β0Hα ) , as follows from energy conservation during real-time evolution. In Figs. 8a and 8b we show the energy per unit of length E β of the Gibbs state ρ(β) as a function of β and the (conserved) energy per unit of length E(t) of the state |Ψ(t) . We subtracted from both quantities the energy per unit of length of |Ψ(0) . The intersection between the curves determines the value of β 0 . Because we simulated the thermal evolution with steps dβ = 0.05 we can only determine β 0 up to 0.05/g. For α = 0.75 we find β 0 g = 1.35(±0.05), α = 1.25 we find β 0 g = 0.85(±0.05) and for α = 1.5 we find β 0 g = 0.70(±0.05).
Appendix B: Scaling to the continuum limit of the real-time results In this paper we consider the following quantities: (a) The electric field: which counts in the weak coupling limit (m/g 1) the number of electrons and positrons per unit of length that are created out of the vacuum or destroyed in the vacuum due to turning on the electric background field α at t = 0. As we will show below, a UV quantity is obtained by considering the renormalized half chain entanglement entropy ∆S(t) = S(t) − S(0).
The fact that these quantities are UV finite is corroborated by Fig. 9 where we show the evolution of the electric field E(t, x), the particle number N (t, x) and the renormalized entropy ∆S(t, x) as a function of time for x = 100, 200, 300, 400. Note that we here explicitly denote the x−dependence of the quantities. We observe that for all these quantities the graphs are almost on top of each other, see Figs. 9 (a) -(c). One can also obtain a continuum estimate for these quantities by a polynomial extrapolation, see Fig. 9d where we perform a polynomial extrapolation for ∆S(t) for tg = 5. It is also clear from this example that we can already expect at x = 100 to be close to the continuum limit. (For the current j 1 (t) this follows from Ampère's law:Ė = −gj 1 .) This justifies that we restrict ourselves to x = 100 for the discussion on the continuum results. .