Solving the QCD non-perturbative flow equation as a partial differential equation and its application to the dynamical chiral symmetry breaking

Non-perturbative renormalization group approach to the dynamical chiral symmetry breaking is an effective method which can accommodate beyond the ladder (mean filed) approximation. The usual method relying on the field operator expansion suffers explosive behaviors of the 4-fermi coupling constant, which prevent us from evaluating the physical quantities in the broken phase. In order to overcome this difficulty, we solve the flow equation directly as a partial differential equation and calculate the dynamical mass and the chiral condensates. Also we formulate a beyond the ladder equation and it gives almost gauge independent results for the chiral condensates.


I. INTRODUCTION
The Wilsonian renormalization group approaches to the continuum quantum field theory have been formulated [1][2][3] and developed in the past decades (see reviews [4][5][6][7]). In these approaches, the renormalization group (RG) flow equation is defined by a functional differential equation, the solution of which gives the partition function defined by the functional integral. Using this framework, we are able to obtain new approximation methods extracting the non-perturbative information of the partition function. Therefore we call this type of the RG "non-perturbative renormalization group" (NPRG).
In order to solve the flow equation approximately, we usually expand the equation in terms of the field operators and their derivatives. The derivative expansion has been applied to evaluation of the universal quantities such as critical exponents and anomalous dimensions.
For example, in the three-dimensional scalar field theories, the expansion with respect to the field operator without the derivatives converges very well [8][9][10][11]. Although it is difficult to confirm the convergence with respect to the order of the derivatives, the result of the expansion up to the 4-th derivatives agrees well with the Monte Carlo simulations [12].
The NPRG has been also applied to the analysis of the dynamical chiral symmetry breaking (DχSB) in the strong coupling gauge theories and its effective theory. As first noticed by Nambu and Jona-Lasinio [13], the scalar 4-fermi operators become the source of the DχSB. From the viewpoint of the NPRG, when we lower the renormalization scale of the effective action, the strong gauge interactions induce the effective 4-fermi operators, which bring about the DχSB at the low energy scale [14]. Unless there are explicit symmetry breaking terms such as mass terms, the running 4-fermi coupling constant diverges at a low energy scale. This explosive behavior is nothing but a signal of the DχSB [4,15].
On the other hand, the (improved) ladder Schwinger-Dyson (SD) equation has been frequently used for the analysis of the DχSB in the strong coupling gauge theories. This ladder approximation has the strong dependence on the gauge-fixing parameter. Moreover it is difficult to improve the ladder approximation systematically. However the framework of NPRG allows us to take account of the effects of the non-ladder diagrams in a systematic fashion [16,17].
Since the running 4-fermi coupling constants diverge at a critical scale as mentioned above, the RG flow can not go beyond the critical scale towards the infrared limit. In the simple framework of NPRG which maintains the chirally symmetric structure of the effective action, we can not evaluate the physical quantities in the broken phase such as the dynamical mass of the fermion and the condensates of the fermion bilinear composite operator.
To overcome this problem, we introduce the bare mass term, which also works as a source term of the chiral condensates [18]. We may expect that the bare mass prevents the divergent behavior of the 4-fermi coupling constants and might allow us to effectively evaluate the order parameter of the DχSB at the infrared limit scale. The field operator expansion of the NPRG equation, however, does not converge well at least in the region of the bare mass as small as the current masses of up and down quarks. Consequently we will take another way of directly solving the NPRG flow equation as a partial differential equation (PDE) without relying on any field operator expansion.
By the way, another method like the Hubbard-Stratonovich transformation has been used in many works to avoid the divergent behavior [19][20][21]. In this method, the composite operators of fermions are partially transformed to scalar fields, one of which obtains the nonzero expectation value as a chiral order parameter. Introduction of those scalar fields has a merit that the meson physics can be argued simultaneously. However it is more complicated to evaluate the convergence of the physical quantities with respect to the operator expansion including the scalar fields.
This article is organized as follows. In Sec. 2 we briefly review the flow equation for the effective average action. We introduce the basic truncation which projects the complete operator space onto the subspace relevant to the DχSB so as to solve the flow equation approximately. In Sec. 3 we examine the truncation method in detail and obtain two types of truncations: one corresponds to the ladder approximation, and the other contains the main parts of the non-ladder corrections. In Sec. 4 we explain how to evaluate the chiral order parameters. In Sec. 5, we solve the flow equation by using the field operator expansion, and examine the convergence of chiral order parameters. In Sec. 6, we directly solve the partial differential equation and show the results of physical quantities. In Sec. 7, we summarize our methods and results, and discuss further issues along the line of thought of this article.

A. Formulation
In order to evaluate non-perturbative effects of the quantum field theory, we introduce the so-called "effective average action" [3] that interpolates between the bare action and the full quantum effective action. For this purpose, we define the generating function with the infrared cutoff supplied by the cutoff term ∆S Λ [Ω] as follows: where Ω represents various fields generically and the functional integral is regularized by the ultraviolet (UV) momentum cutoff Λ 0 . The cutoff term ∆S Λ is defined as the following mass term depending on the momentum: The regulator function R Λ (p) suppresses the quantum fluctuations with the momentum lower than the infrared cutoff Λ. Therefore the regulator function satisfies and it implements the "coarse graining" in the Wilsonian method.
The effective average action is defined by a slightly modified Legendre transformation: It satisfies the following boundary conditions, provided that the regulator function R Λ (p) has the following properties, where we define the dimensionless scale parameter t = log Λ 0 /Λ. The right hand side should be called β functional, and consists of the second order functional derivative of the effective average action, It is considered as a matrix with respect to the momentum and the species of fields. This flow equation is firstly derived by Wetterich [3]. Because of the boundary condition (5), the flow equation interpolates between the bare action S bare and the full quantum effective action Γ.

B. Application to QCD
Hereafter we consider the N f -flavor massless QCD with N c -color. The bare action in Euclidean space is where T a is the generator of the fundamental representation of SU(N c ). This action has the chiral symmetry SU(N f ) L × SU(N f ) R , which is to be broken down to SU(N f ) V dynamically by the strong gauge interactions.
To evaluate the dynamical chiral symmetry breaking (DχSB) by using the nonperturbative renormalization group (NPRG) flow equation, we truncate field operators which are not essential to the DχSB from the complete operator space of the effective average action. Here we define the following truncated effective average action, where we use the covariant gauge with the gauge-fixing parameter ξ, and do not represent the ghost sector for simplicity. This truncated subspace of the complete effective action is spanned by the operators of the bare QCD action and the fermion self-interaction operator V (ψ,ψ; Λ), which we call the fermion potential. The discarded operators including higher derivatives or gluon fields may affect quantitative evaluation, but are not the essential sector to drive the DχSB. Actually we will confirm that the operators of the bare QCD and the fermion potential bring about DχSB using the NPRG flow equation [19].
Here it should be noted that the flow equation induces operators breaking the gauge symmetry, such as the mass operator of the gauge boson, because the cutoff function explicitly breaks the gauge symmetry. However we do not discuss this issue because we truncate out such operators. Eventually the flow equation of the gauge coupling constant agrees with that of the one-loop perturbation theory. We will concentrate on evaluating the flow equation for the fermion potential.
C. Techniques to derive the NPRG flow equation In the rest of this section, we will explain some techniques to explicitly write down the NPRG flow equation (7) for the fermion potential. We denote the degrees of freedom of the fields as a vector Φ t = (A a µ , ψ t ,ψ). Then the dressed inverse propagator including the cutoff function is given by The regulator function is defined by where (D −1 0 ) ab µν (p) ≡ p 2 δ ab (δ µν − qµqν q 2 (1 − 1 ξ )) and δ(p − q) ≡ (2π) 4 δ (4) (p − q). Functions, r(p) and r ψ (p) are defined to satisfy the properties (3) and (6).
Next we explain how to calculate the "super trace" in the NPRG flow equation (7). We transform the flow equation as follows [5], Here the symbol∂ t is defined bỹ where p ≡ 4 . Then we split the inverse propagator matrix (11) into submatrices as follows, where M BB (M FF ) corresponds to the second derivative with respect to bosonic (fermionic) fields, while M BF and M FB correspond to bosonic and fermionic field derivatives. Physically the submatirx M BB corresponds to the inverse free propagator of the gauge field, the submatrices M BF and M FB are the gauge interactions, and the submatrix M FF contains the inverse propagator of the fermion and the fermion self-interactions. Using this notation, we can rewrite the "super-trace log" in the flow equation (13) into the following formula [22]: Note that another expression, has been often used in the NPRG analyses. Here Eq. (16) is more appropriate for our purpose of improving the gauge parameter dependence.
To derive the flow equation for the fermion potential, the matrix M in Eq. (15) should be evaluated by replacing the fields with their zero-momentum components. Replacing

A. Scalar 4-fermi operator
The gauge interactions induce all possible fermion operators respecting the symmetry of QCD, which are enhanced by themselves when we lower the cutoff scale. Even in the truncated subspace of the effective action (10), we can not treat exactly all possible operators and interactions. Therefore, as an approximation of the flow equation for the fermion potential, we project its full operator space onto a specific subspace and restrict the interactions so that we evaluate the DχSB most effectively.
As in the QED with one flavor [16], the central operator for DχSB is undoubtedly the following scalar 4-fermi operator, where λ I (I = 1, · · · , N 2 f − 1) are the generators of the fundamental representation of SU(N f ), and λ 0 = 1 √ 2N f 1 flavor is defined so that they satisfy the proper normalization, It is the only chiral invariant 4-fermi operator which gives corrections to the mass operator, and it becomes the relevant operator in the region of the strong gauge coupling constant. Therefore, for a first-step approximation, we project the operator space of the fermion potential onto the subspace spanned by polynomials in the scalar operator ρ.
To project the flow equation onto the subspace defined above, we determine the coefficient from all possible operators included in the full fermion potential. This is equivalent to count the coefficients of powers of (ψψ) 2 , even though this operator itself is not chirally invariant.
It is due to the fact that (ψψ) 2 operator does not appear in chiral invariant operators other than powers of ρ. Eventually, we may work with a potential function in the simplest scalar operator σ =ψψ: Here we note that the original chiral symmetry is not maintained in this subspace, but the discrete chiral symmetry still remains: the Lagrangian (the fermion potential) is invariant under the following discrete transformation, The discrete chiral symmetry forbids the operators of the odd powers of σ, such as a mass term.
Next we pick up the interactions that are expected to be most important for the DχSB and for improvement of the gauge-fixing parameter (ξ) dependence. In Eq. (18), we further select the large-N c leading interactions in the fermion self-interactions, which leads to the following simplification: Applying the above approximation and the usual field renormalization, ψ → ψ/Z (13) and (16), we obtain the flow equation for the fermion potential, where η ψ is the anomalous dimension of the fermion field. Functions A(p) and B(p) are matrices in the space of the color adjoint representation and the Euclidean Lorentz vector space as follows: Here we should pay attention to the difference between the two types of traces: tr acts the space of the fermion's representation, the Dirac spinor and the color fundamental representation and the flavor fundamental representation, andtr acts the space of the matrix A, B.
The propagators of the fermions and the gauge bosons including the cutoff function in Eq.
(24) and (25) are defined by anomalous dimension η ψ (≡ ∂ t log Z ψ ) appears in Eq. (23), and the gauge coupling constant is renormalized as A . In the next two subsections, we will explain how to extract the scalar operators from the right hand side of Eq. (23).

B. Ladder approximation
In Fig. 1, the flow equation (23) is diagrammatically expressed by using the "corrected vertex" defined in Fig. 2. We regard A and B as the ladder element and the crossed element of the diagrams, respectively. The diagrams consisting only of the ladder element A are the ladder diagrams. Moreover we find that the diagrams consisting only of the crossed elements B are also the ladder diagrams if we rotate all the fermion lines of the diagrams to untangle the crossed gluon lines. Therefore the ladder approximation with only the ladder diagrams should be defined To project Eq. (28) onto the subspace of the scalar operators σ n , we adopt the following approximation rule of picking up σ, where we use T 0 ≡ 1 √ 2Nc 1 color . The right hand side of the above formula is nothing but the scalar part of the general Fierz transformation obtained by using the completeness of the space of the spinor and the color and the flavor. Thus the coefficient F n in the above formula is given by According to this rule, the summation in Eq. (28) is calculated as follows: where C 2 is the second Casimir invariant of the SU(N c ) representation, C 2 = N 2 c −1 a=1 T a T a . Finally we obtain the ladder flow equation, where we rescaled V and σ by common factor As for the regulator function, we adopt the following sharp regulator function, Performing the momentum integration in Eq. (32), we obtain the ladder flow equation as a partial differential equation (PDE), (34) Apart from the anomalous dimension term, this flow equation agrees with the local potential approximated Wegner-Houghton equation, and it was proven that the flow equation gives the results equivalent to the improved ladder Schwinger-Dyson equation [19]. Actually this is the reason why we call this approximated flow equation "the ladder".
Using the momentum scale expansion [23,24], the anomalous dimension of the fermion field η ψ = ∂ t log Z ψ is given by where we define the running mass, As will be seen in the numerical results in Sec. VI, the chiral order parameters given by the ladder flow equation (34) strongly depend on the gauge-fixing parameter ξ.

C. Beyond the ladder approximation
In order to improve the gauge dependence, we have to add the crossed element as well. We evaluate the flow equation (23) using the full corrected vertex in Fig. 2, which is calculated as follows: Here we will ignore the commutator term for simplicity. Note that this ignorance is consistent with our approximation that we do not include diagrams with gluon self couplings. From the following discussion we expect that this ignorance does not induce the strong gauge dependence in the truncated subspace.
We associate the corrected vertex with the gauge independent set of diagrams for the S-matrix in case of the Abelian gauge theory. In the non-Abelian gauge theory, the diagram exchanging one gluon is also needed for the gauge independence of the S-matrix. On the other hand, such corrections can not be added directly due to the one-loop nature of the NPRG β function, and therefore the gauge dependence which would be canceled in the Smatrix appears in the commutator term of Eq. (37). In the NPRG, the correction by the diagrams exchanging one gluon is treated through the effective operators such as ∂ µ F µνψ γ ν ψ.
However such operators are not included in the truncated subspace. Thus we may omit the commutator term without suffering strong gauge dependence.
Then, according to the general Fierz transformation (29), we can pick up the scalar operators as follows: Finally, adopting the sharp regulator function, we obtain the flow equation beyond the ladder approximation as the following partial differential equation:

IV. CHIRAL ORDER PARAMETERS
Now we explain how to evaluate the two chiral order parameters, the dynamical mass of quarks and the chiral condensates ψ ψ , which are generated by the DχSB.
In the framework of the NPRG, calculating the non-zero chiral order parameter is nontrivial because the NPRG flow equation maintains the chiral invariant structure of the effective action, which forbids appearance of the dynamical mass operator. 1 1 In the theories including the scalar fields whose symmetries spontaneously break, we can evaluate the nonzero expectation values of the scalar fields as the order parameters by searching the minimum point of its effective potential. In the theory we consider, this method cannot be used directly because the chiral order parameters are not such expectation values of the fields. Therefore, the scalar fields corresponding to theψψ are introduced by methods like the Hubbard-Stratonovich transformation in many works. In this article, however, we adopt another method as will be seen in Sec. VI. Because of this divergence, the RG flow can not go beyond the critical scale Λ c toward the infrared limit in the chiral invariant operator space. Therefore the divergence seems to imply that the chiral invariant RG flow cannot exist at the cutoff scale lower than Λ c , where the true RG flow might be in the chiral variant operator space including the mass operator.
The relation between the divergence and the DχSB has been discussed in Refs. [4,16,25].
Here, in order to go beyond the critical scale Λ c and to effectively evaluate the chiral order parameters, we introduce the bare mass term, which explicitly breaks the chiral symmetry, in addition to the chiral invariant bare action: [18] When the cutoff scale Λ(t) lowers, the running mass m(m 0 ; t), being m 0 at the initial scale t = 0, is rapidly enhanced around Λ c by the 4-fermi interaction, and consequently the enhanced mass suppresses the β functions due to the decoupling effect. Therefore the 4fermi coupling constant is expected to stay finite at the infrared scale. Taking the zero bare mass limit after solving the flow equation, the infrared limit mass becomes a chiral order parameter, the so-called dynamical mass of quarks: The bare mass m 0 also works as an external source for the composite operatorψψ, and its expectation value is evaluated as follows: whereḠ 0 (m 0 ; t) denotes V (ψ,ψ; t)| ψ=ψ=0 . Here the infrared limit valueḠ 0 (m 0 ; ∞) corresponds to the Helmholtz free energy.

V. FIELD OPERATOR EXPANSION AND ITS CONVERGENCE
Here we attempt to solve the flow equation using the field operator expansion so that we evaluate the chiral order parameters introduced in the previous section. As seen in the end of this section, however, the field operator expansion does not work well in this model. Since such infinitely coupled equations can not be evaluated numerically, we stop the expansion at some maximum power N: Let us see the truncation dependence of the dynamical mass m dyn calculated for the ladder flow equation (34). In App. A we explain the input parameters and the running gauge coupling constant which obeys the one-loop perturbative β function. In Fig. 3, we plot the running mass m(m 0 ; t) at the infrared limit ( t → ∞) for each bare mass, and show its dependence on the truncation order from N = 2 to N = 20. In the bare mass larger than about 0.02 GeV, the truncated solutions converge well with respect to truncation order N. It also shows that the mass is dynamically generated by the DχSB since the infrared running mass is much larger than the corresponding bare mass. In the small bare mass region, however, the truncated solutions do not converge but diverge more badly for larger truncation order. Therefore we can not take the zero bare mass limit straightforwardly in this method with the field operator expansion. It is also difficult to make a reliable extrapolation. The non-ladder flow equation (39)  Here we show the results using another method of avoiding the explosive behavior of 4fermi coupling constant. In Eq. (43), the fermion potential is expanded around the vanishing value (σ = 0). We expand it around the nonvanishing value σ 0 ( = 0), In this expansion, we can go beyond the critical scale Λ c without introducing the bare mass because the nonvanishing value σ 0 plays a similar role as the bare mass. Note that it does not work as an external source for the chiral condensates. The dynamical mass corresponding to Eq. (41) is given by the following limit, As seen in Fig. 4, however, the expansion around the nonvanishing value does not converge in the small value of the expansion point σ 0 .

VI. SOLVING THE FLOW EQUATION AS PDE
In Sec. IV and V, so as to evaluate the chiral order parameters, we have introduced two methods, introduction of the bare mass and expansion around the non-zero point, both of which are expected to avoid the divergent behavior of the 4-fermi coupling constant.
Note that the expansion around non-zero point means that the fermion potential V (σ; t) is evaluated at the non-zero point, σ = 0. In these methods, however, the field operator expansion does not converge in the small value region of the bare mass or the expansion point, and we cannot obtain reliable chiral limit of the dynamical mass.
Obviously the poor convergence of the field operator expansion originates from the divergent behavior of the 4-fermi coupling constant, which corresponds to the second derivative of the fermion potential V (σ; t). On the other hand, if the flow equation is solved as a partial differential equation (PDE), the fermion potential is expected to behave analytically at least except for the origin. 2 Therefore, we will directly solve the flow equation as a PDE without the field operator expansion.
In the practical calculation we solve the flow equation in terms of the mass function, For the numerical calculation, we will set a finite region for x, x L ≤ x ≤ x R . In order to approximate the mass function at the end point x L to be the dynamical mass, we need to choose a small enough x L << −1. Near the boundaries, where we can not take full 6 points, we use the 5-point formula, the 3-point formula, and eventually at the very end we take only the next point. Usually these constraints enhance the numerical error near boundaries. It will be found that actually these low order approximated definition of the finite difference near the boundaries does not induce the instability. In the practical calculation we choose the boundaries such that x L = −19 and x R = 0.
Through the procedure explained above, the coupled ordinary differential equations for the discretized mass function on the grid are obtained, and we numerically solve the equa-tions with respect to the dimensionless scale t using the 4-th order Runge-Kutta method. In Fig. 5 we present the RG evolution of the mass function given by the ladder flow equation (34). Here the Landau gauge ξ = 0 is adopted, and the anomalous dimension is ignored, which is consistent with the local potential approximation (LPA). We can see that the mass function is dramatically increased from the vanishing value when lowering the cutoff scale Λ(t). Particularly the dynamical mass generation is observed below a critical scale, rather rapidly in a short range of scale t. The infrared limit is reliably evaluated, whose size reaches the Λ QCD scale, and thus the DχSB occurs. If the chiral symmetry is not spontaneously broken, the mass function goes to be zero in the limit x → −∞. with the Landau gauge has been proved to give the result equivalent to the improved ladder SD equation. Actually the two chiral order parameter obtained now agree well with the ones obtained by the SD approach in Ref. [29], which assures the total consistency of our method.
We present the numerical results of the two chiral order parameters obtained from the two approximated flow equations, the ladder one (34) and the non-ladder one (39), with the various values of gauge fixing parameter ξ. Table I (Table II) shows the numerical values of the dynamical mass (the chiral condensates) with or without the anomalous dimension.
Here the chiral condensates are the renormalized ones at 1 GeV, ψ ψ 1GeV [29].  We show these results graphically in Fig. 6 Table I.

VII. SUMMARY AND DISCUSSION
In this paper we have derived a new approximated flow equation beyond the ladder approximation so as to improve the dependence on the gauge fixing parameter. We have developed various methods to evaluate the two chiral order parameters, the dynamical mass of quarks and the chiral condensates. These methods are expected to avoid the explosive behavior of the 4-fermi coupling constant in the course of solving the flow equation.
Within these methods, however, the field operator expansion, which is usually applied to solve the flow equations, shows poor convergence induced by the infrared singularity. Then, we stopped the field operator expansion, and solved the flow equation as a partial differential equation by using the grid method. We have obtained the chiral order parameters successfully without any instability or extrapolation.
As for the chiral condensates, we have seen that the gauge dependence of the non-ladder extended result almost disappeared. Therefore our non-ladder extended approximation almost respect the gauge invariance. A next step for further improvement of approximation would be to take into account of higher order operators including the first derivatives. This improvement is implemented by replacing the coefficient of the kinetic term Z ψ (t) with the function of σ, Z ψ (σ; t). Consequently we need to solve coupled partial differential equations in terms of the fermion potential V (σ; t) and the kinetic function Z ψ (σ; t).
Also we have seen that particularly in the Landau gauge the gauge dependent ladder result of the chiral condensates coincides with the almost gauge independent non-ladder extended one. This agreement might be related to a statement that only in the Landau gauge the ladder approximation satisfies the Ward-Takahashi identity. However, at finite temperature and chemical potential, this relation is broken. Therefore we would encounter the new phase structures by applying this non-ladder extended approximation to the hot and dense QCD in the framework of the non-perturbative renormalization group.
where we set a fixed dimensionless scale for t 1 to be t qcd + 1, and t ir is left as an infrared cutoff scale parameter. Obeying Ref. [29], t ir should be parametrized by ∆ ir as t ir = t qcd − 0.5 · (∆ ir + 1), and we take the following parameter: Λ QCD = 484 MeV, ∆ ir = −0.5. (A4) By the analysis using the ladder Schwinger-Dyson equation, it is confirmed that the physical quantities such as the chiral condensates are not sensitive to the choice of the infrared cutoff scale parameter t ir .
Finally, we discuss the initial condition of the fermion potential V (σ; t). The initial cutoff scale Λ 0 has to be large enough so that the obtained RG flow well approximates the renormalized trajectory starting from the ultraviolet limit, Λ 0 → ∞. At the ultraviolet limit, the fermion potential vanishes because the effective average action agrees with the bare QCD action, or if it exists there it would be strongly suppressed soon by its higher dimensionality.
In practical calculation, we take the vanishing fermion potential as the initial condition, and set the initial cutoff scale Λ 0 to be large enough so that the infrared quantities does not depend on the Λ 0 within a given numerical precision. Actually we set Λ 0 to be the Z boson mass scale, M Z = 91.2 GeV.