Lefschetz thimble structure in one-dimensional lattice Thirring model at finite density

We investigate Lefschetz thimble structure of the complexified path-integration in the one-dimensional lattice massive Thirring model with finite chemical potential. The lattice model is formulated with staggered fermions and a compact auxiliary vector boson (a link field), and the whole set of the critical points (the complex saddle points) are sorted out, where each critical point turns out to be in a one-to-one correspondence with a singular point of the effective action (or a zero point of the fermion determinant). For a subset of critical point solutions in the uniform-field subspace, we examine the upward and downward cycles and the Stokes phenomenon with varying the chemical potential, and we identify the intersection numbers to determine the thimbles contributing to the path-integration of the partition function. We show that the original integration path becomes equivalent to a single Lefschetz thimble at small and large chemical potentials, while in the crossover region multi thimbles must contribute to the path integration. Finally, reducing the model to a uniform field space, we study the relative importance of multiple thimble contributions and their behavior toward continuum and low-temperature limits quantitatively, and see how the rapid crossover behavior is recovered by adding the multi thimble contributions at low temperatures. Those findings will be useful for performing Monte-Carlo simulations on the Lefschetz thimbles.


Introduction
The sign problem is the longstanding obstacle which prevents us from applying nonperturbative lattice simulations directly to the physical systems with complex actions, including quantum chromodynamics (QCD) at finite baryon chemical potential µ. The fermion determinant at finite µ becomes complex, which invalidates the importance sampling algorithm. In contrast, the determinant is real at finite temperature (T ) with µ = 0, and lattice simulations of QCD have proved now to be a reliable nonperturbative method to evaluate (e.g.) the equation of state of strongly interacting matter. Nonetheless, studies of QCD-inspired models at finite T and µ have suggested a variety of phase changes from nuclear liquid-vapor transition, to chiral symmetry restoration, and to color-superconducting phase transition, etc. With this situation, in order to unveil the QCD phase diagram from the first principles, many attempts have been made to circumvent the sign problem in lattice QCD simulations, although the complete resolution is still not available [1].
To study the physical systems with complex actions, two alternative approaches have attracted much attention recently -complex Langevin equation [2][3][4] and Lefschetz thimble integration [5][6][7], both of which involve complexification of the dynamical field variables.
Statistical sampling with the complex Langevin equation has been applied to various models , including the massive Thirring model with chemical potential [21,22], as testing grounds, and it is successful in some cases but not in other cases. A formal proof for the correctness of the method has been elaborated under certain conditions [13,17], but full justification of the complex Langevin approach is not established, where logarithm terms such as the ferimon determinant in the action cause a subtlety [44]. Noteworthily, complex Langevin simulations have been applied to full QCD at finite T and µ [26, 30, 34-36, 38, 40, 45, 47], showing consistent results with those obtained by the reweighting method in the parameter region where both methods are stable [47].
Path integration on the Lefschetz thimbles was introduced in the study of analytic property of gauge theories [5], and it was soon recognized as a mathematically sound way to resolve the sign problem [48][49][50]. It can be regarded as a functional generalization of the steepest descent method of complex analysis. In this approach the original integration cycle is deformed to a sum of the curved manifolds, called Lefschetz thimbles, in the complexified field space. On a thimble the imaginary part of the action ImS is constant, and this property allows importance sampling with the weight e −ReS ≥ 0. This advantage was first applied to numerical simulations for 4-dimensional λφ 4 theory with chemical potential with use of Langevin [49] and hybrid Monte Carlo (HMC) [50] algorithms on a single thimble, and successfully reproduced the known results including the so-called Silver Blaze behavior [51] -complete insensitivity of the system to µ below a certain critical value at T = 0. The residual phase problem from the Jacobian due to the curvature is mild and can be efficiently taken into account by reweighting for this theory [50]. The Lefschetz thimble integration has been examined in other models [52][53][54][55] and has been studied from other aspects [56][57][58][59][60][61][62] which involve the sign problem. This approach also shed new light on the complex Langevin sampling method [27,33,62], and vice versa [61].
In this paper, we study the path integration on the Lefschetz thimbles in the (0+1) dimensional massive Thirring model at finite chemical potential µ [63], in order to clarify the effects of the fermion determinant on the structure of the thimbles contributing to the partition function [55]. The lattice model is formulated with the staggered fermions [64,65] and a compact auxiliary vector boson (a link field). This model shows a crossover transition from the low to the high density phase at finite T as a function of µ, and the transition becomes first order in T = 0 limit. Furthermore the exact solution of this model is available on the finite lattice as well as in the continuum limit, and therefore one can assess the validity of the approach precisely by comparing the results with the exact ones.
The fermion determinant has zero points on the complexified field space and actually those zeros form continuous submanifolds on which the effective action S becomes singular. At the same time, the determinant brings in many critical points, each of which a thimble is associated to. We classify the critical points into subsets according to the subspaces they belong to, and identify all the critical points and thimbles in each subspace by noting a one-to-one correspondence between a critical point and a zero point of the determinent. The thimbles whose critical points are located in the uniform-field subspace are shown to dominate the integral toward the continuum limit. Hence we study within the uniform-field subspace how the set of the contributing thimbles to the partition function changes via the Stokes phenomenon as the chemical potential µ varies. We will see that in the crossover region multithimble contributions are inevitable, and become more significant for small inverse coupling and/or in low temperature limit. We study this interplay in more detail by reducing the model degrees of freedom to the uniform-field subspace and show how the crossover behavior is reproduced as adding the multi-thimble contributions to the observables. This paper is organized as follows. In section 2 we introduce the (0+1) dimensional massive Thirring model with chemical potential on the lattice in terms of the staggered fermions and the compact auxiliary vector boson. In section 3, after a briefly review of the Lefschetz thimble approach, we study the critical points and determinant zeros of the lattice model in the complexified field space. The critical points are classified by the subspaces they live, and all the critical points are identified. In section 4, we study the thimble structure of the model at µ = 0, discuss the importance of each thimble by looking at the relative weight at µ = 0. In section 5, we show within the uniform-field subspace the change of the thimble structure with increasing the chemical potential µ via the Stokes phenomenon, and show that the multiple thimbles contribute to the partition function in the crossover region. In section 6, taking the uniform-field subspace, we examine the validity of the single thimble approximation, and discuss the continuum and low temperature limits. Especially in the low temperature limit, the importance of the multi-thimble contributions are clarified. Section 7 is devoted to summary and discussions. The exact solution of the model is derived in Appendix A.

one-dimensional massive Thirring model on the lattice
The (0+1)-dimensional lattice Thirring model we consider in this paper is defined by the following action [21,22,63], where β = (2g 2 a) −1 , ma, µa are the inverse coupling, mass and chemical potential in the lattice unit, and L is the lattice size which defines the inverse temperature as 1/T ≡ La. The fermion field has N f flavors and satisfies the anti-periodic boundary conditions: The partition function of this lattice model is defined by the path-integration, where D denotes the lattice Dirac operator, The functional determinant of D can be evaluated explicitly (see Appendix A for derivation) as This fact can cause the sign problem in Monte Carlo simulations. This lattice model is exactly solvable in the following sense. The path-integration over the field A n can be done explicitly and the exact expression of the partition function is obtained where I 0,1 (β) are the modified Bessel functions of the first kind. The number density and scalar condensate of the fermion field are then obtained as The µ-dependence of these quantities are shown in Fig. 1 for L = 8, ma = 1, and β = 1, 3, and 6. It shows a crossover behavior in the chemical potentialμ (in the lattice unit) around µ m + ln(I 0 (β)/I 1 (β)).
The continuum limit (a → 0) of this lattice model at finite T may be defined as   In this limit the partition function scales as T cosh m T , (2.9) and the continuum limits of n and χχ are obtained as follows: (2.10) From these results, one sees that the model shows a crossover behavior in the chemical potential µ at non-zero temperatures T > 0, while in the zero temperature T = 0 limit it shows a firstorder transition at the critical chemical potential |µ c | = m + g 2 . See Fig. 2.

Preliminaries
Now we consider the complexification of the Thirring model on the lattic and reformulate the defining path-integration of Eq. (2.2) by the integration over Lefschetz thimbles. In the complexification, the field variables A n are extended to complex variables z n (∈ C L ) and the action is extended to a holomorphic function given by S[z] = β L n=1 (1−cos z n )−ln det D[z]. 1 For each critical point z = σ given by the stationary condition,

∂S[z]
∂z n z=σ = 0 (n = 1, · · · , L), the thimble J σ is defined as a union of all the (downward) gradient flow curves determined by The thimble so defined is an L-dimensional real submanifold in C L . Then, according to Picard-Lefschetz theory (complexified Morse theory), the original path-integration region C R ≡ [−π, π] L can be replaced with a set of Lefschetz thimbles 2 , where n σ stands for the intersection number between C R and the dual submanifold K σ , which is another L-dimensional real submanifold associated to the same critical point σ and is defined as a union of all the gradient flow curves s.t. z(+∞) = σ. With denoting the set of the critical points as Σ ≡ {σ}, the partition function and the correlation functions of the lattice model can be expressed by the formulas 3 The functional measure D[z] along the thimble J σ is specified as d L z Jσ = d L (δξ) det U z by the orthonormal basis of tangent vectors {U α z |(α = 1, · · · , L)} which span the tangent space as δz = U α z δξ α (δz ∈ C L , δξ ∈ R L ).
The integration on each Lefschetz thimble is convergent because the real part of the action increases monotonically to ∞ while the imaginary part stays constant along the downward flow, The sign problem remains in the Lefschetz thimble approach in two facts. First, it seems that when we factor out the complex weight e −S[σ] , the integrand of each thimble, e −(S[z]−S[σ]) > 0, is real positive. But a complex phase appears from the Jacobian factor det U z in the integration, which is called residual sign problem. For λφ 4 theory it is demonstrated that the residual sign problem can be treated by the reweighting method safely [50]. Second, the terms Z σ and O[z] σ are actually complex quantities although the total averages Z and O[z] should be real. If there is a certain symmetry in the thimble structure of the system, one can show the cancellation of the phases in the sum [58]. The multi-thimble contributions to the partition function and observables will be more elaborated in this paper.

Critical points and determinant zero points
Given the above mathematical results, however, it is not straightforward to work out for general fermionic models all the critical points Σ = {σ}, the thimbles {J σ }, and their intersection numbers {n σ }. Fortunately in our lattice model, we can find all the critical points determined by the stationary condition Eq. (3.1). The critical point condition for the Thirring model is written as The key observation is that the second term depends on the field configuration only through the sum s, so that all sin z n (n = 0, · · · , L − 1) of a critical point σ must have the same value to cancel the common second term. Let us denote it as sin z, then the field components can be either z n = z or π − z and the sum s is written as where n ± are the numbers of z and π − z in the components {z n } with n + + n − = L. The critical point condition for z is now explicitly written as This can be regarded as the critical point condition for a one-variable model; The case of n − = 0 corresponds to a uniform field configuration, where z n = z (n = 0, · · · , L− 1), and the case n − = 1 means that there is one flipped component π − z, · · · , etc. In the case of n − = L/2, the second term of (3.9) becomes independent of z. The case of n − > L/2 gives the same critical points as in the case L − n − with z ↔ π − z. Hence we need to consider n − = 0, · · · , L/2 − 1.
Thus we have classified the critical points with index n − . By solving the condition Eq. (3.9) of the one-variable model for each n − , we can locate all the critical points of the model. Note that a critical point is L n − -ply degenerated for n − = 0 due to the combination about which components to be flipped.
One of the distinctive features of fermionic theories is the fact that the fermion determinant has many zero points within a compact domain in the complexified space. The real part of the effective action ReS diverges at these zeros, and therefore the downward cycle J σ may flow into one of these zeros, otherwise it must extend outward to the safe exterior region where ReS = +∞. Hence, in addition to the critical points, we need to locate all the zeros of the fermion determinant Thanks to the concise expression of det D[z] in Eq. (2.4), one can easily find all zero points: This only fixes s = L−1 =0 z , and thus defines submanifolds with the complex dimension L−1, embedded in the L dimensional complexified configuration space. Note that these zero points are independent of β, and that nonzeroμ simply shifts the zero points along the imaginary axis. Restricting this submanifold of the zeros in the subspace n − = 0, where s = Lz, we find 2L isolated zeros of while in the subspace (n − = 1) with a single link flipped to π−z (and thus s = (L−2n − )z+π), we have 2(L − 2) zeros of (3.14) Figure 3 shows two sections of the gradient flows in the uniform-field subspace (n − = 0; left) and in the subspace with one link flipped (n − = 1; right) of the model with L = 4, β = 3 − 0.1i and ma = 1 atμ = 0. (The reason for complex β will be explained in the next section.) Globally, the flows are streaming out of the remote points z = ±i∞ and flowing away towards the safe remote points z = ±π ± i∞. We solve Eq. (3.9) numerically, and find ten (eight) critical points 4 for n − = 0 (1), as shown with green dots in Fig. 3. For later convenience, we have numbered the critical points as shown here. We also put the zero points  4 Thimble structure at µ = 0

Bosonic theory
It would be instructive to start our discussion with the bosonic theory without fermions, , whose complexified configuration space is a direct product of (S 1 × R) L . The downward flow is simply given by which is depicted in Fig. 4 for a certain z n with β = 3 − 0.1i. Let us focus on this complex z n plane for a moment. The action is periodic in the direction of the real axis, so that the configuration space is equivalent to S 1 × R, a cylinder. There are two critical points, z n = 0 and ±π (shown in green dots), 5 corresponding to the Gaussian and doubler solutions, respectively. The downward cycle (thimble) J 0 associated to z n = 0, extends to the "safe" exterior regions toward z n = ±π ± i∞ depicted with light-red color at the corners. There is another thimble J −π associated to the doubler solution z n = −π, which connects these two "safe" regions vertically along the imaginary direction. In other words, the two safe regions are connected by two cycles with and without winding around the cylinder. These cycles constitute the base of homology of this restricted space S 1 × R.
We notice here that the original integration path from z n = −π to π is ill-defined as a homology cycle. A well-defined downward cycle should extend to a "safe" region where the Morse function (h = −ReS) approaches −∞ [5]. Actually, the thimble J 0 coincides with this original path only for real β, which is the very parameter for the Stokes phenomenon to occur between z n = 0 and ±π (the action is real at both points; S n = 0 and 2β). Hence in Fig. 4 we have added nonzero imaginary part to the coupling β = 3 − 0.1i 6 to make the thimble J 0 well-defined.
Thanks to the periodicity of the action S n (z), we can exptend the original integration path without changing the value of Z to a U-shaped integration cycle which starts at z n = −π + i∞ and comes down along the imaginary direction to z n = −π then moves along the real axis to z n = π, and goes up to z n = π + i∞ 7 . This U-shaped cycle (which we simply denote with C) is equivalent to the sum of the two thimbles: (4.2) Here we set the orientation of the thimbles so that "+" sign is appropriate here. One can confirm that both the upward cycles K 0 and K −π intersect this integration cycle C.
There are 2 L critical points in the (0+1) dimensional bosonic theory with L lattice sites from combinatorics, and its thimble structure is obtained as a direct product of the thimbles J 0 and J −π . The integration cycle equivalent to the original integration path is symbolically written as The safe exterior region where the real part of the action ReS diverges has the complex dimension (L − 1) because it is characterized by the condition L−1 n=0 (1 − cos z n ) = ∞, i.e., at least one of {z n } is fixed to π ± i∞. We comment on the continuum limit (β → ∞ with fixed β/L). In this limit the contribution to the partition function from each variable becomes Gaussian: The integration along the vertical path, which we have added to make the integration cycle well-defined, becomes irrelevant giving only a contribution which is exponentially suppressed. For example, Thus we see that the doubler contribution J −π is suppressed by e −2β and the free theory with L degrees of freedom is correctly reproduced by the integration on the thimble J L 0 .

Thirring model
We have already idintified the critical points and determinant zeros of the Thirring model in the previous section and shown them in Fig. 3. There we also noticed a certain correlation between a critical point and a determinant zero. Now let us look at the thimble structure of the model atμ = 0.
In the uniform-field (n − = 0) subspace shown in Fig. 3 (left), the thimble J σ 0 extends from one safe remote z = −π − i∞ to another safe remote z = π + i∞ passing through the critical point σ 0 at the origin, and the U-shaped cycle is still equivalent to the sum of the two thimbles, associated to the Gaussian and doubler critical points: where " c ∼" indicates the equivalence as the integration cycles under the constraint n − = 0. In the subspace of n − = 1 (Fig. 3 (right)), on the other hand, the critical point σ 0 contains one doubler component π − z, and the thimble J σ 0 ends at determinant zeros 8 . The U-shaped cycle within n − = 1 space is covered by the sum of four thimbles: with two J σ2 contributions canceling out in the end. The strong correlation between a critical point and a zero point may be expected by noticing the fact that because a zero point z zero is a simple pole of the flow field, one can always find in its vicinity the point on which the first term of Eq. (3.9) can be counter-balanced by the would-be pole contribution, especially when β is large.  Table 1. ReS at the critical points σ i (i = 0, 1, 2, 3, 4,0) with L = 4k, β = 3k, and ma = 1/k (k = 1, 2, 4) at µ = 0. The rightmost column shows the difference of ReS between σ 0 and σ0.
One can also understand the paring between them by considering the thimble structure of the one-variable model assigned by n − , where a thimble J σ becomes a line segment associated to a critical point σ and connects between the zeros and/or safe remote points z = ±π±i∞. In n − = 0 case, for example, two safe remote points are connected by the two thimbles, J σ 0 and J σ0 . Because the thimbles form the basis of independent cycles, no trivial loops are allowed. That is, a set of thimbles must be a connected skeleton graph on S 1 × R subspace assigned by n − . One can add a new critical point, which is accompanied by a new thimble, only when one has a new zero point. Hence the number of thimbles coincides with the number of the critical points, and furthermore with the number of the end-point zeros (including two safe remote points) in our model.
The formulas (3.13) and (3.14) with L = 4 give us eight and four zeros for n − = 0 and 1, respectively, as seen in Fig. 3. Adding two remote zeros, we have ten thimbles, ten critical points, and ten zeros for n − = 0, and six of those for n − = 1. We have just two thimbles in n − = 2 subspace because we have no determinent zeros there.
Note that a thimble J σ is not a simple curve but extends with real dimension L, and its section with the subspace is seen as a curve in Fig. 3. For example, integration on the thimble J σ 0 associated to the Gaussian critical point z = σ 0 in n − = 0 subspace contains the perturbative fluctuations in all the directions around z = σ 0 .
Finally in this subsection, let us look at the real part of the action ReS[σ i ] at these critical points for real β = 3, which is listed in the first row (L = 4) of Table 1. The background brightness of Fig. 3 actually indicates the value of ReS[z] (in arbitrary unit). We only list the values at σ 0,1,2,0 because the critical points which interchange with each other by the reflection about the real and imaginary axes have the same ReS[σ i ] atμ = 0 for real β.
The value ReS[σ0] of the doubler solution is larger than ReS[σ 0 ] by 2βL = 24 for n − = 0 and 2β(L − 2) = 12 for n − = 1. This difference comes from the bosonic part β(1 − cos z) of the action. On the other hand, the action ReS[σ 0 ] at σ 0 in n − = 1 sector is larger than that in n − = 0 sector by a factor of order 2β = 6 because the former point contains one doubler component z n = π. One may notice that ReS[σ 1 ] in n − = 1 sector takes a smaller value than ReS[σ 0 ], indicating the larger weight for it. But K σ 1 has no intersection with C at µ = 0, and the thimble J σ 1 is not a member of the integration cycles for Z.
It is intriguing to check this behavior with changing the lattice size L towards the continuum limit. By increasing L and β with keeping β/L and Lma fixed, there appear more zero points and accordingly the critical points aligned in two rows. We compute ReS[σ i ] and list the results for L = 8 and 16 in the lower part of Table 1. We observe that the contributions from the n − = 1 sector to Z are more suppressed by the factor e −2β for the larger L and β. Within the n − = 0 sector, we can estimate the difference between ReS[σ 1 ] and ReS[σ 2 ] as 4π 2 β/L for larger L, basing on the bosonic part Lβ(1 − cos z) and expanding it with approximation σ k ∼ z zero,k ≡ (2k − 1)π/L − im. This gives us a factor 3π 2 ∼ 30, which is consistent with the numerical result of In summary, we have clarified the thimble structure of the Thirring model in this section. The determinant zeros form submanifolds with complex dimension L − 1, and their sections in the subspace assigned with n − appear as isolated zero points. The critical points of the model are classified with n − , and each of them pairs up with a zero point in the subspace assigned with n − (except for the Gaussian critical point σ 0 and its doubler counterpart σ0). Thus all the thimbles are identified in the (0+1) dimensional Thirring model. Towards the continuum limit (β → ∞), ReS[σ] with nonzero n − , which contains n − "doubler" components, acquire the large values of order 2n − β compared to ReS[σ 0 ] in the n − = 0 subspace. This implies that the relative weights of their contributions to Z are strongly suppressed toward the continuum limit, even when they join the set of the integration cycles as µ increases.

Stokes phenomenon and structure change at finite µ
In this section, with increasing µ, we study the change of the intersection numbers and thimbles which contribute to the partition function Z from the viewpoint of the Stokes phenomenon and jumps. We restrict our discussion in the uniform configuration space n − = 0. Figure 5 shows the downward gradient flows of the model with L = 4, β = 3 and ma = 1 forμ = 0.6, 1.2 and 1.8. The zero points z zero 's and their associated critical points σ's move upward as µ increases. Then the critical points which align on the lower side, cross the real axis at certain values ofμ ∼m (see Eq. (3.13)), and accordingly the intersection numbers of K σ 's with C change on the way. Now one encounters the situation where certain thimbles join and/or leave the set of integration cycles for the partition function Z. For a large enoughμ, as can be inferred from Fig. 5 (c), the single thimble J σ 0 comes to connect the two safe remote points z = ±π + i∞, to become an equivalent cycle to the original U-shaped cycle: C ∼ J σ 0 .
The downward and upward cycles J σ and K σ of a critical point σ generally extend to "safe" and "unsafe" remote regions, respectively, without crossing other cycles J σ and K σ which have different values of ImS. When multiple critical points share the same value of ImS, the cycle associated to one of those critical points may meet another critical point. This is the so-called Stokes phenomenon. Change of the intersection number is achieved only by a jump of one endpoint of a upward cycle K σ from (e.g.) z = −i∞ to z = i∞, and in between the critical point σ must undergo the Stokes phenomenon with another critical point σ in which K σ and J σ just overlap. As has been discussed in the previous section, zeros of the fermion determinant become endpoints of the thimbles. Because the determinant appears as −N f log det D in the action S, the imaginary part ImS changes by −2πN f when we encircle a zero point counterclockwise from one side to the other side of a thimble which terminates at this zero. However this difference is not reflected in the gradient flow. Therefore the necessary condition for the Stokes phenomenon to occur between critical points σ and σ is Incidentally, the imaginary part ImS on the upward cycle (e.g.,) K σ may differ by a multiple of 2π depending on which side of the thimble J σ the cycle starts. Moreover, since the value of ImS changes around a zero point, two thimbles can meet at the zero point making an angle determined by the difference of their ImS(σ i ). Thus, one can read the relative phase of the two thimbles from their relative angle when they meet at the zero point.

Stokes jumps with increasing µ
Let us study the Stokes phenomenon with increasing µ in more details. We set N f = 1. Because the configuration subspace for real β is symmetric under reflection about the imaginary axis as seen in Fig. 5, we discuss the thimble structure on the right-half plane hereafter. Even at finite chemical potential µ = 0, this reflection symmetry z → −z guarantees the realness of Z; the thimbles which interchange under this transformation give the contributions which are complex conjugate to each other and whose sum becomes real [58]. In Fig. 6, we compare the values of the action at the critical points σ i . We first note that ImS = 0 at σ 0 and σ0 independently of the chemical potential µ. Indeed, in Fig. 5 (a), we see the Stokes phenomenon between σ 0 and σ0, where the cycles J σ 0 and K σ0 overlap, and At µ = 0 the critical points σī on the upper side have positive values of ImS(σī) and their associated upward cycles K σī extend to the unsafe region toward z = +i∞. With increasing µ the values of ImS(σī) increase monotonically and K σī continue to have no intersection with C. On the other hand, the critical points σ i on the lower side move upward and the imaginary parts ImS(σ i ) at these points increase from negative to positive values with increasing µ, as seen in Fig. 6 (a). In the enlarged plot in the panel (b), the lines of ImS σ i show three crossings atμ =μ * 1 = 0.7,μ * 2 = 0.735 andμ * 3 = 0.86. We now discuss the Stokes phenomenon and the change of the intersection numbers at eachμ * i . In Fig. 7 we show typical thimble structures at several values ofμ. Atμ <μ * 1 , the cycles K σ 1,2 starting from σ 1,2 extend to the lower unsafe region toward z = −i∞, while K σ1 ,2 extend to the upper unsafe region toward z = +i∞. None of them has nonzero intersection with C, and C c ∼ J σ 0 + J σ0 as was discussed previously. Atμ =μ * 1 , ImS σ 0 = ImS σ 2 is achieved, and the two cycles J σ 0 and K σ 2 overlap. Acrossμ * 1 , one end of the upward cycle K σ 2 jumps from −i∞ to +i∞, to give the intersection number n 2 = 1 with C (see panel (b)). (And one end of the cycle J σ 0 jumps from σ0 to z zero,2 .) At the same value ofμ =μ * 1 , the point σ 2 shows the Stokes phenomenon with another critical point σ0 because ImS σ 0 = ImS σ0 = 0. (This coincidence could be avoided by adding a small imaginary part to β, again.) In this case, the two cycles, J σ 2 and K σ0 overlap, and one end of the cycle J σ 2 jumps from π − i∞ to π + i∞ acrossμ =μ * 1 . Hence, we have the equivalence of the cycles 9 Atμ =μ * 2 (panel (c)), the Stokes phenomenon happens between the critical points, σ 0 and σ 1 . The two cycles J σ 0 and K σ 1 overlap there. Whenμ passesμ * 2 , one end of the cycle K σ 1 jumps from −i∞ to +i∞ and one end of the cycle J σ 0 from z zero,2 to z zero,1 , and therefore the critical point σ 1 now acquires the intersection number n σ 1 = 1. Hence, Atμ =μ * 3 (panel (e)), ImS of σ 1 and σ 2 coincide, which allows the Stokes phenomenon between them. Acrossμ =μ * 3 one end of the cycle K σ 2 flips down from +i∞ to −i∞, while changes from 1 to 0. Thus, we have (μ * 4 introduced below) So far we have discussed only the cases where the Stokes phenomenon occurs between the critical points having the same value of ImS. Forμ larger thanμ * 3 we need to take into account the multivaluedness of the logarithm because the edge of J σ 0 is now going around the zero points. The condition for the Stokes phenomenon to occur is the equality of ImS modulo 2π between the two critical points as announced in Eq. (5.1). For our model parameters, there are three more critical valuesμ * 4,5,6 . Atμ =μ * 4 ( Fig. 8 (g)), the condition ImS σ 0 +2π = ImS σ 1 is fulfilled, and forμ * 4 <μ <μ * 5 (Fig. 8 (h)) the equivalent integration cycle becomes Atμ =μ * 5 (Fig. 8 (i)) the condition ImS σ 0 + 2π = ImS σ 2 is fulfilled, and the equivalent integration cycle changes to Atμ =μ * 6 ( Fig. 8 (k)), the condition ImS σ 0 +4π = ImS σ 2 holds and the equivalent integration cycle now consists of a single thimble C c ∼ J σ 0 forμ * 6 <μ . (5.8)

Multi-thimble contributions and weight factor
We have seen how the original integration cycle C is decomposed equivalently into a set of thimbles with increasingμ. The partition function Z is correctly reproduced only if we evaluate the contributions from all the thimbles in the set, in principle. Especially in the crossover region ofμ, multiple thimbles take part in the set of the integration cycles.
However, importance of their contributions depends on the weight factor exp[−ReS(σ)]. For example, the integration cycle consists of J σ 0 and J σ0 for 0 ≤μ <μ * 1 . But the contribution from J σ0 is numerically negligible because ReS(σ0) is larger than ReS(σ 0 ) by a large amount ∼ 2Lβ as seen in Fig. 6 (c) (see also Table 1 forμ = 0 value).
Forμ * 1 <μ <μ * 4 andμ * 5 <μ <μ * 6 , the thimbles J σ ±1 and/or J σ ±2 are in the set of the integration cycles in addition to J σ 0 . According to the weight factor exp(−ReS(σ)) in Fig. 6 (c), the thimble J σ 0 will give the largest contribution and J σ ±1 will contribute as the second largest. The contributions from J σ ±2 will be strongly suppressed. This behavior is mainly controlled by the bosonic part Lβ(1 − cos z) of the action. (The thimble J σ1 is not a member of the integration cycle, although ReS(σ1) becomes smallest asμ increases.) In Fig. 9 we plot the β dependence of the critical chemical potentialμ * i for L = 4 and ma = 1. Outside of the intervalμ * 1 <μ <μ * 6 the single thimble J σ 0 becomes (almost) equivalent to the original integration cycle C, but within this interval multiple thimbles need to be considered. Especially, the second-dominant thimbles J σ ±1 contribute in the interval µ * 2 <μ <μ * 4 . We notice that the crossover regionμ ∼m is indeed covered by this interval µ * 2 <μ <μ * 4 , which indicates that the multi-thimble contribution is requited to reproduce the crossover behavior correctly. The interval becomes wider (narrower) for smaller (larger) β. From this β-dependence there may be a possibility that the approximate evaluation of Z with the single thimble J σ 0 becomes better for larger β. Note that for larger β the difference in the relative weights among the critical points also becomes more significant and the thimbles whose critical point locates away from σ 0 in the real axis direction is expected to less contribute to Z.
In summary, for L = 4 case, we have clarified the change of the Lefschetz thimble structure and the set of the thimbles contributing to Z asμ increases. At small and large chemical potentials outside of the intervalμ * 2 <μ <μ * 4 , the evaluation of Z with the single thimble J σ 0 is legitimate, provided that J σ ±2 contributions are negligibly small. But in the crossover region J σ ±1 contributions must be taken into account in addition to that of J σ 0 . The approximate evaluation by taking only one thimble J σ 0 is performed in numerical simulations for several models so far [49,50,52,53]. Hence it would be worthwhile to examine the validity of the single thimble approximation across the crossover region with varying β. Furthermore it would be intriguing to study how the crossover behavior is reproduced by the multi-thimble contributions with increasing the lattice size L toward the continuum and/or low temperature limits.

Multi-thimble contributions in uniform-field model
In order to examine the single thimble approximation and to investigate how the crossover behavior is reproduced by contributions from multiple thimbles, we study the Thirring model in the uniform-field subspace. The limitation to uniform-field configurations corresponds to the classical approximation with neglecting the quantum fluctuations. The partition function of this restricted model is analytically evaluated to be and the fermion number density and chiral condensate are obtained by Interestingly, in the T = 0 limit this classical model shows a first order transition at the same value of |µ c | = m 2 + g 2 as the original model.

Single-thimble approximation
We compare the values evaluated on the single thimble J σ 0 to the exact ones by taking their ratios in Fig. 10. We show the results with L = 4 and ma = 1 for β = 1 (left) and 3 (right).
On the other hand, the results deviate from unity in the range ofμ * 2 <μ <μ * 4 , indicating that the contributions from J σ ±1 need to be included to reproduce the original integral quantitatively. The much smaller deviation for β = 3 case can be understood if one recalls the rough estimate for the weight factor exp(−ReS(σ ±1 )) ∼ exp(−βπ 2 /(2L)) as discussed in subsec. 4.2. Furthermore, we notice that the missing J σ ±1 contribution to Z changes the sign from positive to negative, and back to positive again, asμ increases. This is the reflection of the fact that ImS(σ 1 ) increases from 0 atμ =μ * 2 to 2π atμ =μ * 4 . Because ImS(σ 0 ) = 0 for anyμ, the two thimbles J σ 1 and J σ 0 contribute additively just aboveμ =μ * 2 . But when ImS(σ 1 ) = π, they contribute with opposite signs. At this point they are connected at z = z zero,1 with an angle π between their edges as seen in Fig. 7 (f). The J σ ±1 contributions return to be positive aŝ µ approaches the critical valueμ * 4 for the Stokes phenomenon. Regarding n and χχ , their integrands have non-constant imaginary parts on J σ ±1 , and the contributions of J σ ±1 to these densities alternate in different ways in the intervalμ * 2 <μ <μ * 4 .

Toward continuum limit
In Fig. 11 we examine the behavior of the fermion number density n J 0 evaluated only on the single thimble J σ 0 as a function of µ/m for L = 4, 8, 16 toward the continuum limit. The parameters are set to (a) (β/L, Lm) = (1/4, 4) and (b) (β/L, Lm) = (3/4, 4). In Fig. 11 (a), some discrepancy from the exact value (dashed line) is seen betweenμ * 2 <μ <μ * 4 for L = 4, where the thimbles J σ ±1 have the nonzero intersection number and need to be included in the integration. This behavior persists when we increase the lattice size to L = 8, 16 toward the continuum limit (thin black dashed curve). The critical values µ * i /m for the Stokes phenomenon with the thimbles J σ ±1 only slightly shift to largerμ toward the continuum limit. In Fig. 11 (b), The discrepancy from the exact values is practically invisible and again the results are relatively insensitive to the size of the lattice with our parameters. This implies that at finite temperatures Monte Carlo simulations on a single thimble may work well for a certain parameters.

Toward low temperature limit
Next we change L as 4, 8, and 16 with fixed β = 1 and ma = 1, toward the zero temperature limit in Fig. 12. We find that the agreement between n J 0 and n 0 is getting worse as L increases. Even in β = 3 case (Fig. 12 (b)) we see a significant discrepancy from the exact result (dashed line) for larger L. As L increases, the slope of the exact curve becomes steeper in the crossover region and eventually converges to a step function, while the single thimble result n J 0 behaves almost as a linear function between two kink points. The singular points indicate the Stokes jump occurring there, through which the thimbles J σ ±1 join or leave the set of the integration cycles for the partition function Z.

Multi-thimble contributions
We draw the thimble structure on the right-half plane for L = 16 atμ = 0.8, 1.0, 1.35, 1.7 in the crossover region with β = 1, ma = 1 in Fig. 13. Atμ = 0.8 the three thimbles J σ0 and J σ±1 have the nonzero intersectin numbers with the original integration cycle, while at larger µ the thimbles J σ0,±1,±2,±3,±4 (according to our numbering) intersect, and they need to be included as the integration cycles to reproduce the partition function Z 0 .
Based on this observation, we extend the evaluation by including the contributions from J σ ±1 for β = 1, 3 and those from J σ ±2 further for β = 1, as shown with dots and crosses in Fig. 12. Indeed, the agreement between the exact and multi-thimble evaluations becomes systematically improved by taking into acount the multi-thimble contributions.
In Table 2 we listed the contributions to the partition function Z 0 and the fermion density n 0 from each thimble with L = 16, β = 1 and ma = 1. The thimbles J σ ±i give the contributions which are complex conjugate to each other so that their sum becomes always real. Regarding partition function Z 0 , the thimble J σ ±0 gives the largest contribution, but the thimbles J σ ±1 also provide a substantial contribution in this crossover region. Those from J σ ±i (i ≥ 2) decrease rather quickly as i = 2, 3, 4 increases, which will be very favorable for a systematic expansion. But we notice that a cancellation occurs between the J σ 0 and J σ ±1 contributions atμ = 1.35 owing to the negative sign of the J σ ±1 contributions. For the fermion density n 0 the cancellation between the J σ 0 and J σ ±1 contributions becomes more delicate atμ = 0.8 and 1.0, while those come to contribute additively atμ = 1.7. Insensitivity of the observables in small chemical region at low temperatures, especially at zero temperature, is sometimes called Silver Blaze phenomenon. We find here that when multiple thimbles contribute to the partition function they show a delicate cancellation between them. The alternating sign exp(−ImS) of the thimbles atμ = 1.35 manifest in Fig. 13 as the fact that the critical points and zero points are aligned and the thimbles are connected at each zero point with the angle about π. In order to check this alternating pattern, we extend our calculation to L = 32 as listed in the bottom row in Table 2. We find that the thimble-bythimble alternating sign and cancellation become more striking not only for Z 0 but also for n 0 . In this case, we need to include the thimbles up to J σ ±3 to evaluate the observables with a few % accuracy. At larger L more zero points appear near the imaginary axis (Eq. (3.13)), and in between the critical points and associated thimbles are aligned atμ in the crossover region. The weight factor from the bosonic part of the action does not suppress these thimble contributions as far as Re(Lβ(1−cos σ i )) < 1. Therefore we need to treat the neat cancellation in multiple thimble contributions in order to reproduce the sharp rise of the fermion density at low temperature (large L). Implication of this observation to the feasibility of the numerical simulations with large lattice size is left for future study.

Summary and discussions
We have studied the Lefschetz thimble structure of the (0+1) dimensional Thirring model at finite chemical potential, which is formulated on the lattice of size L with the staggered fermions and a compact auxiliary vector field. This model suffers from the sign problem by the complex fermion determinant.
The fermion determinant brings in two important features in the complexified field space: many isolated critical points of the gradient flow and submanifolds of the zero points with complex dimension (L − 1). Those critical points accompany the Lefschetz thimbles and the    submanifolds of the zeros serve the ending points for the thimbles. We have identified all the critical points of this model, and furthermore we have pointed out a one-to-one correspondence between a critical point and a zero point within a projected configuration subspace assigned with n − .
We argued that the thimbles associated with the critical points in n − = 0 subspace become more important toward the continuum limit because the relative weights of the other critical points located in n − = 0 subspaces are suppressed by powers of e −2β . The critical points with nonzero n − actually involve the doubler components and they are expected naturally to decouple from the system in the continuum limit.
Hence, restricting our analysis to the critical points in the n − = 0 subspace, we have shown how the thimble structure changes via the Stokes jumps as the chemical potential µ increases. We found that at small and large chemical potentials the single thimble J σ 0 is sufficient as the integration cycle to reproduce the partition function of the model. However in the crossover region we must include multiple thimbles in the set of the integration cycles for the partition function Z. Their relative weights depend on the lattice size L and the coupling strength β.
Taking the uniform-field model as a concrete example, we have examined the importance of the multi-thimble contributions and how the crossover behavior is generated by them. The single-thimble approximation is justified for large β/L ∼ T /g 2 , even in the continuum limit. But as we increase the lattice size L, i.e., lower the temperature T with β and ma fixed, we have seen the breakdown of the single-thimble approximation, which indicates the necessity of the multi-thimble contributions. The sign of those contributions is alternating, which yields a neat cancellation to reproduce the correct values of Z and observables at large L. We notice that the contributions from the thimbles away from the origin diminish rather quickly. The Silver Blaze behavior and the following abrupt rise of the density n with increasing µ are achieved by the interplay among the multi-thimble contributions in the crossover region.
We have performed HMC simulations for the (0+1) dimensional Thirring model with finite chemical potential on the single thimble J σ 0 in Ref. [66]. We observed scaling behavior of the results to the continuum limit at finite temperature and to the low-temperature limit. The single thimble evaluation in the crossover region is getting worse for smaller β and/or larger L, which is consistent with the results obtained in the uniform-field model. We show one example of the simulation results for L = 16, β = 3 and ma = 1 in Fig. 14. For comparison, we also tried the complex Langevin simulation as yet another approach with complexification and as a possible way to include the "multi-thimble" contributions, which is shown in Fig. 14. We find that the Langevin result also deviates from the exact one in the crossover region, but in a different manner. We observed that the sampling points in the Langevin simulation are distributed around the thimbles J σ 0 and J σ ±1 . The details of the Langevin simulation will be reported elsewhere.
We have seen that an interplay among multi-thimble contributions are necessary and important to describe the rapid crossover behavior of the fermion system. However it is a difficult task to identify all the critical points in generic models. Our analysis suggests that the thimbles whose critical points locate in the uniform-field subspace will give dominant contributions, while those with critical points in non-uniform-field subspace will decouple by the suppressed weight factor toward the continuum limit because they have doubler components. Assuming that we can identify all the relevant thimbles to be integrated over, we will face another challenge -how to add up the multi-thimble contributions in the Monte Carlo simulation. In our model analysis we can sum up them by knowing the partition function values Z σ precisely, but in Monte Carlo simulations we compute only the average of the observables not the partition function. It is, therefore, extremely important to devise the efficient way to perform the multi-thimble integration by extending the Monte Carlo algorithm for practical applications of the Lefschetz thimble integration to fermionic systems with the sign problem. A Exact expression and asymptotics of Z In this appendix, we give the exact expression for the partition function of the Thirring model with the compact action. We assume N f = 1 and L is even.
A useful formula for a matrix determinant is known in [67]: where m ± = m ± √ m 2 + 1. Withm ≡ sinh −1 m and with even L, this can be written as in Eq. (2.4).
Because the A n -odd terms in the determinant vanish after the integration over A n with weight e −β(1−cos An) , we can write the partition function as where I 0 (x) and I 1 (x), respectively, are the zeroth and first order modified Bessel functions of the first kind. The fermion number density and the scalar density can be derived by differentiating ln Z with respect to µ and m, respectively. Using the asymptotic expression of the modified Bessel function I 0,1 (β) for a large β, we find that in the continuum limit at finite T , the partition function Eq. (2.5) scales as It is interesting to observe that in the T → 0 limit both models show a first order transition at the same point |µ c | = m + g 2 .
If we take L large with β fixed, we find In the infinite-L limit these models show a first-order transition at |μ c | =m + ln(I 0 (β)/I 1 (β)) and |μ c | =m + β − η, respectively.