Density of States of the lattice Schwinger model

Using a recently introduced tensor network method, we study the density of states of the lattice Schwinger model, a standard testbench for lattice gauge theory numerical techniques, but also the object of recent experimental quantum simulations. We identify regimes of parameters where the spectrum appears to be symmetric and displays the expected continuum properties even for finite lattice spacing and number of sites. However, we find that for moderate system sizes and lattice spacing of $ga\sim O(1)$, the spectral density can exhibit very different properties with a highly asymmetric form. We also explore how the method can be exploited to extract thermodynamic quantities.


I. INTRODUCTION
Interacting quantum many-body systems represent a challenge for analytical and numerical methods, yet they are key to the understanding of many fundamental physical phenomena. For this reason, and because most of the interesting systems are not exactly solvable, considerable effort is devoted to the development of very different numerical techniques to address these problems. Prominent examples include Monte Carlo algorithms [1] and tensor network (TN) methods [2][3][4].
One of the most significant properties of a quantum many-body problem is the density of states (DOS). Knowing the distribution of energy eigenstates gives access to the partition function and consequently to all the thermodynamic properties of the system. In the context of lattice gauge theories (LGT), approximating the DOS has been proposed [5] as a method to overcome the sign problem [6] that appears, for instance, in the presence of a finite chemical potential. But computing the DOS of a many body problem in general, or a LGT in particular, is a difficult task.
Several approximate methods have been devised to address this question. For classical statistical models broad histogram methods exist, as the Wang-Landau algorithm [7,8], and they have inspired variants that can be used for some quantum problems [9,10]. A different technique is the kernel polynomial method (KPM) [11]. Based on a Chebyshev expansion of the Dirac delta function, it is widely used in the single-particle quantum scenario. Only recently the application of tensor network algorithms to this problem has been explored [12,13]. The advantage of TN methods is that, in principle, they are free from the sign problem, and thus suitable for scenarios out of reach for Monte Carlo based algorithms. The method introduced in [13] precisely adapts the KPM method to an interacting scenario where the Hamiltonian can be expressed as a matrix product operator (MPO) [14][15][16] and can be used to estimate DOS, microcanonical averages and other spectral properties, but also to probe thermalization.
Specific methods have been proposed in the particular context of LGT. They include the Linear Logarithmic Relaxation method (LLR) [5,17], based on the Wang-Landau algorithm, which was introduced to im-prove the precision in the calculation of the DOS, with the specific aim of determining observables at finite density lattice field theories, and the Functional Fit Approach (FFA) [18][19][20] and have been used in a variety of models, such as compact U (1), SU (2) and SU (3) LGT [5], SU (2) gauge theory at finite densities with heavy quarks [21] and the Z 3 spin model at finite density, where the sign problem is present [22]. These methods rely on Monte Carlo techniques and therefore present challenges when applied to systems with dynamical fermions [23,24].
In this work, we explore the performance of the TN method introduced in [13] in a U (1) LGT. More concretely, we consider the lattice Schwinger model, which includes gauge and dynamical fermionic degrees of freedom. The features of the Schwinger model, one of the simplest LGT that nonetheless exhibits non-trivial phenomena common to more complex theories, have established it as a usual testbench for lattice techniques. It has also been the object of the first experimental quantum simulation of a LGT using trapped ions [25], and efficient representations of the model have been suggested for its study in a quantum computer [26]. More recently also a closely related model with U (1) gauge symmetry has been experimentally realized with ultracold atoms [27]. Many properties of the Schwinger model, including its mass spectrum, thermal equilibrium properties and dynamics, have been systematically investigated using TN approximations for the relevant states [28,29]. These studies have demonstrated the suitability of TN to describe the physically relevant states in LGT scenarios, and to attain precise continuum extrapolations. However, the standard techniques used in such studies do not provide easy access to high energy eigenstates, or to the spectral properties at energy densities far from the edges of the spectrum.
Using the above mentioned technique, allows us to directly study the DOS of the lattice model. In particular, we examine the dependence of the DOS with system size and lattice spacing, parameters that need to be varied in the (classical or quantum) simulation in order to approach the continuum limit. We find that the shape of the DOS changes dramatically between a lattice-dominated regime, where it is highly asymmetrical and exhibits sharp features, and a continuum-like regime, in which it resembles a Gaussian distribution. The last observation is in accor-dance with the results of [30], to our knowledge the only previous study of the DOS of the model. As a second goal, we explore the potential of the DOS approximation to calculate different observables in the canonical ensemble, in the spirit of the LLR method, and what are the limitations of this method in comparison to directly approximating thermal states with TN.
The rest of the paper is organised as follows: In Sect. II we briefly introduce Schwinger model on the lattice with the Kogut-Susskind formulation of staggered fermions, that is then mapped to a system of long range spin-spin interactions. We then proceed with presenting the methods used in Sect. III and the results follow in Sect. IV. We close with a discussion in Sect. V.

II. MODEL
In this paper we focus on QED in two space-time dimensions, known as the Schwinger model, and, more concretely, on its discrete lattice version. The Schwinger model is one of the simplest gauge theories, yet shares some of the most interesting features of quantum chromodynamics (QCD) [31,32], such as confinement and a broken chiral symmetry. This justifies its role as a standard testbench for LGT techniques. The model has also constituted a natural first target to benchmark the performance of TNS techniques for LGT. This enterprise started with a successful application of the density matrix renormalization group (DMRG) algorithm in [33], and has led in more recent years to the systematic exploration of spectral properties, thermal equilibrium or dynamics of the lattice model and its continuum extrapolation (see e.g. [28] for a review). Also the quantum simulation of the model has been proposed for various experimental platforms [29].
Specifically, the Kogut-Susskind staggered fermion formulation [34] of the lattice Schwinger model reads: with a the lattice spacing, m the fermion mass and g the coupling constant. Operators Φ n (Φ † n ) annihilate (create) the (single component) fermion mode on each lattice site n, and satisfy canonical anticommutation relations {Φ † n , Φ m } = δ nm and {Φ n , Φ m } = 0. Gauge degrees of freedom residing on the link between sites n and n + 1 are represented by canonically conjugate operators θ n and L n , which satisfy [θ n , L m ] = iδ nm and, in the continuum limit, respectively correspond to the vector potential and the electric field. Additionally, physical states need to satisfy the discrete version of Gauss law [35]. The Hamiltonian of Eq. (1) can be mapped to a spin model via the Jordan-Wigner transformation Φ n = k<n (iσ z k )σ − n with σ ± = 1 2 (σ x ± iσ y ), where σ α , for α = x, y, z, are the Pauli matrices. Additionally, for a system with open boundary conditions, it is possible to explicitly solve Gauss law, which results in a spin chain with long-range interactions [35,36]. Usually, the Hamiltonian is multiplied by a factor 2 g 2 a to result in an adimensional operator which, for a chain of N sites, reads where the relevant (adimensional) parameters are now x = 1 g 2 a 2 and µ = 2m g 2 a , while l represents the background field. Written in this form the model is suitable to be studied numerically with TNS methods.

III. METHODS
In order to compute the density of states of the model (2) we use the TNS techniques introduced in [13], which we summarize here for completeness.
Given a Hamiltonian H and an operator O we can define a generalized density of states function where the sum runs over all energy eigenstates |k of the Hamiltonian, H |k = E k |k . Notice that, up to a normalization factor, the usual density of states corresponds to the function for the identity operator g(E; 1). The functions (3) can be approximated by a finite sum of M Chebyshev polynomials (see appendix A), based in the corresponding approximation of the Dirac delta function δ(x) ≈ δ M (x). Explicitly, with the coefficients γ M n (explicitly shown in appendix A) corresponding to the Jackson kernel in the kernel polynomial method [11].
Specifically we define and by substituting Eq. (4) we get where we have defined the rescaled and shifted energỹ E = αE + η, and correspondingly the HamiltonianH = αH + η, such that the spectrum lies in the interval [−1, 1], and we identify the moments [13] µ n (O) := α tr T n (H)O .
As described in [13] (see also [37][38][39]), if the Hamiltonian is written as a matrix product operator (MPO) [15,16,40], one can also construct MPO approximations to its Chebyshev polynomials, starting from the exact T 0 = 1 and T 1 =H, and sequentially applying the recurrence relation T n+2 (H) = 2HT n+1 (H) − T n (H). The bond dimension of the resulting polynomial increases with n, and truncating it to a fixed value D produces a MPO approximation T (D) n (H). Using the latter, we can estimate the moments in Eq. (7) for operators of interest, and thus approximate the desired functions (see [13] for details).
The above method allows access to thermodynamic observables. The partition function in the canonical ensemble at inverse temperature β can be obtained from the density of states, and, consequently, different thermodynamic quantities can be computed, such as the energy the specific heat or the entropy with F (β) = − 1 β ln Z(β) the free energy. Also for a general observable O, the expectation value in the canonical ensemble can be expressed as Using the expansion and MPO approximations g M (E; O) described above, and performing the one-dimensional integration with the Boltzmann factor e −βE , we can obtain approximations for the partition function Z M (β) and the observables O M (β). For the partition function, in particular, withẼ min andẼ max the rescaled estimates of the edges of the spectrum. The proportionality sign accounts for a factor e βη α /α, which is not explicitly written in the expression above.
We can approximate the functions as an alternative series using the Chebyshev expansion of the exponential [41] where I n (b) are the modified Bessel functions of order n. Substituting the Boltzmann factor by this expansion in Eq. (8) or the numerator of (12), and using the orthogonality relation of the Chebyshev polynomials (shown in Appendix A), we can analytically integrate each term of the sum and express the result as a series of the modified Bessel functions I n (−β/α). Specifically, where N n = (1+δ n0 )π/2 are the norms of the polynomials.
In the same way, we get for the approximation of the numerator of Eq. (12)

IV. RESULTS
The Hamiltonian (2) can be represented exactly as a matrix product operator with small bond dimension (D = 5) [42]. This allows us to employ the method described in Sec. III to approximate the density of states of the model. To the best of our knowledge, the only previous estimate of this quantity was performed in [30], using discrete spectra obtained by the method of discretized light-cone quantization (DLCQ) [43][44][45]. The calculation was limited to finite sizes, but the authors observed that the form of the DOS, after an extrapolation to the continuum, appeared to roughly represent a Gaussian.
It is interesting to notice that, while for local Hamiltonians it is not surprising that the DOS resembles a Gaussian, because it weakly converges to this shape as the system size increases [46,47], this is not necessarily the case for the model in Eq. (2), as it contains long-range interactions. Here we explore in more detail how the DOS of the lattice model changes with the system size and the lattice spacing, and try to understand how the approximately Gaussian behavior arises as the parameters approach a regime close to the continuum. We also evaluate thermodynamic observables like the energy, the entropy, the chiral condensate and the specific heat. In this work we consider exclusively the case of massless fermions µ = 0 and zero background field l = 0. We calculate the DOS of systems with system sizes between 20 and 60 sites and lattice spacing corresponding to x is 1 ≤ x ≤ 300 (see table I for the precise parameters used for each size). In the following, the DOS is normalised to one, dEg(E; 1) = 1, which we denote asg(E; 1).  We observe that the lattice effects dramatically affect the shape of the DOS, as shown in Fig. 1. For small lattice spacing (corresponding to large values of x), the distribution looks similar to a Gaussian, a behaviour which was suggested in [30] for the continuum limit [48]. For a fixed system size N , as we decrease the value of x the peak becomes sharper in accordance with the behaviour suggested in [30]. For a fixed system size N , as we decrease the value of x the peak becomes sharper and closer to the lower edge of the spectrum, and the DOS takes a very asymmetric form, with a fast increment of the density close to the lowest energy and a much slower decrease after the peak. We further appreciate that in the small x regime (see left panels of Fig. 1 for x = 5), this heavier right tail of the distribution is approximately exponential.
In order to investigate quantitatively the deformation of the DOS, we first notice that to the left of the peak the shape of the probability distribution is still close to Gaussian. If the position of the peak is E p , we can then find the Gaussian distribution that best describes the DOS for , i.e. the Gaussian distribution with the peak at the same position and with variance σ 2 fixed by the half-width at half maximum. Specifically, σ = (E p − E 1/2 )/ √ 2 ln 2, where E 1/2 is the value of the energy to the left of the peak satisfyingg(E 1/2 ; 1) =g(E p ; 1)/2. Finally, the normalisation constant N ensures that the maximum value matches the peak of the DOS. The resulting Gaussian function is shown as a purple dashed line for each of the cases illustrated in Fig. 1. For small x, the right tail of the DOS decays much slower than the right tail of this Gaussian fit. Instead, an exponential decay provides a much better fit (shown in Fig. 1 with black dashed lines).
We can quantify the asymmetry of the distribution by computing the difference between the areas under the DOS curve to the right and to the left of the peak, This difference approaches zero as the DOS becomes closer to a symmetric distribution. The actual values of this asymmetry for a range of system sizes N are shown in Fig. 2 as a function of x. We observe that for all system sizes, the asymmetry vanishes as x grows (i.e. towards small values of the lattice spacing), more slowly for larger system sizes. To further analyze this, we fit the last few data points for each system size, in Fig. 2, to an exponential decay A = A 0 e −λx , where the parameters A 0 and λ depend on N . The decay constant λ decreases as the system size grows (see inset of Fig. 2), consistent with a faster approach to a symmetric shape for smaller systems.
As mentioned above, to the left of the peak, the distribution is very close to a Gaussian form. To quantify this observation, we compute the difference between the area of the Gaussian G(E) and the area of the DOS, to the left of the peak, As we can see from the lower panel of Fig. 2, for all system sizes and values of x this value is very small, indicating that, as x increases, the DOS approaches a symmetric distribution that resembles a Gaussian. We note here that in the case of approaching the continuum limit, with x → ∞, the Hamiltonian of Eq. (2) reduces to the exactly solvable XY model multiplied by a factor. In that case, there exists a unitary transformation U such that which implies a symmetric spectrum. This is not necessarily the case when x is finite and the long-range interactions are present.

B. Physical observables
As argued above, the DOS determines the values of thermodynamic properties, because it gives us access to the partition function and all thermal observables in the canonical ensemble through Eq. (12). We can thus investigate how the features of the distribution observed in the previous subsection affect different thermal properties. Specifically, we compute the energy E(β) and the entropy S(β), as well as the value in thermal equilibrium (denoted Σ(β)) of the chiral condensate operator which corresponds, in the continuum limit, to the order parameter of chiral symmetry breaking. In the massless case, the temperature dependence of this parameter in the continuum limit has been solved analytically [50], while for massive fermions, it has been studied numerically using MPS techniques [49,51,52]. As described in section III, the partition function and thermal observables can be directly approximated as finite sums of modified Bessel functions I n (−β/α), where α is the Hamiltonian rescaling factor, proportional to its operator norm. However, because the norm of the Hamiltonian grows with the system size (scaling as fast as N 3 for fixed small x), the argument of the functions involved, and thus the magnitude of the functions, grows fast for fixed β, leading to numerical instabilities, except for very small values of β. [53] Thus, we find it more convenient to evaluate numerically the integrals of Eq. (8) and of the numerator and denominator of Eq. (12).
In Fig. 3 we plot the energy, the entropy and the condensate obtained with this method for system sizes N = 20 and N = 30 and x = 5. It is also possible to compute the observables by finding a MPO approximation to the Gibbs ensemble of the model [2,14,16]. For the Schwinger model this MPO produces very precise results for thermal observables [49,51,52], and we use it here as reference [54]. While we observe that the results obtained with the method of section III (blue line) agree well with the reference (grey line) for the energy, the figure shows that, specially for the condensate, the agreement only holds for small values of β. We attribute the limitations of our calculation to the fact that for larger values of β, the Boltzmann factor e −βE in the integrals enhances the contributions of the lower edge of the spectrum, where the DOS is many orders of magnitude smaller than at the peak and the truncated approximation is less accurate. For more detailed explanation see Sec. IV C.
In [30] it was suggested that the shape of the DOS gave rise to a peak in the specific heat Eq. (10) which appeared to diverge in the continuum limit. These conclusions were contradicted by later, more precise studies [55,56], and attributed to the limited precision of the simulations which would not allow for a reliable continuum extrapolation. It is interesting to investigate whether limiting the study to finite systems and lattice spacing could show similar spurious signals in the specific heat. To this end, and because the precision obtained with our method is limited for large systems, we examine the specific heat using a standard MPS simulation [49,51,52]. We fix a system size N = 80 and systematically increase the parameter x.  17), as a function of the parameter x, for various system sizes (upper panel). As we approach bigger values of x, we notice that the form of the DOS becomes symmetric. For each system size N , the decay of the asymmetry can be fitted to an exponential A ∝ e −λx , with a faster decay for smaller systems. The resulting parameter λ is shown in the inset as a function of the system size. The lower panel shows the deviation with respect to a Gaussian in the area of the DOS to the left of the peak, as defined in Eq. (18).
Notice that in this way we are exploring properties of the (finite) lattice model, since this does not correspond to approaching the continuum (which would require increasing the system size as N ∝ √ x to maintain a consistent physical volume). Figure 4 shows the dependence of the specific heat with the inverse temperature as x increases. We observe a smooth peak in the specific heat, that drifts slowly towards infinite temperature as x increases (as mentioned before, in the limit x → ∞ the Hamiltonian reduces to the exactly solvable XY model with a multiplication factor x). In order to connect the observed quantities to their continuum correspondence, it is important to highlight here that β used in figures 3 and 4 is adimensional since the Hamiltonian of Eq. (2) is also adimensional. It is related to the physical inverse temperature β phys (the inverse temperature in the continuum) in the following way [49] C. Convergence and error analysis The calculated quantities, g M (E; O) and g M (E; 1), have two sources of errors, as discussed in [13]. One type of error comes from truncating the Chebyshev expansion up to order M . As shown in [11], when the Jackson kernel is used, and the function being approximated is continu- FIG. 4. The specific heat for a system size with N = 80, for different values of x. As the lattice spacing decreases and x increases, the peak of the specific heat becomes sharper. The inset shows how the rescaled inverse temperature, βmaxx, with βmax the inverse temperature that corresponds to the peak, behaves as we increase the parameter x. The grey line indicates the exact value for the XY model, to which the Hamiltonian converges in the limit x → ∞. The specific heat was computed using standard MPS methods.
ous in the interval [−1, 1], the errors are of order O(1/M ) and consequently are reduced as we increase the Chebyshev cut-off M . The second source is the truncation of the MPOs used to represent the polynomials to a finite bond dimension T D n (H), with D ≤ D max . The observation in [13] indicates that the bond dimension required to approximate the n−th polynomial with fixed precision grows fast (even exponentially) with the order n.
The effect of these errors is that of effectively limiting the energy resolution attainable with a given bond dimension, since the latter limits the number of moments that can be reliably extracted. When applied to the DOS and related functions, written as sums of Dirac delta terms, we expect that the effect of this limited resolution becomes more significant wherever the DOS exhibits sharp features or where its absolute value is small. In our case, this corresponds to the edges of the spectrum, where eigenstates become more sparse, while the effect is minor where the DOS is sufficiently large, such as near the peak.
In the computation of the spectral quantities g M (E; O) and g M (E; 1), only traces of the polynomials or (of their product with a given operator) appear. These appear more converged in bond dimension than the global error analyzed in [13], which allows us to explore relatively large values of the Chebyshev truncation parameter M . In particular, we found that characterizing the overall shape of the DOS, as discussed in IV A, was possible with moderate computational resources. Concretely, for N = 60, we found M = 500 to be sufficient for large values of x, 100 ≤ x ≤ 300. For smaller x, the peak of the DOS is sharper and requires a larger number of polynomials, growing to M = 1000 for 30 ≤ x ≤ 50 and up to M = 2000 for x = 5, as shown in figure 5. For all these cases, we found bond dimension D = 400 to be enough to accurately capture the shape of the DOS (see right panels in figure 5), and to have convergence of the results shown in figures 1 and 2.
For the same parameters, however, the edges of the spectrum are not converged, as can be appreciated in Fig. 5. Although the same moments are used to compute the DOS at all energies, the final result is the sum of a large number of terms with energy dependent coefficients that can oscillate with large frequencies (increasing with the order). Near the peak the sum is large, and the relative contribution of large order moments (more affected by truncation) is small. However, close to the edges precise cancellations of terms are required to recover the much smaller value of the DOS. Hence, in such areas, the result is much more sensitive to truncation errors in the moments of all orders. In practice, the computed DOS exhibits large fluctuations that depend on the truncation parameters in the regions where the true magnitude of the spectral density is small, up to the spectral edges. Because of the small magnitude of the DOS, these fluctuations do not appear to be present in Fig. 1, but they are visible when we plot the DOS in a logarithmic scale in Fig. 5.
Very close to the edges, we may expect that the discreteness of the energy spectrum for finite systems further contributes to the fluctuations. The typical energy gaps are exponentially small in the system size, but close to the edges, the gaps can be much larger. In order to be able to capture the sparse spectrum in these regions with our method, M would need to take much larger values, which in turn would require even bigger bond dimension D.
The discussion above refers to the determination of the DOS. The thermodynamic observables shown in Fig. 3 are computed using the same approximation to the Chebyshev moments, but combined with different coefficients, to compute Eq. (12). The truncation parameters that allow us to determine the form of the DOS are not enough to explore the whole range of temperatures. Instead, even if the errors at the edges do not affect the form of the DOS, they do constrain the interval of temperatures that is accessible when calculating thermodynamic quantities with our method (see Fig. 3). Specifically, as β increases, because of the Boltzmann factor, so does the contribution from the energies near the lowest edge of the spectrum E min to the quantities g(E; 1) and g(E; O). Because the computation with fixed truncation parameters is less accurate near the edges, attaining a precise calculation of the different thermodynamic quantities becomes more computationally expensive as β increases. The truncation errors related to the bond dimension D at the edges of the spectrum appear to be more significant as we increase the system size N , consistent with a narrower spectral density with much smaller relative magnitude of the DOS in the tails of the distribution, as compared to the peak, and smaller gaps that require better energy resolution. We do not carry out a detailed error analysis for the observables, but instead show in Fig. 3 the best results obtained, as compared to the exact results. These correspond to the biggest value of M that is accessible and a fixed bond di- mension D = 500. We notice that, with the given bond dimension, our results seem to improve significantly as we increase M , for most of the observables of interest.

V. DISCUSSION
In this work, we have applied the recently introduced method [13] to characterize the density of states of the lattice Schwinger model.This technique gives access to the overall shape of the DOS over the full range of energies and can be used to derive thermodynamic quantities that are expressed as integrals of the DOS with Boltzmann factors.
We have shown how the shape of the distribution for the lattice Schwinger model changes with the system size and the lattice spacing parameter, from a very asymmetric spectrum at large spacings, to an approximately Gaussian form for sufficiently small ones. It is worth noticing that, for any finite system size, there appear to be values of the parameters for which the DOS strongly deviates from a Gaussian behavior. Furthermore, the density of states for large lattice spacing exhibits a sharper feature at low energies, with a very fast increase, and a narrow peak, decaying much more slowly towards high energies.
Our observations on the spectral properties may be relevant for the interpretation of experimental realizations of the model, which has already been simulated in some pioneering quantum simulations. Since experiments are necessarily limited to finite (for the moment relatively small) system sizes, it is interesting to understand which ranges of parameters ensure a behavior close to the continuum limit for the spectrum, and the thermal properties. Our analysis can be similarly carried out for other LGT models in one spatial dimension.
Regarding physical observables, our method gives in principle access to all thermodynamic quantities. However, we found that in practice the computation is limited to only relatively high temperatures, and exploring lower temperatures has a much higher computational cost compared to the determination of the DOS.
In conclusion, the method employed here opens another aspect of the quantum many-body problem, and in particular LGT, to the exploration with tensor networks. The DOS method allows us to distinguish the features of the density of states, but directly using it for thermodynamic observables seems to be limited to relatively high temperatures. Since this limitation does not affect other stan-dard MPS techniques, a combination of methods seems to be the most promising strategy to explore other LGT problems in the future.