One-dimensional QCD in thimble regularization

QCD in 0+1 dimensions is numerically solved via thimble regularization. In the context of this toy model, a general formalism is presented for SU(N) theories. The sign problem that the theory displays is a genuine one, stemming from a (quark) chemical potential. Three stationary points are present in the original (real) domain of integration, so that contributions from all the thimbles associated to them are to be taken into account: we show how semiclassical computations can provide hints on the regions of parameter space where this is absolutely crucial. Known analytical results for the chiral condensate and the Polyakov loop are correctly reproduced: this is in particular trivial at high values of the number of flavors N_f. In this regime we notice that the single thimble dominance scenario takes place (the dominant thimble is the one associated to the identity). At low values of N_f computations can be more difficult. It is important to stress that this is not at all a consequence of the original sign problem (not even via the residual phase). The latter is always under control, while accidental, delicate cancelations of contributions coming from different thimbles can be in place in (restricted) regions of the parameter space.


Introduction
The very first proposal of thimble regularization was intended to extend our capabilities to properly define quantum field theories [1]. After a while it was realized that it was a very natural and powerful candidate solution for the sign problem [2,3]. Since then it has attracted quite a lot of attention. By now also other proposals have been put forward which are more or less inspired by the thimble approach, e.g. the holomorphic gradient flow [4] or the idea of combining the latter with the complex Langevin method [5].
All in all, the basic idea underlying thimble regularization amounts to deforming the original domain of integration of a given field theory into a new one, which is made of one or more manifolds. These manifolds live in the complexification of the original domain of integration; in the original formulation (which is the one we adhere to) they are the Lefschetz Thimbles themselves. Thimbles are defined as the union of the Steepest Ascent (SA) paths attached to critical points σ of the (complexified) action and have the same real dimension of the original manifold. Proceeding straight to the field theoretic quantities one is interested in, we denote the thimbles J σ and in a sketchy way we write O = σ n σ e −i S I (pσ) Jσ dz e −S R O e iω σ n σ e −i S I (pσ) Jσ dz e −S R e iω (1) where the z are a shortcut for complexified field configurations and the p σ stand for the configurations which are the stationary points of the action S, the sum formally extending to all of them, even though the n σ can be zero (thus, not all the critical points do contribute). The action is written in terms of a real part S R and of an imaginary part S I . In the previous formula the denominator reconstructs the partition function Z. Notice that a positive measure e −S R is in place and constant phases e −i S I (pσ) have been factored out of the integrals. This is a consequence of the main virtue the thimbles have: the imaginary part of the action stays constant on them. A so-called residual phase e iω is there that accounts for the relative orientation between the canonical complex volume form and the real volume form, characterizing the tangent space of the thimble.
Solving the sign problem via a deformation of the integration domain is conceptually satisfying and the thimble approach is potentially very powerful. However one can not omit difficulties: thimbles are non-trivial manifolds, for which a local characterisation is missing and thus, not surprisingly, devising Monte Carlo methods to sample integrals on thimbles is a delicate issue. Moreover, recent works have stressed how taking into account multiple thimbles can be tricky [4,6,7]. Finally, the final goal of virtually any attempt to solve the sign problem is to eventually attack the study of the QCD phase diagram, and so in the end one struggles to tackle that ultimate goal.
Despite its simplicity, the study at hand addresses virtually all the issues we have just sketched. We will present a numerical study of QCD in 0+1 dimensions, showing that thimble regularization can solve it: analytical results are known (e.g. for the chiral condensate and the Polyakov loop) and those have been obtained via Monte Carlo simulations on thimbles, on a wide range of values for the number of flavor N f , the mass parameter m and the ratio µ/T (results depend on the chemical potential µ via this ratio; T is the temperature). Even if in the end there is no real gauge symmetry in place, the model is a perfect ground to see the thimble formalism for SU (N ) theory at work. Moreover, the sign problem one has to tackle is a genuine one, originating from the (quark) chemical potential via the fermionic determinant. We have already presented preliminary results of this study in [8]; for other, independent work on this subject see also [9].
The paper is organized as follows: in section 2 we discuss 0 + 1 QCD and how to treat it in thimble regularization, in particular enlightening the role of symmetry and discussing the semiclassical approximation; in section 3 we present our results, both by flat, crude Monte Carlo and by importance sampling in the steepest ascents space; finally, in section 4 we present our conclusions.
2 Thimble regularization for 0 + 1-dimensional QCD 2.1 QCD in 0 + 1 dimensions We shall study a lattice formulation of Quantum Chromodynamics in 0+1 dimensions. There is some abuse of terminology in naming the theory at hand QCD. There can be no Yang-Mills action (and no plaquette either) in less than two dimensions. Nevertheless, though much simpler than its 4-dimensional counterpart, this model provides an excellent setting to test the thimble formalism for SU (N ) theories. Moreover, the sign problem is the genuine one as in real QCD, being due to the presence of a (quark) chemical potential. We have staggered fermions on a one-dimensional lattice with (even) N t sites in the temporal direction. The lattice extent is related to the temperature by aN t = 1/T , where a is the lattice spacing. The partition function of the theory for N f degenerate quark flavours of mass m is where D is the lattice staggered Dirac operator andδ ii is the anti-periodic Kronecker delta. Explicitly By an appropriate gauge transformation we can set to 1 all the links {U i } except one. The only remaining link is simply the Polyakov loop U ≡ U Nt , so that we are effectively left with an SU(3) one-link model. We have When N t is even the Dirac determinant can be shown to be equal to the determinant of a 3 × 3 matrix. Namely, the partition function we have to compute is where A = 2 cosh(µ c /T ) and µ c = sinh −1 (m) (from now on, we set a = 1 in all the calculations). As usual, the quark determinant can be turned into an effective action We will be concerned with three main observables. The chiral condensate is the first one The other two are the Polyakov loop Tr U and the anti-Polyakov loop Tr U † = Tr U µ→−µ . The latter two can be related to the quark number density n ≡ T ∂ ∂µ log Z by a relation which takes quite different forms for different values of N f [10]. We will have numbers to compare to, since analytical results for 0 + 1 QCD are available [10,12,13,14]. We also notice that the sign problem has also been solved by means of the so called subset method [10] and complex Langevin [11].
In the following we will be interested in critical points, i.e. stationary points of the action: to each of them a thimble will be attached. In order to write the equations of motion we first of all introduce the Lie derivative Stationary points are now defined as solution of We are also interested in the Hessian There are three critical points {U k = e 2πik/3 1} with k = 0, 1, 2.

Thimble regularization for SU (N ) theory
The basic ingredient we need to construct thimbles for SU (N ) theory were already discussed in [2] and more recently in [15]. Here we recollect the results we need.
First of all, we need to complexify the degrees of freedom. Going to complex fields means SU (N ) U = e ixaT a → e izaT a = e i(xa+iya)T a ∈ SL (N, C) .
We want to stress that In the previous section we found the critical points of the action. We now see how to attach thimbles to them: this is the point at which the theory at hand will reveal itself as simple.
Not having a local gauge symmetry in place, the construction will go on quite smoothly: we simply need to construct the union of the SA paths originating from critical points.
Having defined the Lie derivative (3) we can write the SA equations as It is easy to show that the solutions of these equations display the main properties we expect: the real part of the action is nondecreasing, while the imaginary part stays constant. Namely, since d dτ =∇ a n,μS ∇ a n,μ + ∇ a n,μ S∇ a n,μ we have that ∇ a n,μS ∇ a n,μ S + ∇ a n,μ S∇ a n,μS = ∇S 2 ≥ 0 ∇ a n,μS ∇ a n,μ S − ∇ a n,μ S∇ a n,μS = 0.
Starting from a critical point, we can now reach any point on the thimble by integrating one particular SA equation. This is not the end of the story. In order to construct the tangent space at each point of the thimble, we need to select a basis for the tangent space at the critical point and then transport it along the flow. Lie derivatives obey non-trivial commutation relations from which we can get commutation relations for vectors V ≡ V n,μ,a ∇ a n,μ +V n,μ,a∇ a n,μ Taking V n,μ,c =∇ c n,μS we can derive the equation for transporting a vector V from the critical point to any point along the flow described by (4)

Takagi vectors at critical points
Takagi's factorization provides the characterization of the thimble in the vicinity of the critical point p σ (with coordinates z σ ). We introduce the vector notation and expand the action to second order around z σ where we have assumed z σ = 0 for the sake of simplicity. Takagi's factorization theorem states that, given the complex symmetric matrix H(S; p σ ) (the Hessian, in our case), there exists a unitary n × n matrix W such that W T H(S; p σ )W = Λ, with Λ = diag (λ 1 , · · · , λ n ) and the λ i (called Takagi values) are all real and non-negative. We will find that in the case at hand they are all positive. The columns of W are n normalized Takagi vectors v (i) , that is so that we can rephrase Takagi's theorem as or, equivalently, Takagi's vectors provides a basis for the tangent space at the critical point. This also mean that each SA leaves the critical point along a direction which is a given linear combination of Takagi's vectors. As we will see, our preferred way of singling out one particular point on the thimble is indeed simple: we choose such a direction and then specify when to stop while integrating the SA.
We saw we have three critical points: We thus have one single Takagi value λ k , while the 8 Takagi vectors are thus recognized as v

Reflection symmetry and its consequences
We have already stated that a contribution from each critical point is expected: this is a direct consequence of the fact that they are all sitting on the original domain of integration 1 . Collecting contributions of more than one thimble to solve the theory is in general a delicate issue. We now see that symmetries can play a major role, providing useful checks for results and even making our life simpler. These observations were pointed out in [16].
We will now show that the action of 0+1 QCD fulfills a reflection symmetry: This ensures the reality of the partition function (and of the expectation value of the Polyakov loop as well). This symmetry of the theory is manifestly fulfilled by the decomposition in thimbles and holds at every order in perturbation theory as well, so we shall recover it in the semiclassical expansion. Consider the QCD partition function The action (in our case, the Dirac determinant is the only component) is invariant under charge conjugation C defined by Thus, we can employ charge conjugation to substitute det D(U, µ) → det D(Ū , −µ) leaving the action invariant. We also recall the generalization of γ 5 -hermiticity at finite chemical potential This implies that We have shown that the aforementioned reflection symmetry is fulfilled and thus we expect thimbles to appear in conjugate pairs. This is indeed the case: consider the three critical points {U k }. U 0 = 1 is real and therefore self-conjugate; the consequence of this is that computations on the associated thimble yield real results. As for the other two critical points, being e 4πi/3 = e −2πi/3 , we immediately see that U 2 = U 1 . This implies that U 1 and U 2 form a conjugate pair of critical points and results of integration on U 2 should be the complex conjugate of those on U 1 , yielding an overall real contribution to the partition function (and also to the expectation value of observables). The chiral condensate and the quark number density automatically respect this symmetry, being derivatives of the partition function. Notice that the symmetry holds for the Polyakov loop (it is obvious, since Tr U = TrŪ ) and anti-Polyakov loop as well. All this is well evident in numerical simulations: in figure 1 we show an example in the case of the chiral condensate. In all the numerical results that we present in the following we take advantage of this symmetry: results from thimble 1 and 2 are averaged.

Semiclassical expansion
In this section we will compute semiclassical expansions around thimbles in 0 + 1 QCD. In general if in (6) we change variables according to Z = W η, the action in the η variables becomes After setting Z = W η, so that d n z = det W d n η = e i ωσ d n η, the partition function becomes In our case, we have det Λ k = λ 8 k and e i ω k = e −i ϕ k /2 8 , so that (8) reads It is obvious that, in general, [H(S; p σ ), H O σ ] = 0 and therefore we cannot expect W to "diagonalize" both H(S; p σ ) and H O σ . The expectation value of the observable is given by The first term in the expansion of O comes out of the integral, giving the partition function itself (actually, the contribution of the thimble J σ ). The second term is linear in η i and therefore gives no contribution to the Gaussian integral. For the same reason, the third term contributes only with terms which are quadratic in η i , that is when Therefore we need to compute only the diagonal terms of C O σ and, after performing a Gaussian integral we arrive at The expectation value of the Polyakov loop can be computed starting from expression (9) Being Tr U k = e 2πik/3 Tr 1 = 3 e 2πik/3 , it follows that From the previous considerations on the reflection symmetry featured by 0 + 1 QCD, we can see that reality of Z and Tr U is achieved by setting n 1 = n 2 . This is so since the contribution of J 2 to Z and Tr U is the complex conjugate of the contribution of J 1 . This is manifest in the semiclassical expansion thanks to S(U 2 ) = S(U 1 ), B 2 =B 1 , λ 2 = λ 1 , e i ϕ 2 = e i ϕ 1 = e −i ϕ 1 , all following from e 4πi/3 = e −2πi/3 = e 2πi/3 . Thus we can rephrase Z as The semiclassical expansion on thimbles also provides an easy way to compute an estimate for the relevance of J 1,2 with respect to J 0 in the computation of e.g. the partition function. We define the relative weight r 1,2 and study it at different values of µ T and m. This, as we shall see, provides a reliable estimate which can be compared with the results of numerical simulations. We note that, being B 0 = A + 2 cosh(µ/T ) ∈ R and B 1 = A − cosh(µ/T ) + i √ 3 sinh(µ/T ) for any value of µ T and m (the ratio λ 1,2 /λ 0 is independent on N f ). One thus expects that integrating only over J 0 will give more accurate results at large N f , i.e. there is a regime where the leading thimble dominance scenario actually shows up. This is of course a semiclassical estimate: the reliability of this prediction will be checked against numerical simulations.

Monte Carlo computations on thimbles
Our preferred way of characterising points on a thimble goes through a constructive approach, which we now recall in the formalism which is valid for any generic theory. Given a critical point, we saw in subsection 2.3 how to determine the tangent space. By performing the Takagi factorization of the Hessian we were left with Takagi values λ i > 0 and Takagi vectors v (i) , which provide a basis for the tangent space. The tangent space contains all the directions along which the SA paths defined by 2 leave the critical point. If we impose a normalization condition It is thus quite natural to single out any given point on a thimble by the correspondence with S n−1 R the (n − 1)-sphere of radius √ R. In [17] we made use of this approach to solve a Chiral Random Matrix Model by means of thimble regularization. We now recall how to make use of (12) to rephrase the integrals in the thimble decomposition of the path integral. A change of variables in the integrals in (1) will be involved. Let us first of all define In the previous formula d n y stands for the real volume form on the thimble J σ (i.e. the thimble attached to critical point p σ ). With a slight abuse of terminology we will refer to this expression as a partition function, which can be rewritten 2 We denote by t the time coordinate parametrizing the flow along the SA path.
in terms of the measure over S n−1 R Dn ≡ n k=1 dn k δ | n| 2 − R and the partial partition functions The partition function Z σ has been decomposed in contributions Z where the λ At the same time, the phase of det V σ (n, t) provides the residual phase e iω(n,t) . We can now go back to (1) and rewrite it in terms of the (n, t) variables. Notice that at this point Z σ will be irrelevant: it will reappear later on. All in all we have Before proceeding we make a couple of observations: • The expression in (16) is not the same appearing in [17]: the two are equivalent 3 .
• The notation σ could look generic: it reflects the fact that (12)

Simulations by flat, crude Monte Carlo
We can make use of flat, crude Monte Carlo to compute the integrals in (18). The recipe is very simple • We pick up randomly (with flat distribution) a directionn.
• Since we want to compute the contribution coming from the SA leaving the critical point p σ alongn, we prepare convenient initial conditions both for the field and for the tangent space basis vectors for such a SA. We can do this, since near the critical point solutions of (11) are know as 4 i t which we can compute for t = t 0 0.
• We then integrate the SA equations for the field and the equations for transporting the basis vectors all the way up till we reconstruct the +∞ −∞ dt integrals appearing in (18) 5 .
In Figures 2 and 3 we display results obtained following this procedure for the chiral condensate and the Polyakov loop. We cover a range of values for N f , m and µ/T . We used this very same flat Monte Carlo prescription to solve the Chiral Random Matrix model in [17], but here there is a noticeable difference. In [17] the contribution from one single thimble was needed. Here, as we said, all the three critical points belong to the original domain of integration (SU (3)) and all of them are in principle relevant in the thimble decomposition. However, the semiclassical arguments of 2.5 provide a deeper insight with respect to the actual weight of each contribution entering such decomposition. Figure 4 (first line) depicts r 1,2 0 (defined in (10)) as a function of µ T and m; by studying this quantity one can predict for which values of the parameters ( µ T , m) integration only over J 0 is expected to capture substantially correct results. On the other side (second and third line) figure 4 also shows that there are regions in which taking all the thimbles into account is compelling: for m = 1, N f = 2 and m = 0.1, N f = 6 the semiclassical evaluation of r 1,2 0 predicts that taking only J 0 into account is going to miss the correct result. This is indeed confirmed by the evaluation of the Polyakov loop. We point out that (as predicted by the semiclassical approximation) at large values of N f (e.g. N f = 12) the contribution from the thimble attached to U 0 essentially captures the correct results: a fortiori for even higher values of N f , the single thimble dominance scenario indeed holds true.
The reader will notice that (in particular in figure 2

Simulations by importance sampling
We now discuss a method to perform importance sampling on thimbles. For the sake of simplicity, we start in a simplified setting, i.e. as if only one thimble contributed. This assumption will be later released to make contact with the case at hand: in the meantime, this assumption allows a simplified notation (for the sake of notational simplicity we will often omit in the following the subscript/superscript σ, e.g. in Takagi values: there is no need to distinguish since it is assumed that only one critical point does matter). In the simplified framework of a single thimble contributing, the computation of (1) simply amounts to where a reweighting with respect to the critical phase is in place and we introduced the notation The reader will recognize the expression for Z σ (13) in the denominator. Making use of the representation (12), and thus of the same notation in which we wrote (14) and (15), we can  We proceed as follow. In our Markov chain we start from the current configuration (which is associated to a directionn) and we propose a new one (associated to a directionn ).n is identical ton apart from two randomly chosen components, say (n i , n j ) with i = j. We define C by C ≡ n 2 i + n 2 j = R − k =i,j n 2 k which is fixed by the normalization | n| = √ R and by the values of all {n k } k =i,j . There is a coordinate system in which we can now parametrize the current values of (n i , n j ) by with φ ∈ [0, 2π). Our evolution step now amounts to change φ → φ , which results in (n i , n j ) → (n i , n j ), while for all the other components (k = i, j) n k = n k . φ − φ is extracted flat in a given (tunable) range. We finally accept the proposed configuration with the standard Metropolis test P acc n n = min 1, Zn Zn .
We actually have a more efficient Monte Carlo (which has been preliminary described in [8] and will be further discussed elsewhere [18]): we notice that this is not applicable here, due to the full degeneracy of the gaussian spectrum of the theory at hand.
In our case three contributions should be in principle taken into account. Actually, due to the symmetry of Section 2.
in which notations should be evident (e.g. Z 0 is associated to thimble J 0 and Z 12 is the average of the contributions associated to thimbles J 1 and J 2 ). In order to proceed we now put at work a recipe that we put forward in [17]. Equation (22) can be rewritten where we defined The idea is now to determine the value of α taking a given observable as a normalization point: all the other observables of the theory can then be computed using this input. This can indeed be done. As an example, at N f = 1, m = 0.1, µ/T = 0.5 we get α = 0.2686 (13) and α = 0.2682(8) (taking different observables as normalization points). Figure 5 confirms the effectiveness of the procedure. For N f = 1, m = 0.1 the point µ/T = 2.0 was completely out of reach for flat, crude Monte Carlo, while results are successfully computed with the improve metod (first row). Notice that for these values of parameters the tiny value of the chiral condensate (of order 10 −2 ) results from a delicate cancelation of the contribution coming from the different thimbles. While this could be seen as a sign problem coming back, one should notice that it is a numerical accident occurring for a given observable at a given    value of parameters. On the other side, it has nothing to do with the original sign problem (nor e.g. with the residual phase). Figure 5 also shows (second row) how errorbars can be cut down in the case of Polyakov loop for N f = 1, m = 1, µ/T = 3.0.

Conclusions and prospects
We showed that QCD in 0+1 dimensions can be effectively computed in thimble regularization. The thimbles attached to three critical points are in principle relevant, but a symmetry relates the results coming from two of them. Importance sampling in the space of steepest ascents enabled us to successfully compute observables in regions where a flat, crude Monte Carlo was failing. Nevertheless, the latter has performed reasonably well on a quite extended region of parameter space.
The sign problem which is place in 0+1 dim QCD is a genuine one, stemming from the Dirac determinant. On the other side, the theory at hand does not show the subtleties of gauge symmetry: the task in front of us is now to proceed to tackle gauge theories.