Monte Carlo study of Lefschetz thimble structure in one-dimensional Thirring model at finite density

We consider the one-dimensional massive Thirring model formulated on the lattice with staggered fermions and an auxiliary compact vector (link) field, which is exactly solvable and shows a phase transition with increasing the chemical potential of fermion number: the crossover at a finite temperature and the first order transition at zero temperature. We complexify its path-integration on Lefschetz thimbles and examine its phase transition by hybrid Monte Carlo simulations on the single dominant thimble. We observe a discrepancy between the numerical and exact results in the crossover region for small inverse coupling β and/or large lattice size L, while they are in good agreement in the lower and higher density regions. We also observe that the discrepancy persists in the continuum limit to keep the temperature finite and it becomes more significant toward the low-temperature limit. This numerical result is consistent with our analytical study of the model and implies that the contributions of subdominant thimbles should be summed up in order to reproduce the first order transition in the low-temperature limit.


Introduction
The physics of QCD at finite temperature and density is one of the most important subjects in high energy physics and also in cosmology and astrophysics. To investigate QCD, especially its static and thermodynamic properties, the Monte Carlo simulation of lattice QCD has proved to be a powerful method. However, in the extreme condition of low temperature and high density, the sign problem in lattice QCD, caused by introducing the baryon-number chemical potential, prevents us from the thorough study of the properties of QCD [1]. Recently two alternative approaches to the problem have attracted much attention -complex Langevin dynamics [2][3][4] and Lefschetz thimble method [5][6][7]. Both methods are based on the complexification of dynamical field variables. 1 In our previous work [64], we have applied the Lefschetz thimble method to the onedimensional lattice Thirring model. The model is exactly solvable and shows a phase transition with increasing the chemical potential of fermion number, the crossover at a finite temperature and the first order transition at zero temperature, which is similar to the expected property of QCD. In this model, we have obtained all the critical points and examined the thimble structure by inspecting the solutions of the gradient flow equation, the values of the action at the critical points and the Stokes phenomena. And we have identified the set of the thimbles which contribute to the path-integral and have classified the dominant thimbles for given parameters, L, β, m and µ. Our result there suggests that one should sum up the contributions of subdominant thimbles in order to reproduce the rapid crossover and the first-order transition in the low-temperature limit.

JHEP12(2015)125
In this article, we consider the same one-dimensional Thirring model at finite density and perform Monte Carlo simulations taking the most dominant thimble (referred to as J σ 0 in [64]) with the HMC algorithm proposed in ref. [51]. We will examine to what extent the HMC simulation on the single dominant thimble J σ 0 works for this model by comparing our numerical results with the exact ones.
This paper is organized as follows. In section 2, we introduce the one-dimensional lattice Thirring model and apply the Lefschetz thimble method to the model. In section 3, we describe our HMC simulation details and present our numerical results. Section 4 is devoted to summary and discussion.
2 One-dim. lattice Thirring model complexified on Lefschetz thimbles In this section, first we introduce a lattice formulation of the one-dimensional massive Thirring model [21,65] and discuss its property at finite temperature and density. Next we apply the Lefschetz thimble method to this lattice model. The method is based on the complexification of the field variables and the decomposition of the original path-integration contour into the cycles called Lefschetz thimbles. See refs. [5,49,51] for the detail of the approach and ref. [64] for the detail of the Lefschetz thimble structure of the Thirring model

One-dimensional massive Thirring model on the lattice
The one-dimensional lattice Thirring model we consider in this paper is defined by the following action [21,22,[65][66][67], where β = 1/2g 2 a, 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 (T ≡ 1/La).
The fermion field χ f ,χ f has N f flavors and satisfies the anti-periodic boundary conditions: The auxiliary field A n , which should couple to the vector current of the fermion χ f ,χ f , is introduced as a compact link variables e iAn . The partition function of the lattice model is defined by the path-integration, where D denotes the lattice Dirac operator,

JHEP12(2015)125
This is not real-positive in general when µ = 0 and it has the property (det This fact can cause the sign problem in Monte Carlo simulations. We consider the case of N f = 1 for simplicity in the following sections. 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 with the modified Bessel functions of the first kind as (2.5) The number density and condensate of the fermion field are then obtained as follows: The µ-dependence of these observables are plotted in figure 1 for L = 8, ma = 1, and β = 1, 3, 6. It shows a crossover behavior in the chemical potentialμ (in the lattice unit) aroundμ ≃m + ln(I 0 (β)/I 1 (β)). In the limit L → ∞, these quantities reduce to the following forms, where H 1/2 (x) is the Heaviside step function and µ * (L→∞) is the critical density in this limit given by µ * (L→∞) =m + ln(I 0 (β)/I 1 (β)). Figure 2 shows the β-dependence of the critical density at ma = 1/3, 1/2, and 1.
The continuum limit of the lattice model (a → 0) may be defined at a finite temperature as the limit: β = 1/2g 2 a → ∞, L = 1/T a → ∞, while β/L = T /2g 2 fixed. In this limit, the partition function scales as  From these results, one can see that the model shows a crossover behavior in the chemical potential µ for a non-zero temperature T > 0, while in the zero temperature limit T = 0, it shows a first-order transition at the critical chemical potential µ c = m + g 2 . We note that at the zero temperature T = 0, the number density n vanishes identically for µ ≤ µ c , which is sometimes called as the Silver-Blaze behavior [68].

Thirring model complexified on Lefschetz thimbles
Next we consider the complexification of the above lattice model and reformulate the defining path-integral of eq. (2.2) by the complex integrations 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 the complex function given by Then, for each critical point z(= {z n }) = σ given by the stationary condition, the thimble J σ is defined as the union of all the (downward) flows given by the solutions of the gradient flow equation

JHEP12(2015)125
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 other L-dimensional real submanifold K σ of C L associated to the same critical point σ, defined as the union of all the gradient flows s.t. z(+∞) = σ. Namely, the partition function and the correlation functions of the lattice model can be expressed by the formulae, It is not straightforward in general to find all the critical points {σ} and to work out the intersection numbers {n σ } of the associated Lefschetz thimbles {J σ }. Fortunately, in our lattice model, we can obtain all the solutions of the stationary condition eq. (2.12) and therefore all the critical points. In the separated paper [64], we have shown that the critical points can be classified by an integer n − (= 0, 1, · · · , L/2 − 1) as where n − is defined as the number of the components z n which take the value π − z. Moreover, by inspecting the solutions of the gradient flow equation, the values of the action at the critical points {S[σ]} and the Stokes phenomena, we have identified the set of the thimbles which contribute to the path-integral for given parameters, L, β, m and µ. Especially, we found that the dominant thimbles are associated with the critical points of the type n − = 0, z n = z (n = 1, · · · , L), These critical points are shown in figure 3 for β = 3, ma = 1, L = 8 and µa = 0.6 in the complex plane of z ∈ C which parameterizes the field subspace of z n = z (n = 1, · · · , L) in C L . We denote these critical points by the labels σ i andσ i with i = 0, ±1, · · · , ±L/2 also shown in figure 3 (for the case L = 8).
We note that some of the thimbles terminate at the zeros of the fermion determinant,  Among these thimbles associated with the critical points given by eq. (2.19), the most dominant thimble is the thimble J σ 0 , which is labeled by 0 in the figure. It turns out that its value of the action S[σ 0 ] is closest to that of the classical vacuum of the model. In the following numerical study, we consider this most dominant thimble J σ 0 .
3 Hybrid Monte Carlo study of the Thirring model on the thimble J σ 0 In this section, we describe our numerical simulations of the Thirring model performed on the single thimble J σ 0 . First, we review the Lefschetz thimble HMC method proposed in ref. [51], and discuss a few improvements of the method necessary in applying to the (fermionic) Thirring model. Secondly, we summarize the simulation parameter details. Lastly, we present and discuss our simulation results.

Simulation method: hybrid Monte Carlo on Lefschetz thimbles
The hybrid Monte Carlo (HMC) algorithm on Lefschetz thimbles proposed in [51] is a Monte Carlo method to evaluate the path-integral of an observable O[x] over a given thimble J σ , space as δz = U α z δξ α (δz ∈ C L , δξ ∈ R L ). In this HMC algorithm, a series of field configurations {z (k) } (k = 1, · · · , N conf ) are generated with the real-positive weight e −

(S[z]−S[σ])
Jσ through the Molecular dynamics steps constrained to the thimble and the Metropolis accept/reject procedure, while the residual complex phase factor e iφz = det U z is reweighed to the observable as In the algorithm, any field configuration z on the thimble and the associated tangent vectors {V α z n }(α = 1, · · · , L) are computed by solving the flow equations 3 assuming that the solutions take the asymptotic forms in the sufficient past at t = t 0 (t 0 < 0, |t 0 | ≫ 1) as Here e α (α = 1, · · · , L) is a real vector (e α ∈ R; L α=1 e α e α = L), and v α n (α = 1, · · · , L) are the orthonormal tangent vectors at the critical point σ which factorize the Hesse matrix K nm ≡ ∂ n ∂ m S[z σ ] with the real-positive diagonal elements κ α (α = 1, · · · , L): v α n K nm v β m = κ α δ αβ . By this procedure, one can parameterize any field configuration z on the thimble by the set of the parameters, the flow-direction vector e α and the flow-time t ′ = t − t 0 , defining a map (e α , t ′ ) → z ∈ J σ as z n [e, t ′ ] = z n (t)| t=t ′ +t 0 . (3.5) We employ the 4th-order Runge-Kutta method to solve the flow equations and use Diag package [70] to perform the factorization of the Hesse matrix. The molecular dynamics is then formulated as a constrained dynamical system and solved by the constraint-preserving second-order symmetric integrator as

JHEP12(2015)125
In the lattice Thirring model we are considering, the given thimble J σ can terminate at the zeros of fermion determinant. In such a case, the flow reaches a zero within finite time and the flow time t ′ (= t − t 0 ) is bounded. Moreover, the force terms in the flow equations become quite large in the vicinity of the zero. These points cause problems in solving the flow equations or in solving the Molecular dynamics with finite time steps. To improve this situation and to achieve the necessary precision of the solutions, we implement in this work the adaptive step size in the 4th-order Runge-Kutta method: we simply adjust the step size ∆t depending on the size of the force terms F n [z] as |F n [z]|·∆t = L·const. In this respect, an estimate of the error of the solutions can be obtained by using R = |∂S/∂z n −V α z n κ α e α | 2 /2L, which should vanish for an exact solution. We also introduce and adjust a scale parameter λ as z n → λz n to keep the values of the diagonal elements κ α of the Hesse matrix in a reasonable range, for otherwise the exponential growth of the field configurations could be very rapid with a finite step size, the errors in the solutions of the flow equations eqs. (3.3) could become out of control, and the iterate method to solve the constraints eqs. (3.9) could not converge.

Simulation details
The parameter sets in our simulations are summerized as follows. The base simulations were performed for ma = 1, β = 1, 3, 6 on the lattice L = 4, 8 in order to measure and examine the averages of the residual phase, number density and scalar condensate. A series of simulations for L = 8, 16, 32 with L(ma) = 16 and β(ma) = 2, 3 were done in the study of the continuum limit behavior, and a series of simulations for L = 4, 8, 12, 16, 24, 32 with β = 3, ma = 1 were used for the study of the low-temperature limit behavior. For each parameter sets, the chemical potential was varied in the range µa ∈ [0.0, 2.0] with the increment 0.2.
In solving the flow equations by the Runge-Kutta method, we set t 0 = −4. The initial values of the number of steps and the step size are N t = 20 and ∆t = 0.1, respectively. The scale parameter λ is chosen in the range 0.05 ≤ λ ≤ 0.1. With these parameters, the condition R < 10 −5 was satisfied.
For the Molecular dynamics, the trajectory length and the number of steps are set to τ = 0.5 and N τ = 10, respectively. We generated 1,000 configurations for all the parameter sets and estimated errors using the jackknife method with a bin per 20 configurations.

Simulation results
First of all, we show in figure 4 the result on the averages of the residual phase for ma = 1, β = 1, 3, 6 and L = 4, 8. The average Re exp(iθ) sometimes deviates from unity, but it stays greater than 0.8 almost always. The similar results were observed for the larger lattice sizes L = 12, 16, 24, 32. From these results, we can say that the reweighting should work for this model with our choice of the parameter sets.
We next show the results of the number density and the scalar condensate for L = 4 in figure 5 and for L = 8 in figure 6, respectively. At the larger inverse couplings β = 3, 6, our numerical results are in good agreement with the exact results. But at the smaller inverse coupling β = 1, discrepancies are observed in the crossover region on the both lattice sizes    figure 10 of [64] using the "uniform-field model". In figure 7, on the other hand, we show the lattice size dependence of the number density and scalar condensate at ma = 1 and β = 3. We find that the agreement between the numerical and exact results gets worse as L increases from L = 4. The discrepancies become significant for the larger lattice sizes, L = 16, 24, 32, while the contributions of the thimble J σ 0 seem saturated at about L = 12 as shown in figure 8. These results on the lattice size dependence are quite consistent with the analysis shown in figure 12 of [64] based on the "uniform-field model".
Finally, in figures 9 and 10, we show the results on the continuum limit at a fixed temperature. We find that the discrepancies observed in the crossover region persist in this limit. It seems that the size of the discrepancy scales, too.

Summary and discussion
In this paper, we have applied the Lefschetz thimble method to the one-dimensional lattice Thirring model at finite density and performed HMC simulations on the single thimble J σ 0 , which is expected to dominate the path-integral. We have measured the average residual phase, number density and scalar condensate. The average residual phase almost always stays greater than 0.8 and the reweighting works in this model for our choice of the parameter sets. By comparing our numerical results with the exact ones, we have examined to what extent the HMC method works and the single thimble J σ 0 reproduces the exact result.
The numerical results of the number density and scalar condensate reproduce the exact ones at small L ≃ 4, 8 and large β ≃ 3, 6. We also observed that these numerical results scale toward the continuum limit keeping L(ma) and β(ma) fixed. These results imply that the single-thimble approximation with J σ 0 would work in the weak coupling region of g 2 /m ≤ 1/6 and/or in the high temperature region of T /m ≥ 1/8.

JHEP12(2015)125
However, we observed the discrepancy in the crossover region for smaller β and/or larger L. It persists in the continuum limit at a fixed temperature and becomes more significant toward the large L limit, or the low-temperature limit. These numerical results are quite consistent with our analytical study of the model [64]. Our studies clearly show that the contributions of subdominant thimbles should be summed up in order to reproduce the rapid crossover and the first-order transition in the low-temperature limit.
In the Monte Carlo methods formulated on the Lefschetz thimbles, it is not straightforward to sum up the contributions over the set of the relevant thimbles. This is because one need to obtain the relative (complex) weight factors {e −S[σ] Z σ } (see eqs. (2.15)). However, a general method to compute these quantities is not known so far. It is then highly desirable to devise an efficient way to perform the multi-thimble integration by extending the Monte Carlo algorithms for practical applications of the Lefschetz thimble integration to fermionic systems with the sign problem.