Finite Density $QED_{1+1}$ Near Lefschetz Thimbles

One strategy for reducing the sign problem in finite-density field theories is to deform the path integral contour from real to complex fields. If the deformed manifold is the appropriate combination of Lefschetz thimbles -- or somewhat close to them -- the sign problem is alleviated. Gauge theories lack a well-defined thimble decomposition, and therefore it is unclear how to carry out a generalized thimble method. In this paper we discuss some of the conceptual issues involved by applying this method to $QED_{1+1}$ at finite density, showing that the generalized thimble method yields correct results with less computational effort than standard methods.


I. INTRODUCTION
Observables in non-perturbative field theory are predominantly calculated through stochastic methods and importance sampling. These techniques rely on the integrand of the path integral in Euclidean space, e −S , being real and non-negative, so that it can be interpreted as a probability. In several important cases, among them most theories at finite density, S is complex, resulting in a nonreal Boltzmann factor e −S . Furthermore, e −S typically is a rapidly oscillating function of the field variables and computing the path integral becomes a challenge for any numerical method due to cancellations between different regions in field space. This issue is known as the sign problem and is the main roadblock to the lattice field theory study of dense quark matter (and several important problems in condensed matter physics).
In past years, a program to bypass the sign problem based on complexifying the field variables was set in motion. The path integral is evaluated, not as an integral over all real fields, but as an integral over a chosen integration contour in complex field space. Cauchy's theorem guarantees that, for a wide class of such contours, the value of the path integral is not altered by this procedure. The original suggestion [1][2][3][4] was to use a certain combination, M T , of Lefschetz thimbles, the multidimensional analogue of the stationary phase paths of the complex function of one variable.
A number of algorithms applied the idea of integrating over M T to toy models [2,[5][6][7] and fermionic field theories [6]. Difficulties with this approach soon became apparent. First, it is difficult to determine the correct combination of thimbles which is equivalent to the original path integral, a feat accomplished analytically only * aalexan@gwu.edu † gbasar@uic.edu ‡ bedaque@umd.edu § hlamm@umd.edu ¶ srl@umd.edu in some simple models [8][9][10][11]. Another obstacle is correctly sampling different thimbles with the appropriate weight [12,13]. For these and other reasons, new manifolds have been suggested. In [14], a so-called generalized thimble method (GTM) was proposed, in which a manifold is obtained by evolving the real space by the holomorphic flow. If every point of R N , the original domain of integration, is evolved for a "time" T , a manifold M T is obtained that i) yields equivalent results to the original real space and ii) approaches M T in the limit T → ∞, consequently improving the sign problem. However, for large flow times T , the parameterization of M T by the real plane become ill-behaved. In particular, small regions on the real plane are mapped to large regions of M T , and most of the real plane is mapped to singular points of M T which do not contribute to the integral. A Markov chain on the real plane is then unable to explore M T effectively; this phenomenon is referred to as 'trapping'. The GTM requires choosing T large enough so the sign problem is sufficiently ameliorated but small enough that trapping is not an issue [6,14,15]. Recently, more general manifolds have been used in order to cut down the computational cost of the complexification approach. For instance, a machine-learning technique was used to create a simple parametrization of a manifold L T interpolating points obtained by the computationally expensive holomorphic flow [16]. The holomorphic flow equations are completely avoided in another method proposed in [17][18][19], where M obtained by maximizing the average sign within families of manifolds.
The application of these methods to gauge theories has developed more slowly, and only 0 + 1 dimensional cases and one-plaquette models have been studied [20]. Gauge invariance brings a host of conceptual issues in the application of these methods. In this paper we applied the generalized thimble method to QED 1+1 . Unlike in previously-studied models, the thimble decomposition in this model is not well-defined. This feature of gauge theories calls into question the use of the GTM, which is designed to approximate the thimble decomposition. In this paper, we show that the GTM formally yields correct results in this model. Additionally, we show how the definition of Lefschetz thimbles can be extended to be welldefined in an abelian gauge theory. Finally, we numerically demonstrate that the GTM both yields correct results, and improves the sign problem.
The paper is organized as follows. In Sec. II, we briefly review the generalized thimble method, in Sec. III we present the QED 1+1 model we study and discuss the particular features of it relevant to us. Sec. IV presents the properties of the holomorphic flow in abelian gauge theories that insure correct answers. The results of the numerical calculations are presented in Sec. V and further discussion is found in Sec. VI.

II. GENERALIZED THIMBLE METHOD
We wish to calculate the thermal expectation value of an observable O. In the path-integral formalism this is given by where S[φ] is the Euclidean action. If S is real, then the Boltzmann factor e −S is real and positive, and can be interpreted as a probability distribution. This expectation value can be efficiently approximated by importance sampling according to e −S , and averaging the value of O over the collected samples: In the limit N → ∞, this approximation converges to the correct answer. It follows from the central limit theorem that the standard deviation of this estimator scales as δ O ∼ N −1/2 . In this paper we are concerned with a fermionic field theory at finite chemical potential where the Euclidean action S = S R + iS I is complex. In this case, the complex Boltzmann factor cannot be interpreted as a probability. The expectation value O must now be evaluated by reweighting: where O S R is introduced to denote the expectation value taken with respect to the real part of S. This estimator will still converge to the correct answer; however, the standard deviation of the estimator now scales with the average sign σ ≡ e −iS I For fermionic theories in a finite volume V and at finite density n, the average sign depends exponentially on V . Therefore, estimating a typical O with a fixed precision requires a time exponential in V . This is the sign problem.
In the remainder of this section, we describe a strategy that has been effective for ameliorating the sign problem in other models. Our strategy is to deform R N to M ⊂ C N , where M is an N-real-dimensional manifold embedded in complexified field space on which the path integral remains unchanged but the sign fluctuations are milder. Cauchy's integral theorem guarantees that the integral over R N will be the same as the integral over M, as long as the integrand is holomorphic and the asymptotic behavior of M coincides with R N . As shown in Appendix A, the integrands in the numerator and denominator of Eq. (1) are both holomorphic for our theory; therefore we have (3) How this deformation ameliorates the sign problem can be seen by examining what happens to the average sign: The numerator of σ has a holomorphic integrand, and therefore is independent of M. However, because the integrand of the denominator is not holomorphic, σ will take a different value on different manifolds. Integrating over an appropriately chosen M instead of R N decreases the denominator, thereby alleviating the sign problem. We now introduce a class of manifolds, M T , obtained by continuously deforming R N via the holomorphic flow equation: Letφ t (φ) denote a solution to this equation, with initial conditionsφ 0 (φ) = φ ∈ R N . The manifold M T is defined by the collection of pointsφ T (φ) for all φ ∈ R N . In other words, M T consists of all points in R N after flowing via Eq. (5) for fixed time T . The flow time is a free parameter that we can tune in our computations to optimize total computational expense. We parametrize every pointφ T (φ) on M T by the real coordinates in φ. Using the φ coordinates the expectation value in Eq. 3 can be written as where we introduce the Jacobian matrix given by In order to perform a Monte-Carlo computation of Eq. (6) we need a real, non-negative probability distribution.
The remaining phase, exp(−i Im S eff ), is reweighted as before: where subscript means the importance sampling is performed with the probability distribution e − Re S eff .
We now explain briefly why deforming the original domain, R N , to M T is likely to improve the sign problem. In the limit of T → ∞, most fields acquire a large S R and decouple from the path integral due to their exponentially small Boltzmann weights, e −S R [φ T ] . In terms of the parameterization by the real fields (where the Monte-Carlo sampling is performed), the support of the path integral comes from the points that flow into the critical points where the flow stops, and the small neighborhoods around them. These neighborhoods are stretched by flow and generate N -real-dimensional surfaces around the critical points. In other words, in the large T limit, M T is the union of a set of approximate thimbles attached to these critical points. At the same time, S I is constant with flow, so the variation of S I on a given approximate thimble (in the real field parameterization) is small. Consequently the sampled fields has a small variation of the phase, hence a milder sign problem 1 .
In order to compute det J, the system of equations with J ij (t = 0) = δ ij needs to be solved where the Hessian matrix is This procedure is expensive, so we instead use an estimator W in place of det J, defined by [21] The difference between W and det J is reweighted when computing observables where S eff ≡ S − ln W , ∆S ≡ S eff − Re S eff , and · S eff is the average with respect to exp(− Re S eff ). In this way, the expensive det J needs to be computed only for the configurations used for measurements, rather than for every step of the Monte Carlo chain.

III. QED IN 1+1 DIMENSIONS
As a demonstration of the generalized thimble method in a gauge theory, we will use a three-flavor version of QED in 1+1 dimensions. In the continuum, the Euclidean action is m is the bare fermion mass, g is the gauge coupling constant, the latin index a sums over N F flavors, and Q a denotes the charge assigned to flavor a. Because QED 1+1 provides a long-range force, states with a net charge have a divergent energy in the thermodynamic limit and are excluded of the spectrum. In order to have a finite fermion density while maintaining charge neutrality we then need more than one flavor. The two flavor case with fermions with opposite charge (and equal mass) is not interesting for our purposes as the sign problem disappears. We use a model with three flavors with charge assignments Q 1 = Q + = +2, and Q 2,3 = Q − = −1, so states with a finite fermion number but zero charge are possible (for instance, a three fermion state with one fermion of each flavor). We will refer to each fermion as a "quark" and groups of three quarks, one of each flavor, as a "baryon". In this model, we take the quark chemical potential µ Q to be the same for all flavors.
Discretizing Euclidean spacetime and integrating out the fermions yields a lattice action in which the only degrees of freedom are the bosonic vector field A µ , with two components at each site of the lattice. The action is where D (a) denotes the fermionic matrix for flavor a, and P r denotes the primitive plaquette with r at the lowerright corner: Abovet andx are the unit vectors in time and space direction. We discretize the fermion action using the staggered formulation. The Kogut-Susskind staggered fermion matrix for flavor a is given by Note that one staggered flavor corresponds in the continuum limit to either a pair of two-spinor-component fields or one four-spinor-component field. We take all flavors to have the same mass m a = m. For our specific assignments of flavors and charges, the lattice action may be written where D + (D − ) gives the fermion matrix for the flavor with charge Q + (Q − ). This notation will be used in the following sections.

IV. FLOW AND THIMBLES IN AN ABELIAN GAUGE THEORY
The integration over the manifold M T yields correct physical expectations values: this, we reiterate, is guaranteed by Cauchy's theorem. However, it is not sufficient for the GTM to get the correct answer; in order to be of practical use, it must also improve the sign problem over the real plane calculation. Previous, heuristic arguments for indicating that GTM improves the sign problem relied upon M T approaching the Lefschetz thimbles at T → ∞. Where no such unique thimble decomposition exists, it is not obvious what the large-T limit of M T is, and whether this manifold will have a reduced sign problem.
In QED 1+1 dimensions, there are two separate obstacles to defining a unique thimble decomposition: critical points are degenerate (the action does not change along gauge orbits), and lines of flow may connect one critical point to another (Stokes' phenomenon). In this section we discuss these difficulties, showing that the Lefschetz thimble formalism can be changed to accommodate gauge orbits in general, and neither obstacle substantially adversely affects the performance of the GTM.
We look first at the difficulties presented by gauge orbits in configuration space. A Lefschetz thimble is defined from an isolated critical point z c as the union of all solutions to Eq. (5) which asymptote to the critical point at large negative flow times: lim T →−∞ z(T ) = z c . For a holomorphic function of N complex variables, each isolated critical point is a saddle point with N stable and N unstable directions; therefore, the Lefschetz thimble defined in this way is an N -real-dimensional surface. The Lefschetz thimble decomposition is a union of a certain subset of these thimbles. However, in a gauge theory, all critical points can be continuously connected to infinitely many other critical points by gauge transformations, i.e., critical points form gauge orbits. Therefore, in order to use the Lefschetz thimble decomposition to understand how the GTM behaves when applied to a gauge theory, we must understand how the Lefschetz thimble formalism can be repaired in the presence of gauge orbits.
For QED 1+1 , where the gauge group is abelian, the degeneracies introduced by gauge symmetry can be resolved in a straightforward fashion by gauge-fixing. Before proceeding, we stress that we will ultimately perform calculations without gauge-fixing; gauge-fixing is merely a convenient tool to see how the Lefschetz thimble formalism may be repaired for an abelian gauge theory. In the complexified field space a general gauge transformation is given by where α(x) is any complex-valued function on the lattice points andμ denotes the unit vector along the direction µ ∈ {0, 1}. With V lattice sites, A x,µ can be represented as a 2V dimensional vector, so the complexified field space is C 2V . The gauge orbit of any configuration A x,µ is obtained by adding to it all vectors of the form α(x + µ) − α(x), for every α(x). This is a (V − 1)-dimensional space since a constant α(x) = α does not generate a gauge transformation. So the gauge orbits form a set of parallel (V − 1)-dimensional linear subspaces, and each gauge-fixed slice is a (V + 1)-dimensional subspace. Every A µ (x) can be decomposed with A µ (x) parallel to the gauge orbit and A ⊥ µ (x) orthogonal to it 2 . We can choose A ⊥ µ (x) as the representative of A µ (x) in every gauge orbit. In this way, we decompose the original real configuration space as , and G is a single gauge orbit. The gradient ∂S ∂A is orthogonal to the gauge orbits so a gauge-fixed slice, defined by constant A µ (x), is invariant under the holomorphic flow.
In each gauge-fixed slice, critical points are now isolated, and there is no gauge-related obstruction to the definition of a Lefschetz thimble decomposition. Furthermore, the behavior of flow in this gauge theory can be understood by considering the behavior of flow in each gauge-fixed slice. Let M T be the result of flowing M 0 for some time T . Because the flow commutes with gauge fixing, the result of flowing the entire gauge-free integration domain is M T ⊕ G. In our simulations, we do not gauge fix, so that our simulations perform a random walk in the gauge orbit space G -this has no effect since we only evaluate gauge-invariant quantities.
The second difficulty encountered in the thimble decomposition in QED 1+1 is Stokes' phenomenon, where a flow line connects two critical points. A schematic representation of Stokes' phenomenon is shown in the central panel of Fig. 1. Stokes' phenomenon can be avoided by adding a small i to some parameter of the action. While the resulting thimble decomposition is well-defined, it differs depending on the sign of . The left and right panels of Fig. 1 show the flow assuming such a modification of the action with opposite signs. In each case, the thimble decomposition is well-defined, but the combination of thimbles equivalent to the real domain is different depending on the sign of .
The GTM does not rely on the thimble decomposition, but it is instructive to see how Stokes' phenomenon affects the behavior of the holomorphic flow. We investigate this particularly in the case of QED 1+1 , for which Stokes' phenomenon occurs at all values of chemical potential µ and coupling g. For instance, consider fields with A 0 (x) = 0 and A 1 (x) = A 1 (constant). This particular slice of field space is depicted in Fig. 2. In this slice, there are three critical points: one at A 1 = 0 and two others at non-zero values of A 1 , both of which have flow lines The presence of Stokes' phenomenon does not cause any discontinuity in the flowed manifold. On the other hand, Stokes' phenomenon does produce undesirable "bumps" in the deformed paths, depending on the orientation of the critical points and the starting manifold for the flow. In Fig. 2, several M T are drawn, flowed from a shifted version of the real line M = R + 0.01i. As the amount of flow increases, a bump is created near the origin.
The effect of this bump is best understood by considering the T → ∞ limit, in which it becomes the most pronounced. In this limit, the bump travels up one half of the imaginary axis, and directly back down again, producing a closed contour of integration. The sum of these two segments of the contour cancel and do not contribute to the integral. However, the average sign is decreased by the presence of such closed contours. The effect on the sign problem is bounded because the bump is produced along a path of steepest descent, so that the integral of |e −S | along the bump is finite. This ensures that in the large T limit, σ asymptotes to a non-zero value.

V. RESULTS
In this section we present some numerical results supporting the point that the GTM leads to the correct results -and improves the sign problem -even in gauge theories where the thimble decomposition is not well-defined. We choose the bare parameters g = 0.5 and m = 0.05 so that the renormalized baryon mass lies below the lattice cutoff scale. We estimate the baryon mass m B by observing the onset of non-zero charged fermion density at low temperature (see the Silver Blaze discussion below). The relation between the baryon chemical potential and the quark chemical potential allows us to find am B ≈ 0.6. Since the "baryon" in our model is made up of three fermions we use µ B = 3µ to quote the values of the chemical potential.
In this study, we have undertaken calculations on a fixed spatial lattice size at four different temperatures: n t ×n x = 6×10, 10×10, and 14×10 using both the GTM and a calculation done with real fields. The parameters for the R N and M simulations can be found in Table I.
To separate the effect of the phase fluctuations on M T and the fluctuations induced by approximating ln det J with W , we write the reweighting factor ∆S = i Im S eff + ∆J. The phase factor exp(−i Im S eff ) is a pure phase and converges to the "residual phase" on the thimble as T → ∞. The real factor ∆J, with exp(−∆J) = | det J|/|W |, is the reweighting necessary to correct for using W instead of the Jacobian det J in the Monte-Carlo process.
In order to study the speed-up of the GTM compared to real-plane calculations, one must consider: t config , the wall-clock time required to generate a statistically in- Σ is bounded between 1/n config and 1 and is an estimate of the fraction of configurations that effectively contribute to the statistics of the observable. A small value of Σ indicates that the averages are dominated by a small number of configurations and, consequently, more configurations need to be used for a reliable estimate. The reweighting is only necessary in the calculations on M T , since Σ is 1 when no flowing is done. We find empirically that the statistical power tracks the average sign, decreasing with volume and µ, although as explained below for the average sign, it recovers at saturation densities. Unlike the average sign, Σ tends to decrease with T (although it asymptotes to a finite, non-zero value). As a consequence, the increase of flow time T does not necessarily yield superior performance.
We define a figure of merit h T for a fixed flow time T as A ratio h T1 /h T2 estimates the relative speed-up of flow time T 1 over flow time T 2 , for a fixed desired precision. We expect σ to increase with T because flowing improves the relative sign, but t config also increases because longer flow requires more computational time. Due to the use of W , an approximate Jacobian that is computationally faster, the statistical reweighting plays a nontrivial part in judging the speedup. As we present our results we will focus on the fermion density, as a representative physical observable, and the sign average σ = exp(−i Im S eff ) as a characteristic of the sign problem on the flowed manifold. We will compare our results for the density against the results from a free gas with the same mass as our "baryon" mass. This model is expected to describe well the results at small density but as the number of particles increases the interactions will modify the results significantly. As a first test, we studied a 6 × 10 lattice where σ is large enough on R N that computations can be taken at all µ B /m B up to lattice saturation. A small flow time T = 0.01 was used to show that GTM can improve the sign and that the two methods agree within the small error bars over the entire µ B /m B range as seen in Fig. 3. This should be taken as an empirical confirmation that GTM respects all of the conditions upon complexification that are necessary to produce physically correct results despite the lack of a unique Lefschetz thimble decomposition. On this small lattice, the average sign and statistical power are nearly the same for M T =0.01 and R N , but t config is 6 times larger for the flowed manifold, therefore h 0.01 /h 0 ≈ 1/6, indicating that there is no computational benefit to using the GTM.
In contrast to other models (e.g. the 1+1 Thirring model of Ref. [18]), in QED 1+1 σ appears to approach unity for large µ B /m B , with σ → 1 as µ B /m B → ∞. The absence of a sign problem at saturation densities can be understood directly from the fermionic part of the effective action. In the large-µ B limit, the fermion determinant becomes det D ± → exp [βV (µ + iQ ± A 0 )]. Thus, the product of all three fermion determinants, in this limit, is a positive real number e 3βV µ , and there is no sign problem. This phenomenon has been previously noted in QCD [20], and appears to be a general feature of gauge theories.
On the 10 × 10 lattice (seen in Fig. 4), the average sign obtained with a flow time T = 0.01 is slightly improved compared to R N . For µ B around 2.8m B the average sign on both manifolds become indistinguishable from zero. For those values of µ B /m B , we increase T . For all but µ B /m B ≈ 3, T = 0.02 is sufficient to discriminate the sign from zero, that is σ > 2 σ where σ is the statistical error on σ. For µ B /m B = 3 we increase the flow time T further, taking advantage of the fact that Σ will asymptote to a fixed value (so that flow time can be increased without a penalty to the statistical power). The results for h T /h 0 can be seen in Table II. We find that for this lattice, flowed manifolds can reduce computational time compared to the real plane for some values of µ B /m B .
On our coldest lattice, 14×10, one can see the onset of Silver Blaze phenomenon at small µ B /m B in Fig. 5, as well as development of a plateau at the one baryon threshold. For this lattice, t config for T = 0.01 is 7 times that of the real plane. We find for different µ B /m B , Σ T ranges from 0.65 to 1 and h 0.01 /h 0 ranges from 0.15 to 28. Since the figure of merit for certain values of µ B /m B again exceeds unity, holomorphic flow is seen to reduce the computational time required to achieve a fixed-precision result at µ B /m B ≈ 1.5 .

VI. DISCUSSION AND PROSPECTS
Working with QED 1+1 , we have shown that the generalized thimble method may be applied to gauge theories without encountering any fundamental obstacles. None of the difficulties that plague the definition of Lefschetz thimbles in a gauge theory need to be confronted despite the fact that the original motivation for the GTM explicitly invoked the thimble decomposition. Additionally, despite the holomorphic flow requiring greater computational time to produce a single configuration compared to the real plane, we have shown that the improvement in the average sign can reduce the total time needed to compute observables at a fixed precision to less then an equivalent real time calculation. Although the holomorphic flow gives a relative speed up over the real plane for QED 1+1 , the computational cost is still large. Extending this work to larger lattices, higher dimensions, and non-abelian theories will require the development of faster algorithms to deal with the Jacobian that do not suffer from small statistical power as the approximate Jacobian does at large flow times and larger lattices. For simplicity, we will consider an action involving only a single species of fermions, and therefore only one fermion determinant det D. All arguments given here generalize easily to the cases of two or more differently-charged species. That e −S[φ] is entire is easily established: by definition, this integrand may be written where the 'bosonic' action S B is manifestly holomorphic. Each component D ij of the fermion matrix is a holomorphic function of the gauge fields; since the determinant is a polynomial in those components, det D is a holomorphic function of φ as well.
To see that integrands involving fermionic observables are entire, we write an expectation value in terms of the original, fermionic path integral.
With N sites, the fermionic exponential e −ψDψ may be expanded in 2 2N terms: each Grassman variable may be included in a term or not, and there are 2N Grassman variables. The commuting part of each term is a product of finitely many components of D, and therefore is a holomorphic function of the fields φ. Multiplying by any combination ofψψ and integrating over DψDψ has the effect of selecting one of these coefficients. Therefore, the integral over fermionic fields yields a holomorphic function of φ.