Phase structure of NJL model with weak renormalization group

We analyze the chiral phase structure of the Nambu–Jona-Lasinio model at ﬁnite temperature and density by using the functional renormalization group (FRG). The renormalization group (RG) equation for the fermionic effective potential V (σ ; t) is given as a partial differential equation, where σ := ¯ ψψ and t is a dimensionless RG scale. When the dynamical chiral symmetry breaking (D χ SB) occurs at a certain scale t c , V (σ ; t) has singularities originated from the phase transitions, and then one cannot follow RG ﬂows after t c . In this study, we introduce the weak solution method to the RG equation in order to follow the RG ﬂows after the D χ SB and to evaluate the dynamical mass and the chiral condensate in low energy scales. It is shown that the weak solution of the RG equation correctly captures vacuum structures and critical phenomena within the pure fermionic system. We show the chiral phase diagram on temperature, chemical potential and the four-Fermi coupling constant.


I. INTRODUCTION
Solving Quantum Chromodynamics (QCD) is a challenging problem in elementary particle physics. Due to the strong dynamics of QCD in macroscopic scales, the phenomena such as the quark confinement and the dynamical chiral symmetry breaking (DχSB) are predicted. The exploration of the phase structure and the equation of state in hot and dense QCD are quite important for understanding the evolution of quark matters in the universe.
Numerous studies have been performed and have shown that QCD at finite temperature and density has various phase structures, e.g., the quark-gluon plasma and the color superconductivity. In particular, the study by lattice Monte Carlo simulation has substantially contributed to the understanding of nature of QCD. At present, lattice QCD, however, suffers from the sign problem. Since the Boltzmann weight at finite density has imaginary parts, the generated configurations oscillate and thus the statistical errors become large.
In contrast, the analytical methods, e.g., the mean-field approximation (MFA) and the Schwinger-Dyson equation (SDE) with the ladder approximation, have been applied to QCD and have succeeded to understand its dynamics and phase structures both qualitatively and quantitatively. Furthermore, effective theory approaches have also helped us to understand QCD. The Nambu-Jona-Lasinio (NJL) model [1, 2] and its extended models [3,4], especially, have also played crucial roles as effective models describing the DχSB and the confinement in QCD; see also [5][6][7]. The Lagrangian of the NJL model is given by the following four-Fermi interaction structure: This Lagrangian is invariant under the chiral U(1) transformation: ψ → e iγ 5 θ ψ andψ → ψe iγ 5 θ . Such a four-Fermi interaction would be generated via the interaction between gluon and quark fields within QCD; see e.g., [8]. For a larger coupling constant G 0 than a certain critical coupling constant G c , the non-trivial vacuum, i.e., the chiral condensate ψ ψ , is generated.
Let us here consider the analysis of the DχSB in the NJL model at finite temperature and density with the functional renormalization group (FRG) method [9][10][11][12][13][14][15][16][17][18][19]; see also [20][21][22][23][24][25][26] as the reviews of the FRG. The advantage of the FRG is that one can systematically improve approximations for the FRG analysis. However, as will be explained in the next section, the renormalization group (RG) equation of the NJL model with the strong interaction has singularities originated from the divergence of four-Fermi coupling constant. Therefore, it is difficult to follow the solution of the RG equation after arising the singularities. The method, which have been employed in order to overcome the singularities, is the auxiliary field method based on the Hubbard-Stratonovich transformation [10]. 1 In this method, the four-Fermi vertex is replaced by exchanges of the bosonic auxiliary field φ via the Yukawa interactions yφψψ, namely, G 0 ∼ y 2 /m 2 . 2 Using this method one can not only avoid the singularities but also include the mesonic quantum effects. However, analyses of the system becomes complicated due to the additional bosonic degrees of freedom and ambiguities of approximations within the FRG method. Although the method of bare mass is also used to system [35], the operator expansion for the effective potential (9) does not converge in the small bare mass limit, and then it is difficult to obtain physical values in the chiral limit [36].
Moreover, the most serious problem in the method of bare mass is that one cannot address the chiral phase transition of first-order.
Recently, to overcome the singularities, a novel method called the weak solution method was suggested [37]. The solution of the RG equation with the singularities is mathematically defined as the "weak" solution. Hence, we call the RG equation applied the weak solution method "weak renormalization group equation". It has been shown in [37] that the weak solution method actually can address the chiral phase transitions of both second-and firstorder and evaluate the physical values such as the dynamical mass and the chiral condensate.
In this paper, we analyze the DχSB in the NJL model at finite temperature and density using the FRG with the weak solution method. The chiral phase diagram of the NJL model will be investigated. It should be emphasized that we will obtain the chiral phase diagram of the NJL model using the FRG without the auxiliary field method. Previously, only the auxiliary field method within the FRG has been able to capture the first-order phase 1 Even if the auxiliary field is introduced at the initial scale, the four-Fermi interactions are generated in the middle of scale by the box diagrams made of the Yukawa interaction between the quark and the auxiliary field. To avoid it, the dynamical hadronization method (or rebosonization) has been developed and used to investigate the vacuum structure of QCD [8,24,25,[27][28][29][30][31][32][33][34]. 2 Note that it may seem that the number of free parameters increases from one to two (from G 0 to y and m 2 ). However, one of them should be redundant. It has been shown in [8,10,32] that the low energy physics actually does not depend on the initial value of the Yukawa coupling constant.
transition. The main purpose of this paper is to show that one can obtain the chiral phase diagram with both second-and first-order phase transitions by the weak solution method without the help of the auxiliary field method. This paper is organized as follows: In the section II, we introduce the model and its effective model. The FRG method is briefly explained in order to apply it to the model.
In the section III, we introduce the weak solution method and give the weak RG equation.
The numerical results are exhibited in Sec. IV. We summarize and discuss our approach to the system and the results in Sec. V. In the appendix A, we show how to derive the RG equation. In the appendix B, the convexity and concavity of the beta function is discussed.

GROUP
Our purpose is to investigate the DχSB described by the NJL model in d = 4 and the chiral phase transition at finite temperature and density. In this section, we give the action of the NJL model and the RG equation for the effective potential. We see that the effective potential in the NJL model becomes singular within the RG evolution.

A. NJL model
The Lagrangian of the NJL model with one flavor is given in Eq. (1). Since in this paper we would like to investigate the fundamental behavior of the DχSB by using the weak RG, we consider the following simplified NJL model with one flavor in four dimension Euclidean space: where we assume that the quark fields have the degree of freedom of color, N c . Since the pseudo-scalar operator (ψiγ 5 ψ) 2 is not introduced, this action has no continuous chiral symmetry. Instead, it is invariant under the discrete type transformation, which forbids the operators with negative power ofψψ. Hence, the mass term of fermion mψψ does not appear in the action. This model describes that the chiral symmetry breaks down due to the strong four-Fermi interaction and the non-trivial vacuum, i.e., the non-zero chiral condensate ψ ψ = 0, is realized. In this vacuum, the fermion obtains the dynamical mass which is proportional to the chiral condensate. It should be noted that we do not really work for the discrete chiral symmetric model with only scalar four-fermi interactions.
Starting with the continuous chiral symmetric model with scalar and pseudo-scalar fourfermi interactions together, we define the renormalized fermionic interactions expressed in terms of scalar bilinearψψ and pseudo-scalar bilinearψiγ 5 ψ operators. The renormalized interactions have the same continuous chiral symmetry, since the RG beta functions are set to respect the symmetry.
After integrating the RG equation to the infrared end, we investigate the free energy (vacuum amplitude) as a function of scalar and pseudo-scalar source fields. The free energy is readily obtained making use of the translational invariance feature of our beta functions.
By the Legendre transformation in a usual way, we define the effective potential as a function of σ = ψ ψ and π = ψ iγ 5 ψ , which has the same chiral symmetric form, i.e., it is a function Therefore only a section of the effective potential where π = 0 is sufficient to get the radial profile of it. It is shown straightforwardly that this section does not depend on the coupling constant of elementary pseudo-scalar interactions in the bare Lagrangian. Then we set the vanishing pseudo-scalar coupling constant. This is our strategy of model set up, which gives the shortest way of evaluating the effective potential profile. 3 We should, however, take account of the pseudo-scalar operator to calculate, for example, critical exponents, which are sensitive to the difference of the discrete and continuous symmetries. 3 We comment on the finite number flavor case. In a general flavor case, the chiral symmetry is given as (1) A symmetry is broken by the chiral anomaly. This effective is described by the 't Hooft-Kobayashi-Maskawa term [38][39][40] within the effective model approach. Since this term with N f = 2 is given as the four-Fermi interaction, we would expect no crucial difference from the present setup. In contrast, in the N f = 3 (or 2 + 1) case, 't Hooft-Kobayashi-Maskawa term is given by the six-Fermi interaction. Although the anomaly effect on the chiral phase transition is investigated in several literatures, e.g. [41][42][43], it might depend on chiral effective models.

B. FRG equation
We introduce the FRG to investigate the DχSB in the NJL model. Now, we briefly explain its idea and derivation. Let us consider the bare action S[φ; Λ 0 ] and the path integral, This integration can be interpreted as the summation of the quantum fluctuations with the momentum 0 < |p| < Λ 0 . In the FRG, the path integral is evaluated by integrating out the quantum fluctuation with the higher momentum Λ < |p| < Λ 0 , namely, where we divided the field φ(p) into the higher momentum mode φ > (p) with Λ < |p| < Λ 0 and the lower momentum mode φ < (p) with 0 < |p| < Λ. The action S eff [φ < ; Λ] is called the Wilsonian effective action and is defined by the lower momentum mode. The FRG describes the change from the bare action S[φ; Λ 0 ] to the effective action S eff [φ < ; 0] through integrating the shell momentum Λ − δΛ < |p| < Λ, which is the so-called coarse graining. This step is formulated as a functional differential equation, and solving this equation is equivalent to evaluating the path integral (4). In this paper, we employ the Wegner-Houghton (WH) equation [44], where φ p := φ(p) and we introduced the dimensionless scale parameter, t := log(Λ 0 /Λ) .
It is rewritten as Λ = Λ 0 e −t , and then the increasing t corresponds to the decreasing the cutoff Λ. Since the two point function δφ −p is generally given as the supermatrix in field space, "str" denotes taking the trace for the supermatrix; see [45,46] for the treatments of the supermatrix and its manipulations. For details of derivation of the WH equation, see e.g. [45]. 4 The Wilsonian effective action is generally spanned by an infinite number of effective operators generated by the interactions in the bare action. Note that symmetries 4 The WH equation is also obtained from the sharp cutoff limit of the Plochinski equation [47,48] which the bare action has are not broken through the shell momentum integral, hence, the generated effective operators respect their symmetries.
Although the WH equation (6) itself is exact, we cannot solve it without some approximations. That is, we have to restrict the theory space of the Wilsonian effective action to its subspace. We give the following action as the approximated Wilsonian effective action for the NJL model (2): If the effective potential is expanded into the polynomials ofψψ, we have We call this potential the "fermionic" effective potential. As mentioned above, the odd powers ofψψ is not generated due to the (discrete) chiral symmetry. The potential at the initial scale Λ 0 is given by that is, G Λ 0 ≡ G 0 and other effective coupling constants vanish. Hereafter, we will use the notations σ :=ψψ. Note that the variable σ is not an auxiliary field, accordingly its derivatives with respect to ψ andψ are given as We will discuss the chiral phase transition by the effects of temperature and density in the next section. For this purpose, the density operator d 4 x µψγ 0 ψ should be introduced in the bare action (2) and the Wilsonian effective action (8). 5 Moreover, the momentum in the time direction p 0 is replaced with the Matsubara frequency ω n and its integration becomes the Matsubara summation. It is complicated to evaluate the shell momentum integral for four dimension momentum. To avoid it, we insert the cutoff Λ into the momentum in the spatial direction. Then the shell momentum integral becomes The RG equation for the effective potential is given by where β = 1/T , E := √ Λ 2 + M 2 and M := ∂ σ V . We show the explicit derivation of the RG equation in the appendix A. Note that the large-N approximation (N c → ∞) and the uniformψψ approximation were employed in order to obtain Eq. (13). The function M corresponds to the dynamical mass of quark. We here call it the mass function.
Here, we revisit the NJL model at vanishing temperature and density. The RG equation (13) at zero temperature and zero density becomes Using the expansion of the effective potential given in Eq. (9), we obtain the RG equations of each coupling constant. In particular, the RG equation of the four-Fermi coupling constant is given by If we use the dimensionless rescaled coupling constantG := GΛ 2 /2π 2 , it becomes Obviously, this equation has the fixed pointG c = 1 and yields the blow-up solution, for the initial valueG 0 >G c . We see the critical scale at which the denominator becomes zero, The divergence ofG itself is physically relevant since it corresponds to the chiral fluctuatioñ G ∼ (ψψ) 2 . Thus, the divergence is the signal of the second order phase transition.
However, the divergence disturbs solving the RG equation up to the infrared (IR) scale Λ → 0. Since the four-Fermi coupling constant is lowest order of the effective potential (9), the effective potential becomes non-differentiable at its origin V (σ = 0) due to the singularity of the four-Fermi coupling constant. In other words, the derivative ∂ σ V cannot be defined after the critical scale t c .
We see in the next section that the weak solution method [37] allows us to analyze the DχSB without introducing the auxiliary field and to investigate the chiral first-order phase transition.

III. WEAK SOLUTION METHOD
In this section, the weak solution of the RG equation (13) is defined. To construct it numerically, we introduce the method of characteristics which are given as the system of ordinary differential equations made from the partial differential equation. Moreover, to obtain a single-valued solution from a multi-valued solution, we use the Rankine-Hugoniot Finally, we comment on the RG equation for the four-Fermi coupling constant from the viewpoint of the weak solution method.

A. Definition of weak solution
As discussed in the previous section, the RG equation for the effective potential becomes singular when the DχSB occurs. In [37], it was shown that the weak solution method is valid for overcoming the singularity of this system. Here, we briefly review basic notions of the weak solution method and show the weak RG equation for Eq. (13).
The RG equation (13) is given as the partial differential equation, where F (M, σ; t) is the beta function with a negative sign. Differentiating the both sides of RG equation with respect to σ, we obtain Let us now call this equation the "strong" RG equation in order to compare with the "weak" RG equation introduced below.
Here, we introduce a test function ϕ(σ; t) which is smooth and converges to zero for both infinite t and σ: ϕ(±∞; t) = ϕ(σ; ∞) = 0. Integrating Eq. (20) multiplied by this function, By utilizing the integration by parts, it becomes We call the RG equation (22) "weak RG equation", and its solution is called "weak solution".
The derivatives with respect to t and σ move from the mass function and the beta function to the test function. Hence, M (σ; t) given as a solution of the weak RG equation (22) can have an arbitrary number of singularities at any point on the σ-t plane.

B. Method of characteristics
The RG equation of the mass function (20) is given as a partial differential equation.
Here, we introduce the method of characteristics which reduces from the partial differential equation to the system of ordinary differential equations. The method of characteristics allows us to numerically solve the RG equation using e.g., the Runge-Kutta method and the finite-difference methods.
Let us consider a curved surface M (σ; t) in the σ-t-z space and rewrite the RG equation (20) as where we used since the total derivative of M is given by dM = ∂ σ M dσ + ∂ t M dt. Thus, the vector (∂ σ M ∂ t M − 1) T is orthogonal to (dσ dt dM ), and we find that the vector (∂ M F 1 0) is also the tangent vector for a point on the curved surface M (σ; t). Taking (dσ dt dM ) as the vector being proportional to (∂ M F 1 0), we have where we introduced an infinitesimal parameter ds as a proportionality constant. The solutions (σ(s) , t(s) , M (s)) of the system of ordinary differential equations with each initial conditions are called characteristics. Since the solution of Eq. (26) is simply t = s, we can replace from s to t in Eqs. (25) and (27).
The fermionic effective potential V can be evaluated within the method of characteristics.
To summarize, we have the following coupled equations for the present system: Actually, these equations are independent of each other, therefore, their solutions are given by where σ 0 := σ(0). Note here that F (M, σ; t) given in the right-hand side of Eq. (13) (with negative sign) does not explicitly depend on σ; F (M, σ; t) ≡ F (M ; t). We also note that this simplification happens in the NJL model with a certain approximation. Although for more  (32). Note that the initial mass function is general cases such as QCD, we have to solve the system of ordinary differential equations for a given system, solving it is easier than the case of the partial differential equation.
A schematic figure for the evolution of the mass function described by Eqs. (30)-(32) is shown in Fig. 1. Note that for a fixed M 0 > 0, σ described by Eq. (31)   where σ * (t) is a discontinuity point, and σ * + and σ * − denote closing to the point σ * from the right-hand side (+) and the left-hand side (−), respectively (see Fig. 2). This condition describes the evolution of σ * (t) and implies geometrically the differential-area invariance, namely,  22) with the RH condition is given by the solid line with a discontinuity at σ * . The discontinuity determined by the RH condition is called "shock". 6 As an example, the evolution of the mass function at vanishing temperature and density is shown in Fig. 3, where the dimensionless mass functionM = M/Λ 0 and the dimensionless fieldσ = σ/Λ 0 are defined; see Eqs. (49) and (50) in the section IV. At the critical scale t c , the slope of the mass function atσ = 0 becomes infinity, which corresponds to the divergence of the four-fermi coupling constant. After t c , the solution of the mass function becomes the multi-valued function. Applying the RH condition, the solution is uniquely determined; see 6 The terminology comes from the "shock wave" in gas dynamics. the right-hand panel of Fig. 3. Next, let us consider the case at zero temperature and finite density. The evolution of the mass function in this case is shown in Fig. 4. We see that two singularities arise at the non-zero field values,σ = 0. They move towards the originσ = 0 and eventually merge with lowering the scale t. For such a behavior of the mass function, we observe the first-order phase transition as will be seen in the next section.
The mass function given as the weak solution satisfies since it reflects the chiral transformation σ → −σ. Then, the dynamical mass is given by To summarize so far, the multi-valued mass function as a solution of the characteristics Finally, we comment on the Legendre effective potential and its convexity and concavity.
Let us start with defining the generating functional for the connected Green function, Here, the initial condition for the fermionic potential is given by where the last term is the source term. From the Legendre transformation of Eq. (37), we can define the Legendre effective potential, This potential is a function of the chiral condensate given by Note that the source j is regarded as a function of φ. It is known that the Legendre effective potential V L becomes the convex form [49][50][51][52]. The weak solution method automatically gives us the convex Legendre effective potential as the convex envelope of the nonconvex function obtained from the solution of characteristic equations (30)- (32). See discussions in [37] on its proof.

D. RG flow of four-Fermi coupling constant within weak solution method
We discuss the RG flow of the four-Fermi coupling constant in view of the weak solution.
Since the mass function M (σ; t) is given by Next, we discuss the relation between the inverse four-Fermi coupling constant g := 1/G Λ and the Legendre effective potential V L . We start with the well-known relation, where W (j; t) is given in Eq. (37). The second derivative of W with respect to j corresponds to the susceptibility: Together with Eq. (42), we obtain When the second-order phase transition occurs, the susceptibility diverges: χ → ∞. That is, the left-hand side ∂ 2 V L (φ;t) ∂φ 2 vanishes at the critical scale t = t c . More explicitly, we see At the origin of V L , the mass-squared m 2 , i.e. the curvature of V L , corresponds to the inverse susceptibility 1/χ. The divergence of χ means the vanishing mass-squared m 2 = 0.
As mentioned in the section II, the four-Fermi coupling constant corresponds to the chiral susceptibility: Then, the divergence of G Λ is the signal of the second-order phase transition, and we see that the inverse four-Fermi coupling constant corresponds to the curvature at the origin of the Legendre effective potential: The RG flow of the dimensionless rescaled inverse four-Fermi coupling constantκ =G −1 = 2π 2 /(G Λ Λ 2 ) is given byκ whereκ 0 is the initial value ofκ at the initial scale t = 0. Besides, the first-order phase transition is not captured by only looking for the behavior of the four-Fermi coupling constant.
In the next section, we analyze the weak renormalization group at finite temperature and density and show the chiral phase structures of the NJL model.

IV. NUMERICAL ANALYSIS
In this section, we show the chiral phase transitions within the NJL model described by Eq. (13). To this end, we first give the RG equation written by the dimensionless variables.

B. Phase transitions
Here setting g 0 = 2 as a benchmark point, we show the dependence of the mass function on temperature and density. First, let us consider the zero density case µ = 0. We show the dependence of the mass function in the IR limit t → ∞ (Λ → 0) on temperature in Fig. 5.
With increasingT , the value of the mass function atσ → 0 smoothly becomes smaller. This behavior corresponds to the second-order phase transition. Note that the dynamical mass is given in Eq. (36). The mass function atσ → 0 vanishes atT c 0.589 which is the critical temperature. Next, we consider the finite chemical potential and vanishing temperature case.
In this case, the evolution of the mass function in the IR limit is shown in Fig. 6. We see that the slope of the origin of the mass function changes from negative to positive atμ 0.707; see the left-hand side of Fig. 6. This means that the four-Fermi coupling constant becomes positive via the divergence, and thus the Legendre effective potential has a local minimum at its origin. Here we call this behavior the "second-order" although it does not always correspond to the phase transition. Note that there is still a non-trivial global minimum at µ 0.707. Aboveμ c 0.861, the discontinuities of the mass function in the IR limit are located at non-vanishing σ. Then, when taking the limit σ → 0, the physical dynamical mass suddenly vanishes with increasingμ, i.e., M (σ → 0; t → ∞; This phase transition is of first-order. We see that the weak solution method can capture the phase transitions with both second-and firstorder. In the next subsection, we show the phase diagram on the g 0 -μ-T plane.

C. Phase diagram
In Fig. 7, the chiral phase diagram on the g 0 -μ-T plane are shown. The red line stands for the second-order phase boundary. As discussed in the section II, this boundary can be obtained by only the behavior of the four-Fermi coupling constant. The blue line corresponds to the first-order phase boundary. The biggest and smallest phase boundaries correspond to the cases of g 0 = 2 and g 0 = 1.0101, respectively. We show the phase diagram on thẽ µ-T plane in Fig. 8, which is obtained by the projection from the three dimensional phase diagram of Fig. 7. A schematic figure of a phase boundary with a fixed four-Fermi coupling constant is shown in Fig. 9. Note that the region, which is surrounded by the second-and first-order phase boundaries towards the high density side, is the broken phase. As discussed in the previous subsection, although the four-Fermi coupling constant diverges at the secondorder phase boundary, that is, the curvature of the origin of the effective potential becomes positive, the global minimum is still located at the non-vanishing chiral condensate. In  Here, we note the case where the finite bare mass is taken into account. In this case, the horizontal axis, we obtain the value of the four-Fermi coupling constant at which the critical end-point merges with the second-order phase boundary at T = 0. we could observe the crossover for smaller chemical potential rather than the second-order phase transition. This behavior within the weak renormalization group is shown in [37].

D. Comparisons with other models and methods
Let us compare our results with other models and methods. We first note that as shown in Refs. [9,53], the FRG analysis in the large-N limit is equivalent to the SDE with the ladder approximation. The ladder approximated SDE is also equivalent to the MFA. Therefore, these methods yield quantitatively the same results.
We compare the present work with the other ones in the two-flavor and three-color case (N f = 2, N c = 3). 8 For the mean-field analysis of the NJL model, we employ the following Lagrangian, Note that the critical four-Fermi coupling is g 0c = 1/N f N c = 1/6. In Ref. [54], the UV cutoff, the four-Fermi coupling constant and the bare quark mass are set to Λ 0 = 587.9, g 0 = 0.247   In the other works with the analytic methods [55][56][57][58][59], the free parameters are set such that one obtains the pion decay constant f π = 93 MeV and the pion mass m π = 138 MeV in the IR regime. 9 Refs. [60,61] performed the lattice QCD simulation with imaginary chemical potential µ I = −iµ. 10 By using analytic continuation to real (quark) chemical potential µ, the dependence of the critical temperature on small chemical potential is investigated. In Ref. [62], the Taylor-series expansion method is employed.
In Table. I, the (pseudo-)critical temperature T c (µ = 0) and the CEP obtained from various models and methods are collected. We also show the value of the curvature of the chiral phase boundary κ µ which is defined as Note that κ µ does not depend on the cutoff scale Λ 0 . We evaluate the curvature in a range 0 ≤ µ/(πT c (0)) 0.1 for the present analysis and the mean-field analysis of the NJL model (III).
We see from Table. I that there are no drastic differences in the critical temperature at vanishing chemical potential and the CEP. In contrast, the curvatures obtained from our work and the mean-field analysis of the NJL model are somewhat larger than that of the others. This could be because the mesonic fluctuations are not taken into account in the present FRG and the MFA computations of the NJL model. Fig. 11 exhibits the dependence of the curvature κ µ on the dimensionless four-Fermi coupling g 0 . We see that κ µ tends to be larger with increasing g 0 . Note that the value of the critical four-Fermi coupling constant is g 0c = 1, and as we have seen in the previous subsection, the CEP vanishes below g 0 1.2.

V. SUMMARY AND DISCUSSION
In this paper, we have studied the chiral phase structure of the NJL model at finite temperature and density using the FRG. We have seen that when the DχSB takes place, the solution of the RG equation has to become singular and cannot be evaluated after a certain critical scale. To overcome this situation, we have introduced the weak solution method to the RG equation. It has been shown that the weak solution method can appropriately 9 The work in [55] obtains f π = 87 MeV and M ∼ 300 MeV in the chiral limit. When the effects of the finite and realistic pion mass are taken into account, the decay constant becomes f π = 93 MeV. 10 Ref. [60] uses two-flavor staggered quarks, whereas two-flavor Wilson quarks is used in Ref. [61]. In particular, it was difficult to obtain the chiral phase diagram with the first-order phase transition. We have shown that the weak solution method allows us to evaluate the chiral dynamics within the pure fermionic system.
Using the FRG with the weak solution method, the DχSB in QCD at zero temperature and zero density was studied in [36], where effects beyond the ladder approximation (nonladder effects) on the dynamical mass and the chiral condensate are involved. 11 It was shown that the gauge dependence of the physical values is suppressed thanks to including the non-ladder effects. As a next step, we should investigate the DχSB in QCD at finite temperature and density. The weak solution method also can be applied to systems with the color superconductivity and ones in an external magnetic field. 12 In future works, we will apply the FRG with the weak solution method to such systems.
Here, we discuss several issues of the weak solution method. First, we comment on the cutoff scheme and the convexity and concavity of the beta function. Through this paper, we have analyzed the RG equation obtained by using the WH equation which is formulated 11 The crossed ladder diagrams and the anomalous dimension of the quark field are included. These effects cannot be taken into account by the SDE with the ladder approximation and the MFA. 12 See e.g., [63][64][65] for the analyses of the system in an external magnetic field using the FRG. with the sharp cutoff regularization. When the smooth cutoff scheme, e.g., the optimized cutoff regularization [66], is employed, we obtain the different form of the beta function for the effective potential; see Eq. (A19). As shown in the appendix B, for the sharp cutoff regularization, the beta function for the fermionic effective potential (13) is concave, i.e., ∂ 2 F/∂M 2 < 0 for an arbitrary RG scale, temperature and density. Therefore, the motion of σ(t) described by Eq. (31) behaves as shown in Fig. 1. In contrast, for the case where the optimized cutoff regularization is used, the concavity of the beta function is not guaranteed.
Hence, there are solutions in which the motion of σ(t) comes back to the right-hand side as shown in Fig. 12. 13 In this case, however, we need some additional information about the solution of the mass function removed once by the RH condition, which corresponds to the dotted line in Fig. 12, in order to define the mass function as a function of σ. At present, we do not know how to mathematically define the weak solution for such a "come-back" solution. 14 Second, let us consider improvements of the approximations. In this work, we have applied the large-N approximation. For beyond this approximation, the beta function 13 Such a behavior in terms of the four-Fermi coupling constant was reported in [17,18]. 14 It should be stressed here that the "come-back" solution (physically correct equal area law solution) is a weak solution because it satisfies the RH condition, but it is not a unique solution since there always appear "rarefaction" solutions in this type of situation. We do not even know if we can define some additional condition to assure uniqueness of weak solution in this situation, nor if such condition, if any, may pick up our equal area law solution.
for the effective potential generally involves its second-order derivative terms with respect to the field σ, namely F (M, M , σ; t) where the prime on the mass function denotes the σ-derivative and M = ∂ σ V ; see Eqs. (A4) and (A5). In this case, we cannot naively transfer the t-and σ-derivatives from the mass function to the test function by performing the integration by parts as shown in Eq. (22) and then cannot define the weak solution. The notion of the weak solution that we have introduced in this work corresponds to the so-called "entropy solution" based on distributions. On the other hand, it is known that there is also "viscosity solution" as a different definition of the weak solution. 15 The viscosity solution can be defined for even some second-order nonlinear partial differential equations to which the definition of Eq. (22) cannot be applied. The viscosity solution could be a breakthrough for improving the approximation. Therefore, it is important to study applications of the weak solution method in order to analyze critical phenomena in elementary particle physics and condensed matter physics.

Sharp cutoff scheme
For the given effective action (8), the WH equation reads where φ T p = ψ T (−p) ,ψ(p) . To derive the RG equation, we introduce the mean-field defined as ψ(p) = (2π) 4 δ (4) (p) ψ,ψ(p) = (2π) 4 δ (4) (p)ψ. (A2) We apply them after evaluating the functional derivatives of S eff with respect to φ(p). The inverse two point function is given by where we can calculate Here we employ the following approximation: which corresponds to the so-called large-N leading approximation (N c → ∞).