Towards overcoming the Monte Carlo sign problem with tensor networks

. The study of lattice gauge theories with Monte Carlo simulations is hindered by the infamous sign problem that appears under certain circumstances, in particular at non-zero chemical potential. So far, there is no universal method to overcome this problem. However, recent years brought a new class of non-perturbative Hamiltonian techniques named tensor networks, where the sign problem is absent. In previous work, we have demonstrated that this approach, in particular matrix product states in 1 + 1 dimensions, can be used to perform precise calculations in a lattice gauge theory, the massless and massive Schwinger model. We have computed the mass spectrum of this theory, its thermal properties and real-time dynamics. In this work, we review these results and we extend our calculations to the case of two ﬂavours and non-zero chemical potential. We are able to reliably reproduce known analytical results for this model, thus demonstrating that tensor networks can tackle the sign problem of a lattice gauge theory at ﬁnite density.


Introduction
Lattice methods have proven to be indispensable when addressing questions in the non-perturbative regime of quantum chromodynamics (QCD). Lattice QCD (LQCD), initially introduced in the 1970s by Wilson [1], has succeeded in elucidating many aspects of the strong interaction from first principles, including e.g. the computation of quark and hadron masses, hadron structure related quantities as well as non-zero temperature properties. The standard tool of LQCD are Monte Carlo (MC) simulations, performed on world's most powerful supercomputers. Despite the great progress achieved in the last decades, there are still areas where LQCD fails to give precise quantitative answers. This concerns, in particular, parameter regimes where MC simulations encounter a sign problem. Arguably the most important such case, at least from the point of view of this paper, is at non-zero chemical potential, i.e. at finite baryon density. When the chemical potential is non-zero, the standard Boltzmann factor of LQCD in the Euclidean formulation becomes complex and can no longer a e-mail: kcichy@th.physik.uni-frankfurt.de be interpreted as a probability measure. This undermines the whole principle of MC simulations. Although partial solutions to this problem are known, they only allow to simulate relatively small baryon densities, in particular much smaller than the one corresponding to the conjectured critical endpoint in the density-temperature phase diagram. As a consequence, there is broad search for alternative approaches, including Lefschetz thimbles [2], complex Langevin simulations [3] and density of states methods [4]. Yet another thread of research, one that we concentrate on in this paper, is related to the tensor networks (TN) approach.
The TN approach 1 , originally introduced in the context of condensed matter physics and further developed thanks to quantum information theory, has been successfully applied to the description of quantum many-body systems, including, in last years, also lattice field theory. The latter included computations of spectra [7][8][9][10], thermal states [11][12][13][14][15], phase diagrams [16][17][18][19] and real-time evolution [8,[20][21][22]. In this paper, we review the approach and our results for the Schwinger model. Moreover, we go one step further by simulating the two-flavour Schwinger model in a setup where MC simulations would suffer from the sign problem, at non-zero chemical potential.
The paper is organized as follows. We formulate the Schwinger model for a general number of flavours in Sec. 2. In Sec. 3, we describe the basics of the tensor network method. Sec. 4 summarizes our findings for the chiral condensate in the one-flavour model at zero temperature. In Sec. 5, we review our results at finite temperature, considering again the chiral condensate. In Sec. 6, we discuss the results from the two-flavour Schwinger model at non-zero chemical potential. Finally, Sec. 7 summarizes and a short discussion of the prospects of the TN approach for lattice gauge theories is given.

Multi-flavour Schwinger model
The Schwinger model [23] is quantum electrodynamics in 1+1 dimensions (QED 2 ). It has been widely used as a toy model for testing new lattice methods, before applying them e.g. to QCD. Interestingly, QED 2 shares certain properties with QCD 4 , such as confinement and chiral symmetry breaking.
Our lattice formulation uses Kogut-Susskind staggered fermions [24], leading to the following Hamiltonian for N sites and F flavours, with open boundary conditions (OBC): where: a -lattice spacing, g -coupling, m f -f -th flavour fermion mass, κ f -its chemical potential, φ n, f -fermionic fields describing flavour f on site n, L n -operator that gives the electric flux of link n, θ n ∈ [0, 2π] -conjugate operator such that e ±iθ n acts as the electric flux raising/lowering operator. The Gauss law is expressed by: and allows to integrate out the gauge fields, leaving only the electric field on the left boundary as an independent parameter (taken to be zero in this work). We write the above Hamiltonian in dimensionless form and we use a Jordan-Wigner transformation, φ n = l<n (iσ z l )σ − n , φ † n = l<n (−iσ z l )σ + n , where σ z/± n are Pauli matrices acting on site n and where 1 For a pedagogical introduction, see e.g. Refs. [5,6].

CONF12
we order the fermions on each site such that φ n, f = φ nF+ f . The Hamiltonian in this equivalent spin formulation describes a spin chain of length NF and reads: where Note that the gauge symmetry (manifested in the Gauss law) has effectively generated long-range interactions, of the form σ z n σ z m (with arbitrary site-flavour index n, m) in the spin language, and also a site-dependent magnetic field. The hopping term always involves fermions of the same flavour, i.e. it acts on pairs of spins at a distance F in the effective NF-site chain.
In the following, we consider the 1-flavour and 2-flavour cases. With the former, one can demonstrate, for example, calculations of the spectrum, of ground state expectation values and of thermal properties. However, the chemical potential has no physical effect with only one flavour -therefore, to exhibit that the TN approach can tackle the model at non-zero density, we employ the 2-flavour case.
In this way, we show that in the regime where a MC simulation would suffer from a sign problem, the TN methods can still produce reliable results.

Tensor network approach
Among the tensor network methods, the most successful one for one-dimensional systems (1+1dimensional quantum field theories) is the matrix product states (MPS) [25][26][27][28][29][30]. Below, we use MPS for all our investigations at zero temperature. The MPS ansatz for a system of N sites has the form: where {|i n } d−1 i=0 are basis states for each site n and d is the dimension of the one-site Hilbert space. Each matrix A i n n is D-dimensional and D is called the bond dimension and determines the number of parameters in the MPS ansatz. It has been shown that MPS describe accurately ground states of local gapped Hamiltonians, but their practical application extends to more general models.
The task of finding the MPS approximation to the ground state (GS) boils down to finding the components of the tensors A n that minimize the expectation value of the Hamiltonian, E = Ψ|H|Ψ / Ψ|Ψ . This is performed variationally, by minimizing E with respect to one tensor A n at a time and changing n successively (sweeping from left to right and back) until global convergence is achieved. Such local minimization problem consists in solving the eigenvalue problem of an effective Hamiltonian acting on site n and its neigbouring virtual bonds. This variational algorithm [29] is related to the original density matrix renormalization group (DMRG) formulation [31,32], which was better understood in the language of MPS. Note that we use OBC, which is, in general, more efficient and results in more stable numerical behaviour. The commonly used graphical representation of an MPS and its contractions that give the energy E is shown in Fig. 1 Apart from the GS, the variational algorithm can find excited states. Knowing the ground state MPS, |Ψ 0 , this proceeds by constructing a projector onto the orthogonal subspace and finding the GS in this subspace. Successively projecting out further excited states, one can find any number of them (however, the accumulated error may become large for higher excited states). In this way, one can access the vector and scalar particles in the Schwinger model spectrum. Note, however, that the usage of OBC violates translational invariance and hence identification of particles, which are zero-momentum excitations, is more difficult. In particular, the scalar meson, naively the second excitation of the one-flavour Schwinger model, is above momentum excitations of the vector particle. Nevertheless, unambiguous identification of particles is still possible.
The application of the MPS approach to non-zero temperatures requires its extension allowing for a description of operators, in this case density matrix operators [33][34][35]. One introduces matrix product operators (MPO) of the form: where the tensors M n play an analogous role to the tensors A n in the MPS ansatz. For density matrices appropriate for the description of a system at finite temperature, an MPO approximation can be obtained using imaginary time evolution of the identity operator [33], ρ(β) ∝ e −βH = e − β 2 H e − β 2 H , with β ≡ 1/T being the inverse temperature. For a detailed description of the technical steps for the one-flavour Schwinger model, we refer to our original papers [11,12]. The efficient application of the imaginary time evolution to find thermal states requires specific algorithms, for instance applying discrete time steps in a Suzuki-Trotter expansion [36,37] plus using appropriate approximations for the long-range terms in the Hamiltonian. We start with the identity operator, ρ(0), which can be represented trivially by an MPO with bond dimension one. Then, we evolve the density matrix and approximate each step by an MPO with the desired maximum bond dimension, using the so-called Choi isomorphism [38], |i j| → |i ⊗ | j that vectorizes the density operators, i.e. transforms the MPO into an MPS with physical dimension d 2 . Finally, we minimize the Euclidean distance between the original and final MPS and repeat the procedure until reaching the inverse temperature β.  We start reviewing our results with a zero-temperature calculation of the ground state in the oneflavour Schwinger model. We illustrate the main steps concerning an example of the determination of the ground state expectation value of the chiral condensate, Σ = ψ ψ , which can be written using spin operators as: In the massless case, this expectation value can be computed analytically, yielding the result Σ T =0,m=0 /g = e γ E /2π 3/2 ≈ 0.159929 [39], where γ E ≈ 0.577216 is the Euler-Mascheroni constant.
The analytical computation at non-zero fermion mass is no longer possible. The non-perturbative MPS determination requires subtraction of a logarithmic divergence, which is present already in the free theory and can be removed with: where Σ free (m/g, x) is the free theory condensate at fermion mass m/g and lattice spacing corresponding to the inverse coupling x.  Table 1. Continuum values of the T = 0 chiral condensate Σ/g for different fermion masses. We compare with results of Ref. [9] and with the analytical result in the massless case or the approximated result from Ref. [40] in the massive case.
An MPS simulation is performed for a finite bond dimension D, with a chain of finite length N and at a finite lattice spacing, given by the parameter x. In Fig. 2, we illustrate the sequence of extrapolations needed to obtain the continuum result, which is the aim of the computation. First, we extrapolate to infinite bond dimension, thus removing the effect of approximating the ground state with an MPS of finite D. The typical behaviour, illustrated in the upper left plot of Fig. 2, is saturation of finite bond dimension effects at rather moderate value of D of order 100. As our final value, we take the one corresponding to the maximum computed bond dimension, D = 160 in this case, and we estimate its uncertainty to be of order of the difference between the results at D = 160 and D = 140. Note that, given the approximately exponential approach to the value at 1/D = 0 that we observe in many cases, this error estimate is rather conservative. It should be emphasized that the central value obtained in this way is always compatible with the result of exact diagonalization (ED), in case the latter is feasible. 2 Having the estimates of the 1/D = 0 values of the condensate for a few values of the volume N, one can perform the infinite volume extrapolation, as illustrated in Fig. 2 (upper right). In order that the leading order 1/N behaviour is reliable, one needs to ensure that the ratio N/ √ x is large enough, in this case of order 20-30. Finally, one can subtract the free theory condensate (Eq. (7)) and perform the continuum extrapolation. In the illustrated case (lower plot of Fig. 2), we use Eq. (4.2) of Ref. [12], which is a fitting ansatz quadratic in the lattice spacing 1/ √ x and including logarithmic corrections of O(log(x)/ √ x). Note that we also test sensitivity to higher order corrections of O(1/x 3/2 ) and to the choice of the fitting interval in x. In the end, we obtain a continuum result with all sources of possible uncertainties fully estimated and quantified.
Our final results are given in Tab. 1, together with a comparison to another recent MPS computation [9] and the exact result (m = 0) or an approximate result from Ref. [40] (m > 0). Note that the former comparison is highly non-trivial, since the authors of Ref. [9] used a different approach of infinite MPS, which works directly in the thermodynamic limit, instead of our approach of extrapolating to this limit from selected values of N. In general, the obtained precision is, in both cases, very good, at the level of around five-six significant digits.
Apart from the GS expectation value of the condensate, we also computed the spectrum of the one-flavour Schwinger model, reaching similar precision. For details, we refer to Ref. [7]. 5 Results -1-flavour Schwinger model at T > 0 As described above, the thermal computation is somewhat more complicated. However, the necessary steps in the numerical procedure are very similar, with the addition of another extrapolation that has to be performed. Since the imaginary time interval [0, β] is divided into discrete steps of length δ, one needs to take the limit δ → 0 to eliminate the error necessarily introduced by non-zero δ. This extrapolation is performed as the second one, after the bond dimension extrapolation (see Ref. [12] for explicit plots).
The continuum results for the massless condensate are shown in Fig. 3. The long-range terms in the Hamiltonian require specific approximations to implement the corresponding steps of the evolution. The left plot shows the result when a Taylor expansion is used. We reproduced the analytical curve of Ref. [39] in a wide range of temperatures. Note, however, that the high temperature region with gβ ≤ 0.5 does not agree with the analytical result in this plot. This fact, resulting from enhanced cut-off effects in this region, required further investigation. It is possible to increase the precision of the computation by applying instead a truncation to the evolution operator, which turns out to be equivalent to a truncated model, with maximum flux on the link set to a specified value L cut . We found that L cut = 10 effectively corresponds to the untruncated model, but the presence of the flux cut-off leads to a reduced computational effort. Thus, we could simulate at smaller lattice spacings and obtain more reliable continuum results, illustrated in the right plot of Fig. 3. In addition to reducing the uncertainties in the whole range of temperatures, we obtained agreement with the analytical curve even for small gβ.
In Fig. 4, we also show our continuum result for two selected fermion masses, m/g = 0.125 (left) and m/g = 0.5 (right). We note that, as expected, the thermal curves approach the T = 0 result when gβ is increased. The agreement with the approximate result of Ref. [40] is moderate. However, after our work, a cross-check was provided by the authors of Ref. [15] using the infinite MPS approach and full agreement was observed.  The essential motivation for the TN approach comes from the possibility of overcoming the notorious sign problem that plagues MC simulations of QCD at non-zero chemical potential. The results presented in the previous sections have demonstrated the feasibility of the MPS/MPO approach when applied to lattice gauge theories. However, the analyzed system, the one-flavour Schwinger model both at zero and non-zero temperatures, does not suffer from the sign problem when simulated with MC. Now, we move on to a case where the sign problem indeed appears -the two-flavour model at non-zero chemical potential. For a more detailed description of our results, see Refs. [41,42].

EPJ Web of Conferences
Our aim was to reproduce the analytical result of Narayanan [43], who computed the isospin number, ∆N ≡ N 0 − N 1 , i.e. the difference between the numbers of particles of different flavours, as a function of the isospin chemical potential, µ I , which corresponds in our notation to N(ν 1 − ν 0 )/2x. It was shown that ∆N undergoes jumps that correspond to first-order phase transitions. In our MPS setup, we could calculate the GS energy of phases with different isospin numbers. The upper inset of CONF12 Fig. 5 shows the GS energies of the phases with ∆N = 0 and ∆N = 2 and the crossing of these energies determines the location of the first transition, which corresponds to µ I /2π = 1/2 in the continuum. Our result is, obviously, influenced by effects of finite D, N and x. We extrapolate away these effects in a similar manner as for the chiral condensate. In the end, we obtain transition locations as illustrated in the main plot of Fig. 5, provided the volume is large enough (see the lower inset for effects of too small N/ √ x). Note that the result of Ref. [43] is volume-independent. In our case, residual finite volume effects can occur, related to the fact that N is finite and thus ∆N can become of the order of N. As soon as we satisfy the condition ∆N ≪ N, we no longer observe any dependence on the volume. In this way, our result is in full agreement with the analytical one. We also investigated the massive case, for which no analytical result is available. Contrary to the zero fermion mass case, volume-independence does not hold any more, since the fermion mass provides an additional energy scale in the problem, which affects the locations of the phase transitions.

Conclusion and prospects
Tensor networks are a poweful tool for treatment of general physical systems, especially in lower dimensions. The research of the past few years has demonstrated that they are appropriate for the description of lattice gauge theories and they can provide precise and reliable results both at zero and non-zero temperature. In particular, they are free of the sign problem and hence remain as an interesting direction for future simulations of QCD at non-zero chemical potential. For this, extension of their efficiency for higher dimensions is obviously necessary. Theoretical progress in higher dimensional tensor networks is steadily being made. The natural generalization of MPS to higher dimensions is Projected Entangled Pair States (PEPS) [44,45]. PEPS and other related developments give good prospects of tackling systems in more than one dimension, although they are computationally much more demanding than MPS. Currently, TN methods have already provided one of the most accurate results for certain classes of spin systems in two dimensions, especially for cases where a quantum MC simulation encounters a sign problem. Nevertheless, the application of TN for 2+1-dimensional lattice gauge theories remains difficult. However, it is our aim and the aim of the TN lattice gauge theory community to continue in this direction, which in the future can lead to a full application to 3+1-dimensional QCD.