Fermions at Finite Density in (2+1)d with Sign-Optimized Manifolds

We present Monte Carlo calculations of the thermodynamics of the (2+1) dimensional Thirring model at finite density. We bypass the sign problem by deforming the domain of integration of the path integral into complex space in such a way as to maximize the average sign within a parameterized family of manifolds. We present results for lattice sizes up to $10^3$ and we find that at high densities and/or temperatures the chiral condensate is abruptly reduced.

Monte Carlo methods are critical to the nonperturbative study of strongly interacting quantum field theories and many-body systems.In the lattice field theory approach, one discretizes spacetime and formulates observables as high dimensional lattice path integrals.For systems in thermal equilibrium, such integrals take the form O = Z −1 DA e −S O where Z is the partition function and S is the (Euclidean) action.Path integrals are typically only computable by importance sampling, which relies on interpreting e −S /Z as a probability distribution.However, many theories of interest have complex actions.This sign problem is a major roadblock to the ab-initio study of such systems, including fermions at finite density.
For systems with complex actions S = S R + iS I , a common method is to sample according to the distribution Pr(A) ∼ e −S R (A) , and to express observables as O = O e −iS I R / e −iS I R , where • R means averaging with respect to the S R .This "reweighting" procedure is effective if the average sign σ ≡ e −iS I R is not too small.However, σ typically decreases exponentially with the spatial volume L d , chemical potential µ and inverse temperature β, so for cold dense matter standard reweighting fails [1].In response to this failure, many ideas have been explored: the complex Langevin [2], the density of states method [3], canonical methods [4,5], reweighting methods [6], series expansions in the µ [7], fermion bags [8], and analytic continuation from imaginary µ [9].
In a recently developed family of approaches to taming the sign problem, the original domain of integration M O of the path integral is deformed to a submanifold M of the complexified field space.A multidimensional generalization of Cauchy's integral theorem guarantees, for suitable deformations, that integrals of holomorphic functions (e.g.physical observables) remain unchanged.In contrast, integrals of non-holomorphic functions such as σ depends on M, and therefore a judicious choice of manifold can increase σ and render reweighting feasible.
The first manifolds suggested were sets of multidimensional stationary phase contours called Lefschetz thimbles, M L [10][11][12][13].Analytically, M L have been found in only a handful of cases which include few dimensional integrals and quantum mechanical models [14][15][16][17].Numerous algorithms have been developed to integrate on M L , but these methods have difficulty addressing which set of thimbles reproduce the results on M O [18][19][20][21][22][23][24].To address this, a generalized thimble method was developed.In this approach, one deforms M O via the holomorphic gradient flow for a fixed time T , which yields a manifold M T that approaches M L as T → ∞ [25].The generalized thimble method has been applied to analyze bosonic and fermionic systems at finite density [26][27][28][29][30], real-time linear response [31,32], and gauge theories [33].One drawback to the generalized thimble method is it requires a computationally expensive Jacobian related to the manifold parametrization.This lead to developments in rapidly computable estimators [34] and in applying machine learning to approximate the manifold [35].
To avoid all these difficulties, the sign-optimized manifold method was introduced in [36] wherein one deforms M O to a manifold M S that maximizes σ within a family of manifolds M λ parameterized by a set of real numbers λ i .A similar method is described in [37].To guarantee that the path integral remains invariant under the deformation to M S , it is sufficient that M O is continuously deformable to M λ without crossing any singularity of the integrand.These conditions are satisfied by construct in our family of manifolds because our deformations are smooth and involve only finite shifts of the fields in the imaginary direction.
In this letter, we explore the finite density phase diagram of the two flavor (2 + 1)d Thirring model using the sign-optimized manifold method, extending the range in (T, µ) space beyond what is possible on M O .
We parameterize the manifold M λ by its projection on the real space M O , so that integration on M λ may be achieved by integrating on M O with the inclusion of a Jacobian, which is included into an effective action.Thus, the expectation value of an observable O can be written as:

arXiv:1808.09799v2 [hep-lat] 30 Aug 2018
where Ã(A) is the point on the manifold M λ parametrized by A, S eff ≡ S − ln det J is the effective action, and J is the Jacobian of the parametrization.The average sign on M λ is given by ( The numerator of Eq. ( 2) is independent of λ because it is the integral of a holomorphic function in A, but the denominator depends on λ because e − Re S eff is not holomorphic.We are interested in maximizing this as a function of the manifold parameters λ -this is equivalent This gradient is the phase-quenched expectation value ∇ λ S eff Re S eff , and is therefore free from a sign problem.This allows ∇ λ log | σ λ | to be computed reliably by a short Monte Carlo simulation at each gradient ascent step.To do gradient ascent we use the Adaptive Moment Estimate algorithm [38].We stress that the sign-free nature of the calculation of the gradient is central to the method and allows our calculations to be efficient even when σ is exponentially small.
One potential issue is that the computation of det J is an expensive operation -for general J, this requires time proportional to the cube of the spacetime volume.In [36] it was shown that this computational cost can be avoided by proposing a family of manifolds for which the Jacobian matrix is diagonal.We use a similar family here (details below).A more general ansatz with non-diagonal Jacobian with nearest-neighbor correlations has been shown to improve the sign problem in bosonic theories, with increased computational expense [39].
To integrate on our curved manifolds, we have implemented a modified version of hybrid Monte Carlo (HMC).We define a Hamiltonian and sample according to the distribution P (π, A) ∼ e −H(π,A) .Marginalizing over the momenta yields the distribution P (A) ∼ |detJ(A)|e −S R (A) .Sampling according to P (A) then reweighting with the residual phase e −i Im S eff yields the correct observables.For generic dense Jacobians the derivatives ∂H/∂A x are extremely expensive to compute, but for manifolds with diagonal Jacobians the derivatives are computed analytically and implemented with sparse matrices.Thus, HMC allows for sampling on M λ as fast as sampling on M O .Due to this Jacobian's structure the evolution of Eq. ( 4) can be calculated with implicit and explicit symplectic integrators.
Both were implemented and found to agree.
We now apply the sign-optimized manifold method to the (2 + 1)d Thirring model defined by the lattice action x,y ψa (x)D xy (A)ψ a (y) (5) where −π < A µ (x) ≤ π is a compact bosonic auxiliary field [40][41][42].By virtue of the compact fields, M O = (S 1 ) N and the deformed manifold are submanifolds in the complexified space (S 1 × R) N .The staggered fermion matrix is given by where η ν (x) = (−1) x0+...+xν−1 , the flavor indices a taking values from 1, . . ., N F /2, g is the coupling, and m is the bare mass.There are different lattice actions which naively appear to have as their continuum limit the (2+1)dimensional Thirring model.A substantial literature exists studying different discretizations of the (2 + 1)dimensional Thirring model at zero density, with emphasis on determining the critical N F below which the chiral condenstate ψψ is nonzero when m → 0 [40][41][42][43][44][45].It is however, unclear which discretizations are equivalent in the continuum limit.For our purpose the action in Eq. ( 5) defines what we mean by Thirring model.Integrating out the fermions in Eq. ( 5) gives (6) We presently study the phase diagram in the (T, µ) plane for N F = 2.For µ = 0, the determinant det D(A) is complex and we must address the resulting sign problem.
For insight into a family of manifolds which may increase σ , we look to the µ → ∞ limit of the theory.In this limit, the density matrix is dominated by forward time links, and the path integral becomes βV (7) where only the leading terms in e βµ are included.In this limit, the path integral factorizes, and the sign problem itself comes only from the integral over A 0 .Consequently, we will consider M λ in which A 1 and A 2 remain on M O , and Im Ã0 (x) depends only on A 0 (x), not on any other link.Such factorizable manifolds have the desirable property that J is diagonal.At weak coupling (g 2 → 0) one expects the partition function to be dominated by the saddle point with the smallest action, which is A 0 (x) = iα, A 1 (x) = A 2 (x) = 0 for all x.As found in lower dimensional Thirring models the thimble attached to this critical point can be approximated by a shift of fields in the imaginary direction.This suggests that a shift A 0 (x) → A 0 (x) + iα will improve σ , and this was confirmed in simulations [20,25,28].
Consistent with these observations, we extend the manifolds used in [36] to the following three-parameter family: Every member of the family of manifolds above can be smoothly deformed to (S 1 ) N with the interpolation Ã0 t = A 0 + it (λ 0 + λ 1 cos A 0 + λ 2 cos(2A 0 )) with 0 ≤ t ≤ 1 shows.Moreover, the imaginary shift is bounded, so the condition for the applicability of Cauchy's theorem is satisfied.
The results presented use bare parameters g = 1.08 and m = 0.01.We quote the results of our simulations using lattice units.To demonstrate that we are in the strong coupling regime and to ascertain whether we are not too far from the continuum and thermodynamic limits, we measure the mass of the lowest fermionic and bosonic excitations by fitting the large-time behavior of correlators Using a spatial volume of L 2 = 10 2 we find m f = 0.46(1) and m b = 0.21 (1).The masses depend slightly on L 2 , but in all cases m b /m f 2. This indicates that the system is strongly coupled since the binding energy of the boson is comparable to the 2m f .
In this work, we calculate on six lattice geometries.We perform a series of simulations with fixed volume L 2 = 6 2 and varying temperature β = 4, 6, 8, 10, 12 to scan the (T, µ) plane, and we perform one simulation with L 2 = 10 2 and β = 10 to investigate the finite volume effects.The parameters λ i are typically smooth functions of µ.For the 12 × 6 2 lattice with µ = 0.30 as an example, we found λ 0 = 0.218, λ 1 = −0.126,λ 2 = 0.042.On M S 's, we performed Monte-Carlo calculations generating between 10 2 and 10 8 independent configurations depending on the magnitude of σ .The advantages of using M S over a naive calculation on (S 1 ) N can be ascertained by computing σ .When computed on (S 1 ) N , σ decreases (exponentially) with µ.On M S , σ initially decreases, but near saturation it increases and approaches unity, as can be seen on Fig. 1.This is consistent with expectations due to the discussion of limiting behavior around Eq. (7).
In order to quantify the speedup gained on M S , note that the number of measurements required for a fixed precision scales like σ −2 .Thus the speedup may be estimated by computing σ 2 M S / σ 2 (S 1 ) N .Computing this ratio is difficult however because σ (S 1 ) N is very small at large µ.We therefore estimate the value of σ (S 1 ) N by performing a fit to log σ (see Fig. 1).Using this fit, we can compare the σ at large µ.We find that on a β ×L 2 = 10×6 2 lattice for µ = 0.45, σ 2 M S / σ 2 R N ≈ 10 4 , indicating a sizeable speedup.
Our results for the L 2 = 6 2 lattices are shown in Fig. 1.The distinctive feature is the rapid transition from ψψ 0 to ψψ ≈ 0 as µ increases.As expected on physical grounds, the transition sharpens with lowering T .We present the phase diagram of ψψ in the (T, µ) plane in the right panel of Fig. 2. The heat map is the smooth interpolation of our results based on the fit discussed above.As expected, ψψ ≈ 0 at large values of T or µ.To estimate the location of the transition from a chirally broken to a chirally restored phase we have highlighted the contour at ψψ µ,T = 0.5 ψψ 0 .
A natural question is whether the transition between these two regimes is a true phase transition.Since chiral symmetry is explicitly broken by m f , we do not expect a second order transition line, but a first order transition could exist at small T and large µ.An indication of a true phase transition would be the sharpening of the transition as the volume grows.In the left panel of Fig. 2 we show ψψ as a function of µ for β = 10 and L 2 = 6 2 , 10 2 .The ψψ transition indeed sharpens with L 2 but the data we presently have does not allow a definitive answer on whether this extrapolates to a genuine transition at infinite volume.
In this work, we have extended the sign-optimized manifold method to reduce the finite-density sign problem of a (2 + 1)d field theory.The integration manifold was chosen by maximizing σ over a family of manifolds for which fast hybrid Monte Carlo calculations are possible.The speed at which independent configurations can be collected compensates for the still substantial sign problem on the family of manifolds.Using this method, calculations on lattice sizes up to 10 3 and 12 × 6 2 were feasible.These calculations were enough to outline the broad features of the system's phase diagram.We find a low temperature/density region with a large chiral condensate and a high temperature/density region where the condensate is very small.Investigation of the detailed nature of the phase transition is saved for future work.
It is likely that other manifolds providing a better compromise between speed of calculation and average sign exist and can be found.Greater analytical insight into the geometry of complexified field theories could yield such manifolds.Another direction for future research is the extension of our methods to gauge theories.Although the general idea of changing the domain of integration is shown to be sound [33], suitable manifolds were found only through the computationally expensive method of solving the holomorphic flow equations.

FIG. 1 :
FIG.1: ψψ (left) and σ (right) as a function of µ for β × 6 2 lattices.Notice the increase of σ for large values of µ as expected from the discussion in the text.The black points are σ for simulations on (S 1 ) N on a β = 10 lattice.

25 FIG. 2 :
FIG. 2: Left: ψψ as a function of µ for β = 10 at two different volumes: L 2 = 10 2 and 6 2 .A sharpening of the chiral transition can be seen as the volume is increased.Right: ψψ as a function of T and µ for a spatial volume of size 6 2 .The thick central band indicates the location of ψψ µ,T = 0.5 ψψ 0 and its width represents the statistical error.The peripheral thin lines indicate ψψ µ,T = (0.5 ± 0.05) ψψ 0 to help gauge the sharpness of the transition.