Confinement and string breaking for QED$_2$ in the Hamiltonian picture

The formalism of matrix product states is used to perform a numerical study of 1+1 dimensional QED -- also known as the (massive) Schwinger model -- in the presence of an external static `quark' and `antiquark'. We obtain a detailed picture of the transition from the confining state at short interquark distances to the broken-string `hadronized' state at large distances and this for a wide range of couplings, recovering the predicted behavior both in the weak and strong coupling limit of the continuum theory. In addition to the relevant local observables like charge and electric field, we compute the (bipartite) entanglement entropy and show that subtraction of its vacuum value results in a UV-finite quantity. We find that both string formation and string breaking leave a clear imprint on the resulting entropy profile. Finally, we also study the case of fractional probe charges, simulating for the first time the phenomenon of partial string breaking.


I. INTRODUCTION
The confinement of color charge in quantum chromodynamics (QCD) is one of the beautiful key mechanisms of the Standard Model. Focusing on the static aspect of confinement, one can probe the theory with a heavy quark-antiquark (qq) pair and examine how the modified ground state evolves as a function of the interquark distance [1]. For small distances a color electric flux tube forms between the pair, resulting in a static potential (i.e. the surplus energy of the modified ground state) that grows linearly with the distance. This flux tube can therefore be conveniently modeled by an interquark string with a certain string tension. One can then describe a heavy quarkonium state as a qq pair that is kept together by this confining string. However, there exists a critical distance at which the string breaks. Beyond this distance the flux tube disappears and the potential flattens out to a constant. At this point it has become energetically favorable to excite light particles out of the vacuum that completely screen both the probe quark and antiquark, leading to two isolated color singlets. In a dynamical setting these would then be the two freely propagating jets of hadrons that emerge as the final product of some particle collision.
This phenomenological picture is corroborated both by experiment and theoretical work. At the computational level, the static potential has been studied extensively over the years with lattice QCD. The linearly rising confining interquark potential has been obtained, both in the quenched case [2][3][4][5][6][7][8][9][10][11][12] that excludes dynamical light quark degrees of freedom and in the unquenched case [13][14][15][16][17][18] that includes these degrees of freedom. In the latter case, where the dynamical quarks can screen the heavy probe charges, the phenomenon of string breaking has also been observed [19] as an asymptotic flattening of the calculated potential. Nevertheless our understanding of confinement is incomplete: the Euclidean space-time lattice Monte Carlo simulations cannot access the realtime aspects of the dynamical string formation and string breaking. Furthermore, even in the static case, it is not settled yet if one can fully describe the confinement mechanism -specifically the nonperturbative string formation -in terms of (semi-) local degrees of freedom (e.g., center vortices, magnetic monopoles) [20][21][22].
In this paper, we study how confinement and string breaking show up in the Hamiltonian setup, as opposed to the Euclidean path integral setup of lattice Monte Carlo. We do this for the simplest nontrivial quantum gauge field theory: (1+1)-dimensional quantum electrodynamics (QED 2 ), also known as the Schwinger model [23]. The Schwinger model has a long-standing tradition as toy model for QCD, sharing its confining and chiral symmetry-breaking properties. (We will therefore often refer to 'quark' and 'antiquark' both for the external probe charges and for the light dynamical fermions.) But notice that in the future the significance of QED 2 could go well beyond being a toy model, as QED 2 or QED 2 -like theories might be realized effectively by quantum simulators [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43].
An important difference with QCD is that QED 2 already exhibits confinement at the perturbative level, as the Coulomb potential is linear in (1+1) dimensions. Furthermore, the theory can also be solved in a strongcoupling expansion, via bosonization [44]. We make extensive use of both the strong-and weak-coupling results in the analysis of our numerical results. Our simulations of the lattice Hamiltonian are performed close to the continuum limit, indeed allowing for a quantitative check against these analytic continuum results in the appropriate regimes.
Specifically, we simulate the modified vacuum structure in the presence of two probe charges and this for different distances and values of the charges. As we will show, already in this static case the Hamiltonian simulations give a complementary view on the confining properties of the theory. At the practical level, the direct access to the quantum state allows for a relatively easy calculation of all local observables. In this way we could not only extract the static interquark potential, but also arXiv:1509.00246v3 [hep-lat] 24 Nov 2016 for instance determine the detailed spatial profile of the electric string or the precise charge distribution of the light fermions around the probe charges. At a more fundamental level, our tensor network state simulations (see below) allow for a direct calculation of the entanglement entropy between different regions. In the past decade it has become clear that entanglement entropy is a very useful quantity for the characterization of quantum manybody systems and quantum field theories [45], in particular also for the investigation of the confining properties of gauge theories [46]. In this context the entanglement entropy is typically calculated either from the dual geometry in the AdS/CFT approach [46,47] through the Ryu-Takayanagi conjecture [48], or from lattice Monte Carlo simulations [49][50][51] through the replica trick, allowing for calculation of the discrete Renyi entropies. In contrast, tensor network state simulations give access to the full Schmidt spectrum of the state. The Schmidt spectrum {λ α } follows from the Schmidt decomposition: if |Ψ ∈ H A ⊗ H B is a state belonging to the tensor product of the two Hilbert spaces H A and H B , then one can write ∈ H B orthonormal unit vectors and λ α , called the Schmidt values, non-negative numbers that sum to one. From the Schmidt values one can calculate all Renyi entropies, including the von Neumann entropy. In our simulations we find that subtraction of the vacuum entropy results in a UV-finite entanglement (von Neumann) entropy and that both the string formation and string breaking leave characteristic imprints on this renormalized entropy.
As was mentioned in the previous paragraph, we use the general formalism of tensor network states (TNS) [52,53] for our simulations. Although the TNS formalism has been mainly developed in the context of condensed matter physics, it is actually a universal method in the same way that the Feynman diagrammatic approach has a universal character. The latter applies whenever the interactions are weak, whereas the TNS method applies whenever the interactions are local. It is in fact precisely the entanglement structure of low-energy states for local systems, captured by the so-called area law [54], which lies at the root of the TNS description.
In the next section we discuss the starting point of our simulations, introducing both the relevant lattice Hamiltonian in the presence of probe charges and the appropriate form of MPS that is dictated by gauge invariance. We then first consider the asymptotic case of two (fractional) probe charges at infinity in section III. In section IV we consider finite interquark distances and study how the ground state evolves as a function of this distance. We distinguish three different cases. First we consider the strong-coupling limit. This is a special case, since in this limit the interquark string never forms and all probe charges, fractional or integer, are screened asymptotically. We then go away from the strong-coupling limit, considering first the case of unit probe charges. In this case we clearly observe the transition from a string state at short interquark distances to a broken-string two meson state at large distances. Then, in addition to unit probe charges we also consider fractional probe charges, simulating for the first time the phenomenon of partial string breaking, with probe charges that get only partially screened. Finally in section V we present our conclusions. Technical details on our MPS simulations and on some perturbative weak-coupling calculations can be found in the appendices.

A. Hamiltonian and gauge symmetry
The Schwinger model is (1+1)-dimensional QED with one fermion flavor. We start from the Lagrangian density in the continuum: One then performs a Hamiltonian quantization in the timelike axial gauge (A 0 = 0), which can be turned into a lattice system by the Kogut-Susskind spatial discretization [75]. The two-component fermions are sited on a staggered lattice. These fermionic degrees of freedom can be converted to spin-1/2 degrees of freedom by a Jordan-Wigner transformation with the eigenvectors {|s n n : s n ∈ {−1, 1}} of σ z (n) as basis of the local Hilbert space at site n. The compact gauge fields θ(n) = agA 1 (n), live on the links between the sites. Their conjugate momenta E(n), with [θ(n), E(n )] = igδ n,n correspond to the electric field. The commutation relation determines the spectrum of E(n) up to a constant: E(n)/g = α(n) + p , with α(n) ∈ R corresponding to the background electric field at link n and p ∈ Z.
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 will then correspond to x → ∞. Note the different second (mass) term in the Hamiltonian for even and odd sites which originates from the staggered formulation of the fermions. In this formulation the odd sites are reserved for the charge −g 'quarks', where spinup, s = +1, corresponds to an unoccupied site and spindown, s = −1, corresponds to an occupied site. The even sites are reserved for the charge +g 'antiquarks' where now conversely spin-up corresponds to an occupied site and spin-down to an occupied site.
In the time-like axial gauge the Hamiltonian is still invariant under the residual time-independent local gauge transformations generated by: As a consequence, if we restrict ourselves to physical gauge-invariant operators O, with [O, G(n)] = 0, the Hilbert space decomposes into dynamically disconnected superselection sectors, corresponding to the different eigenvalues of G(n). In the absence of any background charge, the physical sector then corresponds to the G(n) = 0 sector. Imposing this condition (for every n) on the physical states is also referred to as the Gauss law constraint, as this is indeed the discretized version of ∂ z E − ρ = 0, where ρ is the charge density of the dynamical fermions. The other superselection sectors correspond to states with background charges. Specifically, if we want to consider two probe charges, one with charge −gQ at site 0 and one with opposite charge +gQ at site k, we have to restrict ourselves to the sector: Notice that we will consider both integer and noninteger (fractional) charges Q.
As in the continuum case [77], we can absorb the probe charges into a background electric field string that connects the two sites. This amounts to the substitution E(n) = g[L(n) + α(n)] where α(n) is only nonzero in between the sites: α(n) = −QΘ(0 ≤ n < k); and L(n) has an integer spectrum: L(n) = p ∈ Z. In terms of L(n) the Gauss constraint now reads: and we finally find the Hamiltonian [78]: (2.6) in accordance with the continuum result of [77].
In the following sections we will obtain ground-state approximations of this Hamiltonian, for different values of m/g, different values of the probe charge Q and different distances L = k/ √ x (in physical units g = 1) of the charge pair, all this for different lattice spacings 1/ √ x, focusing on the continuum limit x → ∞. An important point regarding the continuum limit is that the ground-state energy of the Schwinger model is UV divergent but that this UV divergence does not depend on the background field α(n). If we write E 0 = 2N 0 (with N = |Z|) for the ground-state energy of (2.6) with zero background field α(n) = 0, we have √ x 0 → −x/π for the energy density in the x → ∞ limit [79]. For the modified ground-state energy in the presence of the probe charge gQ pair at distance L we can then write

B. Gauge-invariant MPS
Consider now the lattice spin-gauge system (2.6) on 2N sites. On site n the matter fields are represented by the spin operators with basis {|s n n : s n ∈ {−1, 1}}. The gauge fields live on the links, and on link n their Hilbert space is spanned by the eigenkets {|p n n : p n ∈ Z} of the angular operator L(n). But notice that for our numerical scheme, we retain only a finite range: p min (n+1) ≤ p n ≤ p max (n + 1). We will address the issue of which values to take for p min (n + 1) and p max (n + 1) later in this subsection. Furthermore, it will be convenient to block site n and link n into one effective site with local Hilbert space spanned by {|s n , p n n }. Writing κ n = (s n , p n ) we introduce the multi-index With these notations we have that the effective site n is spanned by {|κ n n }. Therefore the Hilbert space of the full system of 2N sites and 2N links, which is the tensor product of the local Hilbert spaces, has basis {|κ = |κ 1 1 . . . |κ 2N 2N } and a general state |Ψ is thus a linear combination of these |κ : with basis coefficients C κ1,...,κ 2N ∈ C.
A general MPS |Ψ(A) now assumes a specific form for the basis coefficients [55]: are boundary vectors. The MPS ansatz thus 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 = max n [D(n)] is called the bond dimension.
To better understand the role of the bond dimension in MPS simulations it is useful to consider the Schmidt decomposition (1.1) with respect to the bipartition of the lattice consisting of the two regions A 1 (n) = Z[1, . . . , n] and A 2 (n) = Z[n + 1, . . . , 2N ] [56]: (2.8) Here |Ψ ) are orthonormal unit vectors living in the tensor product of the local Hilbert spaces belonging to the region A 1 (n) (resp. A 2 (n)) and λ α (n), called the Schmidt values, are non-negative numbers that sum to one. One can easily deduce that for a general MPS of the form (2.7) at most D(n + 1) Schmidt values will be nonzero (for the cut at site n (2.8)). We refer to appendix A and C for the computation of the Schmidt values for the specific case of our simulations and to [56,80] for the general case. We thus see that taking a finite bond dimension for the MPS corresponds to a truncation in the Schmidt spectrum of a state. The success of MPS is then explained by the fact that ground states of local gapped Hamiltonians can indeed be approximated very efficiently in D [81] and that the computation time for expectation values of local observables scales only with D 3 , allowing for reliable simulations on an ordinary desktop.
Another advantage of MPS simulations is that one can work directly in the thermodynamic limit N → ∞ [80,82], bypassing any possible finite-size artifacts. In the following, we work in this limit. In section III, where the Hamiltonian is invariant under translations (over two sites), the tensors A κn (n) depend only on the parity of the site n, see eq. (3.1). While in section IV the MPS ansatz is not translational invariant in the bulk, see eq. (4.1). In that case the tensors will be fixed asymptotically (n 1) to their ground-state value, anticipating that we approach the translational invariant ground state of the zero-background Hamiltonian. In both cases the MPS ansatz depends on a finite number of parameters. Finally, we note that, in the thermodynamic limit, the expectation values of local observables are independent of the boundary vectors v L and v R .
As explained in [57], to parametrize 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 (2.5) then imposes the following form on the matrices: [A s,p (n)] (q,αq),(r,βr) = [a s,p (n)] αq,βr δ q+(s+(−1) n )/2,r δ r,p , (2.9) where α q = 1 . . . D q (n), β r = 1 . . . D r (n + 1). The first Kronecker delta is Gauss' law (2.5) on the virtual level while the second Kronecker delta connects the virtual index r with the physical eigenvalue p of L(n). Because the indices q (resp. r) label the eigenvalues of L(n − 1) (resp. L(n)) and we only retain the eigenvalues of L(n−1) in the interval Z[p min (n), p max (n)] (resp. of L(n) in the interval Z[p min (n + 1), p max (n + 1)]) we have that D q (n) = 0 for q > p max (n) and q < p min (n). The formal total bond dimension of this MPS is D(n) = pmax(n) q=pmin(n) D q (n), but notice that, as (2.9) takes a very specific form, the true variational freedom lies within the matrices a s,p (n) ∈ C Dq(n)×Dr(n+1) .
Gauge invariance is, of course, also reflected in the Schmidt decomposition (2.8): for states of the form (2.9) the Schmidt values can be labeled with the same double index α → (q, α q ). More specifically, the Schmidt decomposition (2.8) now reads (see Appendix A and C): (2.10) As before, we observe that taking a finite bond dimension D q (n+1) corresponds to a truncation in the Schmidt spectrum, now of the charge sector q. The choice for the different bond dimensions D q (n + 1) in the different simulations should then be such that the discarded Schmidt values for each charge sector are sufficiently small. For our simulations with zero background, α(n) = 0, in [57] we could take D q = 0 for |q| > 3, i.e. p min = −3 and p max = +3. For the simulations with a nonzero background field we find that for the same accuracy it suffices to consider eight q-sectors. But -not surprisingly given the first term in the Hamiltonian (2.6) -we find the relevant eigenvalues sectors of L(n) to be centered around a dominant sector p 0 that can be shifted away from p 0 = 0 for some sites n. The largest Schmidt value in each q-sector decreases as we move farther away from q = p 0 . When |q −p 0 | 4 all the Schmidt values λ q,αq (n) are sufficiently small and we can safely take D q = 0, i.e. p max p 0 + 4 and p min p 0 − 4. We refer to appendices A and C for more details on the weight of the different sectors for the different simulations, see in particular figs. 12a and 12b and figs. 20c and 20d for some explicit examples.
From the Schmidt spectrum (2.10) one can extract different measures for the entanglement. In the following we will always use the von Neumann entropy [83]. For the half-chain cut at site n, to which we associate the position z = (n + 1/2)a in physical units, we then have: (2.11)

III. ASYMPTOTIC STRING TENSION
We first study the large distance behavior of the potential as captured by the asymptotic string tension σ Q = lim L→+∞ V Q (L)/L. This is the quantity that indicates whether the probe charges are asymptotically confined (σ Q = 0) or not (σ Q = 0). For the Schwinger model σ Q has been computed analytically in the strong-coupling expansion [77,[84][85][86]. At the numerical front the most successful computation to date used finite-lattice scaling methods in a Hamiltonian formulation [79]. An advantage of our MPS simulations is that in contrast to [79] we can directly work in the thermodynamic limit (N → ∞), leaving only the x → ∞ interpolation to extract the continuum results. The challenge of taking this continuum limit now lies in the diverging correlation length ξ/a (in lattice units), as MPS simulations require larger bond dimensions for growing correlation length [53].
To find the asymptotic string tension we put a probe charge −gQ at −∞ and a probe charge gQ at +∞. As we explain in the previous section, a probe charge pair translates to a background electric field α(n) in the Hamiltonian (2.6). In this case the background electric field will be uniform: α(n) = −Q, ∀n. The Hamiltonian is then invariant under T 2 , a translation over two sites. In accordance with this symmetry the appropriate MPS variational ground-state ansatz takes the form , and A κ (n) ∈ C D(n)×D(n+1) takes the form (2.9) (n = 1, 2). This corresponds to a general MPS (2.7) in the thermodynamic limit (N → +∞) where the tensors A κn (n) depend only on the parity of the site n: A κ2n−1 (2n−1) = A κ2n−1 (1) and A κ2n (2n) = A κ2n (2), ∀n. As a consequence D q (n), p min (n) and p max (n) also depend on the parity of n.
As we explain in appendix A we were able to accurately approximate the ground state and its finite energy per site Q = E Q /2N , with E Q the total infrared divergent energy, within the class of states (3.1). Therefore, we perform imaginary time evolution (dτ = idt) of the Schrödinger equation, i∂ t Ψ A(1), A(2) = H Ψ A(1), A(2) , with the time-dependent variational principle (TDVP) [80,82]. In appendix A we also explain how we chose the virtual dimensions {D q (1), D q (2)} and {p min/max (1), p min/max (2)} by investigating the Schmidt spectrum. In [57], we found the energy of the vacuum E 0 = 2N 0 for the zero-background field α(n) = 0. In the same fashion, we now compute the string tension σ Q as the extra energy density, induced by the uniform background electric field: From the numerical point of view it is important to take the convergence of this UV-finite quantity σ Q (x) as criterion for halting the imaginary TDVP time evolution. As we explain in more detail in appendix A we computed values for σ Q (x) in this way, for x = 100, 200, 300, 400, 600, 800 and perform a polynomial extrapolation in 1/ √ x similar to [79]. This indeed allows us to recover a finite value for lim x→∞ σ Q (x), thereby explicitly verifying that the UV divergencies in the energy densities √ x Q and √ x 0 cancel out.
In fig. 1a we plot our result for the continuum string tension σ Q computed for different values of the mass m/g as a function of the charge gQ of the external quark-antiquark pair. Note that we only consider Q-values ∈ [0, 1[ as the string tension is periodic in Q: Q → Q − p upon L(n) → L(n) + p for p ∈ Z in the Hamiltonian (2.6). Note also that one can combine this transformation for p = 1 with a CT -transformation (C = charge conjugation): This transformation gives Q → 1 − Q in the Hamiltonian (2.6) and therefore σ Q = σ 1−Q . So for our calculations we can restrict ourselves to values Q ∈ [0, 1/2]. In practice we consider the explicit values: Q = 0.05, 0.10, 0.15, ..., 0.45, 0.47, 0.48, 0.5 and perform an interpolating fit.
Our considered values for m/g interpolate between the strong-and weak-coupling regime. In the strong-coupling regime m/g 1 the string tension is computed in mass perturbation theory from the bosonized field theory up As one can observe in fig. 1c for m/g → 0 our results indeed converge to this analytic result that is plotted with a dashed line for m/g = 0.125 and for m/g = 0.25.
In the weak-coupling regime g/m 1 we can easily compute the string tension in standard perturbation theory from the continuum Lagrangian (2.1) (see appendix B). Up to order (g/m) 4 we find the string tension for Q ≤ 1/2: with the value for Q > 1/2 following from the identification σ Q = σ 1−Q for the compact formulation of QED 2 that we are considering. In fig. 1d one can again observe the convergence of our numerical results to this analytic result, now for g/m → 0. Notice here that we subtract the leading order term of (3.4).
Comparing the strong-and weak-coupling regime we observe an important difference: in the strong-coupling limit σ Q is differentiable at Q = 1/2 whereas in the weakcoupling limit this is not the case. Therefore there exists a critical mass (m/g) c with the property that σ Q is differentiable at Q = 1/2 for (m/g) < (m/g) c and not differentiable at Q = 1/2 for (m/g) > (m/g) c . This point (m/g) c corresponds to the first-order phase transition for the Hamiltonian H Q (2.6) at Q = 1/2 [79]. H Q=1/2 is symmetric under the CT transformation (3.2) and the point (m/g) c separates the unbroken phase m/g < (m/g) c from the spontaneously broken phase m/g > (m/g) c that was originally predicted by Coleman [44]. This relationship of the breaking of CT -symmetry with the nondifferentiability of σ Q can be made more concrete by noting where . . . Q denotes the expectation values with respect to the ground state of H Q . We now have the relation , which indeed makes it a good order parameter for the CT breaking at Q = 1/2. We perform an independent computation of E Q , again for Q = 0.05, 0.10, 0.15, ..., 0.45, 0.47, 0.48, 0.5, and now using values x = 100, 200, 300, 400 for our continuum extrapolation (see appendix A). Our results are displayed in fig. 1b. At Q → 1/2 we find for m/g = 0.3, E Q /g = 0 up to a numerical error of 4 × 10 −3 while for m/g = 0.35 we find E Q /g = 0.314(2), consistent with the value (m/g) c ≈ 0.33 that was obtained in [60] and also consistent with the behavior of σ Q in fig. 1a.
Finally we also compute the half-chain entropy S (2.11) for different values of Q and m/g, which in this translational-invariant case does not depend on the position of the cut. As such the entropy is a UV divergent quantity, but one expects the divergence to come from the fermion kinetic term in the Hamiltonian (2.6) and therefore be Q-independent. Specifically, the general results of Cardy and Calabrese [87] predict for two fermionic degrees of freedom a UV divergence (with correlation length ξ in physical units): A. There, we also explicitly extract the coefficient −1/6 of the logarithmic term by a logarithmic fit to the data. The errors are of the order 10 −3 for m/g 0.5 and only of order 10 −4 for m/g 0.5, see table IV in appendix A 3.
The universality of the logarithmic UV divergence then allows us to define a UV-finite renormalized entropy ∆S Q ≡ S Q − S 0 , with a finite continuum value that can be obtained by a polynomial extrapolation in 1/ √ x, see inset fig. 2a. Contrary to the string tension and the electric field, we found sometimes that the results at x = 100 and the continuum results differ by a factor of order one or have different sign. We refer the reader to subsection A 3 in appendix A and, in particular, to fig. 18 for the details about the continuum extrapolation. In fig. 2b we show this renormalized entropy ∆S Q as a function of Q for different values of m/g. Most notably we observe an (almost) divergent behavior for m/g = 0.3 at Q → 1/2 close to the critical point Q = 1/2, (m/g) c ≈ 0.33. From (3.6) we indeed expect a growing entropy for growing correlation length. By the same argument one can understand the behavior at small Q-values: there the correlation length (inverse mass gap) increases with growing g/m [88], which is indeed paralleled by the behavior of ∆S Q .

IV. FROM SMALL TO LARGE DISTANCES
Now, we consider the situation where the external quark and antiquark pair are separated over a finite length L. On a lattice with spacing a and interquark distance L = ka, the pair introduces a nonuniform background electric field α(n) = −QΘ(0 ≤ n < k) in the Hamiltonian (2.6).
As ansatz for our MPS trial state |Φ(B) for the ground state we now write: where r L 0 ≤ k r R and A κ (n) = A κ (n mod 2) corresponds to the MPS approximation (3.1) of the ground state of the zero-background Hamiltonian (α(n) = 0) and depends only on the parity of n. This is a MPS (2.7) in the thermodynamic limit (N → +∞) where we take A(n) = B(n) for r l ≤ n ≤ r R − 1 and take the A(n) corresponding to the ground state (3.1) for α(n) = 0 to the left and to the right of the B(n)'s (n < r L and n ≥ r R ).
The idea behind this ansatz is that the nonuniform background electric field changes the vacuum and breaks translation invariance (all B(n) are different) but that asymptotically (|n| 1) it does not af- ). Note that we allow different bond dimensions on different sites. Also a s,p (n) is of the form (2.9) as we impose this to determine the ground state of the zero-background electric field Hamiltonian. Because (4.1) is linear in each of the B n we can use the DMRG-method [89] to obtain the best approximation for the ground state within the manifold of gaugeinvariant states, by optimizing on the UV-and IRfinite quantity V Q (L). By looking at the Schmidt spectrum we are able to fix the values of the virtual dimension D q (n) and the minimum and maximum eigenvalues p min (n + 1) and p max (n + 1) of L(n) we retain in our numerical scheme to obtain an accurate approximation of the ground state. The choices for r L and r R , which vary between −k/2 − 250 ≤ r L ≤ −k/2 − 150 and k/2 + 150 ≤ r R ≤ k/2 + 250 are checked a posteriori by verifying the convergence of local observables at large distances to their zero-background value. We refer the reader to Appendix C for the details. In fig. 3a we plot our results for the potential for m/g = 0 for Q = 1 and Q = 1.5. This can be compared with the exact continuum result [84]: which is indeed of the Yukawa type. We find very good agreement already for x = 100, for both Q = 1 and Q = 1.5. For Q = 1 we also perform a computation for x = 400, in the inset of fig. 3a one can observe the rate of convergence towards the continuum x → ∞ in this case. The charge density ψ (z)γ 0 ψ(z) of the light quarks is, of course, also an interesting quantity to compute, as it explicitly shows the screening of the external probe charges. The analytical result for the probe charge pair put at ±L/2 reads [84]: (4.3) This indeed corresponds to a charge distribution with two 'clouds' of oppositely charged light (in this case massless) quarks, around the external quark and antiquark, that for large distance L have exactly the same total charge ±Q as the external pair. On the lattice the charge density at In fig. 3b, we plot this density for Q = 1.5 where the charges are separated at distances Lg = 5.1 and Lg = 17.3. Here, too, our results for x = 100 are already very close to the continuum result.
In fig. 4 we show the spatial profile of the renormalized half-chain von Neumann entropy ∆S Q (z) = S Q (z) − S 0 (z) for different values of Lg. We compute this quantity for z = (n + 1/2)a with n even and perform an interpolating fit. When the heavy quarks are close to each other, ∆S Q (z) shows a peak in the middle between the charges and falls off very fast with |zg|, see fig.  4a. For larger values of Lg a cloud of light quarks forms around each of the heavy charges which clearly leaves its imprints on the spatial profile of the von Neumann entropy, see fig. 4b. ∆S Q (z) is nonzero around each of the heavy charges and is zero around zg ≈ 0.
The observed spatial profiles of the von Neumann entropy, however, are lattice artifacts and vanish in the continuum limit (x → +∞). Indeed, from the bosonized Hamiltonian for m/g = 0, it can be observed that any position-dependent electric background field can be transformed away [77,91] up to a position-dependent constant. Therefore the von Neumann entropy of the ground state with α(n) = 0 and α(n) = 0 are the same, hence ∆S Q (z) = 0 for any value of Lg. By investigating the scaling towards x → +∞ in figs. 4a and 4b we indeed observe that ∆S Q (z) tends towards a very small value for x → +∞. Note that we here need to perform an interpolation because we can only take Lg to be an integer multiple of 1/ √ x. Specifically, we perform simulations for For m/g = 0, the asymptotic string tension σ Q vanishes only for integer charges Q. This is taken to be an indication for a screeningà la QCD [77], where the potential exhibits a string tension (V Q (L) ∝ L) at short distances, but flattens out completely at large distances, at least for integer charges Q. At these large distances it becomes energetically favorable to materialize light (yet massive) (anti)fermions out of the vacuum that bind to the external quark and antiquark, resulting in two charge neutral mesons.
Historically, for QCD, lattice Monte-Carlo simulations succeeded first to calculate numerically the short distance confining behavior of the potential -both in the quenched and unquenched approximation -via the expectation value of the Wilson [1,22]. The detection of string breaking has posed a larger challenge. A main problem with  the use of the standard Wilson loop is the poor overlap with the broken-string two-meson state. This problem was finally overcome by including light quark propagators in the Wilson loop and analyzing its mixing with the standard Wilson loop [19,92]. For the Schwinger model the string-breaking phe-nomenon has been confirmed in mass perturbation theory [86] and in a semiclassical approximation of the bosonized version of the theory [93,94]. At the numerical level, for Q = 1, lattice Monte Carlo simulations have detected both the confining and string-breaking behavior of the potential [90,95]. In [90] the problem with the Wilson loop was avoided by computing instead the expectation value of the bosonized Hamiltonian, while [95] turned to very high statistics thereby explicitly showing the poor overlap of the Wilson loop with the broken-string ground state.
For the local quantities (charge density =ψ(z)γ 0 ψ(z), electric field=E(z)) and the potential, we restrict ourselves from now on to lattice spacing x = 100(= 1/g 2 a 2 ); from the previous subsection we can expect these results already to be quite close to the continuum. In fig. 5a we display our results for the potential, and this for different values of m/g. We compute explicitly the ground-state energy at Lg = 0.1, 0.3, . . . 15.3 and perform an interpolating fit. We clearly find a transition from the confining behavior, associated with the string state, towards the constant behavior associated with the broken-string twomeson state. This transition happens more suddenly for larger values of m/g, which is in qualitative agreement with the semiclassical results from the bosonized theory [93,94]. This is also what one would expect from the nonrelativistic weak-coupling regime, where the transition can be understood as a level crossing between the zero-particle string state and the two-particle brokenstring meson state (see appendix D). The dashed lines in fig. 5b corresponding to this nonrelativistic result for E string = Lg 2 /2 and E 2meson = 2m + 1.0188 g 4/3 m 1/3 were plotted for comparison. We can indeed observe the convergence towards this result for increasing values of m/g.
We further illustrate this behavior in fig. 5c, where we plot the total charge Q − of the light fermions on the negative z−axis: (4.4) One can observe indeed that the interpolation between Q − = 0 (string state) for small L and Q − = 1 (meson state) for large L becomes more and more discontinuous for growing m/g in accordance with the nonrelativistic level-crossing picture. In figs. 5d and fig. 5e we investigate the interpolation from the string state to the string-broken state in more detail for m/g = 0.75 by plotting the charge density and electric field. For L/g = 0.5 there is only a very small charge cloud around the external quark and antiquark, notice also the very short electric field string displayed at the bottom of fig. 5d. At L/g = 5.1 the clouds start to build up, lowering the electric field value at the center. At L/g = 10.1 the string is completely broken, the electric field at the center has vanished, and we have two clouds of total charge ±1 around the external quark and antiquark. At L/g = 17.3 the two isolated mesons are the ground state of the one-particle problem in a linear potential (appendix D), where Ai is the Airy function [96] and N the normalization factor. As one can observe, the charge distribution from this nonrelativistic picture matches very well our exact (numerical) result. In fig. 5f we compare the charge cloud at the negative z−axis with the nonrelativistic result for other values of m/g. One can again observe the convergence to the nonrelativistic result for growing m/g, notice that already for m/g = 0.5 the match is quite good. For the renormalized von Neumann entropy ∆S Q (z), we also find a characteristic picture, both for the string state and the string-broken state, see fig. 6 for the case m/g = 2. For the string state, Lg 9.5, the entropy shows a constant surplus in between the probe charges, similar to the electric field. But notice that this effect becomes very small in the continuum limit (green line), we find an extrapolated value: ∆S g (z) ≈ 2.0(5) × 10 −3 for zg ∈ [−2.5, 2.5]. For the string-broken case Lg 10, see fig. 6b, we find that the entropy now shows two clouds around the heavy quark and the heavy antiquark, similar to the charge density. But notice that in contrast to the string state the entropy now survives the continuum limit, with the x = 100 value already close to the continuum extrapolation. nomenon of partial string breaking. Indeed, in the nonrelativistic limit m/g → ∞ of string breaking due to meson formation, probe charges Q can only be screened by an integer number: Q →Q = Q − n, where n is the number of light (anti-)quarks that bind to the external charges. For nonzeroQ, i.e. when Q is noninteger, this still leaves a string between the two separated meson configurations. A visualization of this process in the nonrelativistic limit is shown in fig. 7.
Our simulations allow us to verify to what extent this picture is realized for finite m/g. In fig. 8 we plot our results for different values of Q, both fractional and integer. We do indeed recover partial string breaking, largely following the nonrelativistic picture. To our knowledge this is the first successful simulation of partial string breaking in the Schwinger model, a previous Monte Carlo simulation [95] failed to detect the phenomenon.
In fig. 8b we plot, as in the previous section, the evolution of the total dynamical charge Q − at the negative z-axis, for m/g = 1. For all values of Q this charge Q − indeed makes quasidiscrete jumps of ∆Q − ≈ +1 which should correspond to (partial) string breakings. As we see in fig. 8a these jumps indeed correlate with jumps in the string tensions in the different regions of the potentials. For m/g = 0.5 we still find jumps of Q − but they are smoothened out, as can be seen in figs. 8c and 8d. For m/g = 0.25 the jumps are even more smoothened out as can be seen in figs. 8e and 8f. This smoothened behavior, similar to what we obtain in the Q = 1 case, is expected as we go further from the nonrelativistic large m/g regime. But still note the contrast with the behavior in the massless limit m/g = 0 of subsection IV A, where the charge Q − grows continuously to the external value Q, assuring a complete screening.
For L going from 0 to ∞, different partial string break- ings should lead to the asymptotic behavior of the potential that we examine in section III. In table I we show the difference of the slope of the potential around Lg = 15.3 with the asymptotic string tension at x = 100 that we calculate in the previous section. The former is estimated as the mean of the backward differences for Lg = 13.7, 14.1, 14.5, 14.9. The error is computed as the standard deviation of these backward differences.
One observes that for m/g = 1 the string tension has    already converged to the asymptotic result, almost up to the numerical precision, while for m/g = 0.5 we are already very close to the asymptotic result and for m/g = 0.25 there is a slightly larger (but still very small) difference.
For integer values of Q, the asymptotic string tension vanishes, so asymptotically we expect Q − → Q, corresponding to a complete screening. For the values Q = 1 and Q = 5 that we consider, this is already almost satisfied at Lg = 15.3, as can be seen in table II. In the nonrelativistic limit for general Q, the total dynamical charge Q − that is produced asymptotically, will be the integer number that minimizes |Q − Q − |. For finite m/g we expect corrections to the nonrelativistic limit, but as one can see in the table these corrections are still very small for m/g = 1 and m/g = 0.5. Notice also that for the half-integer values Q = 2.5 and 4.5, for which we have spontaneous symmetry breaking in the asymptotic limit (see section III), we find Q − approaching the smallest of the two possible nonrelativistic values Q − ≈ Q − 1/2.
In fig. 9 we show the spatial charge distribution and electric field for different distances of the probe quarks. For Q = 1.75 we have two partial string breakings. The first one, around Lg ≈ 1.7 (see fig. 8) brings the electric field string at the center from E/g ≈ −1.7 to E/g ≈ −0.7. After the second partial string breaking, around Lg ≈ 9, the probe charge is 'overscreened', Q − ≈ 2, leading to a final electric field string with opposite sign E/g ≈ +0.2. Notice that in contrast to the z · g In fig. 10 we show the effect of different partial string breakings on the entropy profile ∆S Q (z), for m/g = 0.5 and Q = 4.5. For the smallest interquark distance Lg = 0.55, the entropy peaks at the center. At Lg 2.55 (after two string breakings, see fig. 8), we observe a profile with two peaks around the positions of the probe charges. At Lg 7.35 and Lg 13.15, after four string breakings, the profile now shows four peaks around the probe quark positions. In addition, we find an entropy surplus in the center, which now seems to be stable under the continuum extrapolation.
In fig. 11 we show that this characteristic imprint on the entropy is generic. We plot ∆S Q (z) for Lg cases two partial string breakings lead to the asymptotic meson state, but in the former case the final electric field is 'overscreened' while in the latter case the final electric field is 'underscreened'. Finally, notice that we can trust these results to be close to their continuum value, as the variation for the different x-values is very small.

V. CONCLUSIONS
In this paper we employed the MPS formalism for a detailed numerical study of the confining mechanism in the static limit for the massive Schwinger model. Our Hamiltonian setup gives us direct access to the modified vacuum state in presence of two probe charges. This allowed us, not only to compute the interquark potential, but also the spatial profile of the electric field between the probe charges and the charge concentration of the light fermions. Even for relatively small m/g the picture that emerged can be understood as a smoothened version of the nonrelativistic limit, with a level crossing between the electric string state that is the ground state at short distances and the broken-string two meson state that is the ground state at large distances. Here the two isolated mesons each consist of a light (anti-)quark cloud around the heavy probe charge, that is well described by the solution to the Schrödinger equation of the appropriate one-particle problem.
In the case of fractional probe charges, we clearly observed the expected partial string breaking. Again in accordance with the nonrelativistic picture we found the screening of the probe charges to happen in jumps ∆Q ≈ 1 of the light fermion charge; with these jumps becoming more and more discrete for growing m/g.
Our tensor network simulations also give us direct access to the full Schmidt spectrum for the different bipartitions on the state. The numerical simulations show that the UV divergence in the corresponding von Neumann entropy is universal, allowing us to define a UV-finite renormalized entropy by subtracting the vacuum value. We have examined the imprint of both the string formation and string breaking on the profile of this renormalized entropy. Most notably we found that string breaking leaves a very distinct imprint on this entropy profile.
We have checked our results not only against the predictions from the one-particle Schrödinger equation (D.3), but also against the weak-coupling results from the original Lagrangian (2.1) and against the strongcoupling results from the bosonized field theory [44]. In the appropriate regimes we found nearly perfect agreement with these continuum analytic results. This not only demonstrates the potential of MPS simulations close to the continuum critical point of a lattice theory. But it also serves as a nice, if not unexpected, cross-check of the consistency of all different descriptions of the Schwinger model.
We have restricted ourselves to the study of the static limit of the confinement mechanism. An obvious future extension of our work is to consider the dynamical problem, simulating the real-time hadronization that takes place in a realistic scattering process. MPS real-time simulations of this type of problem were considered recently for U (1) and SU (2) quantum link lattice models in [58] and [59]. One could also approach this problem by first calculating the scattering eigenstates [97,98]. See also [99,100] for an approach in the semiclassical limit.
Of course it will also be very interesting to bring this type of analysis to higher dimensions. Specifically, present techniques should already allow one to simulate e.g., the static confinement problem for some simple 2+1 dimensional lattice models (see [73] for homogeneous ground-state simulations of a Z 2 model and [74] for simulations of a U (1) lattice model with dynamical fermions). However, the current algorithms for PEPS simulations scale unfavorably with the bond dimension [101] and we therefore expect that the successful simulation of specific microscopic gauge field Hamiltonians in the continuum limit will require new techniques. Still, given the continuous progress of PEPS methods [102][103][104][105] we are hopeful on that front. In any case, in light of the potential for real-time and finite fermion density simulations, and of the new insight that might come from understanding the entanglement structure, it should certainly be worthwhile to further explore this direction, hopefully succeeding one day in the full simulation of the microscopic Hamiltonian of (3+1)-dimensional QCD.
Appendix A: Computation of the asymptotic string tension and electric field In this appendix we discuss the details of the computation of the asymptotic string tension σ Q and the electric field E. In subsection A 1 we discuss how to obtain the electric field E(x) and the energy density (x) for a fixed lattice spacing ga = 1/ √ x and in subsection A 2 we discuss how we to extrapolate these quantities to the continuum limit (x → +∞).
To where |ψ A1(2n) q,αq (resp. |ψ A2(2n) q,αq ) are orthonormal unit vectors in the tensor product of the local Hilbert spaces in the region A 1 (2n) (resp. A 2 (2n)). The Schmidt values λ q,αq , which are non-negative and sum to one, can be obtained as follows: assume A κ1,κ2 is brought in a canonical form such that the matrices r and l corresponding to the right and left eigenvectors of the largest eigenvalue of the transfer matrix [80], are positive definite and diagonal. Here we assume that the largest eigenvalue of the transfer matrix is normalized to one. Then, because A takes the form (A.6), r and l will also be degenerate in the eigenvalues of L(2n): [r] (q,αq);(r,βr) = r q,αq δ q,r δ αq,βr , [l] (q,αq);(r,βr) = l q,αq δ q,r δ αq,βr . The Schmidt values λ q,αq are now obtained by multiplying r and l: λ q,αq = r q,αq l q,αq where q ∈ Z[p min , p max ] labels the eigenvalues of L(2n) and α q = 1 . . . D q labels the variational freedom of the matrices a s1,p1,s2 . As can be observed from eq. (A.7), truncating to a finite bond dimension thus corresponds to an effective truncation in the Schmidt decomposition of the ground state. Ideally one would want a distribution of D q -values such that the smallest retained Schmidt value is more or less equal for each eigenvalue sector of L(2n). Then if we want a reliable MPS approximation for the ground state, these smallest retained Schmidt values should be sufficiently small, which corresponds to taking D q sufficiently large. In practice we do several simulations and adapted D q until the smallest Schmidt value in each eigenvalue sector of L(2n) was of order 10 −17 , i.e. min αq λ q,αq ≈ 10 −17 .
In figs. 12a and 12b we plot the distribution of the Schmidt values among the eigenvalue sectors of L(2n) for the final MPS ground-state approximations for m/g = 0.75, x = 400 and Q = 0.2, 0.45. As in [57], we observe that the sectors corresponding to q = 0, −1, 1 are the most dominant ones which justifies our choice of taking D q = 0 for |q| > 3. This can be understood physically from the term proportional to [L(n) − Q] 2 in (A.1) which punishes large eigenvalues of L(n). We also display the bond dimensions for each sector and for each simulated value of x in figs. 12c and 12d. One can observe that as x increases we need larger D q for the same accuracy. This is explained by the fact that the correlation length diverges as we approach the continuum limit (x → +∞) and it is well known that critical theories require larger bond dimensions for a good MPS approximation. For the same reason we also need larger D q when we are getting closer to the phase transition at m/g = (m/g) c ≈ 0.33 and Q = 1/2 [44,79].

Continuum extrapolation of the string tension and the electric field
In the second part of this appendix we discuss how we obtain an estimate for the continuum value of σ Q and E Q . Note that the string tension at x = 1/g 2 a 2 is obtained from the energy density by: where Q (x) is the ground-state energy per site of the Schwinger Hamiltonian (A.1). As for reduces to the XY -model we have that and it is argued in [79] that Q (x)/ √ x should behave polynomially as a function of 1/ √ x for large x, we have: This means that the energy densities √ x Q (x) and √ x 0 (x) are UV divergent. But as we see, the string tension which is the difference of these quantities is UV-finite and, thus, we should also have C Q = C 0 . However, from the numerical point of view it is clear that small errors in (A.10) or/and in C Q and C 0 would lead to large errors in the extrapolated continuum value lim x→∞ σ Q . To avoid this problem we first calculate 0 and subtract it from the Hamiltonian (A.1): the ground state of the renormalized Hamiltonian. As follows from (A.11), for large x, σ Q (x) should scale as In our simulations we computed σ Q (x) for x = 100, 200, 300, 400, 600, 800. Our estimate σ est Q is obtained by fitting the σ Q (x) corresponding to the five largest x to continuum estimates by fitting all the data to f 1 (x) (see (A.13)) and all our data to (A.14) The error err σ Q is taken to be the maximum of the difference of σ est Q with these two other estimates. In fig. 14a we show the log 10 of err σ Q as a function of Q for m/g = 0.125, 0.3, 0.35, 0.5, 1. It is clear that these errors are quite small. We have the largest error for m/g = 0.3 and Q = 0.5 which is explained by the fact that the gap is very small there as we are in the vicinity of a phase transition [60]. As mentioned above, it is well known that for smaller mass gaps, for a given bond dimension, the error on the ground-state MPS approximation will be larger.
The continuum extrapolation of the electric field, is found in a similar way. Now we use the values computed at x = 100, 200, 300, 400 and perform a linear fit, through the three largest x−values. The fact that we again have analytical behavior as a function of 1/ √ x can be observed from fig. 15 where we display the electric field as function of 1/ √ x. It is also a consequence of the fact that E(x) = −dσ Q (x)/dQ and we already argued that σ Q (x) is analytical as a function of x. To make our estimate more robust against the choice of the interval and the fitting function we compute estimates by a linear fit (A.16) through all the points (x = 100, 200, 300, 400) and a quadratic fit, through all the points. Again, the error err E is taken to be the maximum of the difference with these two estimates. The log 10 of err E is displayed in fig. 14b. The errors are quite small but become larger again around the phase transition at the critical mass (m/g) c ≈ 0.33 when going towards Q = 1/2. At Q = 1/2 we do not display our error because this is a special case. For m/g < (m/g) c the CT symmetry is not broken and, thus, we should have E = 0, and this for all values of x. Therefore a continuum extrapolation of E is useless, see fig. 16a. To obtain an error bound we take the largest value in magnitude of E(x) for x = 100, 200, 300, 400. It is displayed in table III. When m/g > (m/g) c we have two different vacua with opposite sign for the electric field. We will always take the negative sign which comes down to taking the vacuum in the limit Q → 1/2 for Q < 1/2. In this case it is possible to perform a polynomial extrapolation, see fig. 16b. The results are given in Table III. If possible we compare with [60].   as a function of 1/ √ x, see inset figs. 17a and 17b. Therefore, we should be able to fit S 0 (x) to a function of the form and find B 0 = 1/6. Specifically, we fit our data corresponding to the largest five x−values, x = 200, 300, 400, 600, 800, against f 1 to obtain a first estimate for B 0 . To have some robustness against the choice of fitting interval and the fitting function, we also include our result for x = 100 and fitted all our data against f 1 and against This gave us two other estimates for B 0 . In table IV we give the results for B 0 + 1/6. The value that is shown is the largest value for B 0 + 1/6 (in magnitude) from the three fits, i.e. the largest error on the predicted result of [87]. As one observes these errors are at most 2 × 10 − To compute σ Q in the weak-coupling expansion we start from the Lagrangian (2.1) and include a current j µ = g µν ∂ ν Q; with Q constant everywhere in the bulk, and Q → 0 only at the boundaries at infinity (see [77]): where on the last line we perform a partial integration andF µν ≡ µν gQ.
The effective action, obtained by integrating out both the fermion and the gauge fields in the path integral, will then have the general form: where we can exclude derivative terms sinceF µν is constant. At next-to-leading order we find for the first coefficient C 0 : The zero-order term here is the tree-level result while the g 2 /m 2 term follows from the one-loop Feynman diagram on the first line of fig. 19, which can be calculated with standard techniques (see e.g. [108]). Furthermore one can see that all other nonzero diagrams will lead to contributions to the coefficients C i that are at least order g 4 /m 4 . Finally, we can then identify S ef f = d 2 x σ Q , leading to the result (3.4).
To obtain the best approximation within this class of states of the ground state of the Hamiltonian (C.1) we have to minimize with respect to b(r L ), . . . , b(r R − 1). This is a perfect problem to tackle with the DMRG [89]. We briefly sketch how this works in our case. The DMRG first minimizes H(b, b) with respect to b(r L ) while keeping b(r L +1), . . . , b(r R −1) fixed, then minimizes H(b, b) with respect to b(r L + 1) while keeping b(r L ), b(r L + 2), . . . , b(r R − 1) fixed and so on until b(r R − 1). After this sweep, it will sweep back: minimizing H(b, b) with respect to b(r R − 1) while keeping b(r L ), . . . , b(r R − 2) fixed, then minimizing H(b, b) with respect to b(r R − 2) while keeping b(r L ), . . . , b(r R − 3), b(r R − 1) fixed and so on until b(r L ). The algorithm keeps sweeping until convergence of the quantity H(b, b) is reached.
Because we are only interested in the smallest eigenvalue E 0 and its eigenvector we can use the Lanczos iteration [110]. We conclude this appendix by discussing how to fix p min/max (n) and D q (n). As choosing finite values for these quantities means an effective truncation in the Schmidt decomposition (C.13) we need to look at the weight of the Schmidt values λ q,αq (n) over the sectors q corresponding to the eigenvalues of L(n) for any n with r L ≤ n ≤ r R . Assuming that the ground-state approximations for n < r L and n > r R are accurate (see appendix A subsection A 1), we don't have to care about the Schmidt values (C.14) at the boundary. In practice we start with a certain distribution of D q -values for each n, anticipating that the dominant eigenvalue sector of L(n) would shift from q = 0 at large n to q ≈ Q at the center. After a first full DMRG-optimization, the initial D q values are updated: increased in case that the minimal retained Schmidt value in the particular eigenvalue sector is larger than λ min = 10 −18 , decreased in case that the minimal retained Schmidt value is smaller. This is repeated a few times until all retained minimal Schmidt values are smaller or of the same order as λ min . As for the choice of r L and r R , we verify a posteriori that the inhomogeneous interval of the MPS (C.10) is taken to be large enough, by verifying the convergence of local observables at large distances to their value for the homogeneous ground state.
Let us give a specific example. In fig. 20 and fig. 21 we show some details of the simulation of the ground state for m/g = 0.25, Q = 5, x = 100, Lg = 10.1. In our setup with lattice spacing 1/g 2 √ x = 0.1/g this corresponds to a distance of 101 sites between the external antiquark with charge −gQ and the external quark with charge gQ. Specifically, we put the antiquark at site 151 and the quark at site 252. And we reserve 150 sites on the left of the antiquark and 150 sites on the right for the nonuniform part of our MPS ansatz. In total we thus have 151 + 101 + 150 = 402 tensors B n that need to be optimized. By looking at the 10-base logarithm of the expectation value of some local quantities with respect to the Schwinger vacuum, see fig. 20a, we observe that we take the range of the nonuniform part large enough: the errors by taking a finite range for the nonuniform part are of order 10 −6 . In fig. 20b we show the distribution of the minimum charge p min (n) and maximum charge p max (n) we used. For q < p min (n) and q > p max (n) we thus put D q (n) = 0. The p min (n) and p max (n) we took at the boundaries, i.e. n 1 and n 402 correspond to the p min and p max of the Schwinger vacuum, i.e. the vacuum without external charges, that we simulated in [57]. Between the boundaries and the external charges we anticipate the increasing electric field and raised p max (n) to 4 + Q = 9 anticipating the dominant eigenvalue sector p 0 ≈ Q at the center.
In figs. 20c and 20d we plot the distribution of the Schmidt values among the eigenvalues sector q of L(n) at the sites n = 150 (c) and n = 200 (d). As we explained above, we adapted the bond dimensions such that for each site n and at each eigenvalue sector p of L(n): min αq λ q,αq (n) 10 −18 . Comparing with figs. 12a and 12b we observe that the dominant eigenvalue sector is shifted to q = 2 for n = 150 and to q = 5 for n = 200. One can also see that our p min (n) and p max (n) are not entirely optimal: for certain charge sectors the largest Schmidt-value is still well below 10 −18 , and these sectors could have been discarded altogether. As we can see by looking at fig. 21a the most dominant eigenvalue sector of L(n), i.e. the eigenvalue sector q with the largest value for Dq(n) αq=1 λ q,αq (n) shifts from q = 0 to q = 5 as we go from the left boundary to the middle and then decreases to q = 0 as we go to the right boundary. We also show the maximum bond dimension max q D q (n) in fig. 21b. The largest bond dimension is required in the region where the electric background field is applied.  A ground-state solution will therefore be of the form Ψ(x A , x B ) = φ A (x A )φ B (x B ) where now φ A (x A ) and φ B (x B ) are both ground states of the nonrelativistic one-particle problem for a linear potential. All eigenstates for this nonrelativistic Hamiltonian H A (and similar for H B ) can be written in terms of the so-called Airy function Ai [96]: where N is the normalization factor and E n is the (kinetic) eigenenergy of the eigenstate. These energies follow from the continuity requirement on φ A and φ A at x A = −L/2, leading to either even or odd φ A under x A +L/2 → −(x A +L/2). The ground-state wave-function is even and the ground-state energy E 0 is related to the first zero of the first derivative of the Airy function, Ai (x 1 ) = 0, at x 1 ≈ −1.0188: E 0 = − x1 Notice that relativistic corrections to this approximation will necessarily involve quantum field contributions from other particle sectors. The relativistic one-particle Dirac equation has no bound state solutions for a linear (vector) potential [111,112]. Finally, in the nonrelativistic approximation we can then understand the transition from the string state to the broken-string state as a level crossing at the critical length L, where E string = E 2meson .