Lefschetz-thimble analysis of the sign problem in one-site fermion model

The Lefschetz-thimble approach to path integrals is applied to a one-site model of electrons, i.e., the one-site Hubbard model. Since the one-site Hubbard model shows a non-analytic behavior at the zero temperature and its path integral expression has the sign problem, this toy model is a good testing ground for an idea or a technique to attack the sign problem. Semiclassical analysis using complex saddle points unveils the significance of interference among multiple Lefschetz thimbles to reproduce the non-analytic behavior by using the path integral. If the number of Lefschetz thimbles is insufficient, we found not only large discrepancies from the exact result, but also thermodynamic instabilities. Analyzing such singular behaviors semiclassically, we propose a criterion to identify the necessary number of Lefschetz thimbles. We argue that this interference of multiple saddle points is a key issue to understand the sign problem of the finite-density quantum chromodynamics.


Introduction
Understanding strongly-correlated quantum many-body systems has been an ultimate goal in contemporary physics. Numerical simulation formulated on the discretized spacetime, especially lattice quantum Monte Carlo method, is a powerful ab initio tool for that purpose. Exact diagonalization of a Hamiltonian provides us complete information on the system; however, it usually requires the huge computational cost and is limited to small systems. Monte Carlo method, on the other hand, is based on the importance sampling in the phase space of the system, whose computational cost scales algebraically with the system size. Many thermodynamic quantities can be computed for various systems using this method, such as finite-temperature quantum chromodynamics (QCD) in hadron physics [1], and liquid helium [2], ultracold atomic gases [3], Bose-Fermi mixtures [4] in condensed matter physics.
In many quantum systems of great interest, however, Monte Carlo simulation suffers from the so-called sign problem [5,6]. The path-integral weight e −S becomes negative or even complex, and thus the importance sampling cannot be applied. Since the cancellation between large positive and negative noises gives a physically meaningful signal, the computational time to reduce statistical errors grows exponentially with the system size [7]. In condensed matter physics, famous examples of the sign problem are strongly-correlated electrons, such as the Hubbard model away from half-filling [5], frustrated spin systems, such as the XY model on the Kagomé lattice, and so on. In hadron physics, QCD at finite quark densities attracts much attention for understanding the interior structure of neutron stars [8][9][10], but we have no ab initio simulation due to the sign problem [11].
The sign problem of the finite-density QCD is known to become too severe when the quark chemical potential exceeds half of the pion mass [12][13][14][15][16][17]. At low temperatures, there is no phase transition until the quark chemical potential reaches about one third of the nucleon mass; however, the Monte Carlo simulation of two flavor QCD shows a phase transition at half of the pion mass [12,13]. No one can describe the true onset at one third of the nucleon mass by using the path-integral formulation of lattice QCD [14,15]. This is one of the biggest issue in finite-density QCD. In this paper, we address this problem, which is called the "Silver Blaze problem" [14], by studying a toy model.
In this paper, we apply the Lefschetz-thimble approach to the one-site model of electrons. This toy model can be regarded as an extreme limit of strong couplings because it can be obtained by neglecting hopping terms in the Hubbard model. The Hamiltonian of the one-site model can be easily diagonalized and thus we can calculate any expectation value exactly. However, since this model has the severe sign problem in its path-integral expression, it is hard to calculate expectation values by the conventional Monte Carlo method. This toy model provides us a good playground to study theoretical structures of the Lefschetz-thimble approach. In a previous study [34], two-dimensional Hubbard model is studied using the Lefschetz-thimble Monte Carlo method, but one must use the so-called "one-thimble approximation" due to the current limitation of the numerical algorithm. Since our system is exactly solvable, we can understand the complete structure of Lefschetz thimbles, and suggest appropriate approximation schemes.
Based on the semiclassical analysis, we find that interference of complex phases among multiple saddle points is important to reproduce the step-function behavior of the density. The fermion spectrum at a complex saddle point resembles the quark spectrum of phase quenched finite-density QCD, which would suggest the interference of multiple saddle points might also occur at finite-density QCD. We will discuss this point in detail later. This study will help us to understand the sign problem in QCD.
The outline of this paper is as follows. In Sec. 2, the one-site Hubbard model is explained, and its physical properties are calculated by exactly diagonalizing the Hamiltonian. The path integral expression of the toy model is introduced in detail, and the sign problem appearing in its path integral formulation is briefly reviewed. In Sec. 3.1, the Lefschetz-thimble method is introduced, and is applied to the one-site Hubbard model. Complex saddle points play a pivotal role in this approach, and thus we study their property in Sec. 3.2. We discuss structures of Lee-Yang zeros and the fermion spectrum at complex saddle points in Sec. 3.3. In Sec. 3.4, numerical results on the fermion number density are compared with the exact result using several approximations based on the Lefschetz-thimble approach. In Sec. 4, we discuss the resemblance of this model to the finite-density QCD and address the Silver Blaze problem in QCD. We summarize our result in Sec. 5.

Physical properties of the model
The Hubbard model [46][47][48] describes the lattice fermion system, whose dynamics is governed by the (second quantized) Hamiltonian, The summation is taken over all the lattice sites i, and the notation i, j denotes that only the nearest neighbor hopping is considered.ĉ † σ,i andĉ σ,i , obeying ĉ σ,i ,ĉ † τ,j = δ στ δ ij , are creation and annihilation operators of fermions of the spin σ(=↑, ↓) at the site i.n σ,i =ĉ † σ,iĉ σ,i is the number operator. The parameter t is called the hopping parameter, and we will consider only the case t = 0 in this section. The parameter U (> 0) describes the on-site repulsive interaction: The energy increases if two particles get together on the same site. The Hamiltonian is invariant under the global U (1) rotationĉ σ,i , and thus the total number of fermions is conserved. In the Hamiltonian (1), the chemical potential µ is introduced for this conserved quantity.
In the strong coupling limit t = 0, one can solve the Hamiltonian (1) analytically. Let us see the result on the total number density n := (n ↑ +n ↓ ) as a function of the chemical potential µ. Since each site is totally independent from others, one can study the one-site Hamiltonian, In the following, we only consider the one-site case. Since the Hamiltonian (2) commutes with the number density operatorn, we can take the number basis to find the ground state. Let us define the grand partition function by with the inverse temperature β = 1/T . The number basis gives an explicit result: The number density is given as In the zero-temperature limit β → ∞, the number density takes 0, 1, and 2 for µ/U < 0, 0 < µ/U < 1, and µ/U > 1, respectively, and shows the step-function like behavior. In particular in the case µ/U = 1/2, we have n = 1 for any β, which comes from the symmetryn → 2 −n (or,n σ → 1 −n σ ) in the Hamiltonian (2). This point is called half-filling, at which the particle-hole symmetry exists.

Path integral formulation and the sign problem
We have explicitly seen that the one-site Hubbard model can be analytically solved by using the number eigenstates. However, if one changes the basis of the Hilbert space for taking trace, the sign problem emerges and thus the Hubbard model in the strong coupling provides us a good lesson. Let us first derive the path integral expression of the partition function (3). Using two-component complex Grassmannian variables ψ = (ψ ↑ , ψ ↓ ), the fermion coherent state is defined by which satisfiesĉ σ |ψ = ψ σ |ψ . We take the convention in which all the Grassmannian variables anticommute with each other and with fermionic creation/annihilation operators. It is easy to check that and tr(O) = dψ * dψ e −ψ * ψ −ψ|O|ψ .
Let |n ↑ , n ↓ be a fermionic Fock state, and then the inner products with coherent states are Using these relations, one can compute the matrix element of the imaginary-time evolution operator as Here we take an approximation for small ∆τ . In order to circumvent the quartic interaction, let us introduce an auxiliary real field ϕ via the Hubbard-Stratonovich transformation, that is, Inserting Eq. (11) into Eq. (10), we obtain In order to take the naive continuum limit of the path integral expression, we need to exponentiate the auxiliary field in the fermionic bilinear term as Here, ϕ 2 0 means the expectation value of ϕ 2 with the Gaussian weight in Eq. (12), which gives U/∆τ . The path integral expression of the partition function reads with the antiperiodic boundary condition ψ N +1 = −ψ 1 . Equation of motion in terms of ϕ shows The expectation value of ϕ is nothing but that of the total number density. This formula will be used later in order to compute the number density. Since Eq. (14) is quadratic in the fermionic field ψ and ψ * , we can perform the Grassmannian integration explicitly. The spectrum of the fermion bilinear operator in the exponential in Eq. (14) is given by and then, in the continuum limit, the fermion path integral (determinant) becomes Since fermions have two spin degrees of freedom, the overall square is taken. The fermion determinant (17) contains only the Matsubara zero mode of ϕ. Thus the path integral of non-zero Matsubara modes of ϕ gives a trivial Gaussian integration and does not depend on µ. The path integral of our interest is now reduced to an integral of zero Matsubara mode ϕ bg = β 0 dτ ϕ(τ )/β, and we have This integral can be performed analytically to find Eq. (4). Instead, we shall apply the Lefschetz-thimble method because Eq. (18) has the sign problem coming from the complex weight e iβϕ bg . This simple model enables us to study the structure of the sign problem in terms of the Lefschetz-thimble approach.
In order to explicitly show that the path integral (18) contains the sign problem, let us compare it with the phase quenched partition function This integration can be performed analytically as The positive sign in front of U is different from that in Eq. (4). If the last term is negligible, the sign problem is mild. This is the case for small chemical potentials µ/U −1/2 and low temperatures βU 1; the ratio Z/Z is almost one because the dominant contribution comes from the unoccupied state. Therefore the sign problem is sufficiently weak. In the later section, we can interpret this semiclassically because the important configurations accumulate around the origin ϕ bg 0. The sign problem becomes quite severe for µ/U −1/2 and low temperatures βU 1. For example, at µ/U = 0, the ratio becomes This severe sign problem comes from the last term in Eq. (20), which causes the wrong onset of the density at µ/U = −1/2. A similar problem happens in QCD at finite density when the quark chemical potential exceeds the half of the pion mass [12,13,49]. For large chemical potentials, µ/U 3/2, we can switch the occupied and unoccupied fermions by changing the integration variable from ϕ bg to −ϕ bg + 2iU , so that the sign problem is weakened. There is another remarkable condition in which the sign problem disappears again by a simple shift of the integration variable: At the halffilling µ/U = 1/2, a new integration variable ϕ bg = ϕ bg + iU rewrites the original integral (18) into The integrand of Eq. (22) is nonnegative. This is because the particle-hole symmetry at the half-filling becomes manifest in this new integration variable. As a conclusion of this section, the sign problem can be circumvented in these three special cases, µ/U −1/2, µ/U 3/2 and µ = 1/2, even using the path integral expression; however, this does not happen in general cases.

Lefschetz-thimble method
Let us start with a multiple integration that gives the partition function where S(x) is a complex action functional of the real field x = (x 1 , . . . , x n ) ∈ R n . In order to circumvent the oscillatory integral (23), we complexify the integration variables Integration is performed on steepest descent paths, called Lefschetz thimbles, instead of directly evaluating Eq. (23). Each Lefschetz thimble is an n-dimensional space spanned around a saddle point z σ in C n (σ ∈ Σ). Let us consider Morse's flow equation for complexified variables z by introducing a fictitious time t [22][23][24]: Along this flow, the integrand becomes smaller since and thus exp −S(z(t)) → 0 as t → ∞. This flow equation has an important conserved quantity: The imaginary part of S obeys d dt Im S(z(t)) = 0.
Therefore, along each flow line, the integrand has no oscillatory phase, which is a remarkable property to circumvent the oscillatory integral. Now, let us define a new integration cycle called Lefschetz thimble as follows: The Lefschetz thimble J σ is identified as the set of points reached by some flows emanating from z σ : Complex contour integral over J σ definitely converges since Re S diverges at boundaries of J σ . On each integral, we have no sign problem because Im S is constant. The partition function can now be computed by the sum of the nicely converging integrals as ‡ There is another approach to the sign problem also based on the idea of complexification, the complex Langevin method [20,21]. See Refs. [50][51][52][53][54] for recent developments.
The coefficient n σ is given by the intersection number between R n and K σ . The dual thimble K σ is defined as the set of the points reached by flows getting sucked into z σ : There are two possible origins where the sign problem reappears in the Lefschetzthimble method. One is a complex phase coming from the Jacobian of the integration measure d n z on Lefschetz thimbles, which is often called the residual sign problem [37][38][39][40]. Geometrically, the residual sign problem is mild if the Lefschetz thimble is almost straight. Another one comes from the possible cancellation in summing up Lefschetz thimbles in Eq. (28), because the complex phase Im S can be different among several Lefschetz thimbles.

Semiclassical analysis
Let us apply the Picard-Lefschetz theory to the path integral (18). The effective action of this system is The effective action satisfies the reality condition and then the Lefschetz-thimble decomposition manifestly respects this symmetry so as to ensure the reality of physical observables [45,[55][56][57].
The logarithmic function has branch singularities at fermionic Matsubara modes z = i(µ + U/2) + (2 + 1)π/β for ∈ Z. The integrand converges to 0, i.e., S → +∞, at these points. The flow equation of Eq. (30) reads Let us find the set of saddle points, which is an important step not only for the saddle-point approximation but also for the Lefschetz-thimble method. The saddlepoint condition of the effective action (30) is In order to simplify our analysis, we consider a limiting case where T U, |µ|. We can approximately obtain the saddle points as for m ∈ Z by assuming that the second term is much smaller than the first one. The reality condition (31) says that z −m and z m form a pair. If µ > 3U/2 or µ < −U/2, there is another solution Behaviors of the downward flow equation (32) are shown in Fig. 1 with U = 1 and βU = 30. Let us discuss the case when the chemical potential is sufficiently large (µ/U = 2). In this case, Morse's flow is shown in Fig. 1 (a). Only one Lefschetz thimble J * associated with z * contributes in the Lefschetz-thimble decomposition. Other dual thimbles shown with dashed green lines do not intersect with the original integration cycle R as in Fig. 1 (a). Thus intersection numbers in Eq. (28) vanish and we have Within the saddle-point approximation using the complex saddle point z * , the number density is given as and thus the saturation of fermions is well described. The exactly similar thing happens also for µ/U −1/2, and then the saddle-point approximation gives n −iz * /U = 0. This analysis explicitly shows why the sign problem is significantly weakened with an appropriate shift of integration variables for µ/U −1/2 and µ/U 3/2, as we have discussed in Sec. 2.2.
In the following, we concentrate on the case where the sign problem is severe, −1/2 µ/U 3/2, at low temperatures βU 10. In this parameter region, the typical behavior of the flow is shown in Fig. 1 (b). All dual thimbles K m intersect with the original integration cycle R. This shows that all the saddle points z m contribute to the partition function, so that the interference among them may not be negligible. This interference requires a careful treatment of the semiclassical analysis in order to solve the sign problem. The big difference between Figs. 1 (a) and (b) can be explained by the Stokes phenomena [23], which occur around µ/U 3/2 and also around µ/U −1/2. However, this is irrelevant to the non-analytic behavior of the number density, and then we stick to the analysis for −1/2 µ/U 3/2.
Let us denote the classical action at the saddle point z σ as S σ (:= S(z σ )). Substituting the approximate expression (34) into Eq. (30), we have For Eqs. (39) and (40), we have written down the m-dependent leading terms in terms of large βU . Equation (39) shows that subdominant thimbles J m can be comparable with the dominant one J 0 for βU 1 so long as |m|( = 0) is not so large. According to Eq. (40), these different thimbles have different complex phases, and thus a contribution from one Lefschetz thimble is generally canceled by other ones. Exactly at the halffilling, µ/U = 1/2, the complex phase is always an integer multiple of 2πi, and thus such a cancellation is absent. This gives an interpretation on the reason why the sign problem disappears at the half-filling from the viewpoint of Lefschetz-thimble approach.
Let us compute the partition function Z cl only using these semiclassical information: In this formula, Gaussian fluctuation around z m is neglected because it only gives an unimportant overall factor at low temperatures. Indeed, within our approximation, ∂ 2 z S(z m ) is independent of m, and the result of the Gaussian integration slip by the summation. This summation (41) can be computed analytically by using the elliptic theta function as Here we used an approximate expression of the saddle points, Eqs. (39) and (40), to obtain Eq. (42). Using this result, in the limit β → ∞, the number density is given as This good agreement might be accidental to some extent, since we have neglected higher order contributions in the semiclassical analysis. This result still indicates the usefulness of the semiclassical approximation even when the sign problem is severe.

Lee-Yang zeros and the fermion spectrum
Let us analyze the non-analytic behavior (43) from the viewpoint of Lee-Yang zeros [58][59][60][61][62][63][64][65][66] in the semiclassical expression (42). Using the infinite-product expression of the elliptic theta function, We find that the zeros of the partition function in the complex µ plane satisfy for some positive integer . Since Eq. (42) is originally derived as an approximate expression for −1/2 < µ/U < 3/2, we restrict the real part of the complex chemical potential to −U/2 < Re µ < 3U/2. Then, the zeros of the semiclassical partition function are In the limit T → 0, these zeros form straight lines crossing to the real axis at µ = 0, U . This behavior of zeros is a signal of the first-order phase transition in the context of the Lee-Yang zero analysis of a therodynamic system, and indeed we find the non-analytic jumps (43) of n at µ = 0, U in the limit of T = 0 in the one-site Hubbard model. Let us apply the same analysis to the phase quenched partition function. We define the phase quenched semiclassical partition function by Note that this expression does not correspond to the semiclassical limit of Eq. (19) because the saddle points at −1/2 < µ/U < 3/2 do not lie on the original integration cycle R. In the phase quenched approximation, the interference of complex phases e −i ImSm disappears, and then the µ dependence is given only by e −S 0 . Lee-Yang zeros do not exist if −U/2 Re µ 3U/2. Therefore, the only possibility to explain the non-analytic behavior in this phase quenched approximation originates from the Stokes jumps around µ/U −1/2, 3/2, which naturally explains the non-analytic but continuous change of n (see Ref. [44]).
In order to investigate this failure of the phase quenched approximation more deeply, we analyze the fermion spectrum (16) at complex saddle points. For simplicity, let us take the continuum limit at first with an appropriate normalization, The spectral representation of the number density at a background field ϕ is given by If the sign problem is mild enough, µ < −U/2, then the eigenvalue at the saddle point z * = 0 is Since Re λ con (z * , µ) cannot be zero, the chemical potential dependence of the number density is exponentially suppressed as, in the limit T → 0, However, for µ > −U/2, Re λ con can be zero in the vicinity of saddle points, and indeed one finds Within this approximation, the number density becomes The near-zero fermionic mode provokes the fake onset of the number density at µ = −U/2 again. From the viewpoint of the fermion spectrum at complex saddle points, nothing special happens in the vicinity of the correct transition points µ = 0, U . In the language of Lefschetz thimbles, this implies that their topological structure does not change at all around µ = 0, U , as we can see in Fig. 1 (b).

Numerical results
Let us consider one-, three-, and five-thimble approximations in order to analyze the importance of the interference among multiple Lefschetz thimbles. We consider the case where the sign problem is severe, i.e., −1/2 µ/U 3/2 and βU 10. The (2m + 1)-thimble approximation takes into account only J 0 , J ±1 , . . . , J ±m in Fig. 1 (b). The partition function (18) in this approximation reads Each contribution of the summation in the second term comes from a pair J k and J −k , and becomes manifestly real because of the reflection symmetry z → −z. Since the Lefschetz thimble J k is homologically equivalent to the interval i(µ + U/2) + [(2k − 1)πT, (2k + 1)πT ], the above integration (54) is reduced to with x = z − i(µ + U/2). The result on the number density at βU = 30 is shown in Fig. 2, in which one-, three-, and five-thimble approximations are compared with the exact computation. The one-thimble approximation, which is shown with the solid red line in Fig. 2, is an approximation to integrate only over the Lefschetz thimble J 0 . It almost gives the mean-field result, that is, for −1/2 < µ/U < 3/2. Figure 2 shows that the one-thimble approximation is not sufficient to describe the plateaus of the number density in the region −1/2 < µ/U < 3/2. One can easily check that Eq. (56) is also obtained from the phase quenched approximation (53) after shifting the integration variable in Eq. (18) as ϕ bg → ϕ bg + i(µ + U/2). The result of the three-thimble approximation provides us a useful lesson for an application of the Lefschetz-thimble approach to the sign problem. The number density diverges at several chemical potentials around a rapid crossover, although the result is improved around each plateau. Let us analyze this divergence by using the semiclassical analysis. For that purpose, we introduce the semiclassical partition function with the (2m + 1)-thimble approximation by In the three-thimble approximation, the semiclassical partition function approximately behaves as Since βU 1 and thus e −2π 2 /βU ∼ 1, the second term is not suppressed, compared with the first term. This approximate partition function vanishes at some chemical potentials around µ/U = 0 and 1 (Z (3) cl vanishes at µ/U ±0.04, which is roughly consistent with Fig. 2). This example clearly demonstrates that picking up a part of the Lefschetzthimble decomposition can violate the physical requirement, such as the positivity condition on the thermodynamic quantities. In this instance, the incompressibility (∂n/∂µ)/n 2 should not be negative.
According to Fig. 2, we can conclude that five Lefschetz thimbles J 0 ∪J ±1 ∪J ±2 are necessary in order to explain rapid jumps of the number density in terms of the chemical potential at a low temperature βU = 30. In fact, the five-thimble approximation (dotted blue line) well reproduces the exact computation (solid black line) as seen in Fig. 2. How many Lefschetz thimbles are necessary in general at a given lower temperature β? We propose the criterion to neglect the Lefschetz thimbles J m for large |m|: Assuming that the number of Lefschetz thimbles, (2m + 1), is sufficiently large, one can replace Z (2m+1) cl in the denominator by the semiclassical partition function Z cl . Since |Z Here ε is a controlling parameter of an error, and ε 1. Using the approximate results (39) and (42), we can solve this inequality with respect to m: This criterion gives different results depending on µ, and thus let us first derive the strongest restriction. If the sign problem is severe, i.e., µ = 0, the elliptic theta function is exponentially small with respect to βU : Therefore, for reasonable ε such as ε 0.1, the effect of ε is negligible and the criterion (61) gives |m| βU 4π (63) in the limit of βU 1. It means that we need at least (2 βU/4π + 1) thimbles in order to describe the rapid crossover of the number density in the case of the one-site Hubbard model. Let us discuss this behavior from another point of view. Since β becomes large, the relevant integration region becomes smaller as |Re z| U/β. However, Fig. 1 shows that the length of each Lefschetz thimble J k is proportional to 1/β, and thus we need more Lefschetz thimbles in order to cover the relevant integration region. That number is clearly proportional to √ βU , and is much smaller than Eq. (63). Because of the interference of complex phases among Lefschetz thimbles, the semiclassical partition function becomes much smaller than that of the naive expectation if µ/U ∼ 0, 1. This indicates that it is much easier to reproduce the plateau in the vicinity of the half-filling than that of the rapid crossover. The necessary number of Lefschetz thimbles for each phenomenon is also largely different. In fact, for µ/U = 1/2, one can show that and, roughly speaking, the criterion (61) gives for J m being negligible with keeping a good accuracy only around the half filling. This is consistent with the result of the above heuristic argument. Now, we can understand that the large gap of these two estimates (63) and (65) comes from the sign problem in summing up Lefschetz thimbles, and the semiclassical analysis provides a reasonable rough indication.

Speculation on the Silver Blaze problem in finite-density QCD
We studied the one-site Hubbard model as a toy model of the sign problem, and the Silver Blaze problem. We elucidate that the interference of complex phases among multiple saddle points is important to analyze these problems. In this section, we discuss the analogy of the Silver Blaze problem in the one-site Hubbard model with that in finite-density QCD. Let us pay attention to the behavior of the baryon number density n B for finitedensity QCD at T = 0. In general, the quark determinant, Det γ 4 ( / D A + m) − µ qk , will have the µ qk -dependence, which naively means that the baryon number density arises for any chemical potentials. However, since we know that the lightest baryon in the QCD spectrum is nucleon, which is a composite particle of three quarks, the baryon number density n B must be zero for any 3µ qk < m N − B 923 MeV (m N is the nucleon mass and B is the nuclear binding energy). The µ qk -dependence of the QCD partition function must vanish for µ qk m N /3 at T = 0, which empirical fact is highly nontrivial from the viewpoint of the path integral of lattice QCD as discussed above (Fig. 3 illustrates this situation schematically). Understanding this behavior of n B must be a key first step to correctly perform the lattice QCD simulation at arbitrary temperatures and baryon chemical potentials. In Ref. [14], this is named the baryon Silver Blaze problem after the famous story of Sherlock Holmes. So far, the baryon Silver Blaze problem has been well understood only for µ qk < m π /2 (m π 140 MeV is the pion mass) [14,15,65]. In general, the eigenvalues of γ 4 ( / D A + m) − µ qk can be expressed as λ (j,n) (A, µ qk ) = ε j (A) − µ qk − iφ j (A) + i(2n + 1)πT with −πT < φ j ≤ πT [14]. Therefore, Det γ 4 ( / D A + m) − µ qk = j n λ (j,n) . The product over n can be taken in the exactly same way as we usually do for that over Matsubara frequencies, and we find that it is proportional to e β(ε j −iφ j −µ qk )/2 + e −β(ε j −iφ j −µ qk )/2 for each label j. The fermion determinant becomes where N is the normalization factor (see Ref. [15] for its more rigorous derivation using the zeta-function regularization). Its µ qk -dependence at sufficiently low temperatures and at µ qk > 0 becomes crystal clear according to the following expression [14,15]: The quark determinant becomes independent of µ qk at T = 0 if µ qk is smaller than the minimum of ε j (A)(> 0). For statistically significant gauge fields, it is confirmed that min(ε j (A)) = m π /2 [14,15,65], and thus n B = 0 for µ qk < m π /2. This also means that the Lefschetz-thimble decomposition of lattice QCD at µ qk < m π /2 and at sufficiently low temperatures is identical to the original integration cycle up to an exponentially small correction because the quark determinant at µ qk = 0 is positive. In contrast, as computed in Eq. (17), the fermion determinant in the one-site Hubbard model is We can ask the same question: How can we show n = 0 for µ < 0 in the limit β = ∞ by using the path-integral expression when the fermion determinant depends on µ? This Silver Blaze problem for µ < −U/2 can be solved in the same manner as the baryon Silver Blaze problem at µ qk < m π /2. At µ < −U/2, the fermion determinant becomes independent of µ as T → 0. This also explains why the Lefschetz thimble J * at µ < −U/2 is almost identical to the original integration cycle R (The similar discussion shows the Lefschetz thimble J * at µ > 3U/2 is almost identical to the line with Im(z) = 2). For −U/2 < µ < 0, the µ-dependence of the fermion determinant (68) becomes exponentially large, and thus the fact that n = 0 at T = 0 is still veiled. At µ −U/2, zeros of the fermion determinant (star-shape black points in Fig. 1) come close to the real axis, which is almost identical to the Lefschetz thimble J * for µ < −U/2. In the case of the one-site Hubbard model, this triggers Stokes jumps and the original integration cycle is decomposed into multiple Lefschetz thimbles for µ > −U/2. One can no longer replace the expectation value of the number density by the number density at a complex saddle point. The significance of interference among multiple Lefschetz thimbles is identified in order to explain not only the rapid jumps of n but also the µindependence of n (n = 0) for −U/2 < µ < 0 at the zero temperature. Furthermore, if the number of Lefschetz thimbles is insufficient, the thermodynamic stability is violated in the one-site Hubbard model.
We speculate that this decomposition of Lefschetz thimbles and accompanying interference of complex phases also play an important role in finite-density QCD, in particular, in the Baryon Silver Blaze problem and in the sign problem beyond half of the pion mass. For µ qk > m π /2, the quark determinant becomes dependent on µ qk and highly oscillatory on the statistically significant domain of real gauge fields. This means that zeros of the quark determinant come close to that domain. It would be natural to conclude that the statistically dominant region of the original integration cycle starts to be decomposed into multiple Lefschetz thimbles around µ qk m π /2. The complex phase of the fermion determinant becomes ill-defined at its zeros. Thus, the constantphase condition of the integrand on a Lefschetz thimble would be hard to be satisfied without such decomposition, when zeros of the fermion determinant are located in the vicinity of the Lefschetz thimble. Once the original path integral is decomposed into multiple thimbles, we might need to carefully take their interference into account to explain the µ qk -independence of n B until µ qk exceeds (m N − B)/3.
Because of the complexity of QCD, we cannot show these statements rigorously so far. Future study of QCD-like models based on Lefschetz thimbles and justifying our speculation will be crucial to develop our understanding on the baryon Silver Blaze problem and on the sign problem of the finite-density QCD beyond half of the pion mass.

Summary
One-site repulsive Hubbard model shows a non-analytic behavior by changing the chemical potential at the zero temperature β = ∞. This model is obtained by taking the strong coupling limit, and thus our analysis of Lefschetz thimbles forms basics of studying the strongly-coupled repulsive Hubbard model using the path integral. Since the number eigenstate can solve this model exactly and its path integral expression has the sign problem, this is a good testing ground to analyze the structure of the sign problem in various solutions. Furthermore, the non-analytic behavior of the number density looks quite similar to the Silver Blaze problem in finite-density QCD, and thus the understanding of this model from the functional integral itself is an important task.
In this paper, the sign problem in the Lefschetz-thimble approach is carefully studied. Topological structures of the Lefschetz thimbles are analyzed based on the semiclassical analysis, and their relation to the fermion spectrum is investigated. If one picks up only the most relevant Lefschetz thimble, the mean-field approximation is almost recovered. The significance of interference among multiple Lefschetz thimbles is identified in order to explain the non-analytic behavior. As we lower the temperature, there exist regions where the necessary number of Lefschetz thimbles increases linearly in βU . On the other hand, in order to explain the Silver Blaze like behavior, that number only grows as √ βU . If the number of Lefschetz thimbles is insufficient, we found not only large discrepancies from the exact result, but also thermodynamic instabilities due to artifact of the approximation. This means that the Lefschetz-thimble decomposition does not manifestly ensure the thermodynamic stability, or the positivity condition on the partition function, although its reality is always true with a reasonable condition [45]. Developing our understanding on the stability from the viewpoint of Lefschetz thimbles is important.
Since the one-site Hubbard model is quite a simple model, it is not clear whether interference of multiple Lefschetz thimbles is relevant for studying thermodynamic systems with complicated interactions. Nevertheless, we discuss the speculation on the baryon Silver Blaze problem of finite-density QCD based on its mathematical similarity to the one-site Hubbard model. We conjecture that the statistically dominant gaugefield configurations are decomposed into multiple Lefschetz thimbles when the quark chemical potential exceeds half of the pion mass. To justify our conjecture on finitedensity QCD and QCD-like models is significantly important to solve the baryon Silver Blaze problem and to numerically simulate the finite-density QCD beyond the half of the pion mass.
In order to understand strongly-correlated electron systems with the sign problem, applying the Lefschetz-thimble method to the two-site Hubbard model is a nice straightforward extension of this study. This enables us to take into account the effect of the hopping term. If the atomic potential is sufficiently strong, localization of fermions is a natural scenario and thus two-site models will become benchmark. We believe that semiclassical analysis with complex saddle points provides us a deeper understanding of the sign problem, and developing that technique in general cases is an important future study.