Regularization dependence on phase diagram in Nambu-Jona-Lasinio model

We study the regularization dependence on meson properties and the phase diagram of quark matter by using the two flavor Nambu-Jona-Lasinio model. We find that the meson properties and the phase structure do not show drastically difference depending the regularization procedures. We also find that the location or the existence of the critical end point highly depends on the regularization methods and the model parameters. Then we think that regularization and parameters are carefully considered when one investigates the QCD critical end point in the effective model studies.


I. INTRODUCTION
The phase structure of quark matter on finite temperature and density has actively been studied for decades [1]. Under usual condition, meaning low temperature and density, quarks are confined inside hadrons and they never be able to observed as a single particle. On the other hand, due to the nature of the asymptotic freedom [2], quarks and gluons can be free from the confinement at high temperature and density, because the coupling strength becomes weak at high energy. It is, therefore, expected that quark matter undergoes the confined/deconfined phase transition at some temperature and density. This is important subject both in theoretical and experimental studies since it crucially relates to the quark matter properties at relativistically high energy collisions and extremely dense stellar objects such as neutron stars.
The first principle for quarks and gluons is Quantum Chromodynamics (QCD) which is a non-Abelian gauge field theory for fermions. Our goals is to evaluate the phase structure based on this first principle QCD, however, it is difficult to extract theoretical predictions due to the nature of complicated strongly interacting system. One of the most reliable approaches is to use the discretised version of QCD called the Lattice QCD (LQCD) in which theoretical calculation is performed on the discrete spacetime [3]. Although the LQCD works well at finite temperature T for small chemical potential µ ≃ 0, there is the technical difficulty called the "sign problem" when one tries to investigate the system at intermediate chemical potential. There effective models maybe nicely adopted because some models can consistently treat the system at finite temperature and chemical potential.
For the sake of evaluating the phase structure of quark matter at finite temperature and chemical potential, we will employ the Nambu-Jona-Lasinio (NJL) model [4] which is the most frequently used one in this context (there are a lot of nice review papers on the model, see, e.g. [5][6][7][8][9].) The model is constructed by incorporating the four point quark interaction into the model Lagrangian, so it is not renormalizable due to this higher dimensional operator. Therefore, the physical predictions of the model inevitably depend on the regularization procedure and the model parameters chosen. The resulting phase diagram on the T − µ plane is as well affected by the parameters and regularization prescriptions. So it is an important issue to study whether the phase structure obtained in one regularization method is consistent with the ones from different regularization methods.
In this paper, we are going to study the phase structure of quark matter in the NJL model with various regularization ways, which are three dimensional (3D) momentum cutoff, four dimensional (4D) momentum cutoff, Pauli-Villars (PV) regularization, proper-time (PT) regularization, and the dimensional regularization (DR). The 3D cutoff scheme is the most popular method in this model and a lot of works have been done in this way. The 4D cutoff method preserves the Lorentz symmetry in which space and time are treated on equal footing. The Pauli-Villars regularization is based on the subtraction of the amplitude considering the virtually heavy particle to suppress the unphysical high energy contribution coming from loop integrals [10][11][12]. The proper-time regularization makes integrals finite through the exponentially dumping factor [11,13]. The dimensional regularization analytically continues the spacetime dimension in the loop integrals to a non-integer value, then try to obtain finite contribution from the integrals [14]. Beside from the frequently used 3D cutoff way, there have been a lot of works by using the 4D [6,7], PV [6,[15][16][17], PT [6,[18][19][20][21][22][23][24][25][26], and DR [27][28][29][30][31][32]. The physical consequences depend on the regularization [33]. This paper is organized as follows; Section II introduces the model Lagrangian, and show the model treatment on the meson properties and the explicit formalism at finite temperature and chemical potential. In Sec. III, we present various regularization procedures, 3D, 4D, PV, PT and DR prescriptions with explicit equations. We then perform the parameter fitting in Sec. IV. In Sec. V, the numerical results of the meson properties are shown. We then draw the phase diagrams with several parameter sets using various regularization methods in Sec. VI. We also study the phase diagram with the parameters fixed under the condition with the same constituent quark mass in Sec. VII. In Sec. VIII, we give the discussions on the obtained results. Finally, we write the concluding remarks in Sec. IX. Several detailed calculations are shown in Appendix.

II. TWO FLAVOR NJL MODEL
In this paper we consider two light quarks with equal mass. The model has SU L (2) ⊗ SU R (2) flavor symmetry at the massless limit, m → 0.

A. The Lagrangian and gap equation
The Lagrangian of the two flavor NJL model is given by wherem is the diagonal mass matrixm = diag(m u , m d ) and G is the effective coupling strength of the four point interaction. We set m d = m u in this paper. The application of the mean-field approximation leads the following mean-field LagrangianL with the constituent mass m * = m u + σ. The flavor symmetry is broken down, SU L (2) ⊗ SU R (2) → SU L+R (2), by non-vanishing current quark mass, m u , and dynamically generated σ. Thanks to the simple form of the Lagrangian, one can easily evaluate the effective potential, and V is the volume of the system. After some algebra, we see The detailed derivation of the effective potential is presented in [9]. The gap equation is obtained through the extreme condition of the potential with respect to σ, namely, This condition leads the following gap equation with the number of flavors N f and where trace takes the spinor and color indices. This is the key equation in the model because it determines the values of the chiral condensate ψψ and the constituent quark mass m * .

B. Meson properties
The properties of the pion and sigma meson can be studied based on the model with the determined chiral condensate. The interacting Lagrangian of the pion and quarks is written by where τ i are 2 × 2 matrices in the flavor space and π i represent the pion fields. The explicit expression is τ · π = τ − π − + τ + π + + τ 0 π 0 , with τ ± = (τ 1 ± τ 2 )/ √ 2 and τ 0 = τ 3 where τ i are the Pauli matrices. By applying the random phase approximation, we can write the pion propagator as the summation of the geometrical series of the one-loop diagram, which gives where Π π is the following quark loop contribution with the quark propagator The explicit derivation of Eq. (10) is discussed in the review paper [6]. The pion mass is calculated at the pole position of the propagator, so the condition reads It should be noted that the residue at the pole p 2 = m 2 π coincides with the square of the coupling strength g 2 πqq so we have the relation In the similar manner, the sigma meson mass is evaluated at the pole of the propagator, with Therefore the condition which determines the sigma meson mass becomes The pion decay constant is calculated from the following equation The explicit form becomes Thus the equations Eqs. (13), (17) and (19) are the ones which determine the pion mass, sigma mass, and the pion decay constant.

C. Explicit formalism at finite temperature
Since our purpose is to study the phase structure on temperature T and chemical potential µ, we need to extend the equations to finite temperature. According to the imaginary time formalism, the integral region of the temporal component becomes finite due to the periodic or anti-periodic condition of fields as where τ is imaginary time and β is the inverse temperature 1/T . Consequently, continuous integral in the temporal direction is replaced by the following discrete summation, where ω n = 2nπT or (2n + 1)πT depending on the statistical property of field, i.e., for bosons or fermions, and the chemical potential seen in Eq. (20) appears in the way iω n + µ. In this paper, we only treat fermionic quark loop contributions then ω n = (2n + 1)πT is always the case. With the help of the formalism Eq.(21), we see that the gap equation at finite temperature becomes where N c is the number of colors, E = √ k 2 + m * 2 , E ± = E ± µ and f (E) = 1/(1 + e βE ). It is important to note that the contributions can be expressed by the summation of the T independent part (trS 0 ) and T dependent part (trS T ). This characteristic is general if one takes the infinite number of the frequency summation in finite temperature field theory and crucial when we apply the regularization procedures to the appearing integrals.
Since the gap equation is derived by differentiating the effective potential with respect to σ, then the effective potential can be obtained by integrating the gap equation (see, for example, [35]), where we have dropped the suffix in V eff and just written V for notational simplicity. Thereafter the effective potential at finite temperature V = V σ + V 0 + V T is evaluated as It is important to note that, if we apply some regularizations, the results between the direct calculation from Eq. (5) and the one after integrating the gap equation may become different, because regularization essentially means the subtraction and there are several ways of subtractions. Therefore, in this paper, we persistently use the latter way shown in Eq. (25) so that the model treatment becomes consistent. It should be noticed that the finite temperature correction, V T , contains no divergent integral. A finite result can be obtained for the finite temperature correction without applying any regularizations. Next, we carry on the integral in the meson properties; the one loop contribution can be written as with Since trS is already evaluated above, the remaining task is to calculate I(= I 0 + I T ), and it becomes Similarly, the one-loop diagram of the scalar channel can be written as We now have already evaluated all the ingredients of Π σ above in Eqs. (23), (24), (31) and (32), so we do not need further calculations. Finally, let us derive the equation for the pion decay constant. After a bit of algebra we obtain the relation, Here we evaluate f π at p 2 = 0 following [6].

III. REGULARIZATION PROCEDURES
Since the integrals obtained in the previous section include infinities, we need to apply some regularization so that the model leads finite quantities. As mentioned in the introduction, the model is not renormalizable, then the model predictions inevitably depend on regularization procedures chosen. Here, we shall present possible regularization methods in this section.

A. Three dimensional cutoff scheme
The idea of the three dimensional (3D) cutoff is simple; one drops high frequency mode by introducing the cutoff scale Λ 3D into the integrals. We work in the 3-dimensional polar coordinate system and cut the radial coordinate as By performing the integrals, we have for the gap equation σ = 2N f G trS, The effective potential can also be calculated as The quark loop integral in the meson properties I(p 2 ) reads Note that the integral diverges around 4E 2 ≃ p 2 , and we apply the principal integral to avoid this divergence [34]. It may be worth mentioning that the integral can be performed analytically when p 2 = 0 for T = 0, then one has for the pion decay constant, We thus obtain the required quantities in evaluating the phase diagram and meson properties.

B. Four dimensional cutoff scheme
In the four dimensional (4D) cutoff regularization scheme, we introduce the cutoff scale Λ 4D in the Euclidean space after performing the Wick rotation, This is well known four dimensional cutoff method for T = 0 case. As the natural extension to finite temperature, we introduce the cutoff scale by where L 4 is the maximum integer which does not exceed Λ 4D /(2πT ) − 1/2. In the 4D cutoff way, it is difficult to divide the contribution into the temperature independent and dependent parts, since there is also cutoff in the frequency summation.
The explicit form of trS and the effective potential become where ω − n = ω n − iµ. For T = 0, the integral can be performed analytically by using the Feynman parameter method, One should give the special attention in calculating I 4D (p 2 ), because the integral includes divergence to be cured as seen in the 3D cutoff case. The analytic expression of I 0 4D (p 2 ) will be given in appendix A Again we show the explicit form for the pion decay constant at T = 0, (49)

C. Pauli-Villars regularization
In this regularization, the divergences from loop integrals are subtracted by introducing virtually heavy particles as This manipulation induces virtual frictional force so that the contribution from unphysical high frequency mode is suppressed.
In evaluation the gap equation, we apply the following subtraction where By setting the cutoff scales Λ 1 = Λ 2 = Λ PV after the subtraction, we have where E m = √ k 2 + m * 2 and E Λ = k 2 + Λ 2 PV . By integrating the above equation, we obtain the effective potential Since the divergence coming from the integral I(p 2 ) is order of log, one subtraction is enough to make it finite, so we get The pion decay constant at T = 0 becomes

D. Proper-time regularization
The basic idea of the proper-time regularization is based on the following manipulation of the Gamma function, where the lower cut 1/Λ 2 PT induces the dumping factor into the original propagator since, for example with n = 1, 1 Therefore in this regularization high frequency contribution is dumped by the factor e −k 2 E /Λ 2 PT , so the original divergent integral turns out to be finite. For A contains a imaginary part, Eq. (60) is modified as, Under this procedure, the integral of trS in the gap equation becomes where Ei(−x) is the exponential-integral function. For m * 2 ≪ Λ 2 PT , Eq. (63) is expanded as We rotate the contour of the integration in Eq. (64) to the imaginary axis of τ [23][24][25]. For ω 2 0 − µ 2 + m * 2 > 0, the trace becomes and for ω 2 where N = { µ 2 − m * 2 /(πT ) − 1}/2 . Similarly, the effective potential can be calculated through For and for ω 2 I PT can also be calculated by where α is the Feynman integration parameter and ∆ = m * 2 − p 2 /4 + α 2 p 2 . Then the integral can be written where W n (α) = ω 2 n + m * 2 + (α 2 − 1/4)p 2 − µ 2 . The pion decay constant at T = 0 reads the following simple form, For m * 2 ≪ Λ 2 PT , we have

E. Dimensional regularization
In the dimensional regularization method, we obtain finite quantities through analytically continuing the dimension in the loop integral to a non-integer value, D, as where the scale parameter M 0 is inserted so as to adjust the mass dimension of physical quantities. The method is well known since it preserves most of symmetries. Note that this result is the same as the result obtained from the proper-time integral (0 < τ < ∞) and expressed by the poles of the Gamma function.
The trace trS in the gap equation reads where The effective potential becomes In the similar manner, the integral I(p 2 ) appearing in the meson propagator is calculated as Note we need to perform the principal integration for m * 2 < p 2 /4. The pion decay constant at T = 0 reads the following simple form, We show a concrete example of trS 0 DR and f 2 πDR for D ≃ 2, 3, 4 in appendix B.

IV. MODEL PARAMETERS
Having obtained the equation which determines the pion mass and pion decay constant, we are now ready to perform the parameter fitting. In the previous section we suppose that all the cutoff scales are equal in each regularization to reduce the parameters. Thus the model has three parameters: the cutoff scale Λ, the current quark mass m u and the coupling strength G. Whereas in the DR there appears one more parameter, so the total number becomes four: the current mass m u , dimension D, scale parameter M 0 and the coupling G, as discussed in Ref. [36]. In this section, we shall set the model parameters by fitting the pion mass and decay constant. The actual values we use are shown below For the case with the DR, we perform fitting with one more quantity, the neutral pion decay constant to two photons, which will be discussed later.

A. Parameters in various regularizations
Here we align the model parameters in various regularizations in this subsection. Table I, II, III and IV show how the parameters change according to the current quark mass m u , where we first set the value of m u then search the parameters Λ and G which lead m π = 135MeV and f π = 94MeV. One sees the      We also showed the values of the constituent quark mass m * and the chiral condensate ūu 1/3 which are the predicted quantities in the models. We note that the values of m * are about 200 − 300MeV which are comparable to one third of the proton mass. We also note that m * increases with respect to m u , while the absolute value of ūu decreases when m u becomes large. Since the relation m * ∝ G ūu holds, even ūu becomes smaller m * can be larger due to the large value of G, which is actually the case in these regularizations.
We plot the parameters of each regularization in Fig.1. The black circles denote the value which satisfy m * = 311MeV for each regularization. The relation between the cutoff scale Λ and G for the 4D cutoff, Pauli-Villars regularization and proper-time regularization resembles each other. The relation between m u and Λ for these regulatizations also resembles each other. In the case of same value for m * , the m u dependence of Λ is large. However, the values of Λ are close, 660 − 830MeV. In the case of same value for m * , the m u dependence on G is large and the values of G are separate in each regularization. The relation between m u and G is different in each regularization. PV reg.
PT reg.
PT reg.
PT reg.

B. Parameter fitting in the dimensional regularization
For the sake of the parameter fitting in the DR, we present the calculation of the pion to two photon decay rate, Γ π 0 γγ , in this subsection. The decay rate can be evaluated through the following one-loop amplitude, T µν (k 1 , k 2 ), where e is the QED coupling constant and k 1 and k 2 are the external momentum of emitted photons so the square of the total momentum coincides with that of the original pion, namely, (k 1 + k 2 ) 2 = m 2 π . By using T γ , the decay rate, Γ π 0 γγ , is expressed as The detailed derivation is presented in the paper [27]. After some algebra, one obtains with α e = e 2 /(4π), andm 2 π = m 2 π /m * 2 . With the observables, m π = 135MeV, f π = 94MeV and Γ π 0 γγ = 7.8 eV, we perform the parameter fitting in the DR following [27]. Table. V shows the fitted parameters in the DR case. We note that both the constituent quark mass and chiral condensate grow up with increasing m u , which is the characteristic feature in this regularization [36]. We have presented the required equations in the model, then set the parameters for various regularizations. It is now ready for the actual numerical analysis on the model predictions. Here we shall show the thermal meson properties, which are the pion mass, the pion decay constant, and the sigma meson mass.
At finite temperature, there are two ways of the application of each regularization; one is to apply the regularization only for the temperature independent contribution because the temperature dependent contributions are always finite due to the characteristic factor of the Fermi-Dirac distribution, i.e., f (E). The other is to apply the regularization both for the temperature independent and dependent parts, since the regularization essentially relates to the cutoff of the model so the introduction of the same cutoff clearly determines the model scale. On the other hand, the former method retains more symmetry of the model. Then the physical meaning of these prescriptions are that the former one respects the model symmetry, and the latter does the cutoff scale of the model. Since the model is not renormalizable, the predictions depend on the regularization ways and our purpose in this paper is to study the regularization dependence on the model. Therefore we shall study all the cases and compare the results among various regularizations.

A. Results with regularizing T -independent contribution
In this subsection, we show the results of meson properties based on the procedure of applying the regularization only to the temperature independent part. The required integrals are trS and I(p 2 ), and we evaluate the following combination, where trS 0 Reg is trS for T = 0. The lower index indicates each regularization, namely, trS 0 3D , trS 0 4D , trS 0 PV , trS 0 PT and trS 0 DR . For the temperature dependent part trS T , here we use the form shown in Eq. (24). Similarly for I(p 2 ), we use the equivalent expression, with I 0 Reg for each regularization way. Figure 2 shows how the pion mass changes with respect to T and µ for various parameter sets in the previous section. It should be noted that, for some parameter sets, no real solution exists at high temperature as seen in the case with the DR and m u = 3.0MeV [31]. We observe the similar behavior in each regularization; the pion mass remains almost constant for low T and µ, then raises up for higher T and µ. This comes from the fact that the chiral symmetry is broken at low T and µ and restores at high T and µ. The pion has smaller mass when the symmetry is broken due to the Nambu-Goldstone theorem, while the mass becomes large after symmetry restoration. We see that the mass starts to increase around 170MeV which is comparable to the critical temperature for the chiral symmetry breaking. We see that the temperature and chemical potential where the pion mass glows up become larger with respect to m u for the 3D cutoff, 4D cutoff, PV and PT, while they become smaller for the DR. We also see that the discontinuity seen around the transition temperature is considerably larger in the DR case compare to the other regularizations. Then we expect that the tendency of the first order phase transition is strong for the DR comparing the other regularizations. We will discuss it in the next section.
The results of the pion decay are shown in Fig. 3. One sees the similar tendency as well; the decay constant is almost constant for low T and µ, and it decreases when T and µ exceed certain values which are around T ≃ 170MeV and µ ≃ 300−400MeV. It is interesting to note that the decay constant drops discontinuously at high µ for some parameter sets in the 3D, 4D, PV, PT regularizations, while the discontinuity is always the case in the DR. The existence of the gap is the signal of the first order phase transition, and the tendency becomes stronger with increasing m u . This is because the coupling strength is larger for higher m u as seen from the parameter tables, so quarks have stronger correlations when the parameter m u is larger in the 3D, 4D, PV, PT cases.
Having fixed the parameters with the pion mass and decay constant. We will calculate the sigma meson mass. It is one of the predictions of the model. Figure 4 shows the numerical results of the sigma mass. At T = 0 and µ = 0 we find the band around 400 − 900MeV in 3D, 4D, PV, PT regularizations, 900 − 1400MeV in the DR. Then, in the DR case, the predicted values are larger than the experimental value m σ ≃ 500MeV [37], in the leading order of the 1/N c expansion. As is known we can obtain much smaller sigma meson width in this order. We should also check the next to leading order contributions for the sigma meson mass and width. The features of the curves are the similar to that of the pion; the mass decreases with increasing T and µ until some values, then it increases after exceeding the certain values. As seen in the pion mass case, the solution of the sigma mass on the real axis disappears for some parameter set at high temperature.

B. Results with regularizing all contribution
In this subsection, we study the meson properties by applying each regularization procedure to both temperature dependent and independent contributions. It is worth mentioning that the case with 4D cutoff method does not give credible results because enough number of frequency summations is not taken in this method. The cutoff scale of the 4D case is around 1GeV, which means that at T = 100MeV only four terms in the Matsubara mode summation, n = −2, ±1, 0 (ω n = (2n + 1)πT ), are taken into account. It is known that finite temperature field theory does not lead reliable predictions when the number of the frequency summation is small [38]. Therefore, we will not show the results in the case of 4D cutoff scheme here, and consider the other four cases 3D, PV, PT and DR and call these cases 3DRT, PVRT, PTRT and DRRT, respectively. It is also worth mentioning that the calculations technically become impossible at T = 0 in the PTRT case as can be read from Eqs. (67) and (72). Then, we will show the results with T = 10MeV as the representative values on µ dependence at low T . Figure 5 displays the results of the pion mass. One sees that the qualitative feature does not change comparing to the previous case with regularizing only the temperature independent contributions. Quantitatively, we note that the changes become smoother at high T and µ. This can easily be understood because the regularization procedure suppresses the thermal contribution, so the finite temperature term reduces to give smoother curve with respect to T and µ.
We aligned the results of the pion decay constant in Fig. 6. As seen in the pion mass case, the numerical results do not alter qualitatively, the curves become smooth. Note that, although the T dependence becomes considerably smother, the transition chemical potential does not change. This is due to the fact that finite temperature contributions become proportional to the step function θ(µ − m * ) for T = 0, then the transition chemical potential is not affected by the regularization procedure in this model treatment.
The predictions on the sigma meson mass are shown in Fig. 7. The deviations from the results in the previous subsection in Fig. 4 are more or less similar to the deviations on the pion mass; the curves become flatter with respect to T while µ dependence does not indicate much difference for the 3D, PV and PT cases. However, there appears substantial difference between the cases of the DR and DRRT where the values of the sigma meson mass at µ = 500MeV are around 1500 − 2500MeV for the DR and around 600MeV for the DRRT case. This comes from the difference of the integral values between these two cases, which we will discuss in more detail in Sec. VIII.
We find that the solution on the real axis always disappears for high µ and low T (= 10MeV) in the PTRT with m u = 10, 15MeV. This is because for high µ some quantities become pure imaginary number in the calculation due to the complicated counter integral of I PT (p 2 ) as seen in Eq. (72). Then one can not find a real solution in that case. This is the numerical reason why the meson properties behave badly for high µ in PTRT case.

VI. PHASE DIAGRAM
We shall draw the phase diagram in this section. We search the phase transition point by the condition that the maximum change of the chiral condensate with respect to T and µ. In more concrete, we numerically differentiate the condensate with respect to r ≡ T 2 + µ 2 , then the condition can be written We should be careful on the case with the first order transition, because the condensate has the discontinuous point where we need to find the minimum of the thermodynamic potential then determine the chiral condensate. No matter whether the phase transition is the first order or cross over, we can also use the above criterion since at the first order point, d ūu 1/3 /dr = ∞ holds, therefore it is consistent in both cases. By searching the maximum number of the differentiate of the chiral condensate following the condition Eq. (92), we draw a phase boundary for the chiral symmetry. Figure 8 shows the phase diagrams for the various parameter sets and regularizations. One sees that the critical temperatures, T c is found between 150 and 250MeV at µ = 0 if we regularize only the temperature independent parts (left four panels and the bottom panel). On the other hand, a higher critical temperature is observed when the regularization is applied to the temperature dependent and independent parts. The regularizations 3DRT and PTRT give a critical temperature around T c ≃ 200MeV at µ = 0, the PVRT induces a higher critical temperature near T c ≃ 400MeV with m u = 15MeV, and the DRRT indicates it around T c = 300 − 400MeV. One also sees that the critical chemical potential, µ c , has no large dependence on the application of the regularization to finite temperature term, because the terms are dominated by the step function θ(µ − m * ) as mentioned in the previous section. Consequently, the area of phase boundary enlarges in the T direction when we apply the regularization to the temperature dependent term, which is numerically confirmed by the figure.
We think the most interesting comparison from the figure is on the existence of the critical end point where the first order phase transition starts on the phase boundary. For 3D, 4D, PV, PT, no critical end point appears for the parameter sets with small m u . While the critical end point always appears in DR. Then we can numerically conclude that the DR has stronger tendency of the first order phase transition comparing with the other regularizations. The existence of the critical end point in the other four regularizations can be understood by seeing the value of the parameter G. Briefly speaking, the critical end point appears when G is large. This is physically reasonable because G represents the strength of the correlation between quarks, then the larger G makes the condensation stronger.

VII. PHASE DIAGRAM WITH FIXING m *
We have seen the phase diagram and how the location of the critical end point depends on the parameters in the various regularizations. Considering the fact that the transition T and µ is essentially determined by the value of the constituent quark mass since its dependence appears in the thermal distribution, f (E), with E = √ k 2 + m * 2 , we think it may as well be interesting to compare the phase diagram with the parameter sets which lead the same value of m * (= 311MeV) at T = 0 and µ = 0 instead of fixing the current quark mass, m u .
First, we see the results based on the one with regularizing only the temperature independent parts, 3D, 4D, PV, PT and DR. The behaviors of m * for each regularization are shown in Fig. 9. All the results have no large difference, because the gap equations of the temperature independent part for each regularization has the similar behavior. The gap equations contain the following form, in the regularizations, 3D, 4D, PV and PT. While in DR m * and Λ are replaced by M 0 and D in some parts (see appendix B). In Fig. 9, we note that the results almost coincide in three regularizations, 3D, 4D and PT. The behavior in PV shows a smaller and DR gives steeper slope than the others. In Fig. 10 the phase diagram is illustrated in each regularization. The phase boundaries of 3D, 4D and PT show almost equivalent behavior. The area of the chiral symmetry broken phase for PV is larger than the others. This tendency comes from the behavior of m * in Fig. 9. Since the Λ PV is entered in the form of the dynamical mass, the chiral symmetry breaking contribution is enhanced than the other regularizations. The area for the chiral symmetry broken phase for the DR is smaller than the others. The critical end point for the DR locates higher temperature than the one for 3D, 4D and PT. These tendency also comes from the behavior of m * in Fig. 9.
Next, we discuss the results with regularizing both the temperature independent and dependent parts, 3DRT, PVRT, PTRT and DRRT. The behavior of m * for each regularization is shown in Fig. 11. We find that the finite temperature effect becomes smaller or softer than 3D, PV, PT and DR, respectively. The behavior of m * in DRRT is the most closest behavior to the case of its regularizing only temperature independent part.
The chemical potential has the similar contribution for 3DRT, PTRT, and PVRT. The behaviors of m * in DRRT and DR have large difference. This difference is caused by the reduction of the momentum integral dimension, D, from 4 to 3.32.
The phase diagram in each regularization is shown in Fig. 12. To compare with Fig. 10, the region of the broken phase enlarges in PTRT, 3DRT and PVRT for a low chemical potential. However, the critical chemical potentials in PTRT, 3DRT and PVRT for a low temperature are almost equivalent to that in PT, 3D and PV, respectively. From the behavior of m * in Fig. 11, we observe a larger critical temperature and chemical potential for DRRT. We display the location of the critical end points (µ CP , T CP ) for each regularization in Table.VI.

VIII. DISCUSSIONS
We have introduced the regularization methods then studied the meson properties and phase diagram in previous sections. In this section, we are going to present more detailed discussions on the obtained results.
In comparing the left panels v.s. right panels in Fig. 8, one can observe the contribution to apply the regularization procedure to the thermal correction. The phase diagram does not show considerable difference for 3D and PT cases, while in PV the area of the phase boundary becomes larger in the right panel when m u is large, and the areas in DRRT case is larger than DR case for all the parameter sets. These difference can be understood through the following discussion.
The loop integral, I, essentially has the following subtracted forms for 3DRT, PVRT, and PTRT, Where the typical form of F (k, m * ) is given by F (k, m * ) = Cf (E)/E with some constant value, C. Thus the subtracted terms basically relate to the suppression on the high energy contributions which are expected to be small. We numerically confirmed that the subtracted parts are small for almost all the cases, then the phase diagrams does not change drastically. However, in the PV case with large m u , the difference between F (k, m * ) and F (k, Λ) is small since the constituent quark mass becomes comparable to the cutoff scale, e.g., m * = 417MeV and Λ PV = 729MeV for m u = 15MeV in PV. Consequently, the thermal contribution strongly suppressed in the PV case with large m u . We saw that the area of the phase boundary does not alter so much in 3D, PV and PT regularizations. The essential reason is that the infinities appearing from loop integrals are subtracted at high energy. However in DR, the situation is different since the integral is replaced as so this modifies the integral kernel rather than the subtraction of high energy modes. This is the reason why DR shows considerable difference, if we apply the regularization to both temperature independent and dependent contributions. We also saw the location (or existence) of the critical end point is non-trivial. In Fig. 10, the diagram has the critical end point for 3D, 4D, PT and DR. No critical end point appears in the PV. Thus we find that the PV has weaker tendency of the first order phase transition than that of the other regularization methods. Particularly, the temperature of the critical end point in the DR case is higher then others, which may enable us to conclude that the DR has the stronger tendency of the first order phase transition.

IX. CONCLUDING REMARKS
We have studied the regularization dependence on the phase diagram of quark matter on T − µ plane by using the NJL model. We have first presented the regularization procedure at finite temperature and chemical potential, then fitted parameters within various regularization methods. Thereafter, we have studied the meson properties and the phase structure.
We find that the model produces the reliable predictions on the meson properties whose behavior for finite temperature and chemical potential does not alter drastically, which indicates that all the regularizations employed in this paper nicely capture physics on the meson properties. We can conclude that the regularizations are safely adopted. In this context the regularization parameter independent approach is also interesting [39,40].
It is expected that observation of the critical end point can distinguish a suitable regularization for an effective model of QCD by comparing with the resulting phase diagrams. The important difference is the existence of the critical end point. The model predicts that the critical end point appears at intermediate chemical potential around µ ≃ 300 − 400MeV. This density coincides the one in which different quark state such as color superconductivity may occur, there the order of the phase transition might affect crucially on such the dense states. Moreover the color superconductivity may be realized in the dense stellar objects, like quark stars and neutron stars [33]. Therefore the study of the order of the phase transition has important meaning as well in cosmological observations. So we believe that the further and more extensive investigations are necessary on this subject. Since the I 0 (p 2 ) integral can be evaluated analytically, we will present the explicit expression for various regularizations.
One needs special care in performing I(p 2 ) integral since it contains divergent contribution as seen in the 3D cutoff scheme. I 0 4D becomes for m * 2 > p 2 /4, and for m * 2 < p 2 /4, where Concerning on I 0 PV (p 2 ), it is convenient that we divide the integral as where I

0(ΛPV) PV
is subtracted part in the original integral and it becomes with As seen above, we need to separately evaluate the integral I 0(m) PV depending on the values of m * 2 and p 2 . It becomes for m * 2 > p 2 /4, and for m * 2 < p 2 /4,