Anatomy of a strong residual sign problem on the thimbles

Using a simple Gaussian-like Ansatz for the phase distribution of a theory with a complex action, we show how the thimble integration for the average phase factor can be plagued by a strong residual sign problem when the phase of the complex integration measure conspires with the constant phase of the integrand along the thimble. This strong sign problem prohibits the accurate computation of the average phase factor when it becomes exponentially small, and causes a strong sensitivity to the parameters describing the phase distribution.


I. INTRODUCTION
In lattice simulations of Quantum Chromodynamics (QCD) at nonzero chemical potential the action in the partition function is generically complex such that standard importance sampling Monte Carlo methods are unusable. One way to circumvent this problem is to apply the well known density of states method to this particular setting by splitting the complex action in its real and imaginary parts, and considering the density of the phase of the complex weights in the partition function generated by their magnitudes [1].
To apply the density of states method, the phase density p(θ) is measured explicitly in the phase quenched ensemble, and then integrated over to compute the average phase factor e iθ = dθ p(θ) e iθ . ( This can then be used to access thermodynamical observables [2] as the full and phase quenched partition function are related by Z full = e iθ Z pq . It is well known that for growing chemical potential, e iθ is exponentially small in the volume of the simulated system, and its computation is plagued by a strong sign problem [3]. If e iθ is computed from the oscillatory integral (1), its accurate determination requires a very precise knowledge of the phase distribution p(θ), which might be obtained with the LLR method [4]. To further improve the accuracy and stability of the integration, a judicious fit to the measured phase distribution is performed before integrating the phase factor [2]. Herein we present a simple example illustrating that the average phase factor obtained from such a fit is neither necessarily stable under slight variations of the fit parameters, nor can it be computed accurately when the sign problem becomes too strong. To compute this oscillatory integral we opted to use the thimble integration, which is often believed to reduce the sign problem and make it manageable. One of our motivations was to investigate the mechanism by which the thimble integration yields exponentially small values for e iθ . * jacques.bloch@ur.de When integrating along a thimble, the magnitude of the complex integrand falls off in a Gaussian like manner at either side of the saddle point. Nevertheless, there is a potential sign problem due to the residual phase along the thimble, which is caused by the phase of the complex measure along the integration path and the constant phase of the integrand on the thimble. Even though it is often claimed in the literature that the sign problem caused by this residual phase is most probably negligible [5], this is not true for the simple, physically motivated example discussed in this paper, as the phase of the complex measure can conspire with the constant nonzero phase of the integrand along the thimble to cause a sign problem that can even be maximally strong for physically relevant parameter values.
Although we applied the thimble integration to the one-dimensional oscillatory integral (1) occurring in the density of states method, the knowledge about the residual sign problem could have a wider scope as the Lefschetz thimbles, which are paths of steepest descent, are also intensively being explored to resolve the sign problem in the high dimensional integration occurring in the partition function of lattice QCD at nonzero chemical potential [5][6][7][8].
In section II we will describe the choice of the phase distribution p(θ) in (1). In section III we will introduce the basic elements needed to perform the thimble integration, and in section IV we will give the results and discuss the residual sign problem. Finally, we close with some conclusions in V.

II. PHASE DISTRIBUTION
In lattice simulation of QCD at nonzero chemical potential the weights in the partition function are generically complex due to the fermion determinant. The distribution p(θ) represents the probability density of the phase of these complex weights in the phase quenched ensemble [1], which is generated with the magnitude of the complex weights.
As θ is a complex phase it would be natural to assume θ ∈ (−π, +π]. When the sign problem is large, the distribution nears a uniform distribution over this interval, and the exponentially small value for e iθ arises from arXiv:1808.00882v1 [hep-lat] 2 Aug 2018 tiny corrections to this uniform distribution, which cannot be determined accurately enough to allow for a useful determination of e iθ .
To improve upon this, it was suggested to consider an extensive phase, which is defined such that θ is no longer bounded and branch cut discontinuities are avoided [9,10]. Claim is that as the physical volume of the system gets larger, the distribution of the extended phase converges to a Gaussian distribution such that a pure Gaussian Ansatz would suffice to compute e iθ and perform phenomenology at nonzero density [9][10][11][12][13][14]. The argumentation also involves the cumulant expansion for e iθ , which, if it converges, is always real and positive, and whose leading term corresponds to the Gaussian result.
The Gaussian Ansatz was however questioned by the observation that higher order corrections in the cumulant expansion, which involve delicate volume cancellations, could invalidate the Gaussian value of e iθ [15][16][17]. The simple example presented below is very much supporting the latter argument, as we found that the cumulant expansion converges extremely slowly for the simple Gaussian-like distribution (2) when the sign problem becomes strong, and that higher order terms can indeed make e iθ orders of magnitudes smaller than its naive Gaussian value. The detailed discussion of the convergence of the cumulant expansion will be discussed elsewhere [18], as we will presently focus on the analysis of the thimble integration.
The distribution p(θ) of the extended phase is generically described by an exponential of an even polynomial in θ, as is discussed in the context of Monte Carlo simulations of QCD at finite density [2,16,17]. Herein we will consider the simplest extension of the Gaussian distribution within that framework, namely an exponential of a quartic polynomial, with normalization factor N = 2 √ a/(σe κ K 1/4 (κ)), where κ = 1/(16a), K 1/4 is a modified Bessel function of fractional order, and a ≥ 0.
This particular functional form for p(θ) is also suggested by detailed studies of the phase distribution in Monte Carlo simulations of a random matrix theory (RMT) that models some important properties of QCD at nonzero chemical potential [19,20] and which will be reported elsewhere [18]. The quartic term is dictated by the tails of the measured phase distributions, which are slightly narrower than those of a normal distribution. The parameters σ and a can unambiguously be extracted from the second and fourth moments of the phase distribution measured in the Monte Carlo simulations.
Although the analysis performed in this paper is valid for any value of σ in (2), the results will all be given for σ = 4.2. We will investigate the behavior of e iθ for small values of a, where the distribution is very close to normal. This choice of parameter values stems from the  RMT simulations for matrix sizes where the sign problem is strong [18]. The phase distribution (2) is illustrated in Fig. 1 for σ = 4.2 and a = 0.0095. For such small a the distribution is almost undistinguishable from a Gaussian, as can be seen in the figure. In Fig. 2 we show the real part of the oscillating integrand in (1) obtained for this distribution.
Although the main focus of this paper will be the computation of (1) using Lefschetz thimbles, in particular in the region where the sign problem is strong, we can also compute this one-dimensional integral using standard numerical quadrature routines, as long as the sign problem remains amenable to such methods. The results for e iθ as a function of a are presented in Fig. 3 for σ = 4.2. When a = 0 we recover the Gaussian result, where the average phase factor can be computed analytically, As a increases, the average phase factor decreases and has a zero crossing at a 0 ≈ 0.00965632. The region a ∈ [0, a 0 ) is especially important in the context of physical simulations with a complex action as e iθ is known to be positive but exponentially small in the volume. In the remainder of this paper we will investigate how such small values can arise in the thimble framework.

III. THIMBLE ANALYSIS
In the thimble formulation [21], the original integral (1) is rewritten as where Ω is the set of relevant thimbles in the complex plane that contribute to the integral. Thimbles are trajectories of constant phase going through saddle points of the integrand. The integral on a thimble T is where we changed notation from θ to z, to emphasize that the variable has been complexified, and parametrized the thimble by its arc length s. The complex measure along the thimble can be rewritten as dz = ds e iη(s) , such that The phase factor e iη(s) is the Jacobian of the arc length parametrization of the thimble, where η(s) is the angle of the tangential to the thimble. As thimbles are trajectories of constant phase φ, we can write f (z(s)) = r(s) e iφ with r(s) = |f (z(s))| such that In this form the equation will be most useful to investigate the strong residual sign problem. Let us first briefly consider a pure Gaussian distribution in (1). It is well known how the thimble construction trivially solves the sign problem in this case. The single saddle point, given by the zero of the derivative of the action, is located at z 0 = iσ 2 and the thimble is parallel with the real axis, with constant phase φ = 0. The integrand, which was strongly oscillating on the real axis, is now replaced by a Gaussian integrand on the thimble. The integral value becomes exponentially small, while avoiding a sign problem, because the function value in the saddle point becomes exponentially small when the integration contour is pushed up in the complex plain.
A salient feature of the Gaussian distribution is that the average phase factor (3) is always positive and not very sensitive to small changes in the width of the distribution. We will see that this is no longer true when generalizing the phase distribution to (2), as a very different thimble mechanism is at work close to the zero crossing a 0 .
To determine the saddle points and thimble trajectories we rewrite the integrand as f (z) = e −S(z) with complex action S. For the integrand in (1) with phase distribution (2) the complex action is with saddle point equation After rewriting the saddle point solutions as z = it, the equation becomes a cubic equation in t, with real coefficients Depending on the value of the discriminant this equation has either three real solutions if ∆ ≥ 0 or one real and two complex conjugate solutions if ∆ < 0. After substituting (11) in ∆, we find that the discriminant is zero when a = a c with critical value When a ≤ a c the quartic term in (2) is small and the distribution becomes more Gaussian-like. In this case, ∆ ≥ 0 and (10) has three real roots [22] with k = 0, 1, 2. The three saddle points z k = it k are located on the imaginary axis.
Although the solutions (14) can be analytically continued to a > a c , where ∆ < 0, it is more revealing in this case to write the real solution t 0 as [23] and the two complex conjugate solutions t ± as solutions of the remainder quadratic equation, obtained by dividing (10) algebraically by t − t 0 . Its complex conjugate solutions are given by There are thus three saddle points: z 0 = it 0 on the imaginary axis, and a complex pair (z, −z * ) = (it + , it − ), located symmetrically left and right of the imaginary axis. Note that the saddle points are explicit functions of the parameters σ and a of the phase distribution (2).
Once the saddle points are known, a further analysis is performed to determine which thimbles are relevant to the thimble integration.
For a ≤ a c the thimble structure is quite similar to that of the Gaussian distribution and only one thimble, going through the saddle point z 1 = it 1 of (14), contributes to the integral. The constant phase φ along the thimble is zero, as can either be computed explicitly from the action in the imaginary saddle point, or can be deduced from the fact that the integral (1) is known to be real.
On the other hand, for a > a c the thimble through the imaginary saddle point does not contribute to the integral, rather, the thimble integration is now given by the sum of the two thimbles T − and T + that are mirrored about the imaginary axis and go through the (z, −z * ) pair of saddle points corresponding to (17). Crucial is that φ is not constrained to zero on the mirrored thimbles, and its nonzero value will cause the strong residual sign problem on the thimbles, as will be discussed in more detail below.

IV. RESULTS
The aim is to understand how exponentially small integral values arise in the thimble framework, and to investigate the impact of small variations in a on e iθ .  Some thimble properties that will be discussed throughout this section are summarized in Table I. All results in this section were obtained for σ = 4.2.
We first look at the constant phase along the relevant thimbles. For this, we substitute the value of the saddle point on the relevant thimble in the action (8) and set φ = −S I , with iS I the imaginary part of the action.
For a ≤ a c the saddle point is purely imaginary and the constant phase is zero. For a > a c the action for the (z, −z * ) pair of saddle points satisfies S(−z * ) = [S(z)] * , such that the constant phase φ has opposite values on both thimbles. The constant phase φ is shown as a function of a in Fig. 4. The phase is zero for a ≤ a c and becomes nonzero when a > a c , increasing steadily as a function of a in the region of interest.
When studying the thimble integration for various values of a in Figs. 5-7, we will display three panes for each parameter value. From left to right: (A) the path(s) of the relevant thimble(s) in the complex plane with their saddle point(s); (B) the magnitude r(s) of f as a function of the arc length parametrization s along the thimble, see (7). As φ is constant along the thimble, the evolution of r(s) tells us how fast the real and imaginary parts of f fall off when moving away from the saddle point along the thimble. In practice, the thimbles were followed until the function dropped to 10 −12 of its maximal value in the saddle point.
(C) the real part of r(s) e i(φ+η(s)) , which is the actual function to be integrated over after parametrizing the thimble by its arc length s, as in (7). This integrand includes the Jacobian e iη(s) of the parametrization, where η(s) is the angle of the tangential to the thimble.
In panes (B) and (C) the saddle point is chosen as origin of the arc length parametrization. When two thimbles contribute to the integral, these two panes only show the functions along the T + thimble, as those on T − are related by complex conjugation. We now analyze the thimble integration as function of the parameter a. For σ = 4.2 the transition between the single-thimble and double-thimble regions happens at a c ≈ 0.0042, in accordance with (13). For a = 0 the distribution (2) is Gaussian. As a increases, with a ≤ a c , the integral is still represented by a single relevant thimble with saddle point on the imaginary axis and constant phase φ = 0. This is illustrated in Fig. 5, where we show the thimble solutions for a = 0.001 (top) and a = 0.004 (bottom). The content of the three panes is as explained above in items (A)-(C). The thimble path in the complex plane is shown on the left in each row. In the middle we show the magnitude r(s) of f as a function of the arc length s, while on the right we show r(s) cos η(s), which is the real part of the product of f (s) with the Jacobian e iη(s) , as φ is zero below a c . The latter is the actual integrand that yields e iθ after integration over s, see (7). When performing the thimble integration, there is a small variation of the phase of the integrand along the thimble due to the Jacobian e iη(s) . This fluctuation is however insignificant, as cos η(s) remains very close to one. This is confirmed by the observation that panes (B) and (C) of Fig. 5 are almost indistinguishable.
When a nears a c from below, two purely imaginary saddle points move toward each other along the imaginary axis, then merge when a = a c , and move apart again as a > a c , leaving the imaginary axis as a (z, −z * ) pair of saddle points. The integral is then represented by the pair of thimbles T − and T + , which are mirrored about the imaginary axis, as is illustrated in Fig. 6A for a = 0.00425, which is just above a c . From the action (8), it is easy to verify that the contributions of the two thimbles to the integral are complex conjugate. Figure 6 also illustrates how the a ≤ a c region smoothly connects to the a > a c region: as the pair of saddle points emerges, the single-thimble splits into two mirrored thimbles. Let us focus on T + : the thimble is almost vertical for s < 0, such that η(s < 0) ≈ −π/2. As the constant phase φ ≈ 0 (see Table I), the real part of the phase factor cos(φ + η(s)) ≈ 0 for negative s. This can be verified by comparing panels (B) and (C) in Fig. 6, showing the integrand on the thimble T + . Although the magnitude r(s) in pane (B) is fairly symmetric about s = 0, it clearly becomes asymmetric in pane (C) after multiplication with the phase factor cos(φ + η(s)) to compute the real part of the integrand, which almost vanishes for negative s.
If we compare the thimble integration below a c and just above it, we see that the Gaussian-like curve on T in Fig. 5C turns into a sum of two half-Gaussians, one on T + , shown in Fig. 6C, and an identical one on T − . We now further increase a, and zoom in on the region a ∈ [0.009, a 0 ], with a 0 ≈ 0.00965632. In Fig. 7 Table I and the plots of η(s), which fall on top of each other in the left pane of Fig. 8. As expected, the magnitude r(s) in Fig. 7B is Gaussianlike and scarcely changes under small variations in a.
When integrating along the thimble the magnitude r(s) gets modulated with cos(φ + η(s)), and two observations should be made. Firstly, even though η(s) is insensitive to small changes of a, it does evolve significantly as a function of s along the thimble, as is shown in the left pane of Fig. 8. Secondly, when increasing a from a c to a 0 , the constant phase φ on T + increases from φ = 0 to φ ≈ 2.3. More specifically, for a varying from top to bottom in Fig. 7, the constant phase φ varies from 2.1 to 2.3, as can be read off from Table I and Fig. 4. These two facts are gathered in the right pane of Fig. 8 where we show the phase φ + η(s) for the values of a investigated in Fig. 7. We observe that the intersection with π/2, where the cosine vanishes, shifts to values of s well inside the relevant region of the thimble integration, and hence the real part of the integrand in (7) oscillates, which generates a residual sign problem. This can clearly be observed in Fig. 7C, as the positive and negative contributions cancel more and more as the parameter a is varied slightly from top to bottom in the figure and gets closer to a 0 . This is also confirmed by the value of e iθ in Table I, which decreases by seven orders of magnitude compared to the Gaussian result for the values of a considered in the table. In fact, the residual sign problem can become arbitrarily strong as a approaches a 0 : the positive and negative contributions to the thimble integration will cancel to higher and higher precision as e iθ becomes smaller and smaller, and eventually goes to zero.
At first sight, this may look like a rather artificial prob-lem, considering that it is merely caused by the zero crossing of e iθ in Fig. 3, and knowing that the sign problem will go away when increasing a further and the integral turns negative. However, this particular thimble mechanism is physically relevant as we know that e iθ is positive for physical systems with a complex action, as it can be written as a ratio of partition functions. Moreover, the parameter values σ and a of the phase distribution (2), obtained from Monte Carlo simulations of RMT at nonzero chemical potential, turn out to be such that e iθ is many orders of magnitudes smaller than its Gaussian value. In this case the thimble integration does not solve the sign problem as the exponentially small values of e iθ are only obtained at the cost of a strong residual sign problem in the thimble integration itself.

V. CONCLUSIONS
The average phase factor e iθ is central in the density of states method used in simulations of physical systems with a complex action. As its measurement is plagued by a sign problem, it is sometimes improved by integrating over a judicious fit to the measured phase distribution. The hope being that e iθ could then be computed accurately and with little sensitivity to the fit parameters.
Although, the Gaussian distribution is a simple candidate for a fit function that allows for an accurate computation of e iθ and is stable under small variations of its fit parameter, this paper has shown that tiny deviations from the Gaussian distribution can readily ruin these properties.
We extended the Gaussian distribution with a quartic term in the exponential, which was suggested by simulations of random matrix theory at nonzero chemical potential. Higher order corrections were not considered as their contributions in these simulations were consistent with zero. Moreover, the simple extension considered in this paper suffices to show the deficiencies of the fitting Ansatz method.
We performed a detailed analysis of the thimble integration to compute e iθ for parameter values that are physically relevant, i.e., with a quartic term that only slightly perturbs the leading order Gaussian part. We showed that there are fundamentally two regions depending on the strength of the quartic term: one where the integral is represented by a single thimble for which the constant phase of the integrand is zero, and another one where a pair of mirrored thimbles make up the thimble integration. In the latter case, the constant phase is nonzero, and, together with the phase of the Jacobian of the thimble parametrization, a strong residual sign problem arises, allowing e iθ to become exponentially small as a nears a 0 .
In Monte Carlo simulations of physical systems with a complex action, e iθ must be known to good relative accuracy to apply the density of states method. However, the strong residual sign problem means that the precise value of e iθ may not be computable if the fit parameter a is close to a 0 , and moreover, small uncertainties on a will lead to relative uncertainties on e iθ which can be of several orders of magnitude.
As endnote, we observe that, although we have only presented results for σ = 4.2, the strong residual sign problem is a general feature of the computation of e iθ over the distribution (2) for any σ and a → a 0 (σ). As σ is increased, which corresponds to larger volumes of the simulated physical system, the distribution looks more and more normal and the Binder cumulant goes to three, however, the strong residual sign problem remains unmoved. A Gaussian value for the Binder cumulant is thus no guarantee for a Gaussian result for e iθ . This will be the topic of a forthcoming paper discussing the phase distribution in a random matrix model of QCD, and its volume dependence [18].