Schwinger pair production and string breaking in non-Abelian gauge theory from real-time lattice improved Hamiltonians

Far-from-equilibrium dynamics of SU(2) gauge theory with Wilson fermions is studied in 1+1 space-time dimensions using a real-time lattice approach. Lattice improved Hamiltonians are shown to be very efficient in simulating Schwinger pair creation and emergent phenomena such as plasma oscillations. As a consequence, significantly smaller lattices can be employed to approach continuum physics in the infinite-volume limit as compared to unimproved implementations. This allows us to compute also higher-order correlation functions including four fermion fields, which give unprecedented insights into the real-time dynamics of the fragmentation process of strings between fermions and antifermions.


I. INTRODUCTION
Confinement in quantum chromodynamics (QCD) manifests itself, amongst others, in that the energy stored in a gluon string between a quark and an antiquark rises linearly with the string length. A critical distance between the quark and the antiquark exists, at which the string starts to break through fermion-antifermion pair production, reducing the energy stored in the string [1][2][3][4]. The real-time dynamics of the nonperturbative string formation and fragmentation process, however, poses serious obstacles to a computation from first principles. In general, the nonequilibrium dynamics of quantum fields is not amenable to a formulation in Euclidean space-time and the use of importance sampling techniques, such that alternative approaches have to be employed.
Already the dynamics of geometric confinement and string breaking for gauge theories in 1 + 1 space-time dimensions is very rich with a hierarchy of timescales and a close link to the Schwinger pair production mechanism, sharing important aspects with their higher dimensional counterparts. In 1 + 1 dimensions, Abelian string dynamics and the phenomenon of multiple string breaking from supercritical field strength has been studied using classical-statistical reweighting techniques on the lattice [5]. For gauge theories with fermions these methods have also been employed to fermion production in 1 + 1 [6] and 3 + 1 space-time dimensions [7] using staggered fermions, as well as Wilson fermions [8][9][10][11][12]. For real-time Wilson fermions, lattice improvement techniques have been shown to be extremely useful in studies of quantum anomalies in QCD [13,14] and QED [15] in 3 + 1 dimensions.
Much of the recent interest in the low-dimensional dynamics comes from the prospect of studying important aspects of the gauge theory using quantum simulators [16][17][18][19][20][21][22], with a first proof-of-principle implementation on a trapped-ion quantum computer [23]. In this context, the dynamics of Abelian gauge theory has been studied in great detail using classical statistical reweighting techniques [12,24,25]. Important developments concerning string-breaking dynamics in Abelian and non-Abelian gauge theory have also been based on tensor network techniques in 1 + 1 dimensions [26][27][28][29] and recently on a Gaussian variational ansatz [30]. Additionally, the dynamics of flux strings has been studied in the U (1) gauge-Higgs model [31]. Compared to the Schwinger model, (1 + 1)-dimensional SU (N )-theory with massive fermions shares with full QCD the existence of confining strings between color charges as well as mesonic bound states in its spectrum [32][33][34][35]. Non-Abelian Schwinger pair production exhibits several distinct features such as interference and cancellations of currents between the different participating color charges [36][37][38].
In the present work we investigate for the first time the process of non-Abelian Schwinger pair production and the dynamics of string breaking in 1+1 dimensions using classical-statistical reweighting techniques. For Wilson fermions, we analyze the use of lattice improved Hamiltonians in a Kogut-Susskind-like real-time approach. We demonstrate that next-to-leading order (NLO) and nextto-next-to-leading order (NNLO) improvements in the fermion sector of the Hamiltonian can be efficiently employed to significantly reduce the lattice sizes required to approach the continuum and infinite-volume limit of the results. This seems crucial in view of the prospect of implementing these theories on a quantum simulator with a limited number of qubits. Moreover, this allows us to compute higher correlation functions involving four fermion fields. Together with the chromoelectric field strength, color charge and fermion number distributions, these charge-charge correlations provide a very detailed picture of the string dynamics and fragmentation process far from equilibrium.
This publication is structured as follows: In Sec. II the real-time lattice setup is developed, including field equations of motion and formulas for the computation of various observables from the lattice gauge theory. Subsequently, in Sec. III we focus on coherent chromoelectric fields that give rise to a homogeneous production of fermions. The available analytics provide a stringent precision test for the simulated lattice improved fermion number momentum spectra and fermion production rates. We explore thoroughly the occurring non-Abelian plasma oscillations, including an analysis of the effects of lattice improvements up to second order. In Sec. IV the breaking of color strings is investigated. Beginning with a description of corresponding initial conditions, we move on to a phenomenological analysis of the processes that take place in breaking color strings. Observing various dynamical color charge accumulations between the static external charges, we study charge-charge correlations among them and the propagation characteristics of such correlations. In Sec. V we summarize, draw conclusions and give an outlook.

II. REAL-TIME LATTICE GAUGE THEORY WITH WILSON FERMIONS
The Lagrangian density of SU (2) gauge theory with gauge potential A µ (t, x) and one fermion flavor ψ(t, x) transforming in the fundamental representation of the gauge group reads in continuous space-time: In (1 + 1) space-time dimensions there are no colormagnetic fields, and the field-strength tensor F µν is determined by the color-electric field E(t, x) with nontrivial components where the summation over gauge indices a = 1, 2, 3 is implied and T a denote the gauge group generators.
In what follows temporal-axial gauge is employed, i.e. A 0 (t, x) = 0, thus denoting A(t, x) ≡ A 1 (t, x). Since gluon self-couplings vanish in temporal-axial gauge in one spatial dimension, the gauge dynamics simplifies tremendously.

A. Lattice implementation with improved Hamiltonian
To formulate the theory on the lattice we discretize the spatial degree of freedom, such that positions may be labeled by integers n = x/a. The total number of lattice sites is N 1 , the lattice spacing a. The lattice has total length L = N 1 a. We employ spatially periodic boundary conditions. We neither discretize the time direction nor impose temporal periodicity assumptions. In momentum space we find the inverse (6) Spatial link variables are constructed as parallel transporting chromoelectric flux from site n to site n+1. Temporal link variables are identically unity in temporal-axial gauge. Inverse link variables are given by U † n+1 (t), parallel transporting chromoelectric flux from site n + 1 to site n. For fermion fields and link variables the following behavior under a given gauge transformation V n (t) ∈ SU (2) is observed: For later use we construct products of neighboring link variables with integer z > 0, Nontrivial path-ordering of gauge group elements applies here. The gauge sector Hamiltonian of our model reads To implement spinor fields on the lattice without doubler excitations we employ Wilson fermions [8,9,39]. Working in the Hamiltonian formulation of lattice gauge theory, the fermion sector Hamiltonian is given by Here r is the Wilson parameter with 0 < r, and the integer K denotes the order of the lattice improvement for the Hamiltonian: Using an appropriate choice of coefficients C k for k = 1, . . . K it is possible to cancel certain lattice artifacts of order O(a 2K−1 ) in the fermion lattice Hamiltonian [14]. Choosing C 1 = 1 and all other coefficients to vanish, one recovers the unimproved Wilson Hamiltonian, which is accurate to O(a) and which we call leading order (LO). Using C 1 = 4/3 and C 2 = −1/6 we obtain O(a 3 )-accuracy, labeled NLO. Including a third nonvanishing term, C 1 = 3/2, C 2 = −3/10, C 3 = 1/30, we obtain O(a 5 )-accuracy, labeled NNLO. The complete Hamiltonian including gauge and fermion fields reads B. Initial conditions and classical-statistical reweighting To solve the Cauchy problem of the time evolution of field degrees of freedom, we need to specify both fermion and gauge initial conditions at time t = 0. We assume that initially the two sectors decouple, subsequently quenching the system into a coupled state via time evolution with the full Hamiltonian (12). For the physics of fermion production from strong gauge fields or color charges, fermions are initialized as free fermions throughout this work. Details on the fermion initial conditions are given in Appendix A 1.
While the gauge field initial conditions are specified in more detail in Secs. III and IV, we consider strong fields for which the initial color-electric fields are of order of the critical field strength E c = m 2 /g. In this case, well-established classical-statistical reweighting techniques can be employed to replace the full quantum dynamics to very good accuracy by sampling classical gauge field dynamics [40][41][42][43][44]. Observables are then computed as ensemble averages of the results from a sufficiently large number of sampling runs until convergence is observed. For details on the sampling of gauge field quantum initial conditions we refer to Appendix A 2.
To this end, gauge degrees of freedom, i.e. U n and E n , are evolved classically from given gauge field configurations and fermion correlations. To resolve the manifest quantum nature of the fermions, the fermion equation of motion is solved on operator level, applying a modefunction expansion: with time-dependent mode functions φ u n,q,a (t), φ v n,q,a (t). The time-independent creation and annihilation operators b q,a and d † q,a , respectively, satisfy Fermion occupation numbers are determined by We construct the statistical propagator on the lattice, which reads in terms of mode functions, The temporal derivative of an operator O may be calculated using Key steps in the derivation of each of the following timeevolution equations are given in Appendix A 3. For the chromoelectric field one obtains the equation where the trace runs over color and Dirac indices. The time evolution of the link variable follows from with E n and U n acting upon each other via the standard matrix product. Finally, the fermion field operators obey By linear independency of the creation and annihilation operators b q,a , d † q,a the equation of motion for each and every mode function takes this form. To solve the given system of equations, we specify a time-step width a t and employ a four-step Runge-Kutta algorithm.

C. Abelianization for homogeneous fields
Diagonalizing gauge degrees of freedom in color space, which is generally possible by a local gauge transformation, provides a way to simplify the SU (2) gauge group structure. It results in a U (1)×U (1) gauge theory, which may be understood intuitively better compared to the SU (2) theory [36,37]. For later interpretation of lattice results in the context of Schwinger pair production, we consider in the following the Abelianization procedure for homogeneous field configurations.
Let E n (t) be homogeneous with no rotations in color space taking place. We write for all n ∈ Λ, with a constant vector n a , such that n a n a = 1. By Hermiticity of n a T a there exists a unitary matrix U , such that Explicitly, U may be given by Due to Eq. (23) one finds This may be interpreted as effectively decomposing the SU (2) gauge group into U (1)×U (1) with an Abelian coupling constant half as strong as the original non-Abelian one. By homogeneity and constancy of the transformation matrix, the given Abelianization procedure leaves equations of motion invariant.

D. Correlation functions
The time evolution of all fields involved allows the computation of a wide range of bosonic and fermionic observables at any simulation time. Inter alia, the color charge density ρ a n (t) and the fermion color current j a n (t) can be computed with the help of the statistical propagator (16) according to ρ a n (t) = j a n (t) = Here j a n (t) naturally incorporates lattice improvements due to the spatial derivative in the respective continuum expression.
From H (g) one infers the bosonic energy density at site n ∈ Λ, E a n (t) E a n (t) .
Though not indicated explicitly, correlation functions are computed as ensemble averages of the data calculated in individual runs.

E. Abelianized fermion numbers
For the interpretation of results, in particular in situations with homogeneous fields, it is often useful to define fermion pseudoparticle numbers from single-particle energy densities following the lines of Ref. [7]. To this end, gauge degrees of freedom are considered to be spatially uniform at all times, i.e. A a (t) ≡ A a n (t), E a (t) ≡ E a n (t). From the fermion Hamiltonian given in Eq. (11) we read off a Hamiltonian operator with spatial indices, n, m ∈ Λ: From this one can compute a fermion energy density as with a trace over color and Dirac indices. Using the diagonalization matrix U as given in Eq. (24) and E a (t) = n a E(t), A a n (t) ≡ n a A(t) = −n a t 0 dt E(t ), Here we employed We define From this we can compute a diagonalized phase-space energy density ± n,q (t): where φ u/v,± q,q denote the Fourier-transformed Abelianized mode functions, withq ≡ q−N 1 /2. We find an Abelianized single-particle energy density computed from physical masses m ± n,q (t) and physical momenta p ± n,q (t) given by For notational simplicity we introduced Using ω ± n,q (t) we may define an Abelianized fermion pseudoparticle number as follows: From this we can compute marginal distributions via [7] n ± n (t) ≡ and a total pseudoparticle number as [7] n ± (t) = 1 N 1 n,q n ± n,q (t).
Since n ± (t) represents an expectation value obtained from an ensemble average, it gives in general rise to noninteger particle numbers. In what follows we refer to this method of computing fermion numbers as Abelianized fermion numbers, indicating that the single-particle energy ω ± n,p (t) has a well-motivated interpretation if one can diagonalize color degrees of freedom.

F. Gauge-invariant fermion numbers
For comparison, and for the interpretation of the results from inhomogeneous initial conditions, we consider a second definition of fermion number using a Wigner function approach [6]. To this end, first a gauge-invariant statistical propagator is constructed, with a trace over color indices only. W n,m (t) is constructed as the lattice Wilson line along the shortest spatial path that connects the points n, m. It transforms under a gauge transformation Setting δ ≡ n − m and dropping the time argument in the notation, if δ ≥ 0 we find If δ < 0: We account for the periodicity of the lattice by taking the module operation in the second argument of where s, p, v 0 , v are all real. In the free case, (38a) and (38b) reduce to From the decomposition components we calculate various pseudodistributions: corresponding to charge and energy density, respectively. We define a phase-space resolved gauge-invariant fermion quasiparticle number: where n l,q (t), n l,q (t) label particle and antiparticle contributions, respectively. We will refer to this method of computing fermion numbers as gauge-invariant fermion numbers, justified by the included Wilson-line W n,m (t) to construct fermion numbers in a gauge-invariant fashion.

G. Connected charge-charge correlation function
In the study of color string dynamics a measure of connected color charge-color charge correlations is considered. For this purpose we introduce here the connected equal-time charge-charge correlation function In terms of mode functions one finds

H. Gauss's law
To simulate states in the physical Hilbert space, Gauss's law needs to be met. In our model it reads If the color charge ρ a n vanishes, we can iteratively solve this equation for arbitrary n > n ref ∈ Λ: For the G a n commuting with the model's Hamiltonian H, a state that initially fulfils Gauss's law does so for an arbitrary later point in time as well.

III. FERMION PRODUCTION AND IMPROVED HAMILTONIAN BENCHMARKS
In this section we investigate fermion production from homogeneous chromoelectric fields exceeding the critical field strength for Schwinger pair creation. In particular, we establish how improved Hamiltonians can be used to simulate this process and emergent phenomena such as plasma oscillations using significantly smaller lattices as compared to unimproved lattice implementations.

A. Fermion production: Benchmark at early times
To benchmark the lattice setup, we first compare our lattice simulation results with established analytic formulas for the non-Abelian Schwinger mechanism in a static coherent chromoelectric background field [45][46][47]. With this in mind we disregard the backaction of the fermion sector onto the gauge sector and exclude sampling. This approximate description is expected to be accurate at sufficiently early times, for which the analytic results are available. To be able to go to later times, we consider the fully nonlinear dynamics including backaction below.
We introduce the dimensionless field-strength parameter with E c = m 2 /g being the critical field strength. We fix g/m = 0.1 throughout this subsection and apply vacuum initial conditions for fermions, while the homogeneous chromoelectric field is set to the value stated. The simulations correspond to solving the fermion equation of motion, Eq. (21), with a sudden switch-on of the chromoelectric field at initial time t 0 = 0/m. In Fig. 1 the total number of fermions as a function of time for diagonalized color direction + and various background fields = ( 1 , 0, 0) is displayed. The calculations are based on NLO lattice improvements with lattice parameters as given in the figure caption. After a transient regime of enhanced fermion production at small times t tr 1/m, the curves show a linear regime in which fermion-antifermion pairs are produced at a constant rate.
Linearly fitting the data curves for n + (t)/m 2 at times t ≥ t tr , the rates together with results from respective analytical calculations are displayed in Fig. 2. As detailed in Appendix B, one expects a fermion production rate per unit length and per diagonalized color direction ofṅ The results from simulations are seen to agree with the analytics to very good accuracy. An exponential suppression of pair production below E c takes place both analytically and numerically, in accordance with prior studies in the framework of the Schwinger model and QED [6,7,36]. We now study the role of lattice improvements for the approach of the simulation results of the discretized theory in a finite volume to the analytic prediction for the continuum theory in an infinite volume. Figure 3 shows results for fermion production rates for various numbers of lattice sites N 1 in the unimproved LO formulation, which is then compared to NLO and NNLO improvements. Here, the total lattice size L = N 1 a is kept constant such that a is decreasing with increased N 1 to approach the continuum. For the unimproved lattice theory we find that getting close to continuum results for the production rate requires extremely small lattice spacings, in accordance with previous studies in QED 1+1 [6]. We observe that this situation changes dramatically, once improved Hamiltonians are employed. Figure 3 indicates that significantly smaller numbers of lattice sites N 1 lead already to results in the vicinity of the continuum theory. The NLO improved theory is seen to be extremely efficient with only minor differences to NNLO results, both converging very well at least for N 1 250 for the employed parameters. Already around N 1 100 the NLO (NNLO) improvement only slightly overestimates (underestimates) the continuum result, and smaller lattices may be employed if one is willing to accept errors exceeding the few-percent level. Relatively small lattices are crucial, e.g., for any realistic chance to implement this physics in quantum simulators in the not too distant future.
Apart from integrated quantities such as the total particle number, it is instructive to analyze also momentumresolved fermion numbers. In Fig. 4 the momentum spectrum of created fermions, n + p (t), is displayed at time t = 16/m for diagonalized color direction +. Switching to diagonalized color direction -is equivalent to reflecting the graph at p/m = 0. Shown are results without (LO) and with first-order (NLO) lattice improvements for a fixed number of lattice sites as given in the figure caption, together with the analytic result.
The fermion number distribution shows a clear peak around p/m = 0, reflecting the fact that most of the fermion-antifermion pairs are created at rest. Subsequently, they are accelerated in the applied chromoelectric field towards higher momenta. The low-momentum fermions are seen to be rather well described both at LO and NLO. At higher momenta, however, the accelerated fermions are significantly better described using the lattice improved Hamiltonian. These improvements be- come particularly visible for integrated quantities such as the total particle number, which sum over all momentum modes. We note that in analytical computations the initial time is sent to the remote past, t 0 → −∞, such that produced particles occupy arbitrarily high momenta. In contrast, the actual simulations start at t 0 = 0/m and produced particles only occupy finite momenta at finite times in the presence of a constant background field. This initial-time difference is also the reason for the transient regime of enhanced fermion production at small times t tr 1/m visible in the simulation data of Fig. 1, which is not present in the analytic estimates. Of course, only at sufficiently early times restricting to a constant background field is a valid approximation. Since total energy must be conserved, at later times the backaction of the produced fermion pairs on the applied chromoelectric field becomes relevant, which we address in the following.

B. Plasma oscillations: reaching longer times using lattice improved Hamiltonians
We now include the back-coupling of fermion currents onto the gauge sector, correspondingly taking the chromoelectric field equation of motion (19) into account. We keep g/m = 0.3 fixed throughout this subsection. Fig. 5A displays the Abelianized fermion number momentum spectrum n + p (t) + n − p (t) as a function of time. The initial acceleration of the produced fermions is visible along with a subsequent deceleration process, then an acceleration to lower maximum momenta than before and so on. To understand this oscillating behavior, it is helpful to consider also the time evolution of the chromoelectric field displayed in Fig. 5B. A fermion current is induced at times when the gauge field is strong, accompanied by a respective gauge field decay. Once the gauge field decayed fully, via the produced fermions' backaction onto the gauge sector a gauge field builds up again but pointing in the opposite direction. This process occurs again and again, resulting in non-Abelian plasma oscillations. In Fig. 5C different contributions to the total energy density are shown as a function of time. The gauge part reflects well the oscillating behavior of the chromoelectric field. The fermion energy part is seen to be large whenever the chromoelectric field part is small. Their sum stays constant, as required by energy conservation.
All results shown in Fig. 5 have been obtained at NNLO. We now address the role of lattice improvements for the nonlinear dynamics. We have seen that the backcoupling of fermions and gauge fields leads to plasma oscillations with decreasing maximum fermion momenta and multiple zero crossings at later times. In general, finite-volume effects and the associated limits on resolving low-momentum scales are expected to play an enhanced role at late times.
To illustrate the dependence of the results on the volume, we display in Fig. 6 the total fermion number n + (t) + n − (t) as a function of time. The computations are done for two different lattice sizes with N 1 = 512 and N 2 = 768 for fixed lattice spacing such that the volumes differ accordingly. Panel A displays respective results at NLO, for which one observes a practically volume-independent early-time behavior. But already after about the first plasma oscillation, indicated by a plateau in the function n + (t) + n − (t) which is due to a zero-crossing of the chromoelectric field around that time, deviations between the N 1 = 512 and N 2 = 768 setups become visible. This has to be confronted with the corresponding simulations at NNLO given in panel B. In this case, both lattice sizes give very similar results up to about the time of the fourth plateau, and even afterwards the larger lattice seems to give reasonable predictions. The ability to quantitatively describe longer time scales for a given lattice size is a very powerful property of the lattice improvements.
We end this section by noting that the qualitative behavior of the plasma oscillations in QCD 1+1 , including the growth in oscillation frequency with time, coincides well with previous classical-statistical studies in the framework of U (1) gauge theory [6,7]. However, from Fig. 5A a characteristic property of SU (2) theory plasma oscillations can be observed: Against the gauge field background half of the produced fermions are accelerated into positive direction, while simultaneously the other half is accelerated into negative spatial direction. The gauge field changing sign around its first 0, the produced fermions start being accelerated into the respectively opposite direction. The pattern of fermion acceleration closely mimics the oscillating gauge field dynamics at later times, too.
Resulting from temporary fermion production dropouts in the zero-momentum mode, substructures emerge in the time-evolving momentum spectrum of fermions from times t 80/m onwards. In fact, the pattern is a finite-volume lattice artefact, continuously vanishing with an increasing number of lattice sites. For corresponding details we refer to Appendix A 4.
For homogeneous initial conditions we can diagonalize color degrees of freedom as pointed out in Sec. II C. We define the Abelianized fermion current using the already established fermion current (27) and the diagonalization matrix U , given by Eq. (24), U j a n (t)T a U † .
In Fig. 7 the Abelianized fermion currents j ± (t) are shown. One being the negative of the other, they indicate propagation of the produced fermions into opposite spatial directions. The oscillating behavior is a manifestation of the occurring plasma oscillations. The two Abelianized fermion currents sum up to 0, in accordance with both theory and previous studies [36,37].

IV. STRING-BREAKING DYNAMICS AND HIGHER CORRELATION FUNCTIONS
After the benchmark tests for homogeneous fields in the previous section, we now apply the lattice improved Hamiltonian approach to string-breaking dynamics in QCD 1+1 . This involves computations of the nonequilibrium dynamics for inhomogeneous field configurations, constituting a particular strength of the approach. Moreover, we demonstrate that even higher-order correlation functions such as the charge-charge correlator (54) involving four fermion fields are accessible with these techniques.

A. Setup and initial conditions
To motivate our setup and initial conditions, we start by considering a confining gauge string of length d = la between two external color charges, specified classically by the color charge distribution whereas ρ 1 0,n = ρ 2 0,n = 0 due to N 1 0 = N 2 0 = 0. By Gauss's law (55) this results, classically, in a homogeneous chromoelectric field 3-component along the string, The string constructed in such a way from external color charges ±gN 3 0 has an energy content of If we demand that the initial classical chromoelectric field in the string's interior reads by value init = (0, 0, c) for a real constant c, the external color charges need to provide the energy This results in the requirement We apply this argumentation to a string with initial interior chromoelectric field ini = (0, 0, 2.0). Thus obtaining g = m for N 3 0 = 2.0, we keep g = m fixed throughout this section. Simulations to be described in the following demonstrate that to show critical stringbreaking behavior the string length needs to be around d we deduce that the critical SU (2)-string in QCD 1+1 corresponds to an Abelian string in 1 + 1 dimensions with effective external U (1)-charges [5,25] of N U (1) 0 = 2. These initial conditions are supplemented by vacuum fluctuations in the gauge sector. As in previous sections, for fermions we employ vacuum initial conditions. We discuss the impact of initial vacuum fluctuations on the gauge fields below and for details on their implementation we refer to Appendix A 2, again. The number of averaged runs is set to 3 throughout this section, whenever classical-statistical sampling for the generation of the initial vacuum "quantum-half" for gauge fields is applied.

B. String breaking and supercritical color strings
A string between two external charges that breaks completely, in particular with the center chromoelectric field asymptotically approaching 0, is called a critical string. A "supercritical" string is a string that shows multiple string breaking [5], such that the center chromoelectric field oscillates around 0 as it does in plasma oscillations.
Illustrating string-breaking dynamics by means of the gauge field, the color charge distribution and fermion numbers, in Fig. 8   2 0 = 0 all dynamics take place in the 3-components; 1and 2-components stay 0 at all times simulated. If vacuum gauge field fluctuations are taken into account for the present initial conditions, this remains approximately true. The gauge-invariant fermion number definition is employed. For a comparison with Abelianized fermion numbers we refer to Appendix A 5, the essence being that the two agree well with each other.
In row A of Fig. 8 variables are displayed for external charges ±2.0g. At early times, fermion-antifermion pairs are produced via the non-Abelian Schwinger mechanism, as can be seen in fermion numbers (panel A3). The dynamically created fermions and antifermions are initially produced on top of each other, resulting in the absence of color charge inside the string at early times (panel A2). They get separated with nearly the speed of light by the chromoelectric field, gradually screening the external color charges. First, the center chromoelectric field roughly decays linearly, later asymptotically approaching zero (panel A1). Indeed, we find that for a string of length d = 40.5/m complete string breaking occurs. Additionally, part of the produced particles gets accelerated towards the outside of the string and propagates freely with approximately the speed of light beyond the external charges at both string ends, continuously occupying the space that surrounds the string, as can be observed in all three variables (panels A1-A3). Once the string broke and the external charges got screened, the sole dynamical objects that remain -apart from small effects close to and inside the string -are these fermions and antifermions that at constant speed fly away from the string. This behavior matches closely the corresponding Schwinger model string-breaking behavior [5,25].
While in row A of Fig. 8 the color string breaks exactly once, in rows B to E the strings break multiple times. The external charges are set to ±3.0g, ±4.0g and ±8.0g, giving rise to supercritical color strings. Due to the stronger initial chromoelectric field the external charges get screened faster than in the N 3 0 = 2.0-case. The strings break faster. In all four cases, however, the center chromoelectric field does not asymptotically approach 0 within the simulated time interval, but instead shows an oscillating behavior, indicating the occurrence of plasma oscillations and multiple string breaking. In contrast to the homogeneous plasma oscillations described in Sec. III B, the plasma oscillations inside color strings happen inhomogeneously (panels B1-E1). By Gauss's law (55) chromoelectric field inhomogeneities are generated by the presence of color charges. Considering color charge distributions, we find a clear pattern of peaks and falls that oscillates in time (panels B2 to E2), getting more fragmented with increasing time the higher the external charge values get. Within this process, color charge accumulations occur that by sign differ from their immediate neighborhood. As in the critical case, part of the produced fermions and antifermions is accelerated towards the outside of the supercritical string, propagating approximately freely beyond the external charges at both string ends. But now the color charge outside the string varies strongly both spatially and temporally. Moreover, the space-and time-resolved fermion numbers closely follow the chromoelectric field dynamics. A large amount of fermion-antifermion pairs is produced wherever the gauge field is large (panels B3-E3). This is in accordance with the non-Abelian Schwinger mechanism.
We now address the role of vacuum fluctuations in the initial conditions of the gauge field configurations considered. In fact, the data shown in rows A to D are obtained without initializing the vacuum quantum-half for gauge field configurations. The initial vacuum fluctuations translate for each sample run into additional small inhomogeneities, which are averaged in the end. We check in the following whether this could change some of the smaller substructures observed, in particular in the N 3 0 = 8.0-case. The simulations for rows D and E of Fig. 8 are both generated with external charge values of ±8.0g, in row E including vacuum fluctuations in gauge field initial conditions for comparison. While qualitatively both D and E show very similar results, in the chromoelectric field and color charge distributions tiny spatiotemporal oscillations are present if fluctuations are taken into account (panels D1-E2). Fermion numbers are nearly not affected by these oscillations since fermion production by small chromoelectric fields is strongly suppressed in the non-Abelian Schwinger mechanism (panel E3). We expect the tiny oscillations to become smaller with an increasing number of samplings. Furthermore, the behavior of chromoelectric fields in the immediate vicinity of the external charges exhibits some differences: While without initial gauge field quantum fluctuations the chromoelectric field outside the string shows distinctive falls approximately at times 15/m, 33/m and 50/m and peaks at times 22/m, 39/m and 57/m inside the string (panel D1), the last two peaks and falls both inside and outside the string in the external charges' vicinity nearly disappear if fluctuations are taken into account (panel E1). Thus, taking into account the initial vacuum gauge field fluctuations by means of sampling can have an effect on the oscillating chromoelectric field "afterglow" around external charges sitting at color string ends. The color charge distribution resembles this effect (panels D2 and E2); so does the fermion number distribution, though smaller by value (panels D3 and E3).
Fundamentally, the external charges at supercritical color string ends never get screened fully: Always a remnant chromoelectric field inside the string is present. Time-evolving the given initial conditions, supercritical color strings remain confining in our model at all times simulated.

C. Higher correlation functions: charge-charge correlators
In this subsection we study charge-charge correlations involving four fermion fields. This is used to further analyze the color-charge accumulations in the interior of supercritical strings, whose existence is inferred in the previous section from the fermion charge distribution and gauge-invariant fermion numbers of Fig. 8, involving expectation values of two fermion fields. The quantity of interest is the connected anticommutator which is computed using Eq. (54). This quantity is not gauge invariant, but turns out to be very suitable for discussing some characteristic nonequilibrium phenomena in QCD 1+1 associated to four-fields correlations. For all results shown, initial vacuum fluctuations have been taken into account. In the following we consider C ab nm (t) ≡ C ab n (t, x)| x≡ma as a function of time t and position x for given reference site n or position x = n/(10m) for the parameters employed. In Fig. 9, C 33 n (t, x)-data are displayed both for reference site n = 768 and reference site n = 667. The external charges are N 3 0 = 8.0 at position x − 57/m and at x + 97/m, respectively, such that both reference sites are located between the charges.
One observes that charge-charge correlations spread in space as time is increasing, propagating beyond the external charges. The correlations are sharply peaked at equal points x/a = n. In a first temporal regime correlations spread with approximately twice the speed of light and less along the string and beyond, with a comparably large front peak value. The regime of propagation with twice the speed of light lasts until correlations have reached part of the boundary of the external charges' future light cones that lies outside the string. Subsequently, a second temporal regime begins, in which correlations propagate with a maximum front velocity of approximately the speed of light.
Since the correlation (66) represents an anticommutator (as opposed to commutator) expectation value of two bosonic composites, the propagation can exceed the speed of light without being in conflict with any fundamental principle. We have seen above in Fig. 8 that the color charges can move close to the speed of light and in opposite directions. If the charge-charge anticommutator expectation value is approximately a function of the spatial difference, then the spreading of correlations within the medium would exhibit approximately a maximum velocity of twice the speed of light as observed. In fact, once the light cone boundaries of the external charges are reached, since there is no medium outside the light cones, there is no relative motion possible between the fermions inside and outside the light cones and the maximum correlation spreading drops down to the speed of light. We emphasize that this observation is a genuine nonequilibrium phenomenon, which would not be possible in thermal equilibrium where equal-time correlators are constant and anticommutator and commutator expectation values are related by the fluctuation-dissipation relation [48].
Having globally studied connected ρρ-correlations, we now take a closer look at an example Cauchy surface of constant time with the aim of relating the correlation peaks and falls to the color charge accumulations appearing. In Fig. 10A  external charges, not the four additional ones. In Fig. 10B the connected ρρ-correlation function C 33 n (t, x) is displayed at time t = 16.8/m and reference site n = 667. For illustration purposes its negative values are shown. Primarily, we notice that the correlation function peaks match the color charge accumulations. Since the reference site n = 667 is lying inside the color charge accumulation around position x 3 , the pictured correlation function C 33 667 (t, x) can be interpreted to measure the charge-charge correlation at position x with the accumulation around position x 3 at time t = 16.8/m. We clearly find that the peak at x 3 is correlated with the color charge accumulations at positions x 1 , x 2 , . . . , x 5 , but not with the accumulation around x 6 at time t = 16.8/m yet. The reason for this can be seen in Fig. 9B: At the time displayed ρρ-correlations did not propagate yet to the vicinity of x 6 . Indeed, at the respective later times we find also a nonvanishing value for the correlation function at x 6 , such that the color charge accumulations arising inside and around a supercritical color string are correlated with each other. As the data in Sec. IV B already suggest, these multiparticle excitations become increasingly fragmentated with time increasing.

V. CONCLUSIONS
In the present study we investigated real-time fermion production via the non-Abelian Schwinger mechanism as well as the dynamical breaking of color strings between external static color charges in QCD 1+1 with gauge group SU (2) using classical-statistical reweighting techniques. Within this setting we confirmed available analytic results and demonstrated that lattice improvements up to second order can significantly improve convergence towards the continuum limit of nonvanishing fermion numbers in homogeneous configurations. We observed non-Abelian plasma oscillations upon including the backaction of created fermions onto the gauge sector. For SU (2) gauge theory plasma oscillations half of the produced fermions and antifermions is accelerated in positive spatial direction, the other half simultaneously in negative spatial direction. By means of improving the fermion current which feeds back to the chromoelectric field, the fermion number convergence behavior in non-Abelian plasma oscillations at late times can be optimized using higher-order lattice improvement terms, as well.
Additionally, we studied inhomogeneous initial conditions, focusing on configurations with a color string stretched between two external color charges. Being far from equilibrium, non-Abelian string breaking gives rise to a wide variety of involved dynamical processes, including fermion pair production, charge screening, plasma oscillations, fermions and antifermions continuously occupying the surrounding space around the string -all of which emerging from apparently simple initial conditions. In particular, we showed that within supercritical color strings dynamical color charge accumulations can arise. Employing a connected charge-charge correlation function we demonstrated that these charge accumulations are correlated, with correlations propagating initially with almost twice the speed of light.
These phenomena may become accessible in quantum simulators in the not too distant future, if efficient implementations become available. Improved lattice Hamiltonian formulations as those employed in this work can be a crucial ingredient, since significantly smaller lattices may be employed to obtain physical results. This becomes particularly important going beyond the one-dimensional benchmark case, when the real-time confinement dynamics is no longer of geometric origin but arises solely from the non-Abelian character of the gauge theory.
Initially, we implement free fermions with n u q,a (0) = n v q,a (0) = 0 for all q, a, using the free solution to the Dirac equation: the u-and v-eigenspinors reading ω q is computed from m q , p q using Eqs. (48) and (49) via

Bosonic quantum fluctuations
Here we describe the initial bosonic quantum fluctuations, added to the classical gauge string configuration investigated in Sec. IV.
We denote creation and annihilation operators of a gauge boson with lattice momentumq = −N 1 /2, . . . , N 1 /2 − 1 and color index a by a a, † q , a a q , respectively. The quantum fluctuations to construct correspond to the quantum-half contribution to bosonic vacuum occupation numbers, thus requiring In order to stay well below the lattice UV cutoff, we introduced a finite momentum scale Q up to which fluctuations are taken into account, merely. We carefully checked for insensitivity of the obtained results to the precise choice of Q and consistently specified Q = π/(3a).
The quantum dynamics of a a, † q , a a q are in the classicalstatistical approximation imposed by sampling complex numbers α a, * q , α a q from a given distribution function and computing expectation values of a a, † q , a a q according to the latter. As a distribution function W(α a q , α a, † q ) we chose the standard Gaussian distribution,ᾱ a q denoting its average value and σ a q its width, The imposed one-and two-point functions for quantumhalf fluctuations then read a a q = dα a q dα a, * q W α a q , α a, † q α a q = 0, (A6a) 1 2 {a a, † q , a a q } = dα a q dα a, * q W α a q , α a, † q α a, * q α a resulting inᾱ a q = 0 and σ a q = Θ(Q − 2πq/L)Θ(Q + 2πq/L)/ √ 2. α a, * q and α a q being drawn randomly according to W(α a q , α a, † q ), we compute their Fouriertransformed spatial counterparts as Fluctuations of the gauge potential and the chromoelectric field, added to the classical initial conditions, are computed from this as δA a n = √ 2a 2 Re(α a n ), δE a n = √ 2 Im(α a n ).
Since link variables are classically initialized to unity in this work, upon inclusion of classical-statistical sampling they initially read U n (t = 0) = exp iga δA a n T a .
Observables are computed at each time step as ensemble averages of a given number of simulated runs with different fluctuating initial conditions generated this way.

Deriving equations of motion
Crux of deriving equations of motion for the involved field degrees of freedom are commutation and anticommutation relations between the fields to finally employ The chromoelectric field equation of motion follows from using A a n , A b m = E a n , E b m = 0 , in order to obtain by application of Baker-Campbell-Hausdorff E a n , U m = −g δ n,m T a U n , (A11a) E a n , U † m = +g δ n,m U † n T a .
Clearly, we find [E a n , H (g) ] = 0 for all n, a. It remains to compute [E a n , H (f) ], for which we make extensive use of relations (A11a) and (A11b). Expressions such as To obtain equations of motion for the fermion mode functions we evaluate [ψ n , H (f) ]. We make use of the algebraic identity where Dirac and color indices have explicitly been restored. Using this identity, straightforwardly the equation of motion (21) follows.

Volume dependence of substructures in plasma oscillations
In Sec. III B we observed that substructures occur in momentum space-resolved fermion numbers as a result of temporary fermion production dropouts in the zeromomentum mode. In Fig. 11 Abelianized fermion numbers n + p (t) are displayed again, now for both N 1 = 512 and N 1 = 768. We primarily notice that structures are more present at smaller lattice size. Related to finite volume and being smaller at larger lattice size, we take this as a hint that the occurring fermion production dropouts are, essentially, infrared lattice artifacts, disappearing continuously with increasing lattice sizes.

Comparing Abelianized and gauge-invariant fermion numbers
A proof of consistency within our work is to compare both the Abelianized and the gauge-invariant fermion number definitions. In Fig. 12 we find the two fermion number definitions roughly match each other, comparing fermion number momentum spectra for the two in a constant chromoelectric background field. Merely, the gauge-invariant fermion numbers oscillate around the Abelianized ones with amplitudes varying in time. At times, oscillations are larger than in Fig. 12, at times smaller, even nearly vanishing.
Enabling backcoupling of the created fermions onto the gauge sector, we find that both fermion number definitions show the same oscillating pattern. Fermion number momentum spectra and total fermion numbers computed from both definitions agree with each other, though not displayed here.
Furthermore, comparing to the gauge-invariant fermion numbers in Hebenstreit et al. [6], we note that fluctuations in the momentum spectra computed from our simulations are sometimes much larger than those encountered for example in Fig. 3 of [6]. We believe the reason for this to be the different type of lattice fermions employed: While we use Wilson fermions in the present work, Hebenstreit et al. implemented low-cost fermions. We expect that if in our Wilsonian approach we employed classical-statistical sampling also in the homogeneous setting around a coherent initial chromoelectric field, with increasing the number of samples fluctuations in gaugeinvariant fermion number momentum spectra would continuously vanish.  Besides the original computation by Schwinger using the one-loop effective action [49], it is possible to solve the Dirac equation in homogeneous, constant background U (1) gauge fields using quantum kinetic theory [47]. Analogously, we can proceed in SU (2) gauge theory upon diagonalizing color degrees of freedom [46]. Most of the derivation proceeds similarly to [45,47], reducing the number of space-time dimensions to 2. We solely sketch differences here.
We define an inner product for arbitrary mode func-tions φ, χ as with color indices a, b restored. In what follows "in" specifies a free fermionic initial state at time t → −∞; "out" labels the asymptotically free but gauge-rotated fermionic final state at time t → ∞. At t → −∞ we expand the Dirac field color components ψ a in momentum modes, ψ a (x) = dp 2π φ in,+ p,a (x) b in p,a + ψ in,− p,a (x) d in, † p,a . (B2) We expand ψ a as well in out-state momentum modes, ψ a (x) = dp 2π φ out,+ p,a (x)b out p,a +ψ out,− p,a (x)d out, † p,a , (B3) and find the Bogoliubov transformation coefficient in order to arrive with β(p) ≡ (β(p, 1), β(p, 2)) at for the in-vacuum expectation value of the out-particle number density operator. Here, V denotes configuration space volume, finally taken to infinity. We note that if it is possible to Abelianize color degrees of freedom, F(p) is left invariant upon the diagonalization scheme.
Restricting to a uniform, constant chromoelectric field E a n (t) = En a with n a n a = 1 for all n ∈ Λ, we can diagonalize color degrees of freedom as in Sec. II C. We find a diagonalized chromoelectric field diag( /2, − /2) in units of E c and again label diagonalized color components by ±, corresponding to effective electric fields ±E/2. Having applied the diagonalization procedure, we proceed as in U (1) gauge theory [47] to finally find withp ± ≡ ± 4 p ± gEt/2 m , and η ± ≡ ± 2 .
We denote the individual addends ± in Eq. (B6) by F ± (p). F ± (p) approaches a nonvanishing constant value for large p, which translates, approximately, into a constant rate per spatial volume L and per Abelianized color direction ± at which fermion-antifermion pairs are created [50,51], The physically measurable, total production rate of fermions and antifermions per volume L is thus given byṅ