Mean-Field Expansion , Regularization Issue , and Multi-Quark Functions in Nambu – Jona-Lasinio Model †

In this paper, the results of the investigation of multi-quark equations in the Nambu–Jona-Lasinio (NJL) model in the mean-field expansion are presented. The multi-quark functions have been considered up to the third order of expansion. One of the purposes of our computations is the study of corrections of higher orders to parameters of the model. The important problem of the application of the NJL model is regularization. We compare the NJL model with 4-dimensional cutoff regularization and the dimensionally analytical regularization. We also discuss so-called “predictive regularization” in the NJL model, and a modification of this regularization, which is free of the Landau pole, is proposed. To calculate the high-order corrections, we use the Legendre transform method in the framework of bilocal-source formalism, which allows one effectively to take into consideration the symmetry constraints. A generalization of the mean-field expansion for other types of multi-quark sources is also discussed.


Introduction
Quantum chromodynamics as the theory of strong interactions has achieved great successes in the description of various hadronic processes.Due to the asymptotic freedom of this theory it explains processes at short distances.However, for the quantitative description of the region of intermediate and large distances, simplified effective models are necessary.
In 1961, Nambu and Jona-Lasinio [1] constructed a superconducting-type model based on four-fermion interaction for understanding a nucleon mass nature and a mechanism of spontaneous breaking of chiral symmetry.
Eguchi [2] and Kikkawa [3] reformulated the Nambu-Jona-Lasinio (NJL) model for quark fields (see also [4,5]).Hatsuda and Kunihiro [6] used an improved NJL model for description of dense and hot media of hadrons.This model can be applied to study light nuclei as well.The NJL model is currently the most successful model of quantum chromodynamics in non-perturbative region (see [7][8][9][10][11][12] for reviews and further references).
The applications of NJL model in recent years, has intensified in various problems of nuclear physics and a cosmology (see [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29] and references therein).Of great interest is the generalization of the NJL model, which contains the quark confinement mechanism (the so-called Polyakov-loop extended Nambu-Jona-Lasinio model, or PNJL model; see, for example [22] and references therein).Consideration of this model, however, is beyond the scope of our work.Among applications of phenomenological models are investigations of multi-quark functions, which are main topic of the present work.The formalism of multilocal sources are basic idea behind of this approach [30,31].The method of multilocal sources has been used in various field-theoretical studies and may be employed for similar calculations in effective model as well [32,33].
The multi-quark functions appear in the mean-field expansions (MFE) of the NJL model in higher-order terms.To obtain the MFE the iteration scheme of the Schwinger-Dyson solutions with the bilocal fermion source have been used [32,33].This iterative approach is fruitful because equations under discussion have a simple structure, and the solution of these equations at all orders is an algebraic problem.In this work we consider equations for correlation functions of the NJL model up to the third order of MFE.The leading order (LO) and a first step of this scheme include the equations for the quark propagator, the two-quark function and the next-to-leading (NLO) correction to the quark propagator.The second order of MFE includes the equations for the four-quark and the three-quark functions and the equations for the NLO two-quark function and next-to-next-to-leading order (NNLO) quark propagator.
We have solved the four-quark and three-quark equations.The third step of iterations, which includes equations for the six-quark and five-quark functions and the equations for the NLO four-quark and three-quark functions, is also investigated.
The NJL model in the MFE contains quark loops and corresponding ultraviolet divergences.Therefore, regularization procedure is the essential point for this model.Traditionally used regularizations for the NJL model are based on a four-dimensional (4D) cutoff in Euclidean momenta or a three-dimensional momentum cutoff.The Pauli-Villars regularization, or non-local Gauss form-factors, are also used.However, the most popular and very convenient for calculations dimensional regularization is the least used for the NJL model.The reason is related to the next issue: a regularizations parameter for NJL model is involved in physical quantities.Therefore, it is the essential parameter for NJL model.However, the parameter of dimensional regularization is a deviation from the physical dimension of space, which cannot be connected directly with any physical quantity.
An alternative approach is a consideration of the dimensional regularization as a variant of the analytical regularization.In framework of this treatment all computations are made in 4D Euclidean momentum space.For regularization of the divergent integrals, a weight function into integrand is included, and the regularization parameter is a power of this function.This variant of dimensional regularization is consistently developed and applied to the NJL model by Krewald and Nakayama [34].The parameter of such regularization is not at all a deviation from the physical dimension of space.We assume that this parameter can be treated as a "trace of gluons" in the effective four-quark self-action of NJL model.Such treatment in some sense is similar to the non-local versions of the NJL model.
In our paper we analyze the dimensional regularization in the NLO of the MFE in the frame of Krewald and Nakayama approach.(We introduce the new term "dimensionally analytical regularization" (DAR) for avoiding an unnecessary association with the traditional dimensional regularization) An interesting attempt for a deliverance of the NJL model from the regularization dependence made in the work of Battistel, Dallabona and Krein (BDK) [35] (see also [36,37]).Main idea of the method, named by authors as "predictive regularization", consists of a separation of the finite parts of divergent integrals from the divergent ones, and then the integration is done without regularization.A significant result of work [35] is a proof of the fulfillment of symmetry constraints for correlation functions.We discuss some features of the computational scheme of work [35].The main problem of this method is the singularity of meson propagators in the region of Euclidean momenta.The existence of these Landau singularities is an essential problem for any computations beyond the leading approximation.A modification of this scheme is proposed, which is free from Landau poles and conserves the main features of the approach of [35] (see also [38][39][40]) The structure of the paper is following.In Section 2 we construct MFE in the fermion bilocal-source formalism.The DAR of the NJL model is discussed in Section 3. The predictive regularization and the problem of Landau ghost is discussed in Section 4.
A purpose of our computations is to study the corrections of higher orders to parameters of the model, such as the chiral condensate and pion-decay constant.The sigma-meson and pion contributions to the chiral quark condensate are considered in Section 5.The chiral condensate is the principal order parameter of dynamical chiral symmetry breaking (DCSB).As our computations demonstrate, the pion contribution to chiral condensate is inversely proportional to the regularization parameter of DAR and is independent of other parameters of the model.The sigma-meson contribution is small for admissible values of the parameter (see [41][42][43]).
In Section 6 we investigate the corrections to the two-quark function by the method of the Legendre transform with respect to the bilocal source [30,31].As is shown, this method is an effective way to take into account the constraints following from the chiral Ward identity (see also [44]).
In Section 7 we consider the higher orders of the iteration scheme.We calculate the three-quark function, which describes the decay σ → ππ and a correction to pion-decay constant.
In Section 8 the generalization of the MFE for the NJL model in the framework of the diquark and triple-quark sources formalism is shortly discussed.Such generalization can be fruitful, in particular, for the quantum-field description of nucleons.

The Mean-Field Expansion in the Bilocal-Source Formalism
Let us consider the chiral-symmetric NJL model.The model is based on u and d quark fields ψ(x).The two-flavor NJL model Lagrangian is where g is the four-fermion coupling constant with the mass dimension m −2 .Here we omit the bare masses of u and d quarks which are much smaller than the QCD scale.The chiral SU V (2) × SU A (2) invariance prohibits a mass term in the Lagrangian.In Equation (1) ψ ≡ ψ αc j , where α = 1, 2, 3, 4, c = 1, ..., N c and j = 1, 2 are spinor, color and isotopic (flavor) indexes correspondingly, and ∂ ≡ γ µ ∂ µ .The SU(2)-group generators τ a are normalized as tr(τ a τ b ) = 2δ ab .
The generating functional G of vacuum expectation values of T-products of fields (correlation functions) can be written as the following functional integral with bilocal quark-antiquark source η: The n-th functional derivation with respect to source η of generating functional G is the 2n-point (n-particle) function: The translational invariance of the functional integration measure implies relation that can be written as the functional-differential Schwinger-Dyson equation (SDE): To formulate the MFE we shall use the procedure which was proposed in Reference [45].A leading approximation is the SDE (3) with zero r.h.s.: The solution of the leading approximation equation is the functional where Tr is the operator trace and * is the multiplication of operators.A function S (quark propagator) in above formula for G (0) is the solution of following equation The leading approximation produces the linear iteration procedure in which the term with bilocal source η must be considered to be perturbation [32,33].Thus, we can construct the iteration scheme for generating functional This series consists of the step-by-step solutions of the following equations The functional Here P (n) are 2n -th degree polynomials on bilocal source η.
Leading-order correlator S is a quark propagator (a single-particle function).It is where m is the dynamical quark mass that is a solution of the gap equation of the NJL model in the chiral limit (In our work, we include a phase factor in an integration over momentum space d q ≡ d 4 q/(2π) 4 ).The basic order parameter, which defines a degree of DCSB, is a quantity where we take the trace over all discrete indices.It is easy to see that from Equations ( 6) and (7) it follows It is a regularization-independent formula.
The quark chiral condensate c is defined for each flavor separately and in the chiral limit considered here it is Equation ( 7) is the main equation, which describes the principal phenomenon of NJL model-DCSB.The divergent integral in Equation ( 7) should be considered to be a regularization (see Section 3).

First-Step Equations
In this subsection we discuss the structure of first step of iteration scheme.Gap Equation (7) always has the chiral-symmetric trivial solution m = 0.An energetically preferable physical solution m = 0, i.e., a solution with DCSB also exists, and below we shall consider a solution at m = 0 only.
The functional of the first step of iteration scheme can be represented as In this step two-quark function S 2 and NLO correction S (1) to quark propagator emerge.According to Equation ( 5) at n = 1 we obtain for two-particle functions S 2 and NLO correction to quark propagator S (1) following equations: (here and hereafter tr u denotes the upper line trace of functions S 2 ) and 1) ](0).
These equations are reduced in the momentum space to a system of simple algebraic forms.The graphical representations of two-quark function see on Figure 1, where the graphical notations of Figure 2 are used.

Pion and Sigma-Meson
Let us go to the amputated two-quark function and separate the connected part of F 2 As a result of this procedure we get following equation for two-quark function F c 2 : Accounting the color, flavor, and Lorentz structures of the inhomogeneous term of this equation we obtain the general form of solution where A σ and A π are the sigma meson (scalar) and pion (pseudoscalar) amplitudes, correspondingly.These are the solutions of the following equations Going to the momentum space and accounting for gap Equation (7), one gets and Here is the divergent single-loop integral that should be treated as some regularization.
The scalar and pseudoscalar amplitudes, as can be observed in Equations (19) and (20), contain the simple poles, which indicate to the scalar particle with mass 2m (the sigma-meson) and to the massless pseudoscalar particle (pion).Such appearance of the pseudoscalar massless particle (Goldstone boson) in the spectrum has agree with the general Nambu-Goldstone-Bogoliubov (NGB) theorem.
From equation S (1)  (12) for the NLO correction to quark propagator we can calculate the meson corrections to quark mass.This equation is transformed in the momentum space to a system of simple algebraic forms, which contain the NLO mass operator Σ (1) = S −1 * S (1) * S −1 .We obtain from (12):

Regularization Issue: DAR
The NJL model is based on the unrenormalized 4-fermion contact interaction, and the choice of the regularization scheme is an important issue in the effective use of this model.The integral I 0 in frameworks of different regularizations is finite for P 2 = 4m 2 and for P 2 = 0, (see, e.g., [8] for review) (The NJL model as the effective model of QCD is based on the color-current expansion and the Fierz transformation of the one-gluon exchange interaction (see, for example, [9] and references therein).Such approach implies the presence of some parameter which regulates the ultraviolet properties of the model.For this reason, a regularization with momentum cutoff seems to be very natural for this model.Other regularizations, although they do not have such a transparent physical interpretation, in some cases can be more convenient for computations or for the taking into account some special properties of the model), but a most popular gauge-invariant dimensional regularization should be used in this model with some special treatment.
We study in this section the NJL model with dimensional regularization in approach of Krewald and Nakayama [34] (see also [45][46][47][48][49][50][51]).Contrary to renormalized models, a regularization parameter of the NJL model is involved in formulas for observable quantities.This parameter is one of the essentials for the NJL model.However, dimensional regularization parameter, which typically is taken as a deviation of space dimension, does not allow any physical interpretation.Nevertheless, a different approach to the dimensional regularization exists: it can be considered as a variant of an analytical regularization.In this approach all calculations are made in 4D Euclidean momentum space, and the regularization parameter is treated as a power of a weight function, which regularizes divergent integrals.This approach to the dimensional regularization, based on ideas of Wilson and Collins, was consistently developed and applied to the NJL model by Krewald and Nakayama [34] in the mean-field approximation.We stress that in this treatment of dimensional regularization, the regularization parameter is not at all a deviation in the physical dimension of space.
The general characteristics of the approach are: (i) All computations are made in 4D Euclidean space; (ii) Translational invariance is supposed; (iii) The regularization procedure includes of modification of integration measure via the weight function which provides the convergence of integrals.
In this connection we shall use the term "dimensionally analytical" for this regularization to stress its properties which differ from the conventional approach of the dimensional regularization in a formal D-dimensional space.
Let us study of the NJL model gap Equation ( 7) for example.After performing integration by angles, the gap equation in the case m = 0 in Euclidean space is where Ω 4 = 2π 2 .Let us (according the previous comment) include to the integrand the weight function The weight function w Λ,D is the product of two weight functions w Λ and w D .The first function w Λ corresponds to the 4D cutoff regularization, and the second power function w D conforms to DAR.
The integration over dq 2 e results in Here B x (u, v)-the incomplete Beta-function.
(a) Cutoff: By setting D = 4, we have the following result: where we obtain where instead of g we introduce the dimensionless quantity Formula ( 26) corresponds exactly to the result of integration in D-dimensional space with the prescription however, in our situation the computation was performed in the usual 4D space, i.e., in our case D is not a dimension of space.It is a parameter which enables the convergence.In particular, we are not restricted by the limit D → 4 for the analysis of the results.We assume that a possible technique of regularization parameter is a power of some supplementary factor-the measure of gluon influence on the effective local 4-quark self-interaction of the NJL model.We choose the regularization parameter ξ as: (The parameter ξ differs from the typically used parameter ε = (4 − D)/2.They are related as ε = 1 + ξ.Introduction of this notation protect unnecessary associations with the usual definition of the dimensional regularization.) In terms of the parameter ξ the gap Equation ( 26) is The interval of convergence of the integral is 0 < ξ < 1.As we shall see below, this is also the interval of the physical values of the model parameters.
The LO chiral condensate correspondingly is

Two-Quark Amplitude and Model Parameters in Leading Approximation
First-step two-quark amplitude A (a connected part of amputated two-quark function) possesses the following color and flavor structure: Here A σ is the scalar amplitude, and A π is the pseudoscalar amplitude.In momentum space these amplitudes of the NJL model depend on a momentum p only, where p is a sum of quark and antiquark momenta.They have the form [46]: where L S (p) = ig d q trS(p + q)S(q) is the scalar quark loop, and where L P (p) = ig d q trS(p + q)γ 5 S(q)γ 5 is the pseudoscalar quark loop.
Using simple algebra and gap Equation (7), it is easy to obtain for A σ and A π in above-mentioned representations (19) and (20), where integral I 0 (see (21)) can be calculated as above.Going to Euclidean space, using a Feynman parameterization and changing an integration variable (which is possible due to translational invariance of the procedure), we perform the angular integration.According to our prescriptions, then we introduce into the integrand a weight function and calculate the integral over dq 2 e .For DAR we again obtain the result, which exactly corresponds to the result of integration with the formal transition to D-dimensional space: The integral over dq 2 e converges at −1 < ξ < 1. Taking into account the gap Equation ( 29) we obtain: where F(a, b; c; z) is Gauss hypergeometric function.
For 4D cutoff we correspondingly obtain: Formulas for the condensate and the two-particle amplitudes allow to fix values of the model parameters in the leading approximation of MFE.For this purpose, we use regularization-independent Formulas (8), ( 9) and a formula for pion-decay constant in the NJL model (see [8]): For DAR we obtain from (33): and for 4D cutoff (FDC) (see (34)): Correspondingly we obtain for DAR very simple formula: For 4D cutoff the analogous formula is These formulas allow the definition of the values of principal model parameters.We choose the value f π = 93 MeV for the pion-decay constant.Since chiral quark condensate c is not directly measured value, we shall determine sets of parameters for typical values of this quantity.For DAR it is necessary also to fix a value of M ("subtraction point").In work [45] we have used for this purpose a value of decay width π 0 → 2γ.Analysis of results of this work demonstrates that for very large range of condensate values (from −160 to −250 MeV) the value of M is practically permanent and is M ≈ 100 MeV.Since here we shall take M = 100 MeV.
The pseudoscalar amplitude associates with the pion, which in the chiral limit is a massless Goldstone boson.In both regularizations under consideration we can define a pion propagator as a pole term of the pseudoscalar amplitude: where I 0 (0) is defined by Equation ( 36) for DAR and by (37) for 4D cutoff.The situation is different for the scalar amplitude.Function I 0 (p 2 ) possesses in both regularizations a cut which originates in the point p 2 = 4m 2 .For 4D cutoff it is possible to define a scalar sigma-meson propagator as since is a finite quantity.However, for DAR I DAR 0 (4m 2 ) is finite only at ξ < −1/2: To interpret the sigma-meson as a particle in the NJL model with DAR we use the following trick: since in the interval −1 < ξ < −1/2 integral I 0 converges we use the value in the point p 2 = m 2 as a foundation for an analytical continuation of the pole part of the amplitude on parameter ξ to the physical interval 0 < ξ < 1.Then the sigma-meson propagator for DAR will be This formula was used for a computation of the sigma-meson contribution in chiral condensate in work [45].Surely, such procedure of definition of sigma-meson propagator seems to be a somewhat artificial.A more consistent procedure is a separation of a leading singular part of amplitude in the interval of physical values of regularization parameter ξ.
The separation of leading singularity near the point p 2 = 0 results for the pseudoscalar amplitude the same result (40), i.e., the pion in DAR possesses all properties of usual observable particle.For the scalar amplitude it is not so.At p 2 → 4m 2 we have in interval 0 < ξ < 1: , and, correspondingly, the leading singularity is and the leading singularity of scalar amplitude in the model with DAR is of the fundamentally different type in comparison with the cutoff model.Instead of the pole term, which can be naturally interpret as sigma-particle propagator, we obtain for DAR the power behavior which depends on the regularization parameter ξ.Thus, we come to the conclusion that for DAR at physical parameter values the scalar amplitude A σ does not possess a pole term, which can be interpret as a physical scalar meson.

Regularization Issue: Predictive Formulation of the Nambu-Jona-Lasinio Model and Ghost Problem
A very interesting attempt for deliverance the NJL model from the regularization dependence was made in the work [35] (see also [36,37]).The main idea of this approach is to avoid the explicit evaluation of divergent integrals with any specific regularization.The finite parts of integrals are separated of the divergent ones and are integrated without any regularization.Then the NJL model becomes predictive in the sense that its consequences do not depend on the specific regularization of divergent integrals.An important result of work [35] is a proof of the fulfillment of all symmetry constraints.The choice of parameters results in the acceptable values of the quark mass and other parameters of the model.
In this section we analyze some features of this computational scheme We shall name the computational scheme of work [35] as the BDK approach, or the implicit regularization.In particular, we point a connection of the BDK approach with the differential regularization of Gelfand and Shilov [52].The hard problem of the BDK approach is singularity of meson propagators in the Euclidean region of momenta.The presence of this singularity (Landau pole, or Landau ghost [53]) prevents to any computations beyond the leading approximation.We discuss a modification of the scheme (see also [38][39][40]), which is free of Landau poles and maintains the major features of the BDK approach.
The leading approximation of the model is the mean-field approximation, which coincides with the LO of 1/N c -expansion.All correlation functions of the leading approximation are given in terms of single-loop integrals.The problem of computations of these integrals can be transformed to the definition of following five divergent integrals (see [35]): Here m is the dynamical (constituent) quark mass, which is a non-trivial solution of the gap equation for the NJL model.
A basic tool for the definition of integrals ( 44) and ( 45) in the BDK approach is an algebraic identity for the propagator function (see Equation (28) in work [35]).With the identity the divergent parts of integrals ( 44) and ( 45) have been rewritten via three tensorial and two scalar external-momentum-independent integrals which were treated in the sense of some unspecified regularization.Then integrals (44) and (45) have been represented in terms of these five integrals and two standard convergent integrals (see formulas of section III in work [35]).
In this connection we note some general regularization-independent property of integrals ( 44) and ( 45), namely integrals {I 1 ; I µ 1 } and {I 2 ; I µ 2 ; } are connected with following relations and These expressions are regularization-independent and play a significant role in the method.It is not so difficult to verify the validity of these relations for conventional regularization schemes.However, their validity is not evident directly from above-mentioned formulas of work [35].To prove these relations, we must define derivatives of the external-momentum-independent integrals over m 2 .For this purpose, we notice that derivatives of logarithmically divergent integrals are convergent integrals and can be calculated without any regularization.This calculation gives us zero value for derivatives of tensorial logarithmically divergent integrals.Then one can verify that the derivatives of quadratically divergent integrals are the corresponding logarithmically divergent integrals.Taking into accounting the circumstance we can easy to derive Equations ( 46) and (47) for BDK expressions of integrals (44) and (45).Since relations (46) and ( 47) can be used for alternative definitions of quadratically divergent integrals I 1 and I µ 1 without additional regularization, their importance is clear.Using the expressions ( 44) and (45) for integrals in work [35], the symmetry constraints on single-loop correlation functions have been analyzed.Based on these relations all symmetry properties of the theory (such as Furry theorem, Ward identities, etc.) can be fulfilled for all single-loop correlation functions.This point is one of the important results of [35].From the point of view of above discussion the BDK consistency relations simply assert zero values of corresponding integration constants.
The choice of parameters of the NJL model with Lagrangian (1) in the leading approximation is mainly defined by two divergent integrals, namely logarithmically divergent integral I 2 and quadratically divergent integral I 1 .Integral I 1 is a part of the gap equation which defines dynamical (constituent) quark mass m.Integral I 2 determines the structure of meson propagators.Both these integrals can easily be defined with differential regularization without using the above-mentioned algebraic identity for the propagator function.
To define I 2 one can use well-known technique of Gelfand and Shilov [52].Specifically, let us introduce new external variables p = k 1 − k 2 and P = k 1 + k 2 .Then derivatives ∂I 2 /∂p µ , ∂I 2 /∂P µ are convergent integrals, and their computation gives us Basing on these results and using identity ∂p µ , we naturally go to the following definition: where m 0 is the integration constant.This expression coincides with the BDK ones if constant m 0 is related to I 2 (0) as The permutation of two limits-differentiation and regularization removing-is an essential feature of the above computation.Sure, such permutation is implied also in work [35] in the computations of finite parts.This permutation is an essence of Gelfand-Shilov differential regularization [52], which is based on the infinite differentiability of generalized functions.
To define integral I 1 we employ regularization-independent relation (46), which gives where m 2 1 is the integration constant of Equation (46).This definition also corresponds to BDK ones.To calculate m 0 and m 1 and, therefore, to define I 0 and I 1 in full, we take into account two regularization-independent relations of NJL model, namely where f π is the pion-decay constant, and where χ =< 0| ψψ|0 > and c is the quark condensate.We have from Equation ( 53) Then, using Equations ( 52) and ( 54), we obtain for m 1 Equation ( 56), which we shall refer as BDK equation, plays a key role in the BDK approach.
(3) at a > 1 Equation (57) possesses one negative root x 1 < 0 and two positive roots 0 < x 2 < 1 and At c < 0 the second case (a = 1) alone is physically accepted.In work [35], only this case as a solely possible choice for model parameters is taken.The value of the dynamical quark mass, which corresponds to positive root which is also regularization-independent (see, e.g., [8]).Noteworthy, the quark-mass value and the coupling value depend on the quark condensate value only, and do not depend on the pion constant value.The last one defines the value of m 1 , which is m 1 = 879 MeV at f π = 93 MeV.
Let us compare the parameter values of the BDK approach with the parameter values of other regularizations.The value of quark mass (468 MeV) in the BDK approach is noticeably greater compared to the values in the 4D momentum cutoff regularization (236 MeV) and the Pauli-Villars regularization (240 MeV) (at given value of quark condensate c = −250 MeV).Moreover, the quark-mass dependence on the condensate value is quite different.In the BDK approach the quark mass is proportional to the condensate, whereas in four-momentum cutoff and Pauli-Villars regularization the quark-mass increases at decreasing the absolute value of condensate.At c = −210 MeV the quark mass is m = 393 MeV in the BDK approach and m = 423 MeV in the 4D momentum cutoff.At the same time, from the point of view of the phenomenology the parameter values in BDK approach are quite reasonable, and expressions for the correlation functions are much simpler in comparison with any traditional regularization.
Analytical properties of correlation functions in the Euclidean region manifests an essential difference of the BDK approach from other regularizations.Meson propagators in the BDK approach possess a non-physical singularity-a pole at the negative momentum square The existence of similar pole was discovered firstly in quantum electrodynamics [53].The existence of the Landau ghost in a system of fermions coupled to a chiral field has been observed in work [54].It is a characteristic feature of theories without an asymptotic freedom in deep-Euclidean region.This Landau pole is a serious problem of the BDK approach, since any calculations with meson loops become problematic.Though a subject of [35] is a single-loop approximation, and such calculations exceed the framework of this work, its impracticability means a principal impossibility of computations of corrections to leading approximation and cannot be acceptable.
Moreover, due to the smallness of fine structure constant α the Landau pole in quantum electrodynamics is in the very distant asymptotic region ((M 2 L ) QED −m 2 e exp{ 3π α }, where m e is the electron mass and α 1/137), and its presence can be neglected.Indeed, at the much smaller energies the quantum electrodynamics becomes a part of an asymptotically free grand unification theory with self-consistent asymptotic behavior.In the NJL model with the implicit regularization the situation is different.The Landau pole is placed near the physical area of the model, and above reasoning is impossible in principle, since Landau mass value M 2 L = −p 2 L only twice larger in comparison with the quark mass: M L 2m.The conventional regularizations, such as the 4D cutoff or the Pauli-Villars regularization, are free on this problem-the meson propagators in these regularizations do not have the Landau poles.Therefore, BDK approach, with a certain appeal and simplicity, contains the serious shortcoming as the nearby Landau pole.
Hence, a compromise approach is needed, which maintains the major features of implicit regularization and at the same time solves the problem of the Landau pole.Such compromise can be reached with the Feynman regularization for logarithmically divergent integral I 2 : where M r is a regulator mass (M 2 r > m 2 ), and the definition of quadratically divergent integral I 1 as before is made by Equation (46).At p 2 < 0 integral (59) can be rewritten as where One easily can prove the absence of zeroes of this function at x > 0. Evaluating elementary integrals in Equation ( 60) and introducing variables v = 1 2 ( 1 + 4 x − 1) and v r = 1 2 ( 1 + 4M 2 r xm 2 − 1) we obtain the following expression: Since M 2 r > m 2 , then v r > v, and from elementary inequalities log 1+v r 1+v > 0 and v r log(1 ) it follows that F > 0, i.e., the meson propagators do not possess the Landau pole with this definition.Other features of this modified regularization are similar to implicit regularization.In particular, for this modified regularization and I 1 is defined by same formula (52) with substitution m 0 → M r , i.e., integration constant m 0 everywhere is substituted by regulator mass M r .This re-definition of I 1 still implies an existence of additional parameter m 2 1 , which is the integration constant of equation (46).Furthermore, Equation (56) for the quark mass m, has the same form for the modified regularization, and, consequently, all parameter values are the same.Note that the condition M 2 r > m 2 is the direct consequence of formula (53).Therefore, this modification conserves main features of the implicit regularization and simultaneously solves the problem of Landau pole.(An alternative method ridding the Landau ghost is the method of [55], based on the Källen-Lehmann representation.This method has been used in work [56] to the chiral σ-model.However, in the framework of the implicit regularization the suggested technique is likely to be preferable since this regularization deals with a set of integrals while the method of [55][56][57] should be used to the proper two-point functions) To define divergent integrals (44) and ( 45) one can use the following rules for integrals (45): (The upper index BDK means that each integral with mass m or M r is given by formulas of [35]).Then integrals (44) are defined by Equations ( 46) and (47).With such definitions integrals I 2 and I µ 2 are convergent without additional regularization.A consistency relation is the coincidence of integration constants, and It should be noted that due to linearity of the consistency constraints considered in work [35], the analysis of symmetry conservation, which was made in this paper, can be carried to the proposed modification too.(For the more details of the modified implicit regularization see work [38].Further developments of these ideas see in works [39,40])

Meson Contributions in Chiral Condensate and in Quark Propagator
Although for most phenomenological applications of the NJL model enough to consider the LO of MFE, or the LO of 1/N c -expansion (In the LO the results of both approaches coincide), an analysis of the structure of the NJL model in NLO of the expansions also of considerable theoretical interest (see [58][59][60][61][62][63][64][65] and refs.therein).This analysis is necessary for clarification of the area of applicability of results and their stability under a variation of parameters and quantum fluctuations due to higher-order effects (Please note that the destabilization induced by the quantum fluctuations do not necessarily have the direct physical meaning due to absence of the confinement mechanism in the NJL model).In this section, we will discuss such corrections in the framework of MFE.
First-order equations of the MFE (see [41][42][43]45]) define corrections to quark propagator.First-order mass operator Σ (1) = S −1 * S (1) * S −1 , where S (1) is a first-order correction to quark propagator, is defined in x-space by equation Introducing dimensionless first-order mass functions a (1) and b (1) : and defining the first-order condensate and a ratio of the first-order condensate to the LO condensate ρ ≡ χ (1) /χ (0) , we obtain from (64) the expressions for a (1) and b (1) in momentum space: It follows from Equations ( 67) and ( 68) that the corrections to quark propagator consist of two parts: pion corrections (due to pseudoscalar amplitude A π ) and contributions due to scalar amplitude For the ratio of the first-order condensate (66) to the LO condensate (8) we obtain ρ = −gχ (1) After the computation, the condensate corrections one can to define the corrections to quark mass.Inverse quark propagator is If the propagator has a pole in point p 2 = m 2 r , then b(m 2 r ) = m r a(m 2 r ), and expanding a (1) (m 2 r ) and b (1) (m 2 r ) near the point m one obtains the formula for the quark-mass correction δm ≡ m r − m:

Pion Contribution
The contributions of pion to the quark propagator are (1) We shall use the pole approximation (40) for the calculation of these contributions.The pion contribution to first-order condensate is calculated by formula (69).The integral over dp can be calculated for DAR in closed form, and the result is the very simple expression: (See also [45], where this result has been found by a slightly different method.)For the calculation of the pion contribution 4D-cutoff regularization we use Equations ( 40) and (37).Further the pion contribution in condensate is calculated with formula (69).For 4D cutoff the result for ρ π is not described by a simple formula, as for DAR.Nevertheless, the computation is no any troubles.Note that whereas in DAR ρ π is a function of regularization parameter ξ, in 4D cutoff this quantity is a function of ratio x ≡ Λ 2 /m 2 : As examples we give values of (ρ π ) FDC for two characteristic values of this ratio.At x = 3 (this value corresponds to value c (0) = −210 MeV of the LO condensate) the computation gives (ρ π ) FDC = −0.272.At x = 19 (this value corresponds to value c (0) = −250 MeV of the LO condensate) the computation gives (ρ π ) FDC = −0.183.
To calculate the pion contribution in quark mass we apply Equation (70).As a result, we obtain: i.e., for the DAR, the pion correction to quark mass equal zero.For 4D cutoff the pion correction to quark mass is where .
Signs of (ρ π ) FDC (x) and h π (x) are opposite, and their contributions in δm are mutually cancelled.As for the DAR, the pion correction to quark mass is zero.We obtain this result, unlike to the exact result (74) of DAR, by computations in a framework of given accuracy of inputs.Such coincidence of results in both regularizations suggests an idea that the zero-pion correction to quark mass is the regularization-independent fact of NJL model.However, we cannot to proof this suggestion analytically.

Scalar Contribution
In this subsection we consider contributions of scalar amplitude A σ in the chiral quark condensate and in the mass of quark.In correspondence with (68) and (67) we have To calculate these contributions, we use the leading-singularity approximation: For the DAR, this approximation is described by Equation (43).From Equation ( 69) we obtain the quantity ρ σ .A computation gives us the following values for sigma-contribution: at ξ = 0.25 we obtain (ρ σ ) DAR = −0.033;at ξ = 0.4 we get (ρ σ ) DAR = −0.01.As one can see the contribution from sigma is small and possesses the opposite sign in comparison of the pion contribution, i.e., it decreases the common contribution.
In the case of the 4D cutoff the leading-singularity approximation for A σ coincides with the pole approximation (41).The quantity of ρ σ , as ρ π , for the 4D cutoff is a function of x ≡ Λ 2 /m 2 : At x = 3 we find (ρ σ ) FDC = −0.007.At x = 19 we obtain (ρ σ ) FDC = −0.116.In contrast to the DAR, the sign of contribution for sigma is the same as for pion contribution.
A correction from sigma to quark mass for DAR is and at ξ = 0.25 : δm DAR (σ) = −0.086m,at ξ = 0.4 : δm DAR (σ) = −0.056m.Since a pion correction to quark mass in this regularization equals zero, these values are full corrections to quark mass for DAR.
For the 4D cutoff the correction from sigma to quark-mass is where Let consider an issue on accuracy of above computations.A main approximation is the leading-singularity approximation.Although the leading singularity gives a main contribution in an integral, nevertheless let consider a part of other terms.To estimate their part for DAR let us apply the expressions of amplitudes at ξ = 1/2.A computation with formulas (67)- (69) demonstrates that the contributions of non-pole terms in chiral condensate equal zero.Since the values of parameter ξ are near this point, we can maintain that at ξ = 1/2 their contributions are also small in comparison with the main pole-pion contribution.
In the case of 4D cutoff the computations with exact formulas (19) and (20) for the amplitudes also demonstrate, that pole approximation in the case gives the main contribution in condensate.Therefore, at x = 3 the computation with the exact formulas (19) and (20) gives for the pion contribution ρ π = −0.267,i.e., differs from the result of pole approximation less than on 2%.For contribution from sigma the difference is more significant: the computation with the exact formulas gives ρ σ = −0.031,but since this contribution is much less in comparison with the pion contribution, this difference again practically does not affect to final result.

Improved Model Parameters
In this subsection we discuss the modification of the model parameters with the calculated corrections.We modify an expression for the condensate as follows: The formula for f π (see (35)) stays the same, since corrections to amplitudes generate in the next (second) order of MFE.The quark mass: where δm is defined by Equation (70).For 4D cutoff the system of Equations ( 80), ( 39) and ( 23), which determines these parameters, has no solution at f π = 93 MeV and at |c| ≤ 230 MeV.There is very important circumstance-for 4D cutoff the meson contributions can destabilize the NJL model.Though these contributions are relatively small (they do not exceed 25% from the leading contribution), but their opposite sign results in a non-stability of all the system.At that for DAR, the situation is principally different: due to the positivity of the meson contribution in condensate for this regularization a stabilization of the model takes place.The quark mass m r almost does not depend on c, and values of regularization parameter ξ increase in comparison with corresponding LO values, i.e., shift to an area of stability of model, where meson contributions decrease.
As it follows from results of this section, the NJL model with DAR is stable with respect to quantum fluctuations caused by scalar-amplitude contributions in chiral condensate, whereas for the NJL model with 4D cutoff the meson contributions can lead to destabilization.Surely, several physical applications of the NJL model are connected exclusively with the LO of MFE (mean-field approximation), and the possibility of such destabilization can be ignored.On the other hand, some physical applications of the NJL model exist that connected with multi-quark functions (see below), for which a neglecting by the meson contributions in quark propagator is certainly non-correct from the point of view of the MFE and, consequently, the stability of basic model parameters with respect to these contributions becomes essential.

The Corrections to the Two-Quark Function and the Legendre Transform
In this section, we apply the method of the Legendre transform with respect to a bilocal source [30,31] to find the corrections the two-quark function in the NJL model.The principal problem in the computation of the corrections to the two-particle amplitude is the constraints imposed by chiral symmetry.An effective way to solve this problem is the Legendre transform method.
We go from generating functional G to the logarithm Z = 1 i log G and define the quark propagator as a functional of source η Then we consider relation (81) as an equation for η = η[S] and define the generating functional of the Legendre transform (effective action) as a functional of S: From Equations ( 81) and (82) we obtain. where is two-quark function.This function as a functional of S is defined by the relation which is, in fact, the equation for the two-quark function in the formalism of the Legendre transform with respect to the bilocal source η.The kernel K of this equation is defined as the connected part of the second derivative of the generating functional of the Legendre transform and is given by the relation Going to amputated two-quark function (13) and taking its connected part (14) we obtain the equation for connected part F c 2 of the amputated two-quark function (the two-quark amplitude): One of principal things in a theory with spontaneous chiral symmetry breaking is the chiral Ward identity.This identity relates the axial vector part of the two-quark function to the propagator and has the following form: The most important consequence of this relation is that the function S 2 should have a simple pole of the form 1/P 2 in the quark-antiquark channel (here P = p x + p y is the total momentum) in the case of DCSB.This pole means that a massless pseudoscalar particle (Goldstone boson) appears in the spectrum of the theory in full accordance with the NGB theorem.
We stress that this identity is a model-independent relation for any local theory with chiral symmetry and should hold in both the NJL model and QCD-type fundamental theory.In terms of the generating functional of the Legendre transform, this identity reads τ a γ 5 i ∂z (δ(x − z)δ(z − y)) = tr u δη(z, x 1 ) δS(y, x) S(x 1 , z)γ 5 τ a + τ a γ 5 S(z, x 1 ) δη(x 1 , z) δS(y, x) dx 1 . (90) This equation indicates that to check the chiral Ward identity in the Legendre transform formalism, there is no need to calculate the two-quark function, and it suffices to check it for kernel K.The mean-field approximation for the generating functional of the Legendre transform is SDE (84) without the last two terms: Equation ( 91) with the switched-off source is the equation for the quark propagator in the leading approximation (see Equation ( 4)), whose solution is Equation (6).
By differentiation with respect to S we obtain from Equation (91): The connected part of this expression (i.e., the second term in the r.h.s.) is kernel K in the LO of the MFE for the NJL model.With the direct calculation, we can easily see that chiral Ward identity (90) holds in the leading approximation.Prior the computation the two-quark amplitude, we can thus be sure that the mean-field approximation in the Legendre transform formalism satisfies to the major symmetry requirement of the NJL model in the chiral limit, namely the NGB theorem.
The equation for the generating functional of the Legendre transform with corrections has the form, η(x, y) = S −1 (x, y) + i ∂δ(x − y) + igδ(x − y){tr[S(x, y)] + iγ 5 τ a tr[iγ 5 τ a S(x, y)]} + 1 i where and F 0 is the LO two-quark function.Using the equation for F 0 , we can obtain the more compact expression 1 i The differentiation of Equation ( 95) with respect to the functional variable S yields the correction to the kernel of the Bethe-Salpeter equation (BSE) for the two-particle function: The result is This equation together with (96) yields the correction to the kernel of the BSE for the two-particle function.As can be seen, this correction consists of the one-meson and two-meson contributions.
The principal problem of the computation of corrections to the two-quark amplitude is the requirement that these corrections correspond to the NGB theorem.The chiral Ward identity holds true in the determinations of these corrections by the Legendre transform method.Therefore, this method ensures the validity of the NGB theorem.Checking the chiral Ward identity for the kernel of the BSE including corrections (96) and ( 97) is a less trivial procedure than similarly verifying the leading approximation.Nevertheless, the result is positive for this case as well.The obtained result shows that the considered approximation is physically reasonable, i.e., it is a symmetry-preserving approximation for calculating the corrections to the two-particle function (for the more details of this computation see [44]).

The Equations for Multi-Quark Functions
In this section, we will focus on multi-quark functions in NJL model and discuss their field-theoretical status.The multi-quark functions give us a fundamental opportunity to significantly expand the field of application of the model.Therefore, the three-quark function describes the decay σ → ππ.The solution of the equation for three-particle function permits also to close the equation for NLO two-particle function and to calculate a correction to pion-decay constant.The four-quark function S (1) 4 gives us possibility to describe the pion-pion scattering in quark fields context.

Four-Quark and Three-Quark Functions
The generating functional of second step is The equations for four-quark and three-quark functions see on Figures 3 and 4.
From the second step of the iteration scheme we obtain four-quark function S 4 , three-quark function S 3 , and also NLO two-quark function S 2 and NNLO quark propagator S (2) .For these functions we have a system of four integral equations.These equations and all equations of following steps of the iteration scheme possess the structure, which is similar to the structure of Equation ( 11) (We shall use operator notations for the bulky formulas of this section.For full content see [66][67][68]): and differ from each other by the structure of inhomogeneous terms S 0 n .The inhomogeneous term in the equation for four-quark function S 4 is where S 2 is defined by Equation ( 11), i.e., the equation for S 4 does not include another three second-step functions.The inhomogeneous term in the equation for three-quark function S 3 includes four-quark function S 4 , two-quark function S 2 , NLO quark propagator S (1) and LO quark propagator S.This inhomogeneous term is as follows Inhomogeneous terms of other second-step equations ( (S 2 ) 0 and (S (2) ) 0 ) see in Appendix A of work [68].
The solution of the system should be started with the equation for four-quark function S 4 , then should be solved the equation for three-quark function S 3 , etc.
For four-quark correlation function S 4 with inhomogeneous term (99) we obtain simple disconnected solution: The graphical representation for the solution of equation for four-quark function see on Figure 5.The disconnection means the absence of physical effects due to four-particle functions in the given order of MFE, and the pion-pion scattering cannot describe in the given order.It will be appearing only in next step of iteration scheme (third order of MFE).
The expression for the four-quark function gives us the closed equation of type (98) for three-quark function S 3 .The inhomogeneous term of this equation is The solution can be obtained likewise to solving of Equation ( 11) for the two-quark function.We can go to the amputated function    With these definitions we obtain for vertex function W ab the following equation: where l S (x) ≡ tr[S(x)S(−x)] is the scalar quark loop.Inhomogeneous term W ab 0 is: −ig dx 1 (V σ (xx 1 x 1 ) − 4V a 1 a 1 π (xx 1 x 1 ))S ab π (x − x ).
Using above definitions, we have: +∆ a 1 b π (x 2 − x )∆ a 1 a π (x 1 − x )].The Equation (104) can be easy solved in the momentum space and its solution is (See also Figure 6).The above calculation is, of course, purely illustrative in view of the very uncertain experimental status of the sigma meson.In the NLO formula (107) defines correction ( f (1) π ) 2 to expression (53).Surely, to calculate this correction, it is no need in a complete solution of the NLO two-particle equation.In correspondence with (53) it is quite enough to calculate the pion-pole part only.
The presence of the new diquark-source results the connection condition for derivatives of generating functional: δ 2 G δ ῡ(x 2 , x 1 )δη(y, x) = − δ 2 G δ ῡ(x 1 , x)δη(y, x 2 ) . (115) Consequently, SDE (114) can be rewritten in the alternative forms.Being fully equivalent from the point of view of an exact solution of SDE's, these alternative forms, can lead to different approximations in the MFE.The choice of the suitable forms for the construction of MFE in the case should be made following of corresponding physical reasons.
In a similar way one can to enter into consideration three-quark sources.Such sources can be useful for the direct description of nucleons and other baryons omitting the intermediate diquark modelling.
For these systems of equations one can to elaborate the MFE by the method, which is similar to that of Section 2. An elaboration of this construction with solutions of corresponding equations is the object of future investigations.

Results and Conclusions
In the present work we have used an iteration scheme of solution of the SDE with the fermion bilocal source to formulate the MFE for NJL model (Section 2).The equations of any order in this iterative scheme have a simple analytic structure, and solving the equations of any order is actually an algebraic problem.
According to the results obtained in following sections, the NJL model with DAR gave us simple closed formulas not only for the scalar amplitudes and the pion-decay constant, but also for the pion contribution to the chiral condensate.As it follows from our results, in the NJL model with DAR, this pion contribution is significant and should be taken into account with a choice of physical values of the model parameters.
Obtained results demonstrated that the NJL model with DAR essentially differs from the NJL model with 4D cutoff at least in two problems.

Figure 1 .
Figure 1.Graphical representation of equation for two-quark function.

) 3 + F three−meson 3 ,
Then we should separate the functions with the same algebraic structures, which are reproduced by iterations.The connected part of the amputated three-quark function possesses two-meson and three-meson contributions: F c 3 = F two−meson where F two−meson 3 is a quadratic functional of A σ and A π , and F three−meson 3 is a functional of third degree.

Figure 3 .
Figure 3. Graphical representation of equation for four-quark function.

7. 3 .ψγ µ γ 5 τ b 2 ψ 2 .
NLO Two-Particle Function and Correction to Pion-Decay ConstantThe solution of the equation for three-particle function permits to close the equation for NLO two-particle function which enables calculation of a correction to pion-decay constant f π .The pion-decay constant is defined asi f π δ bb P µ =< 0|J b µ5 |P, b >,where |P, b > is the pion state with momentum P and isospin b, and J b µ5 = is the axial current.If the two-particle function has a pole term S pole 2 , which corresponds to the pion, then, taking into account these definitions, this pole term is connected with the pion-decay constant by relation (2π) 4 δ(P − P ) f 2 π = i dxdx e i(Px−P x ) tr u,d [γ µ γ 5 tr u,d denotes the traces over up and down lines of function S pole In the LO we obtain from (107) the well-known expression for the pion-decay constant ( f (0) π 2 (see Equation (53)) ab σππ (pp p ) = i∆ σ (p)W ab 0 (pp p ), (105) where p is σ−mesons momentum, and p , p are pion momenta: p = p + p .The connected part of W ab is the decay amplitude σ → ππ.It has the following form: [W ab σππ (pp p )] con = 2N c i ∆ σ (p)[v P (pp p ) + v P (pp p )]∆ aa 1 π (p )∆ a 1 b π (p ).