Functional Integral Approaches to the Bosonization of Effective Multi-Quark Interactions with $U_\A(1)$ Breaking

Low energy hadron phenomenology involving the (u,d,s) quarks is often approached through effective multi-quark Lagrangians with the symmetries of QCD. A very successful approach consists in taking the four-quark Nambu--Jona-Lasinio Lagrangian with the chiral $U_L(3)\times U_R(3)$ symmetry in the massless limit, combined with the $U_\A(1)$ breaking six-quark flavour determinant interaction of 't Hooft. We review the present status and some very recent developments related to the functional integration over the cubic term in auxiliary mesonic variables that one introduces to bosonize the system. Various approaches for handling this functional, which cannot be integrated exactly, are discussed: the stationary phase approximation, the perturbative expansion, the loop expansion, their interrelation and importance for the evaluation of the effective action. The intricate group structure rules out the method of Airy's integral. The problem of the instability of the vacuum is stated and a solution given by including eight-quark interactions.


Introduction
Non-perturbative methods to deal with the low energy regime of QCD involve lattice calculations, QCD sum rules [1], effective field theories, which are renormalizable in the modern sense and gave birth to Chiral Perturbation Theory [2,3], large N c (number of colours) expansions [4,5,6,7], and effective models with solid symmetry content [8]. The object of our studies is placed in the latter category, consisting in a model Lagrangian which combines the four-quark interactions of Nambu and Jona-Lasinio (NJL) [9], extended to the U L (3)×U R (3) chiral symmetry of massless QCD [10,11,12], and the six-quark interaction of 't Hooft [13], which breaks the unwanted U A (1) symmetry of QCD. We call it the NJLH Lagrangian. The 2N f (N f denoting the number of flavours) multi-quark interactions of 't Hooft were derived after eliminating the gluonic degrees of freedom of QCD in the semiclassical instanton approximation, assuming dominance of the zero modes in the scattering process of quarks and anti-quarks on the same instanton. Recent reviews providing further evidence on multi-quark vertices can be found in [14].
In the last decade some new developments brought further relevance to the NJLH model, as it may help to shed some light on the possible existence of an hierarchy of multi-quark interactions, with prevalence of lower ones. There is some conflict on two fronts: (i) Going beyond the zero mode approximation, an infinite number of multi-quark interactions, starting from the four-quark ones, has been found in the instanton-gas model, all of them being of the same importance [15]. (ii) On the other hand, accurate lattice measurements for the realistic QCD vacuum show a hierarchy between the gluon field correlators with a dominance of the lowest one [16]. A similar hierarchy of the multi-quark interactions can be triggered after averaging over gluon fields. If this is true, there is an apparent contradiction with the instanton-gas model (see also [17]). This point has been stressed in [18].
The hierarchy problem of multi-quark interactions can be addressed on pure phenomenological grounds. For that we suggest to simplify the task considering only the four-and sixquark interactions of the NJLH model. Once one knows the Lagrangian the obvious question arises: does the system possess a stable vacuum state and does this state correspond to our phenomenological expectations? If hierarchy takes place this question is pertinent for the leading four-quark interaction, because in this case the effective quark Lagrangian can be studied step by step in the hierarchy with the assumption that four-quark vertices are the most important ones. In the opposite case we must study the system as a whole to answer the question. For the best known and simplest example, to which we dedicate most of our attention here, the solution may be found analytically. As a result one can obtain definite answers to the above questions with a convincing indication in favour of a hierarchy for the considered example.
Pioneering work on the original NJLH Lagrangian has been done in [19], using Bethe-Salpeter techniques to evaluate quark -anti-quark scattering in one-loop approximation and obtain the meson masses from the poles of the scattering matrix, and in [20], where the quark Lagrangian has been bosonized using functional integral methods and the mesonic spectrum calculated at the leading order of the stationary phase approach (SPA) for one critical point. A good description of the pseudoscalar nonet, especially the η and η ′ masses and mixing has been achieved. Within the respective approximations, the effective potentials for both methods coincide [21]. The model has subsequently been widely and successfully explored at the mean field level, see e.g. [22,23,24].
In the present work we present a detailed picture and mathematically careful study of the functional integration over the auxiliary bosonic fields of cubic order, inherited from the 't Hooft six-quark interaction. They cannot be integrated out exactly. We came across one startling fact and several interesting properties, which we hope help in the understanding of the approximations used and in the identification of the shortcomings and potentialities of multi-quark models. Before turning to the calculations, let us summarize what we have learned from them [25,26]: (1) We show that the stationary phase equations have several roots (critical points), of which one is regular and the others singular. A rigorous SPA treatment, taking into account all critical points, leads to an unstable vacuum for the theory. This is the startling point, since this fatal flaw of the model is in crass contradiction with its phenomenological success.
(2) The results in [20] are obtained if one considers only the regular critical point of the SPA, which of course is accompanied with a different asymptotic behaviour.
(3) This result can be understood in the perturbative sense, where the cubic interaction is taken as a perturbation around the stable NJL vacuum. The situation is in many respects analogous to the problem of a harmonic oscillator perturbed by an x 3 term. This system has no ground state, but perturbation theory around a local minimum does not know this.
(4) By performing a resummation of the perturbative series in powers of the 't Hooft interaction one gets a loop expansion [27], of which the leading term contains all tree graphs present in all powers of the perturbative series. This corresponds exactly to the result of [20]. (5) The loop expansion can be classified in terms which contribute to the measure (all odd numbers of loops) and to the phase, relevant for the effective action (all even number of loops). (6) The global instability of the vacuum can be removed within the multi-quark interaction picture: the addition of a chiral invariant OZI-violating [28] eight-quark interaction to the Lagrangian renders it globally stable, for certain values of the coupling strengths of the several interactions. (7) A comparison between the effective potentials for the NJLH and for the eight-quark enlarged Lagrangian reveals that the global stabilization happens around the local minimum which appears if only the regular critical point is considerd in the NJLH Lagrangian. This is probably at the heart of its success.
Points (6) and (7) will not be discussed in this contribution. In fact, at the time the conference took place, they had not been obtained yet. We added them to the list, as they represent a kind of "happy-end" for the model. We refer to [26] for a detailed presentation of this issue.
The paper is organized as follows. In Section 2 we write out the NJLH Lagrangian and use the functional integral representation to express it in terms of bosonic fields. In Section 3 we show that the chiral symmetry group imposes constraints which are only compatible either with the perturbative approach, or the expansion in a parameter that multiplies the total Lagrangian density (the loop expansion). Otherwise, as the consistent stationary phase treatment shows, the model is unstable. In Section 4.1 we discuss the perturbative treatment; the integration over auxiliary variables leads to a special problem with δ(0) -singularities, to which we give a physical meaning through Feynman diagrams in 4.2. The loop expansion is considered in Section 5. We obtain in closed form the two-loop contributions to the functional integral Z and give arguments to justify this result. Since we have forestalled the conclusions in the Introduction, we shall omit concluding remarks.

The Lagrangian and bosonization
On lines suggested by multicolour chromodynamics it can be argued [29] that the U A (1) anomaly vanishes in the large N c limit, so that mesons come degenerate in mass nonets. Hence the leading order (in N c counting) mesonic Lagrangian and the corresponding underlying quark Lagrangian must inherit the U L (3) × U R (3) chiral symmetry of massless QCD. In accordance with these expectations the U L (3) × U R (3) symmetric NJL interactions, can be used to specify the corresponding local part of the effective quark Lagrangian in channels where the matrices P L,R = (1 ∓ γ 5 )/2 are projectors and the determinant is over flavour indices. The coupling constant κ is a dimensional parameter ([κ] = GeV −5 ) with the large N c asymptotics κ ∼ 1/N N f c . The coupling G, [G] = GeV −2 , counts as G ∼ 1/N c and, therefore, the Lagrangian (1) dominates over L H at large N c . It differs from the counting G ∼ 1/N 2 c , which one obtains in the instanton-gas vacuum [15].
It is assumed here for simplicity that interactions between quarks can be taken in the long wavelength limit where they are effectively local. The 't Hooft-type ansatz (2) is a frequently used approximation. Even in this essentially simplified form the determinantal interaction has all basic ingredients to describe the dynamical symmetry breaking of the hadronic vacuum and explicitly breaks the axial U A (1) symmetry [30,31]. The effective mesonic Lagrangian in leading order SPA, corresponding to the non-local determinantal interaction, has been found in [32,33].
Anticipating our result, we would like to note that if the hierarchy of multi-quark interactions really occurs in nature, the perturbative treatment seems adequate. The NJL interaction alone has a stable vacuum state corresponding to spontaneously broken chiral symmetry. But, as we shall show, the effective quark theory based on the Lagrangian 1 has a fatal flaw: if L H is comparable with L NJL , it has no stable ground state. In equation (3) the current quark mass,m, stands for the diagonal matrix diag (m u ,m d ,m s ), which explicitly breaks the global chiral SU L (3) × SU R (3) symmetry of the Lagrangian.
To bosonize the theory one introduces auxiliary bosonic variables to render fermionic vertices bilinear in the quark fields. This procedure requires twice more bosonic degrees of freedom than necessary [20]. Redundant variables must be integrated out and this integration is problematic as soon as one goes beyond the lowest order stationary phase approximation [34,21]: the lowest order result is simply the value of the integrand taken at one definite stationary point [30]. We shall show how to extract systematically the higher order corrections which contribute to the effective mesonic Lagrangian, and what to do with infinities contained in these corrections.
The many-fermion vertices of Lagrangian (3) can be linearized by introducing the functional unity [20] in the vacuum-to-vacuum amplitude We consider the theory of quark fields in four-dimensional Minkowski space. It is assumed that the quark fields have colour (N c = 3) and flavour (N f = 3). The auxiliary bosonic fields, σ a , and, φ a , (a = 0, 1, . . . , 8) become the composite scalar and pseudoscalar mesons and the auxiliary fields, s a , and, p a , must be integrated out. By means of the simple trick (4), it is easy to write down the amplitude (5) as We assume here that σ = σ a λ a , and so on for all auxiliary fields σ, φ, s, p.
The totally symmetric constants A abc are related to the flavour determinant, and equal to We use the standard definitions for antisymmetric f abc and symmetric d abc structure constants of U (3) flavour symmetry. One can find, for instance, the following useful relations At this stage it is easy to rewrite equation (6), by changing the order of integrations, in a form appropriate to accomplish the bosonization, i.e., to calculate the integrals over quark fields and integrate out from Z the unphysical part associated with the auxiliary bosonic variables (s a , p a ) where The Fermi fields enter the action bilinearly, thus one can always integrate over them, since one deals with a Gaussian integral. One should also shift the scalar fields σ a (x) → σ a (x) + m a by demanding that the vacuum expectation values of the shifted fields vanish 0|σ a (x)|0 = 0. In other words, all tadpole graphs in the end should sum to zero, giving us the gap equation to fix the constituent quark masses m a corresponding to the physical vacuum state. The functional integrals over s a and p a are the main subject of our study. We put here ∆ a = m a −m a , and N is chosen so that Z[0, 0; ∆] = 1. Let us join the auxiliary bosonic variables in one 18-component object R A = (R a , Rȧ) where we identify R a ≡ s a and Rȧ ≡ p a ; a,ȧ run from 0 to 8 independently. It is clear then, that R 2 A = s 2 a + p 2 a . Analogously, we will use Π A = (σ a , φ a ) for external fields and ∆ A = (∆ a , 0).
we find after some algebra with the following important property to be fulfilled Now it is easy to see that the functional integral (8) can be written in a compact way where We have arrived at a functional integral with a cubic polynomial in the exponent.

The integration over the auxiliary variables R A
Given the cubic structure in the functional integral, one might be tempted to solve it using Airy's integral methods. Invoking the existence of a large expansion parameter, such as N c , one uses the well known asymptotics of the Airy's function on the real axis. This would require [25] which cannot be fulfilled. On the contrary, the system of equations based on the first order derivatives is self-consistent and can be solved [21]. Therefore, one can obtain the semi-classical asymptotics through the stationary phase method.

Solving equation (13)
We need to recall shortly the solutions of equation (13). Up to some order in the external mesonic fields, Π A , we may write them as a polynomial and so on. Putting this expansion in equation (13) One can always find the trivial solution H A = 0, corresponding to the unbroken vacuum ∆ A = 0. There are also non-trivial ones for the scalar component, i.e., The number of possible solutions, i, depends on the symmetry group. The coefficients h This system is equivalent to a fifth order equation for a one-type variable which can be solved numerically. For two particular cases, whenm u =m d =m s andm u =m d =m s , equations (15) can be solved analytically [21]. The simplest example:m u =m d =m s (or, equivalently, h s ) corresponds to SU (3) flavour symmetry. In this case equation (15) has two solutions If 4G 2 > κ∆, they are real and will contribute to the stationary phase trajectory.

The lowest order semiclassical asymptotics and instability
Since the system of equations (13) can be solved, we may replace variables in the functional integral (10) to obtain the semi-classical asymptotics where n is the number of solutions, R A , of equation (13), L ′′ AB has been defined in equation (12), and Here we used equation (13) to eliminate the term proportional to κ. Let us stress that L a . This dependence is singular at κ → 0. One can see this, for instance, from equation (16) where h (2) u → ∞ for small κ. This behaviour reflects the fact that we are far from the perturbative regime, meaning that the interactions L NJL and L H are equally weighted.
The linear term in the σ field is written explicitly. This part of the Lagrangian is responsible for the dynamical symmetry breaking in the multi-quark system and taken together with the corresponding part from the Gaussian integration over quark fields in equation (7) leads us to the gap equation.
At leading order, k = 0, and for the terms linear in the fields contributing to the phase of Z we have the estimate for equation (17) where A (i) is real and proportional to 2 Therefore, if one considers the case withm u =m d =m s , Z is given by Let us recall that the quark loop contribution to the gap equation is well known (see, for instance, [35]). Combining this known result with the estimate (18), one can obtain the corresponding effective potential U (m) as a function of the constituent quark mass m where we consider the casem = 0 for simplicity and Λ q denotes the cutoff of quark loop integrals. This system has at most a metastable vacuum state (for G/κ > 0); for G/κ < 0, the effective potential does not have extrema in the region m > 0. We must conclude that the model considered has a fatal flaw and can be used only in the framework of the perturbative approach, which does assume the hierarchy of multi-quark interactions.

Perturbative expansion of Z
We shall restrict ourselves in this section to the perturbative treatment of the functional integral (10). The loop expansion will be considered in the next section.

The perturbative series
Let us divide the Lagrangian (11) in two parts. The free part, L 0 , is given by The 't Hooft interaction is considered as a perturbation L I Thus the perturbative representation for the functional integral (10) can be written as whereX Since the boson fields appear quadratically in equation (20), they may be integrated out, yielding whereΠ A = Π A + ∆ A . The overall factor N = (−2πi/G) 9 N ′ is unimportant in the following. We want to calculate the effective action Γ eff = d 4 xL eff , which by definition is the phase of Z and A(Π A ) is a real function. Comparing (21) and (22), one gets Here Γ 0 represents the leading order result for Γ eff while the second logarithm in equation (23) is a source of U A (1) breaking corrections which arise as a series in powers of the functional derivatives operator To make this statement more explicit let us consider the expansion Taking into account the symmetry properties of the coefficients Φ ABC and our previous result (9), we find represents the effective action (23) as a perturbative series in powers of κ. For instance, up to and including the second order in κ we have where The real factor A(Π) is always chosen such as to cancel the imaginary part of the effective action. To the approximation considered we have, for instance, It contributes to the measure of the functional integral over σ a , φ a .

Giving physical meaning to the singularities
The terms with the δ(0) function require further explanation. The fact that auxiliary fields can lead to special problems with infinities is well-known [36,37,38,39]. To see in our case the origin of the encountered singularities and endow them with physical meaning we shall use the language of Feynman diagrams. The graphs contributing to lowest order in κ are shown in Fig. 1. Figure 1. The lowest order graphs contributing to δ 1 .

+
In these diagrams, a line segment (straight or curved) stands for a "propagator" extracted from the last exponent in equation (21). A filled circle at one end of a line segment corresponds to the external field, i d 4 xΠ A (x), and a vertex joining three line segments is used for iκΦ ABC d 4 x.
The first diagram represents the δ 1 term in equation (26). The contribution of the second tadpole diagram is equal to zero. Indeed, the vertex contains the group factor Φ ABC . The contraction of any two indices in this factor by δ AB from the propagator (situation occurring for the tadpole graph) reduces it to zero, according to equation (9). We thus find that tadpole diagrams do not contribute due to the flavour structure of the 't Hooft interaction.
To next to leading order in κ we have the four graphs shown in Fig. 2. Figure 2. The diagrams of the κ 2 order contributing to δ 2 .
As we just learned, the third diagram does not contribute. The other ones correspond exactly to the three terms of equation (27). The first tree diagram is finite. The second one-loop diagram has a divergent factor δ(0). The last two-loop diagram contributes as δ(0) 2 . These singularities were caused by the local structure of the last exponent in equation (21) or, which is the same, by the δ(x − y) term in the propagator (28). We believe that if one would start from non-local NJL interactions, the singularities could be weaker or would even disappear.
The factor δ(0) requires a regularization. This is an expected trouble in the NJL model which is nonrenormalizable and, as a consequence, the fundamental interactions must be cut off. The cutoff is an effective, if crude, implementation of the known short distance behaviour of QCD within the model. The problem with δ(0) singularities can also be analysed using the spectral representation method to evaluate the integral (8) [25]. For a one-dimensional field theory it takes the form where L refers to the box size in which the system is put. Therefore the second term does not contribute in the limit L → ∞. The cutoff Λ cuts the density of Fourier harmonics related to the auxiliary bosonic variables, its finite value must be fixed by confronting the model with experiment; it may differ from the value of Λ q in equation (19), which is associated with the integration over the fermionic degrees of freedom. In this interpretation the model is effectively finite, including the higher order corrections 3 . The diagrams are a very convenient language to understand another feature related to the multi-loop contributions: the phase factor corresponding to the diagram can be simply calculated. Indeed, it is easy to see that for any diagram the formula is fulfilled. Here E is the number of external fields, I stands for the number of internal lines, and V is the number of vertices. On the other hand, the number of loops, L, is given by This is because the number of loops in a diagram is equal to the number of δ functions surviving after all integrations over coordinates are performed, except for one over-all integration related to the effective action. Every internal line contributes one δ function, but every vertex or external field carries an integration over the corresponding coordinate, and thus reduces the number of δ functions by one.
This result shows that the overall phase factor of a given graph i I−E−V = i L−1 is entirely determined by the number of loops. In particular, diagrams with an even number of loops contribute to the effective action (the argument of Z), while the diagrams with odd number of loops contribute to the measure (the modulus of Z), see equation (22).

The loop expansion of Z
The perturbative series, considered in the previous section, can be resummed. This is done by introducing a parameter t in the perturbative representation (20) of the functional integral according to the substitution: Then equation (25) to second order in κ 2 becomes The groups in terms of inverse powers of t of Γ ′ eff (t) define the new series: first all diagrams with no closed loops (tree graphs) will carry the weight 1. The tree graphs yield the result of [20]. All the one-loop diagrams (and higher odd numbers of loop diagrams) contribute to the imaginary part of the action, as discussed in the previous section; as mentioned before, they will be cancelled by the appropriate choice of the real quantity 1 t A(Π A ), with A being of order 1. They contribute therefore to the measure of the integral over Π A and were obtained in [21], see also equations (31), (32) below. The set of graphs with 2n loops have the factor t −2n and contribute to the effective action; in general they include, as a subset, all graphs of k n+1 th order or higher in this coupling constant. The resummed series is equivalent to the well-known loop expansion [27]. The latter can be written formally as a SPA integral which includes only the regular critical point (SPAr), and constitutes the most economic way to obtain the loop corrections. A detailed derivation of the equivalence of the resummation of the perturbative series, the loop expansion and the SPAr is given in Section 3 of [25] for a one-dimensional analogue of the considered functional.
We proceed therefore to obtain the NLO to the effective action through the SPAr method. Consider the effective mesonic action generated by the functional which in comparison with equation (17) has only one critical point, related to the stable configuration, which is the solution with i = 1 in equations (16). We shall identify L (i=1) st = L st , and L ′′ AB (R (i=1) ) = L ′′ AB in the following. By replacing the continuum of spacetime positions with a discrete lattice of points surrounded by separate regions of small spacetime volume Ω, the functional integral (29) may be reexpressed as a Gaussian multiple integral over a finite number of real variables R A (x) for a fixed spacetime point x. We think of DR A as the infinite product After evaluating the Gaussian integrals one obtains Here M = tr L ′′−1 3 + 6 tr L ′′−1 tr L ′′−1 2 + 8 tr L ′′−1 3 , contains the complete information on the two-loop term and the dots mean the terms corresponding to the three-loop contribution and higher. The I 0 is the one-loop contribution where λ j are eigenvalues of the N × N matrix L ′′ AB . In our case N = 18, related to the 18-component object R A , equation (12), and L ′′−1 denotes the inverse matrix of L ′′ . The totally symmetric symbol δ ABCDEF generalizes an ordinary Kronecker delta symbol δ AB by the recurrent relation Up to the given accuracy we have in equation (30) x with the quantity I 0 containing the one-loop corrections to the measure. It shows that in the continuum limit the two-loop correction contributes to the effective Lagrangian as This is our final expression for the effective Lagrangian in the two-loop approximation.
The field-dependent factor M contains all possible mesonic vertices, including the σ-tadpole contribution to the gap equation, contributions to the masses of scalar and pseudoscalar nonets, as well as interaction terms. Its dependence on the parameters κ, G, m,m enter through the expansion coefficients H  (14), (15). Let us do some estimates to justify the result. For this purpose let us simplify the integral (29). After neglecting the symmetry group and discretizing the spacetime it takes the form To justify the stationary phase approximation for the integral (34) we assume that The dominating role of the Gaussian integral is reflected in the fact that essential values for R x in the integral have the order R 2 x ∼ 1/(ΩL ′′ st ). For the cubic term it follows then where we have used that in the model considered here, L ′′′ st ∼ κ, L ′′ st ∼ G. If the parameters of the model can be chosen in such a way that the inequality ζ ≪ 1 is fulfilled, the cubic power of R x yields terms that go to zero relative to the Gaussian term as ζ → 0, and the stationary phase approximation will be justified. Note that Ω may be written as an ultraviolet divergent integral regularized by introducing a cutoff Λ Therefore, the inequality restricts the value of Λ from above.
Meanwhile, it is interesting to see that the inequality ζ ≪ 1 is an exact equivalent of (35). Indeed, in the essential region, i.e., around a sharp minimum, one has L ′ st ∼ R x L st , L ′′ st ∼ R 2 x L st , and so on, thus Therefore, the asymptotical series (30) with the ultraviolet cutoff imposed in the continuum limit is sensible. One deals here, actually, with a series in powers of the dimensionless parameter ζ. The expansion is formally justified for ζ ≪ 1.

The NLO effective potential
As an example let us obtain the NLO corrections due to (33) to the effective potential for the SU (3) limit and with the current quark masses set to zero. Here only the regular critical point was taken into account. For the case of the full SPA, one has to consider the same expression also for the singular point. The term linear in the scalar fields coming from the two-loop correction is calculated with the help of "Mathematica" [40] and it must be added to the leading order term equation (18). Here the coefficient and ω i = ω u,i = ω d,i = ω s,i = kh i 16G with the index i denoting the two possible solutions for h i to the SPA equations in the SU (3) case, equation (16). One must integrate c i with respect to the quark mass to obtain the correction to the effective potential, since where for the last equality the equations (15) were used to obtain Here m = m u = m d = m s and m(i) = m for both solutions. With x i = ω 2 i one gets As a result we obtain Therefore the complete expression for the SU (3) effective potential is There will be poles at x i = { 1 4 , 1}. In terms of m there can be in principle up to 4 non-degenerate poles altogether, two for each critical point. They are due to the zero eigenvalues of L ′′ , which appear for certain values of the coupling constants [21] in subsets of the 18-field configurations. In the neighborhood of these points the SPA method is unreliable. This puts constraints on the coupling constants, which must be chosen such that the poles appear far from the region of physical interest (the quark mass value at the local minimum of the effective potential). One sees however that the solution h 2 , related with the singular critical point, has always a pole at m = 0, since for this case x 2 = 1. It is one more indication that the contribution of the singular point must be excluded from the analysis. Only the perturbative regime is sensible. In this case the effective potential reads GeV, Λ q = 1.64 GeV one has the two poles at m ∼ 1. GeV amd m ∼ 2.8 GeV. Both sets of parameters lead to good meson mass spectra and weak decay constants in the realistic case with SU (2) × U A (1) symmetry [41]. In conclusion, from the point of view of the complete SPA (inclusion of regular and singular points), there is no improvement, because of the bad pole at zero. From the point of view of the perturbative (regular) solution, the poles are far away on the positive m axis and the present approach can be used to study systematically NLO corrections.