Weak solution of the non-perturbative renormalization group equation to describe the dynamical chiral symmetry breaking

We analyze the dynamical chiral symmetry breaking (D$\chi$SB) in the Nambu-Jona-Lasinio (NJL) model by using the non-perturbative renormalization group (NPRG) equation. The equation takes a form of two-dimensional partial differential equation for the multi-fermion effective interactions $V(x,t)$ where $x$ is $\bar\psi\psi$ operator and $t$ is the logarithm of the renormalization scale. The D$\chi$SB occurs due to the quantum corrections, which means it emerges at some finite $t_{\rm c}$ in the mid of integrating the equation with respect to $t$. At $t_{\rm c}$ some singularities suddenly appear in $V$ which is compulsory in the spontaneous symmetry breakdown. Therefore there is no solution of the equation beyond $t_{\rm c}$. We newly introduce the notion of weak solution to get the global solution including the infrared limit $t\rightarrow \infty$ and investigate its properties. The obtained weak solution is global and unique, and it perfectly describes the physically correct vacuum even in case of the first order phase transition appearing in finite density medium. The key logic of deduction is that the weak solution we defined automatically convexifies the effective potential when treating the singularities.

Now we state the central subject of this paper. In the NPRG method, we integrate the path integral of the theory from the micro scale to the macro scale, slice by slice, and obtain the differential equation of the effective action with respect to the logarithmic renormalization scale t, which is defined by where 0 is the bare cutoff scale and is the renormalization scale.
In the simplest approximation, the effective action is expressed in terms of the field polynomials without derivatives in addition to the normalized kinetic terms. Then the coupling constants C i , the coefficients of these operator polynomials, are functions of t. The change of these coupling constants are given by functions of the coupling constants themselves, where the right-hand side is called the β-function. This differential equation is the renormalization group equation (RGE).
The β-function is evaluated by one-loop diagrams only and the propagators constituting the loop are limited to having sliced momentum. Therefore its evaluation is so simple that only angular integration need be done to leave the total solid angle factor, which is called shell mode integration.
Here we consider the four-Fermi coupling constant G. The four-Fermi interactions play an essential role in the total framework of the Dχ SB analysis. The β-function for the dimensionless four-Fermi coupling constant rescaled by the renormalization scale ,G(≡ 2 G), is given by the following Feynman diagrams, where the first term in the right-hand side represents the canonical scaling due to the operator dimension, the solid lines are quarks, the wavy lines are gluons, and the above loops represent shell mode integration. It is evaluated as follows [7][8][9]: where α s is the gauge coupling constant squared as usual, C 2 is the second Casimir invariant of the quark representation and ξ is the gauge fixing parameter in the standard R ξ gauge. Due to the chiral symmetry, the RGE for the four-Fermi operator is closed by itself; no contribution from higherdimensional operators is allowed. Consider first the Nambu-Jona-Lasinio (NJL) model as an effective infrared model of QCD. Then the β-function is simplified to give and it is depicted in Fig. 1, where the arrows show the direction of renormalization group flows towards the infrared. Note that there are two fixed points, zeros of the β-function, at 0 andG c = 4π 2 , and they are the infrared (stable) fixed point and ultraviolet (unstable) fixed point, respectively. It is readily seen that there are two phases, the strong phase and the weak phase, and the critical coupling constant dividing the phase isG c . If the initial coupling constant is larger thanG c , flows go to the positive infinite, otherwise they approach the origin. We are sure that this is the quickest   argument ever to derive the critical coupling constant in the NJL model and the result is equal to the mean field approximation.
The weak phase has no problem here. The flows approaching the origin do not mean that the four-Fermi interactions vanish at the infrared. Since the variable here is the dimensionless coupling constant, the four-Fermi interactions look vanishing due to the multiplication factor of exp(−2t). As for the dimensional variable, it just stops running to give some constant size infrared interactions determined by the initial coupling constant.
On the other hand, the strong phase has a serious problem. Equation (5) is easily integrated to give the solutionG whereG 0 is the initial coupling constant (Fig. 2). This solution is a blow-up solution and it diverges at a finite t c defined by Because of this blow-up behavior, we cannot obtain the solution after t c ; that is, there is no global solution in Eq. (5) forG 0 >G c . This is not a fake, nor an insufficiency of our treatment. This must happen necessarily [10,11]. The reason is the following. The NJL model is expected to have the Dχ SB with strong coupling constantG 0 >G c . This is the quantum effect, that is, it is caused by the quantum loop corrections.   The schematic view of the chiral condensates or the dynamical mass as a function of the renormalization scale t is drawn in Fig. 3, where at t c the Dχ SB occurs. This figure resembles the temperature dependence of the spontaneous magnetization if the low and high t are reversed. Then consider what happens if t approaches t c from the origin. The correlation length of spins (ψψ operators) must diverge and also the magnetization susceptibility diverges. The susceptibility is proportional to the total spin fluctuation (ψψ) 2 . This is nothing but the four-Fermi interactions. Therefore the diverging behavior of the four-Fermi interactions itself is a normal and naive phenomenon that just tells us we are approaching the second-order phase transition point.
The same story holds for the QCD case. The β-function of the four-Fermi coupling constant now depends on the gauge coupling constant, as shown in Fig. 4 (ξ = 0, C 2 = 1 case). Switching on the gauge interactions, the gauge coupling constant becomes large towards the infrared, the parabola function moves up, and finally after α s > α c s = π/3, it has no zero at all. Accordingly, the stable and unstable fixed points approach each other, and they are merged to pair-annihilate. After this annihilation, the total space belongs to the strong phase and the four-Fermi coupling constant moves to positive infinite. This is our NPRG view of how QCD breaks the chiral symmetry in the infrared, and it also explains the peculiar phase diagram of the so-called strong QED (or the gauged NJL model) as is depicted in Fig. 5 from the RGE point of view [7][8][9]  Now the problem is how to calculate the infrared quantities, such as the chiral condensates and the dynamical mass. They all reside beyond the blow-up point. We stress here that this blow-up is irrelevant to the approximation scheme. Any theory, approximated or not, as far as it exhibits the Dχ SB, must encounter the spontaneous breaking during renormalization, since it is produced by the sufficient quantum corrections. The Dχ SB must be expressed by the singular effective interactions of fermions. Therefore, any model providing the Dχ SB must meet singularities during renormalization. This is inevitable.
Many NPRG analyses of the Dχ SB have been performed by introducing the bosonization of the multi-Fermi interactions [7,9,12,13], which is called the auxiliary field method or the Hubbard-Stratonovich transformation. This type of approach may avoid the blow-up singularity and has indeed worked well. However, it must contain additional difficulties or ambiguities in how to treat and approximate the system of bosonic degrees of freedom.
We take the other way around. We do not introduce the bosonic degrees of freedom and instead we extend the notion of the solution of the NPRG equation. First we go back to the original NPRG equation without expansion in polynomials of field operators. Then it is a partial differential equation (PDE) of at least two variables, the operator and the renormalization scale t. The blow-up behavior is reexamined in the next section and it is made clear that even the solution of the PDE must encounter the corresponding singularity at t = t c . Thus the origin of the singularity is irrelevant to the Taylor expansion itself of PDE. Even if the PDE type of RGE is adopted, we cannot have any global solution.
Then we introduce the notion of a weak solution of a PDE [14]. This allows a solution with singularities and we may define the global solution including the infrared limit. This paper is organized as follows. In Sect. 2, taking the NJL model as an example, we briefly review the Wegner-Houghton equation, an approximation of the NPRG equation. In Sect. 3, the notion of a weak solution is introduced, and a practical method to construct the weak solution is explained in Sect. 4. In Sect. 5, we solve the NJL model to explicitly construct the weak solution. In Sect. 6, we discuss the effect of the bare mass of quarks, which is necessary to define the chiral order parameters, and we construct the Legendre effective potential. In Sect. 7, the method of weak solution is applied to the first-order phase transition case in the finite chemical potential NJL model. The convexity of the effective potential given by the weak solution is fully discussed. Finally, we summarize the paper in Sect. 8. To conclude, we will prove that the weak solution of the NPRG equation is perfectly successful, giving the global solution to calculate the infrared physical quantities. Most impressive is the case of the first-order phase transition. There the weak solution correctly picks up the lowest free energy vacuum among the multiple locally stable vacua, through the procedure that the effective potential is automatically convexified by the weak solution.

Non-perturbative renormalization group
In this section we introduce the NPRG by using the Wegner-Houghton (WH) equation. For simplicity, we take the Nambu-Jona-Lasinio (NJL) model with a simplified discrete chiral symmetry, which is regarded as a low energy effective model of QCD explaining the Dχ SB.
The Lagrangian density is given by where ψ andψ are the quark field and the antiquark field, respectively. There is no continuous chiral symmetry. However, the Lagrangian is invariant under the following discrete chiral (γ 5 ) transformation, This discrete chiral symmetry forbids the mass term mψψ and the chiral condensate ψ ψ as well as the usual continuous chiral symmetry. It is first clarified by Nambu and Jona-Lasinio [1,2] that for strong enough four-Fermi coupling constant G 0 , it exhibits Dχ SB. Using the mean field approximation, the critical coupling constant G c is found to be 4π 2 / 2 0 , where 0 is the ultraviolet cutoff scale. Note that the mean field approximation is equivalent to the self-consistency equation limited to the large-N leading diagrams, where N is the number of quark flavors.
In the NPRG approach, a central object is the Wilsonian effective action S eff [φ; ] defined by integrating out the microscopic degrees of freedom φ H with momenta higher than the renormalization scale , where S 0 [φ; 0 ] is the bare action with the ultraviolet cutoff scale 0 and φ generically denotes all relevant fields. The renormalization scale is parametrized by dimensionless variable t as defined in Eq. (1), The t-dependence of the effective action S eff [φ; t] is given by the NPRG equation, which is the following functional differential equation, which is also called the Wegner-Houghton (WH) equation [15] (see Ref. [16] for the detailed form of β WH ). The right-hand side called the β-function is actually evaluated by the one-loop diagram exactly. No higher loops contribute. This functional differential equation should be solved as the initial value problem, where the initial condition refers to the bare action, Solving the equation towards the infrared (t → ∞) and we get the macro effective action from which physical quantities such as the chiral condensate and the dynamical quark mass are obtained. The WH equation itself is exact, but it cannot be solved exactly, since it has infinite degrees of freedom.
Here we apply the WH equation to the NJL model. As an approximation, we restrict the full interaction space of the effective action S eff [ψ,ψ; t] into the subspace most relevant to Dχ SB. We adopt the so-called local potential approximation where any quantum corrections to the derivative interactions are ignored. Furthermore, we set the local potential as a function of only the scalar fermion bilinear field, x =ψψ. Then our effective action takes the form The potential term V (x; t), where x denotesψψ (integrated over space time), is called the fermion potential here, whose initial condition is set to be according to the NJL Lagrangian (8).
When evaluating the β-function we adopt an additional approximation of ignoring the large-N non-leading part in β WH . Then, the NPRG equation for the fermion potential in the large-N leading is given by the following partial differential equation: Hereafter we sometimes use a quick notation to save space: ∂ t = ∂/∂t, etc. The renormalization scale , the momentum cutoff, is defined by referring to the length of four-dimensional Euclidean momentum p μ , that is, 4 μ=1 p 2 μ ≤ 2 . This is called the four-dimensional cutoff. The right-hand side of the equation is nothing but the simple trace-log formula of the Gaussian integral of the shell modes and it corresponds to the one-loop corrections, though it is exact. Note that the approximate NPRG equation above is known to be equivalent to the mean field calculation [7][8][9]. Now we introduce the mass function M(x; t), the first derivative of the fermion potential, The value of the mass function at the origin is the coefficient of mass operatorψψ in the effective action, and this is why we call it the mass function. Note that it is still a function of x, that is, a function of an operator, or more rigorously, a function of the bilinear Grassmannian variable x.
The chiral symmetry transformation defined in Eq. (9) is represented by the reflection The NPRG equation (16) preserves this reflection symmetry and the initial condition also respects it. Accordingly, the fermion potential at any time t satisfies the reflection symmetry and therefore the mass function satisfies, and it is an odd function of x. Then the mass function must vanish at the origin, assuming it is well-defined at least. This means there is no spontaneous mass generation and the phase is chiral symmetric.
Then what actually happens in the case when Dχ SB and the spontaneous mass generation occurs? We have to abandon the well-definedness of the mass function at the origin, that is, the fermion potential V must lose its differentiability at the origin. Then we can have the following situation, that is, the mass function has a jump at the origin, still keeping its odd function property. This is the spontaneous mass generation. To explain clearly that this does indeed give us the dynamically generated mass, we need the notion of the effective potential, which will be introduced in Sect. 6.
Now we understand what we should deduce. At some finite t c , the fermion potential becomes nonanalytic at the origin, and it generates a finite gap in the mass function. The gap continues to increase and finally gives the dynamical quark mass at the infrared. Just before t c , the second derivative of V must become large and increases towards infinity at t c . The second derivative of V at the origin gives the four-fermion interactions at the scale t, and this divergent behavior is nothing but the blowup solution we find in the expanded NPRG equation in Eq. (5) in the case of G 0 > G c [10,11].
See the plot (d) in Fig. 12, which is exactly what we would like to realize as a solution for the mass function. Before t c the total function is a smooth analytic function, whereas after t c there is a finite jump at the origin which shows the spontaneous mass generation. Such a solution, however, is not allowed as a solution of the PDE in Eq. (16), since the solution of the PDE must satisfy the equation at every point. Therefore, the differentiability at any point (x, t) is mandatory, while our expected solution must break it, though only at the origin.
To cure this situation, we have to modify the original equation so that differentiability of function V might not be necessary. This is the notion of a weak solution, which is defined by the solution of the weak equation, and it will be the subject of the next section.

Weak solution and the Rankine-Hugoniot condition
In this section, we first derive the partial differential equation for the mass function and define the weak solution of it [17,18]. Then we construct the weak solution explicitly using the standard method of characteristics [14] in the next section.
We start with the following partial differential equation for V (x; t), This is an extended system of the NJL-model NPRG equation (16), where the β-function F contains explicit dependence on x. The reason why we include this extension is that the gauge theory analysis requires this additional x dependence, and also the total system becomes more symmetric and better interpreted in terms of classical mechanics. Differentiating the PDE (23) with respect to x, we obtain the following PDE for the mass function, Here it should be noted that we have used somewhat subtle notations that two forms of the partial derivative of F with respect to x have different meanings for the treatment of variable M. This equation belongs to a class of PDE called the conservation law type since it has a form of conservation law in two-dimensional space-time. This class contains the famous Burgers' equation for a non-linear wave without viscosity [19]. We integrate the above equation convoluted with a smooth and bounded test function ϕ(x; t), and then by the partial integration we have  (24) at every point is in fact weak solution. The inverse does not hold. In this sense, the weak solution is a wider notion than the strong solution.
Then even in the case when the strong solution does not have a global solution, the weak solution can be global and then we can calculate infrared physical quantities using the weak solution. In general, the weak solution may not be unique. It also may depend on the initial condition. We will prove that the weak solution of our NPRG equation is in fact unique and global, and we can successfully obtain the infrared physics without any ambiguity.
Here we state what the weak solution is without explicit proof (see, e.g. Ref. [14] for a detailed argument). We suppose a weak solution has some jump singularities at some finite number of points and otherwise it is smooth, that is, it is piecewise differentiable. Then almost everywhere other than the singularity, the weak solution must satisfy the strong equation (24). This is the first property of the weak solution and thus it is not so different from the strong solution.
Second, what is required at the singularity? As is shown in Fig. 6, we consider the neighborhood of a singularity where there is a jump discontinuity on a curve x = S(t), and values of M at the left-hand side and right-hand side of the singularity are denoted by M L and M R respectively. Remember that the original equation (24) is regarded as the conservation law where the charge density is M(x, t) and the current is F(M(x; t), t). Consider a small boxed region d S × dt in Fig. 6 and evaluate the total charge in the section of length d S. The conservation law requires that the change of total charge must be balanced with the difference between total inflow and outflow. We have, This is nothing but the integrated form of the conservation law.
The above condition (26) is called the Rankine-Hugoniot condition [20][21][22], and all discontinuities must satisfy it. This is the second property which the weak solution must satisfy, and no more. This condition actually determines the equation of motion of the singularity position by the first-order differential equation. It needs one initial condition to determine the motion of the discontinuity point, and it is actually given by the spontaneous generation of the singularity, as is seen later.  The Rankine-Hugoniot condition for the mass function M is converted to the condition for the fermion potential V . In Fig. 6, we denote the fermion potential at each side of the singularity by V (L) and V (R) respectively. Then we evaluate its difference between the bottom-left point and top-right point of the dt × d S box as follows: We have the following relation given by the RH condition (26), Therefore the fermion potential develops equally between the left and right sides of the singularity. Taking account of the fact that there is no singularity in the beginning in the case of our physical initial condition, that is, any singularity starts later during t-development, we conclude that Thus the fermion potential function V (x; t) is continuous everywhere as far as it is the weak solution of the PDE. This form of the RH condition helps us to obtain the weak solution in the next section.

Method of characteristics
In this section, we construct the weak solution using the standard textbook method of characteristics. The characteristic is a curve on a two-dimensional world (x, t) and is denoted by x =x(t). Consider the mass function on this curve,M(t) = M(x; t), and we calculate its derivative: We take the characteristics curve to satisfy the differential equation then the mass function on it satisfies, This coupled pair of ordinary differential equations (32) and (33) are equivalent to the original twodimensional partial differential equation (24). The initial conditions for these equations are set bȳ where x 0 is a parameter discriminating characteristics and M 0 is the initial value of M at x = x 0 . When explicitly indicating the initial value of the characteristic, we usex(t; x 0 ). The characteristics will give us the strong solution of the PDE. Taking a characteristic starting at x 0 , we solve the coupled ODE, and then we get the characteristicx(t; x 0 ) as shown in Fig. 7 (a), and the mass function on it,M(x(t); t).
Here we change our viewpoint. The coupled ODE in equations (32) and (33) is regarded as the canonical equation of motion of a kinematical system with coordinatex, momentumM, and timedependent Hamiltonian F(M,x; t). We plot the two variablesx(t) andM(t) in Fig. 7 (b) and this is nothing but the phase space orbit of the Hamiltonian system.
The fermion potential V also has its counterpart in mechanics. Consider the fermion potential on the characteristics,V (t) = V (x(t); t), and we calculate its derivative, The right-hand side is the Lagrangian of this mechanical system. Then the fermion potential is obtained by integration, Thus this is the action of the system as a function of the final time t and the coordinate variable at the final timex(t), plus the initial value. Then the original PDE in Eq. (23) is what is called the Hamilton-Jacobi equation in this kinematical system. The starting point of this particle is controlled by the parameter x 0 , and its initial momentum is the initial mass function M(x 0 ; t = 0). Then we consider all orbits at once as shown in Fig. 8, which is now the motion of a string on the two-dimensional world x × M. The string never crosses itself, since the equation of particle motion uniquely determines the phase space orbit, thus no crossing of orbits is allowed. The shape of the string at t gives the mass function M(x; t), the strong solution of the original PDE.
There is nothing singular, and no explosive blowup so far. The string motion gives the strong solution, locally. The point here is that it may not give the global solution as it is, when the string is folded and cannot define a function of x. The string itself is parametrized by x 0 , and M(x(t; x 0 ); t) is a unique function of x 0 . However, when the characteristics are crossing each other, we have multiple 11   x 0 for a single value ofx(t). This situation is intuitively stated as that we have a "multi-valued We note that the string motion here with the proper initial condition has no singularity at all up to the infinite time t = ∞. The string neither breaks up nor shrinks. The string M(x(t; x 0 ); t) as a function of x 0 is always continuous. This fact will ensure that the weak solution we define now is unique and satisfies the so-called entropy condition. Normally the multi-folded structure is born as three-folded leaves, and in the example in Sect. 7, two three-folded structures meet together to become five-folded.
Let us consider a typical plot for the NJL model (super-critical situation, G 0 = 1.005G c ) in Fig. 9, where the string motion and characteristics for 0 < t < 5 are shown to make a surface. The initial string is a straight line given by the initial mass function, G 0 x, and it develops to generate the selffolding structure. After a finite t c the surface is three-folded. Simultaneously the characteristics pass over the origin. All leaves of this surface are the strong solution, at least locally. Now we construct the weak solution by patchworking leaves to define a single-valued function M(x; t) of x at any time t. After the patchworking, discontinuity singularities appear. How to cut and patch is determined by the RH condition.
We pick up a three-folded sector and draw V (x; t) and M(x; t) as multi-valued functions of x in Figs. 10 (a) and (b) respectively. The fermion potential has the characteristic shape of a multi-folded structure, which is called swallowtail in mechanics or spinodal decomposition in thermodynamics. The mass function M(x; t) is a fixed-t section of the surface in Fig. 9.  The weak solution must satisfy the continuity of V (x; t) according to the RH condition (29), and it is uniquely determined to give the solid line as drawn in Fig. 10 (a). Note that this selection can also be stated as that we take the maximum possible branch among the candidates.
This weak solution gives a solution for M(x; t) as is drawn in Fig. 10 (b). The points A, B, L, and R in (a) and (b) correspond to each other, where L and R in (a) are the same. The vertical line is the cut we take for patchworking and the mass function has a discontinuity from point L to R. We integrate this multi-valued function along the string path, and we obtain where we have used the RH condition in the form of Eq. (29). This assures that the total (signed) area surrounded by the string and the cut (vertical line) must vanish. In other words the two areas of the left-and right-hand side of the cut are equal to each other [23]. These geometrical representations of the RH condition are convenient to determine the weak solution from the string motion diagram M(x; t). When folding is generated, we just cut and patch so that we might keep the equal area rule.
We should stress here that the equal area rule is not a condition for the weak solution but just a resultant feature referring to the strong solutions in the backyard. These multi-leaf strong solutions are not seen by the weak solution, and thus it is impossible to refer to the backyard solutions in order to get the next time step. The target solution in the renormalization group equation is the effective interactions at a scale t and there is absolutely no information for the backyard strong solutions. In fact, the original RH condition (26) is equivalent to the conservation of the difference V (R) − V (L) in Eq. (28). After adding the fact that any singularity is generated during t development, we have the continuity V (R) = V (L) and consequently the equal area rule. The original condition (26) or (28) is enough to determine the weak solution uniquely. Therefore we do not need access to the backyard solutions to solve the renormalization group equation. The renormalization group equation determines its weak solution by itself, without recourse to the multi-leaf solutions.
The above equal area rule may remind the reader of Maxwell's rule in treating a coexistence situation, for example, a liquid and a gas, where the vertical line in Fig. 10 (b) looks like the coexistence isothermal line in the p(= x)-V (= M) phase diagram. However, Maxwell's rule has to refer to the backyard thermodynamic functions which are usually not well established in equilibrium thermodynamics. In this sense Maxwell's rule is to be classified as a phenomenological rule. Nevertheless we may assign such a jump line between L and R the coexisting states of the L and R states with the bulk domain structure. We do not discuss this viewpoint more since it is outside the scope of this paper. 13  As is proved in Eq. (37), the fermion potential is the action up to the initial constant in the mechanics analogy, and in fact the selection of leaves satisfying the RH condition will be done so that the maximum action principle holds, which is directly related to the minimum free energy principle, as is seen later.

Weak solution of the NJL model
We construct the weak solution for the NJL model explicitly. The Hamiltonian F has no explicit dependence on x, that is, the system is translationally invariant, Therefore the momentum M is conserved in the time development and the characteristicsx(t) are the contour of M. Then the string motion is driven by a horizontal move of each point as shown in Fig. 11. The initial position of the string is given by the bare four-Fermi interactions, which is a straight line. The velocity of each point of the string is given by, and it has the opposite sign to the momentum M. The velocity is suppressed for the large M region and it will vanish for large t (small ). The schematic motion of the string is drawn in Fig. 11 in the case of the super-critical phase, where the string starts being folded after a finite t c . Now we solve the characteristics by integrating Eq. (41) with the initial conditions and we get,x If t is small enough, the function C t is a monotonic function and then we can express the initial value x 0 by a function of the final value x(t) using the inverse of function C t , Finally, we obtain the solution by using this inverse function as follows, The solution in Eq. (45) is the strong solution as far as the inverse function C −1 t can be defined. This situation will be broken at finite t for a strong enough coupling constant G 0 . The monotonicity of function C t breaks at the origin first. We evaluate the derivative at the origin: This becomes negative at a finite t c , for strong bare coupling constant G 0 > G c , This result is equal to the argument in Sect. 1 using the blowup solution. Finally, we calculate the fermion potential exactly by integrating Eq. (36): where x 0 and M 0 = G 0 x 0 should be understood as a function of x through the inverse function C −1 t in Eq. (44). The fermion potential is also multi-valued where the inverse function C −1 t is multi-valued. We show the numerical plots of these results in Fig. 12 in the case of the super-critical interaction, G 0 = 1.005G c and t c = 2.65, and explain how to do the patchwork to construct the weak solution. Hereafter all the dimensional variables are rescaled by 0 to be dimensionless. The characteristics are plotted in Fig. 12 (a), and we find that they will start crossing each other after t > t c . The overpopulated region enlarges with a three-folded structure. The string motion is drawn in Fig. 12 (b), where we see that it starts folding after t c .
Then we have to cut and patch after t c . First, note that the original chiral symmetry  area. Then the only way to satisfy the RH condition is to put the discontinuity at the origin and adopt the up and down leaves, without using the middle leaf. Thus we obtain the weak solution, drawn in Fig. 12 (c)/(d). This weak solution is uniquely obtained for any initial four-Fermi coupling constant and it is always global up to the infrared limit t → ∞. From the conservation law point of view, the symmetry breaking and the RH solution after that are understood as follows. The current F is always negative and it is larger for larger |M|. Then the charge M accumulates towards the origin with increasing slope ∂ x M| 0 since the current at the origin is vanishing. At t c , the slope becomes infinite, the singularity appears and we start defining a discontinuous solution. Now the charge flows in to the backyard and flows out from the backyard at the origin, ignoring the vanishing current condition there. This is the RH solution.
We answer here the possible question of the necessity of the weak solution. If we agree that there is a singularity in M in fact, but only at the origin, then according to the total symmetry M(−x) = M(x), we may solve the system only for the half space (0, ∞). This is true. We demonstrated that the numerical analysis of the PDE in (0, ∞) with the free end boundary condition at 0+ could actually work and M(0+) is obtained to grow upwards after t > t c [5]. However, this method of getting the dynamical mass is really tedious and its reliability is not proved. Moreover, most important is that such a trick does not work for the first-order phase transition case in Sect. 7, where the singularity appears off the origin, pairwise, and then they move. In such a case the weak solution is the only method to reliably calculate the physical quantities. 16

Bare mass and the Legendre effective potential
In the previous section, we have obtained the Dχ SB weak solution of the mass function M(x; t) which has a finite jump at the origin. In this section we argue how to calculate the chiral order parameters, the dynamical mass of the quark, and the chiral condensate, in the infrared limit.
The physical mass of the quark should be defined by M(0; ∞) and it is actually not well defined even by the weak solution itself. This is a common issue of the spontaneous symmetry breaking and we take the standard method of adding the bare mass term m 0ψ ψ, explicitly breaking the chiral symmetry, to the bare action and investigate what happens.
The addition of the bare mass term does not modify our system of NPRG at all. It just changes the initial condition of the fermion potential as and the initial mass function takes the form As is stressed in Eq. (39), the NJL system is translationally invariant with respect to x. We shift the origin of the x coordinate as and completely the same calculation holds for the massless case in the previous section. Therefore the weak solution with bare mass is obtained from the massless case by shifting the x argument as follows: As for the fermion potential, we have to take account of the global shift of the bare potential in addition to the shift of x, This looks miraculous considering that the RGE system with bare mass is completely different from the massless case; all the operators are coupled together. Now the mass function at the origin is well defined and we can define the physical mass by the following chiral limit, (55) Figure 13 shows the RG evolution of the physical mass for m 0 = 0 and its chiral limit. The physical mass in the chiral limit shows the second-order phase transition due to the singular behavior of the mass function at the origin, while the physical mass at m 0 = 0 shows the cross-over behavior.
The reader may consider that the weak-solution method is not necessary if we consider m 0 = 0 from the very beginning. This statement is absolutely wrong. First of all, note that the mass function is a function of operator x totally describing the effective interactions. Even with the bare mass included, the mass function still has completely the same singularity. The singularity does not disappear; it is just moved away from the origin. Actually, Ref. [5] shows that the Taylor expansion method to solve the PDE (16) does not work with the small bare mass.
Second, if we are interested only in the origin of the mass function, then there seems to be no singularity as just drawn in Fig. 13. However, this is the case of the second-order phase transition. 17  Such a disappearance of singularity does not hold for the first-order phase transition, which will be explicitly argued in the next section. The difference here is that the second-order phase transition is directly related to the symmetry breakdown and therefore is sensitively affected by the explicit breaking term, while the first-order phase transition is not related at all to the symmetry breaking and therefore is not altered much by the explicit breaking term. These points will be made clear in the next section.
In order to evaluate the chiral condensate and the Legendre effective potential for it, we first introduce the source term jψψ in the bare action. Then the initial condition of the fermion potential is The free energy as a function of j is given by the value of V at the origin, and the chiral condensate is defined by Note that these quantities are all defined at renormalization scale t. We define the Legendre effective potential of the chiral condensate as where j in the right-hand side is considered as a function of φ through the inverse of Eq. (58). We have the standard relation, As seen previously, because of the translational invariance of the system with respect to x, the j-dependence of V is determined by applying Eq. (54): Thus the free energy and the chiral condensate are calculated through the quantities with j = 0 as follows: Therefore, we conclude that the chiral condensate function φ( j; t) has the same structure of multivaluedness as the mass function, as shown in Fig. 14 (a). We have obtained the weak solution of the mass function M(x; t), which determines the weak solution of the chiral condensate φ( j; t). The weak solution of the mass function is easily obtained by the equal area rule. This rule does work as well for the weak solution of the chiral condensate, since the difference between functions of M and G 0 φ is a well-defined single-valued function, m 0 + j, and it does not contribute to the total area. We denote the discontinuity of φ at the singularity by φ (L) and φ (R) , which are the values on each side of the singularity.
The behavior of the Legendre effective potential in the corresponding multi-valued sector is shown in Fig. 14 (b). Note that between φ (L) and φ (R) , the function j (φ) is not monotonic and therefore the Legendre effective potential in that region breaks the convexity (monotonicity of the derivative). We evaluate the jump of the Legendre effective potential at the singularity ( j = j s ), where we have used the continuity of the fermion potential V proved in Eq. (29), which ensures the continuity of the free energy W at the singularity. Taking account of the fact that j is equal to the derivative of the Legendre effective potential at both φ (L) and φ (R) , the condition in Eq. (64) means that there is a common tangent line between these two points. Therefore the weak solution of the Legendre effective potential is the envelope of it. In other words the weak solution condition automatically convexifies the Legendre effective potential. Any patchworking of function φ( j) in Fig. 14 (a) does ensure the monotonic change of j as a function of φ. However, it is not enough to have the convexified effective potential. If we do a wrong patchwork (wrong placing of j s ), violating the RH condition, then the fermion potential V has a discontinuity, and accordingly the Legendre effective potential V L is not even uniquely determined as a function of φ.
In this paper we do not demonstrate how to calculate these potentials for the zero-density NJL model; instead, we make plots for the finite density model in the next section. By examining the 19  plots there, readers may easily understand what happens for the second-order phase transition case of the zero density NJL, since it is simpler and more straightforward.
The automatic convexification proves the greatest feature of the weak solution, although it is not so obvious in the case of the second-order phase transition, where the convexified effective potential just shows up as that with a flat bottom connecting the symmetry-breaking pair of vacua. In the case of the first-order phase transition, however, this is crucial, since the convexified effective potential correctly selects the globally lowest free energy vacuum even when there are meta-stable vacua around. We will demonstrate this in the next section.

First-order phase transition at finite chemical potential
In this section we investigate the first-order phase transition observed at finite chemical potential (μ = 0). The first-order phase transition is highly non-trivial compared to the second-order phase transition because the RG evolution of the physical mass has a finite jump even with the bare mass m 0 = 0 (as shown in Fig. 15). Therefore the bare mass does not help at all. Moreover, the effective potential has multi-local minima, and most important is whether the NPRG solution can pick up the physically correct, globally lowest free energy vacuum or not. We show that the weak solution at finite chemical potential is uniquely constructed, and the effective potential obtained is automatically "convexified," which ensures that the solution is the physically correct lowest free energy vacuum.
We consider the finite density NJL model. The chemical potential μ is introduced by adding the term μψγ 0 ψ to the bare Lagrangian (8). To make the NPRG equation more appropriate for the future finite temperature calculation, we use the spacial three-dimensional momentum cutoff, to define the renormalization scale. Calculation of the β-function for the Wilsonian effective potential proceeds as follows. The β-function is given by the shell mode path integration where the inverse propagator matrix O is given by We omit the large-N non-leading terms (diagonal elements in O) and we get where (x) is the Heaviside step function and C is a divergent constant independent of μ, , and M, and thus we ignore it here. The RGE for the mass function is obtained as (69) The translational invariance with respect to x still holds. Note that due to the change of the cutoff scheme the critical coupling constant at vanishing chemical potential changes, and we now have G c = 2π 2 / 2 0 . The characteristic curve is given by the following equation of motion, The particle motion stops suddenly at some renormalization scale due to the (x) factor, and the smaller the momentum M, the earlier the particle stops. Now we solve the system using the characteristics and construct the weak solution in order. In Fig. 16, we show the characteristics for the finite density NJL model with m 0 = 0, μ = 0.7, and G 0 = 1.7G c . The characteristics start folding at t c = 0.41 at some finite x place, and it occurs at two places simultaneously. This is clear in Fig. 17, where the bare mass is included to be m 0 = 0.01 and we plot the mass function M in the leftmost column, the fermion potential V in the middle column and the Legendre effective potential V L in the rightmost column, respectively. Note that we set a finite bare mass, and accordingly the mass function and the fermion potential are translated to the negative x side as derived in Eq. (61). If there is no bare mass at all, the chiral symmetry ensures these two singularities start symmetrically, keeping the odd function property of the mass function.
This situation is quite different from the second-order phase transition at vanishing density. There, the appearance of the singularity means Dχ SB, and it occurs at the origin. Now this first emergence 21  of the folding does not mean chiral symmetry breaking, and singularities are generated pairwise. It just means that the convexity of the Legendre effective potential is broken at some finite x positions. It should be noted that this is not the emergence of the meta-stable state yet. The pair of singularities grows up and we define the weak solution of the mass function M by using the equal area rule just as before, which is depicted by the dashed line in the mass function. As for 22/27 the fermion potential V , the weak solution, which must be continuous, just takes the upper envelope. It is proved easily that there is no other way of making the continuous function. This selection is equivalent to the maximum potential principle, and it is related to the maximal action principle in the viewpoint of the kinematical system. This variational property of solution can be directly related to another type of weak solution, the viscosity solution [24,25], which will be treated elsewhere.
As time goes by, the singularity points move towards the symmetry center at x c = −m 0 /G 0 , which is the origin in the case m 0 = 0. In the midst of this movement, between t = 0.4 and 0.5, the three folded branches at the right-hand side cross the origin x = 0. This happens exactly when there appears a meta-stable state, a local minimum, if we faithfully adopt all the multi-valued branches. This is seen in Eq. (63) where, in the above situation, the source j can reach zero on a local (unused) strong solution M, and it means the derivative of the Legendre effective potential calculated from the multi-valued M vanishes there.
The pair of singularities reaches the symmetry center x c at t = 0.72. Then the singularities are merged into one singularity on the symmetry center, and it does stop and will not move forever. This is the only way of satisfying the RH condition, which is understandable by the equal area rule. Note here that there are four regions contributing to the area. After this merging of pairs of singularities, all the functions M, V , V L given by the weak solution are very much like those of the second-order phase transition after the breakdown. There are no effects of meta-stable state near the origin.
In Fig. 17, the special time t = 0.5615 is picked up. This is the time of the first-order phase transition for these specific parameters of the model. It is characterized by the time that the right singularity crosses the origin. The vanishing derivative point of the Legendre effective potential jumps from the central well to the right-hand side well as seen in the corresponding figure (4c). At this phase transition point the chiral condensate does jump, although the bare mass is included to explicitly break the chiral symmetry, which is drawn in Fig. 15.
Thus the bare mass does not help us to avoid the jump of the physical quantities in t-development in the first-order phase transition case, and even in that case the weak solution successfully describes this jump behavior. We may in principle regularize this jump by limiting the total degrees of freedom to be finite. Then there cannot be any singularity in all functions and the jump becomes just a quick change of the physical value. In this case not only the second-order phase transition but also the first-order phase transition are smoothed out to be the crossover.
Note that most essential point of these procedures is that the phase transition is correctly described so that the lowest free energy vacuum is automatically selected. This selection is done by the weak solution, not by hand. The weak solution is unique, that is, starting with the regular and smooth function as the initial condition, developing by time, encountering singularities, then it automatically treats the multi-folding structure to pick up a single-valued function, and the resultant solution is physically correct.
The moving of singularities are drawn in Fig. 16, where note that the symmetry center is the origin (m 0 = 0). Pairwise-generated singularities move towards the origin and merge into one, which stays at the origin. The characteristics are all moving into the discontinuity, and thus it satisfies the so-called entropy condition.

Summary
In this paper, we have introduced the weak solution to define the singular Dχ SB solution of the NPRG equation that can predict physical quantities such as the physical quark mass and the chiral 23/27 condensate. The weak solution satisfies the integral form of the partial differential equation. Specifically, we have evaluated the weak solution of the large-N NPRG equation of the NJL model for the mass function which is the first derivative of the fermion potential with respect to the scalar bilinearfermion fieldψψ. We also applied our method to the finite-density NJL model where the first-order phase transition appears to give a highly non-trivial situation compared to the second-order phase transition case of vanishing density.
There is no global solution for our NPRG equation unless we do not consider the weak solution, as far as it describes the Dχ SB due to the quantum corrections, that is, the mass function must encounter singularities at a finite scale. The weak solution helped us perfectly to define the global solution up to the infrared limit of the infinite time where we can evaluate physical quantities. The weak solution is uniquely obtained, and there is no ambiguity or approximation to get it.
The Rankine-Hugoniot condition required by the weak solution is converted into the equal area rule, by which we can easily obtain the weak solution from the multi-folded local strong solutions. However, some questions remain. Since our partial differential equation is the renormalization group equation it must be solved with the initial condition only, that is, the effective action at a time. After it has discontinuity, it does not know the discarded strong solutions behind. The effective interactions must determine its development with the information given by the weak solution only. Therefore there can be no notion of "area" referring to all multi-valued leaves.
This is true. We should note again that the equal area rule is a result given by the weak solution, not an assumption to determine the weak solution. We just use it in order to easily obtain the weak solution. In fact, the Rankine-Hugoniot condition is written down as continuity of V or conservation law for M, it is the local condition without referring to the strong solutions behind, and it can determine the solution without recourse to the discarded strong solutions. Thus the weak solution with singularity can develop by itself without any additional information.
The weak solution defined in this paper has a perfect feature to describe the physically correct vacuum even when there are multiple meta-stable local minima given by all the local strong solutions in the Legendre effective potential. The basic logic to achieve this is that the weak solution successfully convexifies the Legendre effective potential, and it ensures that the lowest global minimum determines the vacuum.
In this paper we have worked out only the NJL-type models with a large-N leading approximation. We can go ahead to define the weak solution of the NPRG equation for gauge theories even with nonladder (large-N non-leading) effects. For example, the large-N leading NPRG partial differential equation is defined by which has the same type of PDE analyzed in this paper with translationally non-invariant Hamiltonian F. Also, the extension including large-N non-leading diagrams has been done to give where B = ∂ x V + 2C 2 πα s x/ 2 and G = 2C 2 πα s x/ 2 [5,6]. This is also manageable by our weak solution method. Applications to gauge theories are reported in forthcoming paper. 24  (73) Therefore this particle is a usual classical particle but with negative and time-varying mass of 2π 2 e 2t / 2 0 [5]. The particle motion becomes gradually slower due to the exponential increase of mass. The string near the origin is approximated by a straight line. The velocity is proportional to the momentum, though with a negative sign, and the string near the origin just turns counter-clockwise until it stops (Fig. 18). This is the central engine for generating the singularity and creating the Dχ SB.
Note that the Dχ SB occurs exactly when the string at the origin becomes vertical. Then the derivative of the mass function M(x; t c ) with respect to x becomes divergent at the origin. We denote this slope byG(t) and it is evaluated asG The RG equation forG(t) reads This is nothing but the RG equation of the four-Fermi interactions obtained in Eq. (5), whose blowup solution encounters the singularity at t c . However, if we work with the motion of string, the verticality of the string does not mean anything singular; it is just vertical. Some authors define the inverse four-Fermi coupling constant g(t) ≡ 1/G (t), which satisfies the following RG equation [26]: This equation has the global solution In the strong phase (g < 1/4π 2 ), g(t) crosses the origin at t c and goes to the negative region without any singularity. It has a global solution up to t = ∞. This tricky procedure is not authorized by itself. However, it is perfectly right if we consider the string motion and just use the inverse slope instead of the slope to describe the string near the origin. Sometimes, in the case of finite-density media, 25 the inverse coupling constant g(t) goes to the negative region and returns back to the positive region again. Even in such a case, it is understood to have given the right result, that is, the chiral symmetry is not broken in the macro effective theory. In gauge theories the same engine does exist but there is a little difference. The gluons are massless and have no intrinsic scale, and the engine does not stop at the infrared. Therefore, the string rotates infinitely around the origin, which corresponds to the infinitely many unphysical solutions encountered in the SD-type analyses.
This singularity generation mechanism quite resembles the generation of a shock and its expansion procedure described by the non-linear wave equation of Burgers [19]. Also, our NJL weak solution satisfies the so-called entropy condition and is understandable as a physical shock. We have done this argument using the NJL (μ = 0) model near the origin. However, it can also be applied to any neighborhood of the singularity birth place. Thus we demonstrated that the Dχ SB is the spontaneous generation and growth of a shock describable by the weak solution of a PDE.
The generation and growth of a non-moving shock describes the "continuous" behavior in the second-order phase transition. On the other hand, a moving shock realizes the "jump" behavior of the first-order phase transition.
We stress again that this is the first demonstration and successful application of the weak solution for the non-perturbative renormalization group equation. Why does it work so well? One intuitive plausible argument follows. The definition of the weak solution is given by some integration equality without recourse to the singularities of the target function. Our target functions are effective interactions which are to be path-integrated with operators. Therefore such functions do not have to be regular as they are, and it may be enough that they satisfy necessary equations as a form of integration convoluted with smooth and bounded test functions.
In the introduction we claimed that encountering the singularity in the midst of the renormalization group solution is inevitable in any, approximated or not, model of the Dχ SB. If the PDE system here is the continuum approximation of the originally discrete many-body system, then we may add higher-derivative terms, for example, the dissipation term, as an improvement of the approximation, and such dissipation may regulate the singularity. This is well known in the case of Burgers' equation. In our PDE of the non-perturbative renormalization group equation, though, there may appear higher-derivative contribution in improving the approximation, but it cannot help at all, because our singularity is intrinsic and is not a matter of approximation.
There may be another way out by changing the system drastically. If we invent some regularization method to limit the number of degrees of freedom N to be finite, the spontaneous symmetry breaking cannot occur and all effective interactions and thermodynamic functions are singularity free. Then, after all calculations, we take the infinite N limit to get the physical results. However, no good method is known of realizing finite degrees of freedom regularization in formulating the non-perturbative renormalization group equation. We have found a way of limiting the depth of quantum corrections, and the Dχ SB is regularized. After obtaining the thermodynamic potentials we take the infinite depth limit, and finally reach the right functions with singularities. We will report it in a separate article.