Confinement Mechanism in the Field Correlator Method

Confinement in QCD is caused by vacuum fluctuations of gluon fields. There are two numerically different scales, characterizing nonperturbative QCD vacuum dynamics: a"small"scale, corresponding to gluon condensate, critical temperature etc, which is about 0.2-0.3 GeV, and a"large"one, given by inverse confining string width, glueball and gluelump masses etc, which is about 1.5-2.5 GeV. We discuss this mismatch in a picture where confinement is ensured by two-point gauge-invariant field strength correlator of special type. It supports most of perturbative and nonperturbative dynamics, while contribution of other terms is argued to be small. This object, on the other hand, can be expressed via gluelump Green's function, whose dynamics is defined in terms of the same correlator. In this way one obtains a self-consistent scheme of mean-field type leading to permanent confinement at temperatures less than the critical one.


Introduction
Confinement of color is the most important property of quantum chromodynamics (QCD), ensuring stability of matter in the Universe. Attempts to understand physical mechanism of confinement are incessant since advent of constituent quark model and QCD, for reviews see e.g. [1,2,3,4]. Among many suggestions one can distinguish three major approaches: • confinement is due to classical field lumps like instantons or dyons • confinement can be understood as a kind of Abelian-like phenomenon according to seminal suggestion by 't Hooft and Mandelstam [5] of dual Meissner scenario. This approach has gained much popularity in the lattice community and studies of various projected objects such as Abelian monopoles and center vortices are going on (see, e.g. [6]).
• confinement results from properties of quantum stochastic ensemble of nonperturbative (np) fields filling QCD vacuum. The development of this approach in a systematic way started in 1987 [7] (see, e.g. [8] as example of earlier investigations concerning stochasticity of QCD vacuum). The nontrivial structure of np vacuum can be described by a set of nonlocal gauge-invariant field strength correlators (FC). We discuss this scenario in the present paper.
QCD sum rules [9,10] were suggested as an independent approach to npQCD dynamics, not addressing directly confinement mechanism. The key role is played by gluon condensate G 2 , which is defined as np average of the following type: In the sum rule framework it is a universal quantity characterizing QCD vacuum as it is, while it enters power expansion of current-current correlators with channel-dependent coefficients. This is the cornerstone of QCD sum rules ideology. It is worth stressing from the very beginning that G 2 and other condensates are usually considered as finite physical quantities. Naively in perturbation theory one would get G 2 ∼ a −4 , where a is a space-time ultraviolet cutoff (e.g. lattice spacing). It is always assumed that this "hard" contribution is somehow subtracted from (1) and the remaining finite quantity results from "soft" np fields, in some analogy with Casimir effect where modification of the vacuum by boundaries of typical size L yields nonzero shift of energy-momentum tensor by the amount of order L −4 . In QCD the role of L is played by dynamical scale Λ −1 QCD . The idea of np gluon condensate has proved to be very fruitful. However the relation of G 2 to confinement is rather tricky. Let us stress that since there can be no local gauge-invariant order parameter for confinementdeconfinement transition, G 2 is not an order parameter and, in particular, neither perturbative nor np contributions to this quantity vanish in either phase. 1 On the other hand one can show (see details in [3]) that the scale of deconfinement temperature is set by the condensate, or, more precisely, by its electric part: T c ∼ G 1/4 2 . Let us remind the original estimate for the condensate [9] given by G 2 = 0.012 GeV 4 , with large uncertainties. One clearly sees some tension between this "small" scale and a "large" mass scale in QCD given by, e.g. the mass of the lightest 0 ++ glueball, which is about 1.5 GeV.
The question about the origin of this mismatch turned out to be rather deep. It was suggested in [7], see also [11,12] that there is another important dimensionfull parameter characterizing np dynamics of vacuum fields: correlation length λ (also denoted as T g in some papers). The crucial feature of stochastic picture found in [13] is that the lowest, quadratic, nonlocal FC, F µν (x)F λσ (y) describes all np dynamics with very good accuracy. It was also shown that a simple exponential form of quadratic correlators found on the lattice [14] allows to calculate all properties of lowest mesons, glueballs, hybrids and baryons, including Regge trajectories, lepton widths etc, see reviews [2,12] and references therein. Since the correlation length λ entering these exponents, is small, potential relativistic quantum-mechanical picture is applicable and all QCD spectrum is defined mostly by string tension σ (which is an integral characteristic of the nonlocal correlator, see below) and not by its exact profile. This is discussed in details in the next section.
However to establish the confinement mechanism unambiguously, one should be able to calculate vacuum field distributions, i.e. field correlators, self-consistently. In the long run it means that one is to demonstrate that it is essential property of QCD vacuum fields ensemble to be characterized by correlators, which support confinement for temperatures T below some critical value T c , and deconfinement at T > T c .
Attempts to achieve this goal have been undertaken in [15], however the resulting chain of equations is too complicated to use in practice.
Another step in this direction is done in [16], where FC are calculated via gluelump Green's functions and self-consistency of this procedure was demonstrated for the first time. These results were further studied and confirmed in [17].
The main aim of this paper is to present this set of equations as a selfconsistent mechanism of confinement and to clarify qualitative details of the FC -gluelump connection. In particular we demonstrate how the equivalence of colormagnetic and colorelectric FC for T = 0 ensures the Gromes relation [18]. We also show, that self-consistency condition for FC as gluelump Green's function allows to connect the mass scale to α s , and in this way to express Λ QCD via string tension σ.
Lattice computations play important role as a source of independent knowledge about vacuum field distributions. Recently a consistency check of this picture was done on the lattice [19] by measurement of spin-dependent potentials. The resulting FC are calculated and compared with gluelump predictions in [20] demonstrating good agreement with analytic results, in particular small vacuum correlation length λ ∼ 0.1 Fm is shown to correspond to large gluelump mass M 0 ≈ 2 GeV.
The paper is organized as follows. In the next section we remind the basic expressions for the Green's functions of qq and gg systems in terms of Wegner-Wilson loops and qualitative picture of FC dynamics is discussed. Section 3 is devoted to the expressions of spin-dependent potentials in terms of FC. We also shortly discuss the lattice computations of FC. In section 4 FC are expressed in terms of gluelump Green's functions, while in section 5 the self-consistency of resulting relations is studied. Our conclusions are presented in section 6.

Wegner-Wilson loop, Field Correlators and Green's Functions
Since our aim is to study confinement for quarks (both light and heavy) and also for massless gluons, we start with the most general Green's functions for these objects in a proper physical background using Fock-Feynman-Schwinger representation. The formalism is built in such a way, that gauge invariance is manifest at all steps. A reader familiar with this technique can skip this section and go directly to the section 3. First, we consider the case of mesons made of quarks while gluon bound states are discussed below. For quarks, one forms initial and final state operators where ψ † , ψ are quark operators, Γ is a product of Dirac matrices, i.e. 1, γ µ , γ 5 , (γ µ γ 5 ), ..., and Φ(x, y) = P exp(ig x y A µ dz µ ) is parallel transporter (also known as Schwinger line or phase factor). The meson Green's function can be written (in quenched case) as and the quark Green's function S q is given in Euclidean space-time (see [21] and references therein): where K is kinetic energy term, m is the pole mass of quark and quark trajectories z µ (τ ) with the end points x and y are being integrated over in (Dz) xy . The factor P σ (x, y, s) in (4) is generated by the quark spin (color-magnetic moment) and is equal to where σ µν = 1 4i (γ µ γ ν −γ ν γ µ ), and P F (P A ) in (6), (4) are, respectively, ordering operators of matrices F µν (A µ ) along the path z µ (τ ).
With (3), (4) one easily gets: where Tr means trace operation both in Dirac and color indices, while In (8) the closed contour C(x, y) is formed by the trajectories of quark z µ (τ ) and antiquark z ′ ν (τ ′ ), and the ordering P A and P F in P σ , P ′ σ is universal, i.e., W σ (x, y) is the Wegner-Wilson loop (W-loop) for spinor particle.
The factors (m −D) and (m ′ −D ′ ) in (7) need a special treatment when being averaged over fields. As shown in Appendix 1 of [22], one can use simple replacement, The representation Eq. (7) is exact in the limit N c → ∞, when internal quark loops can be neglected and is a functional of gluonic fields A µ , F µν , which contains both perturbative and np contributions, not specified at this level.
The next step is averaging over gluon fields, which yields the physical Green's function G qq : The averaging is done with the usual Euclidean weight exp(−Action/ ), containing all necessary gauge-fixing and ghost terms.
To proceed it is convenient to use nonabelian Stokes theorem (see [2] for references and discussion) for the first factor on the r.h.s. in (8) (corresponding to the W-loop for scalar particle) and to rewrite it as an integral over the surface S spanned by the contour C = C(x, y): (11) Here F (u i )ds i = Φ(x 0 , u i )F a µν (u i )t a Φ(u i , x 0 )ds µν (u i ) and u i , x 0 are the points on the surface S bound by the contour C = ∂S. The double brackets ... stay for irreducible correlators proportional to the unit matrix in the colour space (and therefore only spacial ordering P x enters (11)). Since (11) is gauge-invariant, one can make use of generalized contour gauge [23], which is defined by the condition Φ(x 0 , u i ) ≡ 1. Notice that throughout the paper we normalize trace over color indices as Tr1 = 1 in any given representation.
Since (11) is an identity, the r.h.s. does not depend on the choice of the surface, which is integrated over in ds µν (u). On the other hand it is clear that each irreducible n-point correlator F...F integrated over S (these are the functions ∆ (n) [S] in (11)) depends, in general case, on the choice one has made for S. Therefore it is natural to ask the following question: is there any hierarchy of ∆ (n) [S] as functions of n for a given surface S?
To get the physical idea behind this question, let us take the limit of small contour C. In this case one has for the np part of the W-loop in fundamental representation: where S is the minimal area inside the loop C. The np short-distance dynamics is governed by the vacuum gluon condensate G 2 , see (1). Higher condensates and other possible vacuum averages of local operators, in line with the Wilson expansion, have been introduced and phenomenologically estimated in [10], for a review see [24]. Suppose now that the size of W-loop is not small, i.e. it is larger that some typical dynamical scale λ = O(1 GeV) to be specified below. Let us now ask the following question: if one still wants to expand the loop formally in terms of local condensates, what would the structure of such expansion be? It is easy to see from (11) that there are two subseries in this expansion. The first one is just an expansion in n, i.e. in the number of field strength operators. It is just The second series is an expansion of each ∆ (n) [S] in terms of local condensates, i.e. expansion in powers of derivatives. Taking for simplicity the lowest n = 2 term, this expansion reads It is worth mentioning (but usually skipped in QCD sum rules analysis) that the latter series contain ambiguity related to the choice of contours in parallel transporters. The expression (14) is written for the simplest choice of straight contour connecting the points x and y (and not x and x 0 or y and x 0 ). This is the choice we adopt in what follows. In general case each condensate in (14) is multiplied by some contour-dependent function of x, y and x 0 . It is legitimate to ask whether this approximation, i.e. replacement of affects our final results. Since the only difference between the above two expressions is profiles of the transporters, this is the question about contour dependence of FCs. We will not discuss this question in details in the present paper and refer an interested reader to [7]. Here we only mention that the result (weak dependence on transporters' choice) is very natural in Abelian dominance picture since for Abelian fields the transporters exactly cancel. Let us stress that both expansions (13) and (14) are formal since the loop is assumed to be large. Moreover, from dimensional point of view, the terms of these two expansions mix with each other. For example, one has nonzero v.e.v. of two operators of dimension eight: with no a priori reason 2 to drop any of them. Moreover, even if one assumes that all these condensates of high orders can be self-consistently defined analogously to G 2 , it remains to be seen whether the whole series converges or not.
Here the physical picture behind FC method comes into play. It is assumed (and confirmed a posteriori in several independent ways) that if the surface S is chosen as the minimal one, the dominant contribution to (11) results from v.e.v.s of the operators F D k F and this subseries can be summed up as (14). In this language nonlocal two-point FC can be understood as a representation for the sum of infinite sequence of local terms (14).
Physically the dominance of F D n F terms (and hence of two-point FC) corresponds to the fact that the correlation length is small yielding a small expansion parameter ξ =F λ 2 ≪ 1 (F is average modulus of np vacuum fields) or, in other words, typical inverse correlation length λ −1 characterizing ensemble of QCD vacuum fields is parametrically larger than fields themselves. The dimensionless parameter ξ in terms of condensate is given by (G 2 λ 4 ) 1/2 , and it is indeed small: of the order of 0.05 according to lattice estimates. One can say that often repeated statement "QCD vacuum is filled by strong and strongly fluctuating chromoelectric and chromomagnetic fields" is only partly correct -namely, in reality the fields are not so strong as compared to λ −1 . Technically one can say that we sum up the leading subseries in (11) in the same spirit as one does in covariant perturbation theory [25] for weak but strongly varying potentials.
Thus, for understanding of confining properties of QCD vacuum dynamics not only scale of np fields given by G 2 is important but also another quantity -the vacuum correlation length λ, which defines the nonlocality of gluonic excitations. It is worth repeating here that physically λ encodes information about properties of the series (14) and has the same theoretical status as the corresponding np condensates in the following sense: knowledge of all terms in (14) would allow to reconstruct λ. Of course, it is impossible in practice and one has to use other methods to study large distance asymptotics of FC. On more phenomenologic level it was discussed in [26], while rigorous definition was given in the framework of the FC method 3 (see [7], [27] and [11] for review).
It is important to find also the world-line representation for gluon Green's function, which is done in the framework of background perturbation theory (see [28] for details). In this way one finds (where a, b are adjoint color indices and dash sign denotes adjoint representation) where 3 Often referred to as Stochastic Vacuum Model All the above reasoning about asymptotic expansions of the corresponding W-loops is valid here as well.
Thus the central role in the discussed method is played by quadratic (Gaussian) FC of the form where F µν is the field strength and Φ(x, y) is the parallel transporter. Correlation lengths λ i for different channels are defined in terms of asymptotics of (17) at large distances: exp (−|x − y|/λ i ). The physical role of λ i is very important since it distinguishes two regimes: one expects validity of potentialtype approach describing the structure of hadrons of spatial size R and at temporal scale T q for λ i ≪ R, T q , while in the opposite case, when λ i ≫ R, T q the description in terms of spatially homogeneous condensates can be applied. At zero temperature the O(4) invariance of Euclidean space-time holds 4 and FC (17) One has to distinguish from the very beginning perturbative and np parts of the correlators D(z), D 1 (z). Beginning with the former, one easily finds at tree level D p,0 (z) = 0 ; D p,0 where At higher orders situation becomes more complicated. Namely, one has at n-loop order where the gauge-invariant functions G (n) (z), G 1 (z) has the following general structure: where c n is numerical coefficient, µ is renormalization scale and we have omitted subleading logarithms and constant terms in the r.h.s. Explicit expressions for the case n = 1 can be found in [29,30]. Naively one could take perturbative functions D p,n (z), D p,n 1 (z) written above and use them for computations of W-loop or static potentials in Gaussian approximation. This however would be incorrect. The reason is that at any given order in perturbation expansion over α s one should take into account perturbative terms of the given order coming from all FCs, and not only from the Gaussian one. For example, at one loop level the function D p,1 (z) is nonzero which naively would correspond to area law at perturbative level. Certainly this cannot be the case, since self-consistent renormalization program for W-loops is known not to admit any terms of this sort [31]. Technically in FC language the correct result is restored by cancelation of contributions proportional to D p,1 (z) to all observables by the terms coming from triple FC (see details in [32]).
In view of this general property of cluster expansion it is more natural just to take relation D p,n (z) ≡ 0 as valid at arbitrary n, having in mind that proper number of terms from (13) has to be included to get this result. Thus we assume that the following decomposition takes place: and D(z) has a smooth limit when z → 0. As for D 1 (z) its general asymptotic at z → 0 reads: where c and a 2 weakly (logarithmically) depend on z. Notice that by D 1 (0) we always understand np "condensate" part in what follows: The term 1/z 2 deserves special consideration. One may argue, that at large distances the simple additivity of perturbative and np contributions to V QQ potential is violated [33]. It is interesting whether this additivity holds at small distances and, in particular, mixed condensate a 2 was suggested in [34] and argued to be phenomenologically desirable.
As we demonstrate below our approach is self-consistent with a 2 = 0, i.e. when perturbative-np additivity holds at small distances and there is no dimension-two condensate. However, at intermediate distances a complicated interrelation occurs, which does not exclude (22) being relevant in this regime.
We proceed with general analysis of (17). As is proved in [7], D(z), D 1 (z) do not depend on relative orientation of the plane (µν), plane (λσ), and the vector z α . This orientation can be of the following 3 types: a) planes (µν) and (λσ) are perpendicular to the vector z α ; b) planes (µν), (λσ) are parallel or intersecting along one direction and z α lies in one of planes; c) planes (µν) and (λσ) are perpendicular and z α lies necessarily in one of them.
It is easy to understand that this classification is O(4) invariant, and therefore one can assign indices a, b, c to D µν,λσ , resulting in general case in three different functions. Physically speaking the case a) refers to e.g. color magnetic fields F ik , F lm with z α in the 4th direction, and the corresponding FCs will be denoted as D ⊥ (z). The case b) refers to e.g. the color electric fields E i (x), E k (y), connected by the same temporal links along 4th axis, and the corresponding FCs are denoted as D (z). Finally, in case c) only D 1 part survives in (18) and it will be given by the subscript EH, D EH 1 . In case of zero temperature, when O(4) invariance holds, these three , (note that just D ⊥ , D were measured on the lattice in [14]). For nonzero temperature the correlators of colorelectric and colormagnetic fields can be different and in general there are five FCs: D E (z), D E 1 (z), D H (z), D H 1 (z), and D EH 1 (z). As a result one obtains the full set of five independent quadratic FCs which can be defined as follows It is interesting that the structure (24)-(26) survives for nonzero temperature, where O(4) invariance is violated: it tells that the equations (24)- (26) give the most general form of the FCs. Note that at zero temperature, D H = D E , D H 1 = D E 1 at coinciding point; the O(4) invariance requires that colorelectric and colormagnetic condensates coincide: Here we assume the np part of all FCs is finite.
As is shown below, all correlation lengths λ j appear to be just inverse masses of the corresponding gluelumps λ j = 1/M j . Having analytic expressions for FCs one might ask how to check them versus experimental and lattice data. On experimental side in hadron spectroscopy one measures masses and transition matrix elements, which are defined by dynamical equation and the latter can be used of potential type due to smallness of λ j . In static potential only integrals over distance enter and the spin-independent static potential can be written as [7,35] V QQ (r) = 2 At this point one should define how perturbative and np contributions combine in D E , D E 1 , and this analysis was done in [32]. Making use of (19) together with definition of the string tension [7], one obtains from (29) the standard form of the static potential for N = 3 at As is argued below, λ E ≈ 0.1 Fm, so that the form (31) is applicable in the whole range of distances r > 0.1 fm provided asymptotic freedom is taken into account in α s (r) in (31). At the smallest values of r, r < ∼ λ E , one would obtain a softening of confining term, σr → cr 2 [7], which is not seen in accurate lattice data at r > ∼ 0.2 fm [36,37], imposing a stringent limit on the value of λ E ; λ E < ∼ 0.1 Fm, see discussion in [37].

Spin-Dependent Potentials and FC
We now turn to the spin-dependent interactions to demonstrate that FC can be extracted from them [35,40,41] and one can extract the properties of the latter from the spin-dependent potentials. This procedure is used recently for comparison of analytic predictions [20] with the lattice data [19].
To express spin-dependent potentials in terms of FC, one can start with the W-loop (8), entering as a kernel in the qq Green's function, where quarks q,q can be light or heavy. Since both kernels P σ , P ′ σ contain the matrix σ µν F µν , the terms of spin-spin and spin-orbit types appear. As before, we keep Gaussian approximation, and in this case spin-dependent potentials can be computed not only for heavy, but also for light quarks. They have the Eichten-Feinberg form [42]: Spin-spin interaction appear in W σ , Eq.(8) which can be written as where z 1,2 = z(τ 1,2 ) and spin-orbit interaction arises in (17) from the products σ µν F µν ds ρλ F ρλ . It is clear, that the resulting interaction will be of matrix form (2 × 2) × (2 × 2) (not accounting for Pauli matrices). If one keeps only diagonal terms in σ µν F µν (as the leading terms for large µ i ≈ M) one can write for the spin-dependent potentials the representation of the Eichten-Feinberg form (32).
At this point one should note that the term with dε dR in (32) was obtained from the diagonal part of the matrix (m −D)σ µν F µν , namely as product iσ k D k · σ i E i , see [35,40] for details of derivation, while all other potentials V i , i = 1, 2, 3, 4 are proportional to FCs H i ΦH k Φ . One can relate FCs of colorelectric and colormagnetic fields to D E , D E 1 , D H , D H 1 defined by (24)-(26) with the following result: One can check, that the Gromes relation [18] acquires the form [41]: For T = 0, when D E = D H , D E 1 = D H 1 , the Gromes relations are satisfied identically, however for T > 0 electric and magnetic correlators are certainly different and Gromes relation is violated, as one could tell beforehand, since for T > 0 the Euclidean O(4) invariance is violated.

Gluelumps and FC
In this section we establish a connection of FC D(z) and D 1 (z) with the gluelump Green's functions.
To proceed one can use the background field formalism [28] [43], where the notions of valence gluon field a µ and background field B µ are introduced, so that total gluonic field A µ is written as The main idea we are going to adopt here is suggested in the second paper of [15]. To be self-contained, let us explain it here in simple form. First, we single out some color index a and fix it at a given number. Then in color components we have, by definition, and there is neither summation over a in the first term nor over c in the second. As a result the integration measure factorizes: DB c µ so one first averages over the fields B c µ , and after that one integrates over the fields a a µ . The essential point is that the integration over DB c µ provides colorless adjoint string attached to the gluon a a µ , which keeps the color index "a" unchanged in the course of propagation of gluon a in background made of B. This is another way of saying that gluon field of each particular color enters QCD Lagrangian quadratically. Despite such separation of the measure is a kind of trick, the formation of adjoint string is the basic physical mechanism behind this background technic and it is related to the properties of gluon ensemble: even for N c = 3 one has one color degree of freedom a a µ and 7 fields B c µ , also confining string is a colorless object and therefore the singled-out adjoint index "a" can be preserved during interaction process of the valence gluon a a µ with the background (B c µ ). These remarks to be used in what follows make explicit the notions of the valence gluon and background field.
We assume the background Feynman gauge D µ a µ = 0 and assign field transformations as follows As a result the parallel transporter Φ(x, y) = P exp ig x y B µ (z)dz keeps its transformation property and every insertion of a µ between Φ transforms gauge covariantly. In what follows we shall assume that Φ is made entirely of the field B µ . Note that in the original cluster expansion of the W-loop W (B + a) in FC [43].
it does not matter whether (B µ + a µ ) or B µ enters Φ ′ s since those factors cancel in the sum. Thus we assume from the beginning that Φ ′ s are not renormalized, since such renormalization is unphysical anyway according to what have been said above. As for renormalization having physical meaning, it is known [31] to reduce to the charge renormalization and renormalization of perimeter divergences on the contour C. For static quarks the latter reduces to the mass renormalization.
In background field formalism one has an important simplification: the combination gB µ is renorminvariant [28,31], so is any expression made of B µ only. This point is important for comparison of FC with the lattice correlators in [19,20]. Renormalization constants for field strengths Z b used there account for finite size of plaquettes and they are similar to the lattice tadpole terms. Now we turn to the analytic calculation of FC in terms of gluelump Green's function. To this end we insert (39) into (17) and have for F µν (x): Here, the term F µν can be omitted, if summing over all indices a is presumed to be done at the end of calculation.
As a result D µν,λσ can be written as where the superscript 0,1,2 refers to the power of g coming from the term ig[a µ , a ν ].
We can address now an important question about the relation between FCs and gluelump Green's functions. We begin with 1-gluon gluelump, whose Green's function reads According to (40) µν (x, y) is a gauge invariant function. As shown in [16], the first term in (43) is connected to the functions D E 1 , D H 1 and it can be written as follows where ∆ (0) µν,λσ contains contribution of higher FCs, which we systematically discard.
On the other hand one can find G (1g) µν (x, y) from the expression [21] written as where the spin-dependent W-loop is and the gluon spin factor is exp F ≡ exp(2ig s 0 dτF B (z(τ ))) withF B made of the background field B µ only.
Analogous expression can be constructed for Green's function of 2-gluon gluelump. It is given by the following expression At small distances G (2gl) (x, y) is dominated by perturbative expansion terms, however, as we already discussed in details, all these perturbative terms are canceled by those from higher FCs (triple, quartic, etc...), therefore expansion in fact starts with np terms of dimension four.
To identify the np contribution to G (2gl) (z) we rewrite it as follows: where the two-gluon gluelump W-loop W Σ (C 1 , C 2 ) is depicted in Fig. 6 and can be written as (in the Gaussian approximation) and the total surface S consists of 3 pieces, as shown in Fig.7 F where (F ) has Lorentz tensor and adjoint color indices F ab ij and lives on gluon trajectories, i.e. on the boundaries of S i . In full analogy with (45) we have The crucial point is that the Green's function (45), (50) can be calculated in terms of the same FCs D(z), D 1 (z). Indeed one has where the index w stays for 1-gluon or 2-gluon gluelump Hamiltonians. The latter are expressed via the same FCs D(z), D 1 (z) (see [45] and references therein): where the last term H string [D, ν] depends on D(z) via adjoint string tension and Hamiltonian depends on einbein fields µ and ν. The term corresponding to perimeter Coulomb-like interaction ∆H Coul [D 1 ] depends on D 1 (z) (compare with (29), (32)). The self-consistent regimes correspond to different asymptotics of the solutions to these equations. In Coulomb phase of a gauge theory both H 1g , H 2g exhibit no mass gap, i.e. large z asymptotic of (54) is power-like. The function D(z) vanishes in this phase and W-loop obeys perimeter law. In the confinement phase realized in Yang-Mills theory at low temperatures a typical large-z pattern is given by i.e. there is a mass gap for both H 1g , H 2g .
Confining solutions are characterized by W-loops obeying area law. In other words, Hamiltonians expressed in terms of interaction kernels depending on D(z), D 1 (z) exhibit mass gap if these kernels are confining. On the other hand the same mass gap plays a role of inverse correlation length of the vacuum. This should be compared with well known mean-field technique. In our case the role of mean-field is played by quadratic FC, which develops non-trivial Gaussian term D(z). Notice that the exponential form of its large distance asymptotics exp(−|z|/λ) (and not, let's say, exp(−z 2 /λ 2 ) ) is dictated by spectral expansion of the corresponding Green's function at large distances.
Full solution of the above equations is a formidable task not addressed by us here. Instead, as a necessary prerequisite we check below different asymptotic regimes and demonstrate self-consistency of the whole picture.
We begin with small distance region. The np part of contribution to D 1 (z) at small distances comes from two possible sources: the area law term (first exponent and the exp F term in (47)). Both contributions are depicted in Fig.1 and Fig.2 respectively. We shall disregard the term exp F in this case, since for the one-gluon gluelump it does not produce hyperfine interaction and only gives rise to the np shift of the gluon mass, which anyhow is eliminated by the renormalization (see Appendix for more detail). This is in contrast to the two-gluon gluelump Green's function generating D E , D H , where the hyperfine interaction between two gluons is dominating at small distances.
As a result Eq. (44) without the (exp F ) term yields It is remarkable that the sign of the np correction is positive. At large distances one can use the gluelump Hamiltonian for one-gluon gluelump from [47] to derive the asymptotics [16] where M (1) GeV for σ f = 0.18 GeV 2 [47,48]. We now turn to the FC D(z) as was studied in [16]. The relation (53) connects D(z) to the two-gluon gluelump Green's function, studied in [47] at large distances. Here we need its small -z behaviour and we shall write it in the form where G (2gl) p (z) contains purely perturbative contributions which are subtracted by higher-order FCs, while G (2gl) np (z) contains np and possible perturbativenp interference terms. We are interested in the contribution of the FC F F to G (2gl) (z), when z tends to zero.
One can envisage three types of contributions: a) due to product of surface elements ds µν ds λσ , which gives for small surface the factor exp − , similar to the situation discussed for D 1 . This is depicted in Fig.7.
b) contribution of the type dτ 1 dτ 2 F F , where dτ and dτ 2 belong to different gluon trajectories, 1 and 2, as depicted in Fig.8. This is the hyperfine gluon-gluon interaction, taken into account in [47] in the course of the gluelump mass calculations (regime of large z), however in that case mostly the perturbative part of F F contributes (due to D 1 ). Here we keep in F F only the NP part and consider the case of small z. c) Again the term dτ i dτ j , but now i = j. This is actually a part of the gluon selfenergy correction, which should be renormalized to zero, when all (divergent) perturbative contributions are added. As we argued in Appendix of [16], we disregard this contribution here, as well as in the case of D 1 (z) (one-gluon gluelump case). Now we treat both contributions a) and b). a) In line with the treatment of D 1 , Eq. (57), one can write the term, representing two-gluon gluelump as two nearby one-gluon gluelumps which yields for D(z), ) in this case one should consider the diagram given in Fig. 8, which yields the answer (for details see Appendix).
which contributes to D(z) as where at small z, h(z) ≈ D(λ 0 ) 64π 4 log 2 λ 0 √ e z , λ 0 is of the order of correlation length λ.
At large distances one use the two-gluon gluelump Hamiltonian as in [48] and find the corresponding spectrum and wave functions, see [48] and appendix 5 of [16] for details. As a result one obtains in this approximation where M (2) 0 is the lowest two-gluon gluelump mass found in [47] to be about M (2) 0 = (2.5 ÷ 2.6) GeV We shall discuss the resulting properties of D(x) and D 1 (x) in the next section.

Discussion of Consistency
We start with the check consistency for D(z). As is shown above, D(z) has the following behavior at small z Since α s (µ(z)) ∼ 2π/β 0 log(Λz) −1 the first term is subleading at z → 0 and the last term on the r.h.s. tends to a constant: Since from (65) one can infer, that is an increasing function of z at small z, z < ∼ λ 0 , and for z ≫ λ one observes exponential falloff. The qualitative picture illustrating this solution for D(z) is shown in Fig.5.
This pattern may solve qualitatively the contradiction between the values of D(0) estimated from the string tension D σ (0) ≃ σ πλ 2 ≈ 0.35 GeV 4 and the value obtained in naive way from the gluon condensate D G 2 (0) = π 2 18 G 2 ≈ (0.007÷0.012) GeV 4 . One can see that D σ (0) ≈ (30÷54)D G 2 (0). This seems to be a reasonable explanation of the mismatch discussed in the introduction.
As was shown in [16], the large distance exponential behavior is selfconsistent, since (assuming that it persists for all z, while small z region contributes very little) from the equality σ = πλ 2 D σ (0), comparing with (63) one obtains where in α s (µ) the scale µ corresponds roughly to the gluelump average momentum (inverse radius) µ 0 ≈ 1 GeV. Thus (66) yields α s (µ 0 ) ≈ 0.4 which is in reasonable agreement with α s from other systems [49].
We end this section with discussion of three points: 1) D(z) and D 1 (z) have been obtained here in the leading approximation, when gluelumps of minimal number of gluons contribute: 2 for D(z) and 1 for D 1 (z). In the higher orders of O(α s ) one has an expansion of the type Hence the asymptotic behavior for D(z) will contain exponent of M (1) 0 |z| too, but with a small preexponent coefficient.
2) The behavior of D(z), D (np) 1 (z) at small z is defined by NP terms of dimension four, which are condensate G 2 as in (12), and the similar term from the expansion of exp F , namely, therefore one does not encounter mixed terms like O( m 2 z 2 ), however if one assumes that expansion of F (0)F (z) starts with terms of this sort as it was suggested in [34], then one will have a selfconsistent condition for the coefficient in front of this term.
3) To study the difference between D E 1 , D H 1 at T ≤ T c one should look at Eq. (45) and compare the situation, when µ = λ = 4, ν = σ = i and take z = x − y along the 4-th axis (for D E 1 ), while for D H 1 one takes µ = λ = i, ν = σ = k and the same z. One can see that in both cases one ends up with the one-gluon gluelump Green's function G µν (46) which in the lowest (Gaussian) approximation is the same G µν = δ µν f (z), and f (z) is the standard lowest mass (and lowest angular momentum) Green's function.
Hence D E 1 = D H 1 in this approximation and for T = 0 this is an exact relation, as was discussed above. Therefore The value M 0 in (69) is taken from the calculations in [47]. The same is true for D E , D H , as it is seen from (48), where G (2gl) corresponds to two-gluon subsystem angular momentum L = 0 independently of µν, λσ.

Conclusions
We have derived, following the method of [28] the expressions for FC D(z), D 1 (z) in terms of gluelump Green's functions. This is done in Gaussian approximation. The latter are calculated using Hamiltonian where np dynamics is given by D np (z), D np 1 (z). In this way one obtains selfcoupled equations for these functions, which allow two types of solutions: 1) D np (z) = 0, D np 1 (z) = 0, i.e. no np effects at all; 2) D np (z), D np 1 (z) are nonzero and defined by the only scale, which should be given in QCD, e.g. string tension σ or Λ QCD . All other quantities are defined in terms of these basic ones. We have checked consistency of selfcoupled equations at large and small distances and found that to the order O(α s ) no mixed perturbative-np terms appear. The function D 1 (z) can be represented as a sum of perturbative and np terms, while D(z) contains only np contributions.
We have found a possible way to explain the discrepancy between the average values of field strength taken from G 2 and from σ by showing that D(z) has a local minimum at z = 0 and grows at z ∼ λ. Small value of λ and large value of the gluelump mass M gl = 1/λ ≈ 2.5 GeV explains the lattice data for λ ≈ 0.1 Fm. Thus the present paper argues that relevant degrees of freedom ensuring confinement are gluelumps, described self-consistently in the language of FCs.
Taking into account that ∆D 1 (x) = − 2g 2 N 2 c d dx 2 ∆f (x), one obtains for x → 0, ∆D 1b (x) = 2g 2 G 2 I(x) where we have defined The integral I(0) diverges at small v and large u, v. The latter divergence is removed since at large arguments G(x, y) is damped by confining force, keeping the propagating gluon nearby the 4-th axis, x 4 = y 4 = 0 where the static gluonic source resides. One can estimate I(x → 0) ∼ O((ln λ 2 x 2 ) α )). Since G 2 is connected to D, D 1 one can see in (74), that coefficients of D E 1 on both sides of (74) have the same order of magnitude and the same sign, suggesting a selfconsistency on this preliminary level. However, as we argued before and in [16] this contribution is actually gluon mass renormalization, which is zero. This is especially clear when F F in Fig.2 is replaced by its perturbative part.
We turn now to the most important case of the confining FC D E (x). The corresponding two-gluon gluelump amplitude is depicted in Fig.3, and the np contributions to D E (x) are given in Fig.4 and Fig.5.
We start with the amplitude of Fig.4, which can be represented as a doubled diagram of Fig.1 and therefore the contribution to the two-gluon gluelump Green's function G (2gl) (x) will be .
Now using relation (12) one obtains for the np contribution of Fig.4 ∆D a (x) = g 4 (N 2 Note, that this contribution is finite at x → 0, however with the negative sign. We now turn to the np contribution of Fig.5, which can be obtained from (48) inserting there the product of operatorsF f g (w)d 4 wF f ′ g ′ (w ′ )d 4 w ′ , where f g(f ′ g ′ ) are adjoint color indices,F f g = F a T a f g = F a (−i)f af g . Using the formulas one arrives at the expression where G(y) is the gluon Green's function in the two-gluon gluelump; the gluon is connected at large distances by two strings to another gluon and to the static gluon source, see Fig.5. At small x (we take it along axis 4 for convenience) the integral in h(x) grows when G(y) is close to the 4th axis and becomes the free gluon, G(y) → G 0 (y) = 1 4π 2 y 2 . As a result one obtains and one should have in mind that the integral for x → 0 diverges both at small and at large w, w ′ . The small w, w ′ region will give the terms O(ln λ 2 0 x 2 ), and at large w, w ′ the integral is protected by the fall-off of G(y) at large y due to confinement, therefore we must imply in the integrals in (81) the upper limits for w, w ′ at some λ 0 > ∼ 1 √ σ adj . Then introducing for w(w, w 4 ) polar coordinates |w| = ρ sin θ, w 4 = ρ cos θ (and the same for w ′ ), one arrives at the integral π 0 sin 2 θdθ ρ 2 − 2xρ cos θ + x 2 = π 2ρ 2 (ρ ≥ x) or π 2x 2 (x ≥ ρ), and finally one has the estimate of contribution to (80) from the region of small w, w ′ (w|w|, |w ′ | < ∼ λ 0 ) and finally for D(x) one obtains (omitting the perturbative term in (77), which cancels,