Heavy-dense QCD, sign optimization and Lefschetz thimbles

We study the heavy-dense limit of QCD on the lattice with heavy quarks at high density. The effective three dimensional theory has a sign problem which is alleviated by sign optimization where the path integration domain is deformed in complex space in a way that minimizes the phase oscillations. We simulate the theory via a Hybrid-Monte-Carlo, for different volumes, both to leading order and next-to-next-to leading order in the hopping expansion, and show that sign optimization successfully mitigates the sign problem at large enough volumes where usual re-weighting methods fail. Finally we show that there is a significant overlap between the complex manifold generated by sign optimization and the Lefschetz thimbles associated with the theory.


I. INTRODUCTION
Mapping the phase diagram of Quantum Chromo Dynamics (QCD), the theory that governs the strong nuclear force, has been an outstanding problem for many decades.The main challenge in a nutshell is that such a task requires non-perturbative methods which almost always involve a numerical stochastic, "Monte-Carlo" component where the path integral is sampled statistically.Even though such non-perturbative methods have been very successful in computing the thermodynamic properties of QCD at small densities [1,2], for most theories at finite density, including QCD, their applicability is severely limited by the "signproblem" which arises when path integral measure is not positive definite and therefore cannot be interpreted as a probability measure.In many cases such as QCD at finite density, the underlying action is complex and the Boltzmann weight, e −S , leads to severe phase oscillations in the path integral which become exponentially rapid in the large volume / low temperature limit [3][4][5][6].These oscillations make the numerical computation of the path integral practically impossible.
Over the years, there has been many attempts to tackle the sign problem.One set of ideas stem from complexification of fields where instead of sampling the configurations on the original domain of the path integral (for example SU (N C ) for QCD), one samples on a complexified domain (a subset of SL(N C ) for QCD) where the phase oscillations are milder and therefore can be dealt with using conventional methods.Of course, the deformation to the complexified fields must not change the value of the path integral.An example of a complex domain over which the path integral has milder phase oscillations is the Lefschetz thimble decomposition, which is a multi-dimensional generalization of expressing a one-dimensional integral as a linear combination of steepest descent contours [7,8].Just like steepest descent contours, the phase over each Lefschetz thimble is stationary, rendering the phase oscillations milder [9][10][11].However finding the correct linear combination of thimbles for a given set of parameters in the theory and sampling is generally challenging.To overcome these difficulties various methods have been proposed where instead of the thimbles, one samples over other complex domains which still have milder phase oscillations than the original domain but easier to construct and sample compared to thimbles [12].
A more direct approach to mitigating the sign problem by complexification is to convert it into an optimization the problem where one minimizes the sign oscillations within a family of complex domains parameterized by a set of auxiliary variables.In a nutshell, an observable, such as the average phase, that measures the strength of the sign oscillations is minimized within that set of parameters that describe the complexified domain.The set of ideas that stem from this approach goes under the name "sign optimization" and has been applied to various QCD-like models [13][14][15].More broadly, this idea has also been used to mitigate signal-to-noise problems in gauge theories even in the absence of a sign problem [16] In this paper, we implement sign optimization in the heavy-dense limit of QCD which is captured by an effective theory whose degrees of freedom describe heavy quarks at very high density.Even though this is a somewhat academic limit of QCD, it does inherit the sign problem from QCD which becomes exponentially severe with the volume/inverse temperature.Furthermore the degrees of freedom are elements of SU (N C ) the complexification of which can be generalized to QCD.Furthermore, the simplified mathematical structure of the theory in this limit allows analytical expressions which can be directly compared with the Monte-Carlo results obtained from sign optimized lattice simulations.Finally, because of the same reason, the simulations are not as demanding as QCD in terms of computational resources.All these factors combined make heavy dense QCD an ideal test ground for implementing sign optimization with the goal of applying it to QCD in the future.Along similar lines to this work, it has been studied via Lefschetz thimbles [17] and Complex Langevin [18,19] in the context of the sign problem.
In Ref. [20] we studied the one-dimensional limit of QCD in the context of sign optimization.This paper is a natural continuation of the work along those lines.One crucial difference, however, is that heavy-dense QCD is a three-dimensional effective theory where the sign problem grows exponentially with the volume, allowing us to analyze the performance of sign optimization as a function of volume; which was not possible in one-dimensional QCD.At the same time, as we discuss further in the paper, the existence of complex saddle points and the relations between the Lefschetz thimbles attached to these saddles and the sign optimized manifolds is similar to what we have observed in one-dimensional QCD.
The rest of the paper is organized as follows.In Sections II and III we briefly recapitulate the heavy-dense effective theory of QCD, and sign optimization respectively.Our results are presented in Section IV.We then compare the complexified domains obtained by sign optimization to Lefschetz Thimbles in Section V. Finally we discuss our results and present our conclusions in Section VI.

II. HEAVY DENSE QCD
We consider the cold and dense limit of QCD with heavy quarks [21][22][23].In this limit, the quarks cannot move due their large mass, but do not decouple entirely either, due to their high density.The resulting effective theory is three dimensional whose degrees of freedom are Polyakov loops, P (x), that describe these stationary quarks.More precisely, to leading order in this heavy-dense limit, the Dirac determinant (for Wilson fermions on the lattice) limit reduces to [21] det where P x = Nt x0=1 U 0 (x 0 , x) is the Polyakov loop (in the temporal gauge) that is an element of SU (3) and are the quark and anti-quark fugacities with a, µ and m being the lattice spacing, quark chemical potential and constituent quark mass.It is also useful to define the "hopping parameter" κ = e −ma /2.In the heavy-dense limit ma → ∞, µa → ∞ such that κe µa = finite .
At high densities the quarks have much larger fugacity than anti-quarks, h ≫ h and the second term in Eq. (1) can be neglected, which simply means that the anti-quark contribution to the effective action is negligible.We further take the lattice strong coupling limit where the gauge field dynamics can be neglected as well.This turns out to be a fairly good approximation for g 2 Y M ≳ 1 because as explained in Ref. [22], the effective gauge coupling, λ, vanishes rapidly as λ ∼ (2N c /g 2 Y M ) Nt for an SU (N c ) gauge group with coupling g 2 Y M .Finally, for N c = 3, the Dirac determinant can be expressed in terms of TrP as Putting everything together, to leading order in the heavy-dense limit, the partition function of the effective theory is obtained as Apart from h = 0 (zero density) or h = 1 (half-filling), P x and P † x have different weights, therefore the Dirac determinant is complex and the theory has a sign problem.Furthermore the sign problem grows exponentially with volume, which follows from the fact only non-vanishing group integrals are [DP ] = [DP ]P † P = 1, and the path integral factorizes for each x.
Beyond leading order, the dynamics of fermions can be systematically incorporated in the theory as interaction between the Polyakov loops at different spatial points and can be organized in a hopping expansion controlled by the hopping parameter κ [23].These interactions further simplify in the limit N t ≫ 1, which we assume, and lead to the following next-to-next-to-leading order (NNLO) effective action [23] S Here î, ĵ denote the nearest neighbors of x and the sums over î and ĵ run over each direction (x + 1, x − 1, etc...).The effective coupling that controls the hopping expansion is We have simulated the leading order as well as NNLO theory by using sign optimization.

III. COMPLEXIFICATION AND SIGN OPTIMIZATION
In this section we summarize the idea of sign optimization [13-15, 24, 25], and detail how we implement in heavy dense QCD.As mentioned in the Introduction, our strategy is to complexify the path integral, without changing its value, in a way that the sign oscillations in the complexified domain are milder than the original domain.Let us denote a generic complex space whose shape depends on a set of parameters, ⃗ λ, as M ⃗ λ .Our first requirement is that complexification does not change the value of the path integral: This is ensured by Cauchy's theorem as long as the deformation from SU (3) to M ⃗ λ does not cross any singularities.The strength of the sign oscillations are captured by the average phase, When the sign oscillations are strong |σ ⃗ λ | ≈ 0 due to the phase cancellations and when they are mild |σ ⃗ λ | ≈ 1.
Notice that even though the path integral remains unchanged with the deformation, the average sign does change, since the integrand in the denominator involves a non-holomorphic term ReS.Therefore our goal is to find a set of ⃗ λs which maximizes |σ ⃗ λ |.To do this, we follow a gradient descent trajectory where we start from SU (3) (which corresponds to λ i = 0) and update ⃗ λ according to where τ denotes the gradient ascent step.With an appropriate choice of the step size δ that is determined empirically, this procedure converges to a local maximum of |σ ⃗ λ |.
It is clear that the complexified path integration domain, M ≡ M ⃗ λ V , must have the same dimensions as SU (3) V , meaning it is a middle dimensional complex manifold embedded in SL(3) V .To parameterize M ⃗ λ we first express P ∈ SU (3), with 8 angles a lá Bronzan [26]: The path integral is evaluated on the domain θ i ∈ [0, π/2] for i = 1, 2, 3 and θ i ∈ [0, 2π] for i = 4, . . ., 8. In this representation, the group measure is explicitly written as where H(θ) is the Haar measure, identified with the invariant measure over the SU (3) manifold Here g is the invariant SU (3) metric defined as g ij = Tr P −1 ∂P ∂θi P −1 ∂P ∂θj up to an arbitrary normalization that can be fixed by demanding d 8 θH = 1.The path integral in terms of the θ variables can be written as where We can now express M ⃗ λ in terms of complex θs, which we denote as θ.It is convenient to express the θ in terms of its real part, which we simply call θ, as where f i (θ) are non-singular, real functions that generically depends on a set of parameters ⃗ λ.By construction P ( θ(θ)) is an element of a middle-dimensional manifold M ⃗ λ ⊂ SL(3) since it is still parameterized by eight variables.Furthermore the path integral over this manifold M ⃗ λ is equal to the original path integral over SU (3).This is because θ can be smoothly deformed into θ say by considering a set of intermediate surfaces parameterized by θ i + isf i with s ∈ [0, 1] without crossing any singularities in the integrand.By Cauchy's theorem this deformation does not change the value of Z [24].One advantage of the ansatz Eq. ( 17) is that we can express the path integral over M ⃗ λ using the original real variables θ where is the Jacobian of the transformation, and the effective action is given as We shall consider a class of complex fields generated by the so-called "mixing" ansatz : m,n cos(mθ 4 + nθ 5 ) for i = 4, 5 λ (i) for i = 6, 7, 8 (21) This ansatz was introduced in the study of the one-dimensional QCD and it was observed to perform better than the "diagonal" ansatz where Im θi is only a function of θ i , i.e. f i (θ) ≡ f i (θ i ) [20].The heuristic reason for the better performance is that the mixing between θ i s capture the fluctuations around complex saddle points more accurately than the diagonal ansatz.We will elaborate more on this in Section V when we compare the manifold created via sign optimization to Lefschetz thimbles.
We note that since the only degrees of freedom of the heavy-dense effective theory are Polyakov loops, TrP , we could have used the two independent eigenvalues of TrP to parameterize the path integration domain instead of using the full SU (3) parameterized by eight variables.Instead we choose to work with the full SU (3), even though it is redundant to do so, in order to keep the discussion general and pave the way for more realistic cases where the degrees of freedom are the full SU (3).That said, in this case, one would likely need to work with a more general ansatz for the complexified fields than Eq. ( 21) , which can easily be realized by, for instance, taking into account the Fourier coefficient of all the θ i s.A more detailed analysis of the computational performance of different choices of the anstaze is left for future work.
Being equipped with the representation of the path integral over M ⃗ λ in terms of real θ i given in Eq. ( 18) , it is straightforward to simulate the theory using standard Monte-Carlo techniques by using the effective action given in Eq. (20) .The Markov chain on the complex field space, M ⃗ λ is generated by using real θs that parameterize M ⃗ λ with respect to the probability distribution e −Re S , and the remaining phase is reweighted.In other words the expectation value of an operator is computed as where ⟨. ..⟩Re S denotes the average with respect to the Boltzmann factor e −Re S , and the phase quenched partition function and the average phase are given as We kept the subscript ⃗ λ in the average phase to emphasize that it explicitly depends on ⃗ λ, even though the partition function, and therefore the expectation values of physical quantities do not.We now come back to the problem of choosing the optimal manifold within our ansatz that maximizes σ ⃗ λ and hence alleviates the sign problem.As mentioned earlier, the optimal choice of ⃗ λ is determined via gradient ascent.We start from SU (3) (λ i = 0) and update ⃗ λ according to Eq. (11) .The gradient in Eq. ( 11) can be explicitly calculated as [24] ∇ Therefore each gradient ascent step includes the Monte-Carlo computation of the expectation value above which, notably, can be carried out without a sign problem.We have simulated the theory for various parameters, both to leading order and next-to-next-to leading order using sign optimization as outlined above.In the next section we present our results.

IV. RESULTS
As explained in the previous section, the configurations on the sign optimized manifold can be sampled by using the effective action Eq. ( 20) whose argument are real θ i s, and reweighing the remaining phase, Eq. ( 22) , as usual.We simulated the theory using a Hybrid-Monte-Carlo (HMC) algorithm [27] where we used a standard leapfrog integrator for the evolution of the Hamiltonian h(θ, p) = S(θ) + p 2 /2, and used reflective boundary conditions for θ 1,2,3 [28].We fixed the number of points in the Euclidean time direction to be N t = 100 which corresponds to the cold limit.We simulated the theory both to leading order and to next-to-next-to leading order in hopping expansion.In the former case, we fixed the hopping parameter to  κ = 0.01 and in the latter we varied it between 0 and 0.05.Note that the range of the chemical potential where the fugacity, h, varies between 0 and 1 is roughly |1 − µ| ∼ 1/N t that can be seen from Eq. (2) .
In Fig. 1 (left) we show the average sign as a function of volume for µ = 0.998, κ = 0.01 and N t = 100 (which corresponds to fugacity h = 0.45) to leading order in the hopping expansion.Especially for larger volumes V = 5 3 , and 6 3 the original sign problem is bad enough where reweighting is practically not possible anymore.With sign optimization, it is lifted to values where reweighting is possible.To illustrate this, we plotted the Polyakov loop for the same values on the right.The result after sign optimization agrees with that on the real manifold, albeit with significantly smaller error bars.The improvement in the sign problem as a function of gradient ascent is shown in Fig. 2 where we used the same parameters as above with volume V = 4 3 .Similarly, the value of the average Polyakov loop as a function of the ascent step is shown on the right.Notice that initially the statistical uncertainty is quite substantial due to due to the severe sign problem.Gradient ascent indeed stabilizes the sign which in turn decreases the statistical uncertainty in the physical observable, ⟨P ⟩.
The average sign, equation of state and the Polyakov loop as a function of the chemical potential for volumes 5 3 and 6 3 are shown in Fig 3 .For both volumes without sign optimization, the sign problem is too severe as seen from both directly from the top figures, or the error bars in the physical observables below, density and Polyakov loop.For both volumes sign optimization works as expected by stabilizing the sign and therefore significantly reducing the statistical uncertainty for all values of µ.Note that the Silver Blaze phenomenon [29], where the density sharply rises when the (baryon) chemical potential reaches the baryon mass, is correctly reproduced via sign optimization, in line with previous work on heavy-dense QCD [18,23].In our case this happens at µ = 1, where recall that that the chemical potential is measured in units of constituent quark mass.Finally, in figures 4 and 5 we show the average sign, density, and Polyakov loop as a function of κ, before (right) and after (left) sign optimization to next-to-next-to leading order in hopping expansion.For these simulations, we used the same form of the effective action given in Eq. ( 20) and Eq. ( 16) where in this case S[P (θ)] in Eq. ( 16) is identified with the NNLO action, S N N LO [P (θ)] given in Eq. (7) .Note that for different values of the chemical potential, the sign problem becomes severe for different values of the hopping parameter κ as seen in Fig. 4. As seen from the same figure, sign optimization alleviates the sign problem as expected even in the presence of the hopping terms.This is in contrast with the case where the sign problem can be solved by the worm algorithm only in to leading order (free) in hopping expansion and in parallel with Complex Langevin which also works in either case.
Finally in Fig. 5 we show the dependence of the density and Polyakov loop to the hopping parameter for the same values of the chemical potential as above.Similar with the free case, the statistical uncertainties that arise from the sign problem are significantly reduced by sign optimization.In order to check the validity of our results we compared them with second order perturbation theory.The dashed/solid lines show analytic results to first/second order perturbation theory.These analytic results were calculated by expanding e −S N N LO to first (second) order in η and evaluating the integrals over L x = TrP x and L * x = TrP † x by the help of the relations Tr As expected the agreement between the sign optimized Monte-Carlo computations and the perturbative calculations is fairly good for small values of κ.

V. CONNECTION WITH LEFSCHETZ THIMBLES
Although being computationally very effective, by construction sign optimization does not provide much physical insight in how it alleviates the sign problem.In this section we provide a physical interpretation of the results we presented by making a connection between the sign optimized manifold, M ⃗ λ , and the Lefschetz thimbles associated with the theory.In order to construct the Lefschetz thimbles associated with the path integral, Eq. ( 5) , to compare and contrast with M ⃗ λ , we first express the path integral in terms of the eigenvalues of the Polyakov loop, α 1 , α 2 .We choose to work with αs instead of θs in order to visualize the connection with the thimbles easily.The Polyakov loop in the diagonal form is given as P = diag{e iα1 , e iα2 , e −i(α1+α2) }.The group measure for the α variables is a Vandermonde determinant: such that the partition function is expressed as Here the effective action is defined as and the Polyakov loop takes the form TrP (⃗ α) = e iα1 + e iα2 + e −i(α1+α2) .The Lefschetz thimble decomposition of Z can be constructed by evaluating the holomorphic gradient flow starting from α i,x ∈ [0, 2π] and taking the formal limit τ → ∞.Finite values of τ are associated with a family of complex manifolds that interpolate between the original domain, [0, 2π] 2V , and the thimble decomposition.Following [30], we shall call these manifolds "generalized thimbles".Since there are no interactions between different spatial points in the leading order in hopping expansion given in Eq. ( 5) , the flow equation, Eq. ( 31) can be solved for each spatial point, x, independently.Therefore each generalized thimble is V dimensional direct product of a two dimensional complex manifold.We solve the flow equation to generate this two dimensional manifold for different flow times.
In order to visualize the generalized thimble, we then take a projection on the subspace defined by α 1 +α * 2 = 0. Notably, this quantity is "conserved" under the flow, namely This can be seen by observing that the action, S e , given in Eq. ( 30) is a real function of x and y where α 1 = −α * 2 = x + iy.The conservation equation, Eq. ( 32) , then follows directly from the flow equation, Eq. (31) .We now compare the sign optimized manifold and the generalized thimbles on the projection on the subspace α 1 + α * 2 = 0 denoted by the dashed black line in Fig. 6.Projections of the generalized thimbles as well as the sign optimized manifold on this subspace are one-dimensional, which we show in Fig. 7.
Firstly, it can be seen that the generalized thimbles converge into the the thimble decomposition with increasing flow time as expected.Secondly, there are multiple complex saddle points, hence multiple thimbles, which contribute to the thimble decomposition.There are two saddle points on the subspace that we focus on (those that intersect the dashed line in Fig. 6), but the other saddle points contribute equally as well since Re S e (hence the path integral weight) on each saddle point is equal.Thirdly, the sign optimized manifold approximately reconstructs the thimble, especially around the saddle points where most of the contribution to the path integral comes from.
We would like to point out that in general, sampling multi-model distributions like this could be present challenges as the Markov chain could get stuck in a local minimum.The generalized thimble approach is susceptible to this phenomenon as the relevant support of each thimble in the sampled field space, Re α, shrinks with increasing flow time as pictured in Fig. 8.As a result transitioning from one such region to another becomes more difficult.The sign optimization, however, does not.This follows directly from the form of the parameterization given in Eq. ( 21) as Im α i is expressed as a function of Re α i , so that the size of the relevant region in the sampled field space remains the same [12].Finally, the Jacobian associated with the sign optimization ansatz, Eq. ( 19), is fairly simple and can be analytically calculated, as opposed to the Jacobian associated with holomorphic gradient flow which has to be numerically evaluated, introducing extra computational cost.

VI. CONCLUSIONS
We simulated the heavy-dense effective theory of QCD at finite density by using sign optimization to alleviate the sign problem.The optimal complex domain for the path integral is chosen within a family of anstaze via a gradient ascent algorithm.Building up on our previous observation from one-dimensional QCD, we used the "mixing" ansatz for the complexified manifold which is constructed to capture the fluctuations of fields around the complex saddle points.We have simulated the theory for different values of the hopping parameter, volume and chemical potential both to leading order and next-to-next-to leading order in the hopping expansion.Notably the sign problem in the latter case cannot be solved via worm algorithm as opposed to the former case [31].We observed that the sign optimization method can handle both cases with comparable amount of computational resources.We also demonstrated that sign optimization works successfully for volumes 4 3 , 5 3 , and 6 3 where the sign problem is too severe to overcome with usual reweighing.
Even though sign optimization is demonstrated to be an efficient way to alleviate the sign problem, it is not transparent how it does it; i.e. it does not provide much of physical insight.This is, of course, a generic property of optimization algorithms.We showed that sign optimization in fact in a way reconstructs the Lefschetz thimbles (a multi-dimensional generalization of stationary phase contours) around the complex saddle points of the theory.Similar phenomena has been observed before [14].In this paper we explicitly compared the Lefschetz thimble and the sign optimized manifold and showed that they indeed overlap especially in the vicinity of the saddle points.This comparison further solidifies the observation of the similarity between the Lefschetz thimbles emanating from the complex saddles of the theory and the sign optimized manifold made in one-dimensional QCD [20].As a result, one might view sign optimization is a method to construct an approximation to the Lefschetz thimbles.However, the there are advantages of using sign optimization instead of directly sampling the (generalized) Lefschetz thimbles.Firstly, the parameterization of the sign optimized manifold does not require a numerical computation of the Jacobian as opposed to thimbles as it can be performed analytically, reducing the computational cost.Secondly sampling the sign optimized manifold does not lead to severe multi-modal distributions, which are difficult to sample as opposed to the thimbles.We performed sign optimization within a family of complex manifolds which are parameterized by a fairly generic ansatz that can in principle be applied to an arbitrary theory whose degrees of freedom take value in SU (3).This generality comes with a cost, however.Finding a generic solution to the sign problem is an NP hard problem [32], therefore any generic ansatz, such as we used, will still lead to an exponentially hard problem.Of course, even reducing the exponent lead be practically useful results as we demonstrated here.At the same time, by using certain properties of the underlying theory such as symmetries or perhaps some other knowledge such as complex saddle points, one could consider a more specific ansatz tailored for the specific theory in mind.In this case, there is no a-priori reason for the sign problem to be exponentially hard.Here, we took the first step in relating the sign optimized manifolds to the complex saddle points and the fluctuations around them.The natural next step is to improve the complexification ansatz by using some analytical information based on these observations, which is left for future work.

FIG. 1 .
FIG. 1.The average sign (left) and the Polyakov loop (right) as a function of volume computed on the real plane sign optimized manifold, M λ , with 250 and 1000 gradient ascent steps.The dashed gray line on the left denotes the exact value of the Polyakov loop.

FIG. 5 .
FIG. 5.The density and Polyakov loop as function of the hopping parameter, next-to-next-to leading order in the hopping expansion with (left) and without (right) sign optimization for different values of µ, V = 4 3 and Nt = 100.

FIG. 6 .FIG. 7 .
FIG. 6.The histogram of the eigenvalues of the Polyakov loop averaged over the volume.The green stars denote the saddle points and the dashed line represents the subspace defined by α1 + α * 2 = 0 which we plot the Lefschetz thimble on.

3 -FIG. 8 .
FIG.8.The parameterization of the generalized thimbles (left) and the sign optimized manifold (right) in terms of the real fields (bars).For thimbles, with increasing flow time the relevant region in the real filed space that is sampled shrinks into small, disconnected regions leading to multimodal distribution in contrast to sign optimization.