Control the model sign problem via path optimization method: Monte-Carlo approach to QCD effective model with Polyakov loop

We apply the path optimization method to a QCD effective model with the Polyakov loop at finite density to circumvent the model sign problem. The Polyakov-loop extended Nambu--Jona-Lasinio model is employed as the typical QCD effective model and then the hybrid Monte-Carlo method is used to perform the path integration. To control the sign problem, the path optimization method is used with complexification of temporal gluon fields to modify the integral path in the complex space. We show that the average phase factor is well improved on the modified integral-path compared with that on the original one. This indicates that the complexification of temporal gluon fields may be enough to control the sign problem of QCD in the path optimization method.


I. INTRODUCTION
Investigation of the phase structure of Quantum Chromodynamics (QCD) is one of the important subjects in the study of the hot and dense QCD matter. If we directly obtain the phase diagram at finite temperature (T ) and density from the first-principles calculation such as the lattice QCD simulation, there are no fog in the exploration. Lattice QCD, however, has the sign problem at finite real chemical potential (µ) and then we cannot obtain the reliable results at finite density. To circumvent the sign problem, several methods have been proposed such as the Taylor expansion method [1][2][3], the reweighting method [4,5], the analytic continuation method [6][7][8], the canonical approach [9][10][11][12][13] and so on. However, we cannot go beyond the µ/T ∼ 1 line by using those methods; for example, see Ref. [14].
Recently, new ideas have been applied to attack the sign problem such as the complex Langevin method [15,16] and the Lefschetz-thimble method [17][18][19]. Both methods are based on the complexification of dynamical variables. The complex Langevin method is based on the stochastic quantization and thus it does not have the sign problem, in principle. In the Lefschetz-thimble method, one should solve flow equations to construct the new integral path which is corresponding to the steepest descent trajectory; this trajectory is so called the Lefschetz thimble. Unfortunately, these methods have some serious problems and thus it is difficult to apply it to QCD at high density: In the complex Langevin method, there is the possibility that it is converged to wrong results due to the excursion and singular problems [20,21]. In comparison, the Lefschetz-thimble method has the global sign problem when multi Lefschetz-thimbles contribute to the path integral and then there is the serious cancellation between them; for example, see Ref. [19]. In addition, we should evaluate the Jacobian induced by the modification of the integral path and it leads the serious increase of the numerical cost. Furthermore, the Jacobian induces the residual sign problem because the oscillation of the Boltzmann weight arises. Also, we may face the problem to draw thimbles when the classical or the effective action has singular points [22].
In Ref. [23], we have proposed the new method which we call the path optimization method (POM) motivated by the Lefschetz-thimble method. The path optimization method is strongly improved in Ref. [24] by introducing the feedforward neural network to optimize the modified integral path. In this method, we first complexify the variables of integration as in the Lefschetz-thimble method. Then, the modified integral path is constructed to minimize the cost function which reflects the seriousness of the sign problem. Therefore, we can treat the sign problem as the optimization problem. The sign problem in the simple one-dimensional integration [23,25] and the complex λφ 4 theory [24] are found to be under controlled. Unfortunately, the numerical cost of the path optimization method is still heavy and thus we cannot apply it to QCD yet, but this method has large extensibility compared with the Lefschetz-thimble method and thus we can still dream to apply it to QCD in the future.
Because of several difficulties in QCD as mentioned above, QCD effective models have been widely used to investigate the QCD phase structure because such effective model is much easier than original QCD. The Polyakovloop extended Nambu-Jona-Lasinio (PNJL) model [26][27][28] is one of the famous and powerful effective models. Unfortunately, the PNJL model has the model sign problem [21,29] even in the mean-field approximation because one particular global minimum perfectly dominate the path integration and then the thermodynamic potential can be complex. It should be noted that this complex nature of the thermodynamic potential is not related with instabilities; the correct thermodynamic potential should be real.
In this study, we consider the PNJL model and the path optimization method to circumvent the model sign problem. This article is organized as follows. In Sec. II and III, we explain the formulation of the PNJL model and the path optimization method, respectively. The numerical results are shown in Sec. IV. Section V is devoted to summary.

II. FORMULATION
In this work, we employ the PNJL model as the QCD effective model. The PNJL model can describe the chiral symmetry breaking/restoration and the approximate confinement-deconfinement transition. Also, it can well reproduce QCD properties at finite imaginary chemical potential which is a big advantage from the viewpoint of the topologically determined confinement-deconfinement transition [30][31][32]. The following formulation and the computation of the Monte-Carlo PNJL model is based on Ref. [33].

A. Polyakov-loop extended Nambu-Jona-Lasinio model
The Euclidean action of the PNJL model is where β = 1/T , V is the three-dimensional spatial volume, q denotes the two-flavor quark-fields, m 0 does the current quark mass, D ν = ∂ ν − iA ν δ ν4 , G is the coupling constant of the four-fermi interaction, Φ (Φ) means the Polyakov loop (its conjugate) and V Φ expresses the gluonic contribution.
With the homogeneous auxiliary-field ansatz after the Hubbard-Stratonovich transformation, the effective action is simplified as Γ = βV V. Then, the effective potential is expressed as where V NJL is the contributions of the NJL part. In the actual calculation, we employ the Polyakov gauge, The actual form of V NJL becomes where N f (N c = 3) is the number of flavor (color), Λ is the three-dimensional momentum cutoff and The auxiliary fields are redefined as with π 0 = π 3 and π ± = (π 1 ± iπ 2 )/ √ 2. The functional form of the Polyakov loop and its conjugate are where and then A 4 are diagonalized because we use the Polyakov gauge; . In this paper, we choose the logarithmic Polyakov-loop potential [34] as the gluonic contribution; where with The parameters should be set to reproduce the lattice QCD data in the pure gauge limit.

B. System volume and parameters
In this study, we follow the lattice formalism and thus V can be expressed [33] as where N s (N t ) are the number of spatial (temporal) lattice sites and a is the lattice spacing. Then, V should depend on the temperature. In this article, we use the homogeneous ansatz and thus our Monte-Carlo simula-tion reaches the mean-field results in the infinite volume limit. This fixed k treatment with homogeneous ansatz leads the inconsistent results, in principle. However, such inconsistency becomes smaller and smaller when the system volume becomes larger and larger and thus it is a minor problem in this study. Full simulation of the lattice PNJL model is our future work. The present PNJL model has three parameters, G, m 0 and Λ in the NJL part. The actual values of the parameters are taken from Ref. [35]; m 0 = 5.5 MeV, G = 5.498 GeV −2 and Λ = 631.5 MeV. The parameters in the Polyakov-loop potential is taken from Ref. [34]; with T 0 = 270 MeV.

III. PATH OPTIMIZATION METHOD
To deal with the model sign problem appearing at finite density, we here introduce the path optimization method.

A. Introduction to POM
In the path optimization method, we start from the complexification of the variables of the integration, x i ∈ R → z i ∈ C where i = 1, · · · , n with the dimension of integration n. To construct the new integral path in the complex plane, we prepare the cost function which should be related with the seriousness of the sign problem. The functional form of the new integral path is represented by using the feedforward neural network with some parameters which are optimized via the minimization of the cost function. Since the neural network even with the mono hidden-layer can approximate any kind of continuous function on the compact subset as long as we can use sufficient number of units in the layer [36,37], the neural network seems to be suitable for the path optimization method. Details are shown in Ref. [23].
In the path optimization method, we represent z i by using the parametric quantity (t) as where w and b are parameters. Parameters w, b and also parameters in f are determined by using the backpropagation algorithm [38]. Figure 1 shows the schematic figure of the neural network used in this study. In the back-propagation, we choose tanh for the activation function. In the following, we represents parameters in the neural network as c = {c i }. Details are shown in Refs. [23,24]. It should be noted that the path optimization method leads the same results with the original theory because of the Cauchy(-Poincare)'s theorem as long as the integral path does not go across singular points, Je −Γ = 0. In the complex Langevin method, singular points induce the problem because it leads the singular drift term in the Langevin-time evolution, the path optimization method can care for it.
In this study, we complexify A 3 and A 8 and then σ and π are still treated as real variables. Since it is known that the model sign problem can be resolved by the complexification of the temporal gluon fields in the Lefschetzthimble method [21], our treatment should work. In addition, such treatment is consistent with the calculation which imposes the CK symmetry on the fermion determinant at finite density [39,40] where C and K denote the charge conjugation and the complex conjugation, respectively. Unfortunately, the Lefschetz-thimble method is difficult to apply to the PNJL model because we should solve flow equations and it sometimes encounter singular pints of the effective potential. The NJL model with the repulsive vector-type interaction is a typical example [22]. In comparison, the path optimization method can avoid the problem and thus it is suitable for the PNJL model analysis.
It should be noted that the present path optimization with the machine learning is unsupervised learning because we do not need teacher data to obtain optimized parameters in the neural network. We try to increase the average phase factor compared with that in the previous optimization step. The one attempt to introduce the supervised learning to the study of the sign problem has been done in Ref [41] based on the generalized Lefschetz-thimble method [42].

B. Optimization process
We use the following cost function in the calculation; where θ(t) = arg(J(t)e −Γ(z(t)) ), θ 0 = arg(Z), with the partition function Z. In the equation, z i represent complexfied A 3 and A 8 and also the real σ and π. Thus, we have eight dynamical variables. If we wish to take care of the periodicity of the effective potential for Re A 3 and Re A 8 , we should use periodic functional form for those as inputs of the neural network like as Ref. [41]. In this study, configurations are well localized in the range −π ≤ Re A 3 /T ≤ π and −π ≤ Re A 8 /(T √ 3) ≤ π by using the simple complexification of A 3 and A 8 . Thus, we do not introduce the periodic form of inputs in this study. This cost function is proportional to e iθ and thus it reflects the seriousness of the sign problem [24]. In the physical system, θ 0 should be 0. If we consider θ 0 = 0, we can apply the path optimization method to other systems such as the complex chemical potential.
In the hybrid Monte-Carlo method, we should evaluate the expectation value of the cost function. To do this, we replace the cost function as where P is the appropriate probability distribution.
In the back-propagation procedure, we need the derivative of the cost function by c i . We represent it as dF i below. After straightforward calculations, we finally reach the expression; In the hybrid Monte-Carlo method, we rewrite Eq. (17) as similar to Eq. (16) with P. This cost function is responsible to the alignment of the Boltzmann weight with each other on the modified integral path if the points are relevant to the path-integral. To make our optimization easier, we employ the mini-batch training. The configurations are divided as N config = nN batch where N batch is the batch size and then the learning is performed batch by batch. To include all updation of each batch, the parameters in the feedforward neural network is updated by replacing dF i by its mean-value as In one optimization step, we update n-times with N batch configurations.
In this study, we use the simple feedforward neural network with the hidden mono-layer. Input is the original integral path and output is its imaginary part. For the optimizer, we employ Adam algorithm [43]; where j is the fictitious time step, dF i means ∂F /∂c i and β 1 and β 2 are the smoothing factors of the exponential moving average. This algorithm is based on the AdaGrad algorithm [44] and the momentum method with preventing the learning weight decay.

C. Simulation setup
The number of units in the hidden layer is set to N unit = 4. For Adam algorithm, we use η = 0.001, α = 0.999, β = 0.9 and ǫ = 10 −8 . In the mini-batch training, we set to N batch = 10. These parameters are so called hyper parameters in the machine learning. Initial values of parameters in the neural network are prepared based on Xavier initialization [45].
In the calculation of the expectation values of operators, we have generated 80000 configurations analyzed each 50 trajectories for each optimization step. Then, the expectation values are estimated after 2 times optimization. Statistic errors are obtained by using the Jack-knife method where the bin number is set to 10. The expectation value of an operator (O) is obtained via the phase reweighting as In this study, we calculate the chiral condensate and the Polyakov loop. The T -dependence of σ and Φ at µ = 0 where σ is normalized by that at T = µ = 0 in the infinite volume limit. The dotted lines are mean-field results in the infinite volume limit as the eye guide. The imaginary part of the effective potential in the A3-A8 plane with σ 2 = π 2 = 0. The effective action is normalized as βV V with k = 64.

IV. NUMERICAL RESULTS
The T -dependence of the chiral condensate and the Polyakov loop at µ = 0 is shown in Fig. 2. The meanfield results in the infinite volume limit is also shown as the eye guide. Because of the finite size effect, σ is deviated from the mean-field results in the infinite volume limit and this result is consistent with results obtained in Ref [33]. In the present calculation, V depends on T and thus the finite size effect becomes serious when T increases. The expectation values of π are consistent with zero in 2σ error.
To discuss the model sign problem, we consider the finite µ below. Figure 3 shows the imaginary part of the effective action. We can clearly see that the effective potential can become complex if we pick up a certain point as the solution of the mean-field approximation.
The µ-dependence of the average phase factor on the original integral-path at T = 100 and µ = 300 is shown in  It can be clearly seen that e iθ becomes small around the chiral transition region, T ∼ 335 MeV, and then the model sign problem seriously appears. In the present computation, we use k = 64 and thus e iθ does not become small so much, but it will be worse when we consider larger k because e iθ is exponentially suppressed by k.
To improve the average phase factor, we use the path optimization method. Figure 5 shows the average phase factor at each optimization step with T = 100 and µ = 300 MeV as an example. After one optimization step (N opt = 1), the average phase factor becomes worse. Since the starting configurations used in the first optimization are generated on the original integral-path and thus those are expected to be located far from the relevant points of the optimized integral-path and then the average phase factor becomes worse temporally. How- ever, the average phase factor is significantly improved after few optimization steps. In the case at T = 100 and µ = 300 MeV with k = 64, N opt = 2 are enough to make the average phase factor close to 1. It should be noted that the present tendency of the improvement depends on several initial conditions of parameters and thus we should try with several set of parameters when the optimization is not well worked. Also, an increase of units in the hidden layer may lead to better convergence [46]. The µ-dependence of the chiral condensate and the Polyakov-loop after performing the path optimization are shown in Fig. 6. The mean-field results in the infinite volume with imposing the CK symmetry to the fermion determinant is also shown as the eye guide. We can correctly reproduce the relation Φ = Φ with Φ , Φ ∈ R in the PNJL model at finite µ by using the path optimization method. These results mean that the path optimization method can well work in the QCD effective model which has the model sign problem. Also, we can expect that the complexification of the temporal gluon field may be sufficient to control the sign problem in the lattice QCD with the path optimization method since the lattice QCD and the PNJL model share similar properties about the sign problem.

V. SUMMARY
In this article, we apply the path optimization method to the Polyakov-loop extended Nambu-Jona-Lasinio (PNJL) model to circumvent the model sign problem. This study is the first attempt to apply the path optimization method to the effective model with dynamical quarks based on QCD. The PNJL model can describe the chiral phase transition and also approximated confinement-deconfinement transition and thus it is a good starting point to investigate the QCD phase struc-ture at finite real chemical potential (µ). Therefore, we choose it as the typical QCD effective model in this article.
The temporal components of the gluon field, A 3 and A 8 , are complexified and the modified integral path in the complex space is expressed by using the feedforward neural network by minimizing the cost function which reflects the seriousness of the sign problem. Then, the sign problem becomes the optimization problem. Parameters in the feedforward neural network are optimized via the back-propagation method. The neural network tries to increase the average phase factor compared with the previous optimization step and thus it is nothing but the unsupervised learning. The scalar and pseudo-scalar auxiliary fields are treated as real valuables. This treatment is motivated by the Lefschetz-thimble analysis done in Ref. [21]. In the actual optimization process, we use the mini-batch training with Adam algorithm.
We have shown that our treatment of variables of integration works well; the average phase factor is sufficiently improved after the optimization at finite µ. This means that the path optimization method can resolve the model sign problem based on the hybrid Monte-Carlo method. In this study, we use the homogeneous ansatz for the auxiliary fields and thus we cannot go beyond the mean-field approximation, but it is a first step to correctly treat the sign problem in the QCD effective models. By considering the straightforward extension of the present formulation, we will go beyond the mean-field approximation of the QCD effective models. We leave the actual simulation as our future work. From these results, we may expect that the complexification of the temporal gluon fields is the sufficient way to control the sign problem in the lattice QCD with the path optimization method because QCD and the PNJL model share similar properties about the sign problem.