World-Line Formalism: Non-Perturbative Applications

Abstract

This review addresses the impact on various physical observables which is produced by confinement of virtual quarks and gluons at the level of the one-loop QCD diagrams. These observables include the quark condensate for various heavy flavors, the Yang-Mills running coupling with an infra-red stable fixed point, and the correlation lengths of the stochastic Yang-Mills fields. Other non-perturbative applications of the world-line formalism presented in the review are devoted to the determination of the electroweak phase-transition critical temperature, to the derivation of a semi-classical analogue of the relation between the chiral and the gluon QCD condensates, and to the calculation of the free energy of the gluon plasma in the high-temperature limit. As a complementary result, we demonstrate Casimir scaling of k-string tensions in the Gaussian ensemble of the stochastic Yang-Mills fields.

1. Quark Condensate for Various Heavy Flavors

As is well known, because of confinement in QCD, quarks and gluons do not exist as individual particles, but appear only in the form of bound states (for recent reviews, see [1,2]). The latter include mesons, baryons (as well as other possible quark bound states, such as tetra- and pentaquarks), glueballs, and the so-called hybrids consisting of a quark, an antiquark, and one or several gluons. Each bound state can be represented as an average of a certain operator over the Euclidean QCD vacuum. Since the vacuum state is gauge-invariant, only the gauge-invariant operators yield non-vanishing results once averaged over it. (This statement can be viewed as a corollary of the so-called Elitzur theorem [3,4], which states that a local gauge symmetry cannot be broken spontaneously.) For local composite operators, such as $\overline{\psi }\psi$ or ${\left(F_{\mu \nu }^a\right)}^2$, gauge invariance holds automatically, and the mean values of these operators yield the QCD condensates [5,6]. Nevertheless, the quark and the gluon bound states are in general represented by the field operators which are defined at different points of the Euclidean space, so that a phase factor (also called a parallel transporter or a Schwinger string) interpolating between these points is required in order to provide gauge invariance of the full non-local operator. That is, the non-local gauge-invariant operators have, e.g., the form

$${\overline{\psi }}^i\left(x\right){\mathsf{\Phi }}_{x{x}^{\prime }}^{ij}{\psi }^j\left(x^{\prime }\right),F_{\mu \nu }^a\left(x\right){\mathsf{\Phi }}_{x{x}^{\prime }}^{ab}{F}_{\lambda \rho }^b\left(x^{\prime }\right),$$

(1)

where ${\mathsf{\Phi }}_{x{x}^{\prime }}^{ij}$ and ${\mathsf{\Phi }}_{x{x}^{\prime }}^{ab}$ are the phase factors transforming under the fundamental and the adjoint representations of SU(N), respectively, with $i,j=1,\dots ,N$ and $a,b=1,\dots ,N^2-1$. (Throughout this review, we denote the number of colors by N.) Thus, the phase factors provide long-range correlations of color between individual quark operators ψ and $\overline{\psi }$, as well as the gluonic field strengths $F_{\mu \nu }^a$, thereby making the non-local operators (1) gauge invariant. (Note, however, that, in spite of their gauge invariance, the operators (1) still depend on the shape of a contour that enters the phase factor ${\mathsf{\Phi }}_{x{x}^{\prime }}^{ij}$ or ${\mathsf{\Phi }}_{x{x}^{\prime }}^{ab}$.) Such long-range correlations of color are associated with confining interactions between the constituents of the QCD bound states. By accounting for confinement of gauge particles which propagate in the loops, one can arrive at the results inaccessible within the conventional perturbation theory, such as the infra-red finiteness of the running strong coupling. The world-line representation for the Green function of a particle, which propagates along a closed Euclidean trajectory and interacts with quantum gauge fields, contains a Wilson-loop average. In general, the corresponding world-line integral cannot be calculated analytically, because finding the minimal surface for an arbitrary contour $z_{\mu }\left(\tau \right)$ in $D>2$ dimensions is too complicated. For this reason, effective parametrizations of the minimal surface have been invented in the literature [7,8,9] in order to construct the minimal-area functional in terms of $z_{\mu }\left(\tau \right)$. Owing to the fact that the Wilson-loop average is fully defined in terms of the geometric characteristics of the contour $z_{\mu }\left(\tau \right)$, this approach has an advantageous feature of being able to reduce the calculation of vacuum amplitudes in QCD to the calculation of the world-line integrals in an effective Abelian gauge theory, or even in quantum mechanics.

We start demonstrating these ideas at an example of the calculation of the heavy-quark condensate [10]. Numerically, the wavelengths of the heavy, i.e., c-, b-, and t-, quarks are smaller than the Yang-Mills vacuum correlation length λ (see the details below), so that the typical Euclidean trajectories of these quarks fit into a circle of diameter λ. To calculate the corresponding small-sized Wilson-loop average, one can use the non-Abelian Stokes’ theorem and the cumulant expansion to represent this object as [11,12,13,14,15,16]

$$\langle W\left(C\right)\rangle \simeq \frac{1}{N}\mathrm{tr}exp\left[-\frac{1}{2!}\frac{g^2}{4}{\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x\right){\int }_{S_{\min }}d{\sigma }_{\lambda \rho }\left(x^{\prime }\right)\langle F_{\mu \nu }^a\left(x\right)T^a{\mathsf{\Phi }}_{x{x}^{\prime }}{F}_{\lambda \rho }^b\left(x^{\prime }\right)T^b{\mathsf{\Phi }}_{x^{\prime }x}\rangle \right].$$

(2)

In this expression, $T^a$ is a generator of the fundamental representation of SU(N), and the integrations go over the minimal surface $S_{\min }$ bounded by the quark trajectory C. Such a minimal surface is unique for a given contour C, i.e., it is defined by C in an unambiguous way. Furthermore,

$${\mathsf{\Phi }}_{x{x}^{\prime }}\equiv \frac{1}{N}\mathcal{P}exp\left[ig{\int }_{x^{\prime }}^x{T}^a{A}_{\mu }^a\left(u\right)d{u}_{\mu }\right]$$

(3)

in Equation (2) is a phase factor along the straight line, and the factor of $1/2!$ stems from the cumulant expansion, while the factor of $1/4$ is due to the non-Abelian Stokes’ theorem. For surfaces at issue, which fit into a circle of size λ, the field-strength tensor $F_{\mu \nu }^a\left(x\right)$ can be treated as constant, so that the correlation function in Equation (2) can be approximated as

$$\langle F_{\mu \nu }^a\left(x\right)T^a{\mathsf{\Phi }}_{x{x}^{\prime }}{F}_{\lambda \rho }^b\left(x^{\prime }\right)T^b{\mathsf{\Phi }}_{x^{\prime }x}\rangle \simeq \frac{{\widehat{1}}_{N\times N}}{N}·\frac{1}{12}\left({\delta }_{\mu \lambda }{\delta }_{\nu \rho }-{\delta }_{\mu \rho }{\delta }_{\nu \lambda }\right)·\mathrm{tr}\left(T^a{T}^b\right)·\langle F_{\alpha \beta }^a\left(x\right)F_{\alpha \beta }^b\left(x\right)\rangle .$$

(4)

Using this parametrization in Equation (2), and noticing that $\mathrm{tr}{T}^a{T}^b=\frac{1}{2}{\delta }^{ab}$, one obtains

$$\langle W\left(C\right)\rangle \simeq exp\left[-\frac{\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle }{96N}{\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x\right){\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x^{\prime }\right)\right].$$

(5)

In particular, for a flat non-selfintersecting contour C, the double surface integral in this expression becomes ${\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x\right){\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x^{\prime }\right)=2{\mathsf{\Sigma }}^2$, where Σ is the area of the flat surface bounded by C. Thus, one gets the “area-squared” law for the contribution of soft stochastic Yang-Mills fields to a small-sized Wilson-loop average [11,12,13,14,15,16]

$$\langle W\left(C\right)\rangle {|}_{\mathrm{flat},\mathrm{non}-\mathrm{selfintersecting}C}\simeq exp\left[-\frac{\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle }{48N}·{\mathsf{\Sigma }}^2\right].$$

(6)

Notice that the proportionality of the coefficient at ${\mathsf{\Sigma }}^2$ to $\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle$ in this formula has for the first time been found in Refs. [17,18].

Rewriting the double surface integral in this formula as

$$\int d{\sigma }_{\mu \nu }\left(x\right)\int d{\sigma }_{\mu \nu }\left(x^{\prime }\right)=-\frac{1}{2}\int d{\sigma }_{\mu \nu }\left(x\right)\int d{\sigma }_{\mu \rho }\left(x^{\prime }\right){\partial }_{\nu }^x{\partial }_{\rho }^{x^{\prime }}{(x-x^{\prime })}^2,$$

and using the (Abelian) Stokes’ theorem, we obtain

$$\int d{\sigma }_{\mu \nu }\left(x\right)\int d{\sigma }_{\mu \nu }\left(x^{\prime }\right)=-\frac{1}{2}\oint d{z}_{\mu }\oint d{z}_{\mu }^{\prime }{(z-z^{\prime })}^2.$$

One can notice that only the $\left(z{z}^{\prime }\right)$-term from ${(z-z^{\prime })}^2$ yields a non-vanishing contribution to the last integral, so that

$$\int d{\sigma }_{\mu \nu }\left(x\right)\int d{\sigma }_{\mu \nu }\left(x^{\prime }\right)=-\frac{1}{2}\oint d{z}_{\nu }\oint d{z}_{\nu }^{\prime }(-2z{z}^{\prime })={\left(\oint d{z}_{\nu }{z}_{\mu }\right)}^2.$$

(7)

Inserting this expression back into Equation (5), we can further apply the Hubbard-Stratonovich trick to represent the Wilson-loop average as

$$\langle W\left[z_{\mu }\right]\rangle =\frac{1}{{\left(8\pi G\right)}^3}·\left(\prod _{\mu <\nu }^{}{\int }_{-\infty }^{+\infty }d{B}_{\mu \nu }\right){\mathrm{e}}^{-\frac{B_{\mu \nu }^2}{16G}-\frac{i}{2}{B}_{\mu \nu }\oint d{z}_{\nu }{z}_{\mu }},$$

(8)

where we have denoted for brevity $G\equiv \frac{\langle {\left(g{F}_{\mu \nu }\right)}^2\rangle }{96N}$. Since the auxiliary antisymmetric field $B_{\mu \nu }$ does not depend on a point of the Euclidean space, the Wilson-loop average in the form of Equation (8) can be readily applied to the calculation of the one-loop effective action. The latter follows from the QCD partition function upon the integration over the quark fields, and reads [19,20,21,22,23,24,25]

$$\langle \mathsf{\Gamma }\left[A_{\mu }^a\right]\rangle =-2N_{\mathrm{f}}{\int }_0^{\infty }\frac{ds}{s}{\mathrm{e}}^{-m^{2}s}{\int }_P^{}\mathcal{D}{z}_{\mu }{\int }_A^{}\mathcal{D}{\psi }_{\mu }{\mathrm{e}}^{-{\int }_0^{s}d\tau \left(\frac{1}{4}{\dot{z}}_{\mu }^2+\frac{1}{2}{\psi }_{\mu }{\dot{\psi }}_{\mu }\right)}\times$$

$$\times \left\{〈\mathrm{tr}\mathcal{P}exp\left[ig{\int }_0^{s}d\tau T^a\left(A_{\mu }^a{\dot{z}}_{\mu }-{\psi }_{\mu }{\psi }_{\nu }{F}_{\mu \nu }^a\right)\right]〉-N\right\}.$$

(9)

In this formula, P and A stand, respectively, for the periodic and the antiperiodic boundary conditions, so that ${\int }_P\equiv {\int }_{z_{\mu }\left(s\right)=z_{\mu }\left(0\right)}$, ${\int }_A\equiv {\int }_{{\psi }_{\mu }\left(s\right)=-{\psi }_{\mu }\left(0\right)}$. The trajectories $z_{\mu }\left(\tau \right)$ obey the equation

$${\int }_0^{s}d\tau z_{\mu }\left(\tau \right)=0.$$

(10)

Equation (10) means that the center of a trajectory is located at the origin, i.e., the overall factor of volume, which is associated with the translation of a trajectory as a whole, is already divided out. In Equation (9), $T^a$ is a generator of the SU(N)-group in the fundamental representation, and $N_{\mathrm{f}}$ is the number of quark flavors. Since the quark condensate is always associated with a given flavor, we set $N_{\mathrm{f}}=1$. Furthermore, since the quark condensation can occur only due to the gauge fields, we have subtracted the free part of the effective action, so that $\langle \mathsf{\Gamma }\left[0\right]\rangle =0$.

Equation (9) can be simplified by using the fact that the product $F_{\mu \nu }^a{T}^a$, which enters the quark spin term in the world-line action, can be recovered through the area-derivative operator $\frac{\delta }{\delta {\sigma }_{\mu \nu }}$ acting on the Wilson loop [26,27]. This observation allows one to reduce the gauge-field dependence of Equation (9) to that of the Wilson-loop average as follows:

$$\langle \mathsf{\Gamma }\left[A_{\mu }^a\right]\rangle =-2N{\int }_0^{\infty }\frac{ds}{s}{\mathrm{e}}^{-m^{2}s}{\int }_P^{}\mathcal{D}{z}_{\mu }{\int }_A^{}\mathcal{D}{\psi }_{\mu }{\mathrm{e}}^{-{\int }_0^{s}d\tau \left(\frac{1}{4}{\dot{z}}_{\mu }^2+\frac{1}{2}{\psi }_{\mu }{\dot{\psi }}_{\mu }\right)}\times$$

$$\times \left\{exp\left[-2{\int }_0^{s}d\tau {\psi }_{\mu }{\psi }_{\nu }\frac{\delta }{\delta {\sigma }_{\mu \nu }\left(z\right)}\right]\langle W\left[z_{\mu }\right]\rangle -1\right\}.$$

(11)

Substituting into this expression Equation (8) for $\langle W\left[z_{\mu }\right]\rangle$, we have

$$\langle \mathsf{\Gamma }\left[A_{\mu }^a\right]\rangle =-\frac{2N}{{\left(8\pi G\right)}^3}{\int }_0^{\infty }\frac{ds}{s}{\mathrm{e}}^{-m^{2}s}\left(\prod _{\mu <\nu }^{}{\int }_{-\infty }^{+\infty }d{B}_{\mu \nu }\right){\mathrm{e}}^{-\frac{B_{\mu \nu }^2}{16G}}·\left\{{\int }_P^{}\mathcal{D}{z}_{\mu }{\int }_A^{}\mathcal{D}{\psi }_{\mu }\times \right.$$

$$\left.\times exp\left[-{\int }_0^{s}d\tau \left(\frac{1}{4}{\dot{z}}_{\mu }^2+\frac{1}{2}{\psi }_{\mu }{\dot{\psi }}_{\mu }+\frac{i}{2}{B}_{\mu \nu }{z}_{\mu }{\dot{z}}_{\nu }-i{B}_{\mu \nu }{\psi }_{\mu }{\psi }_{\nu }\right)\right]-\frac{1}{{\left(4\pi s\right)}^2}\right\}.$$

The world-line integral in the latter formula is given by the so-called Euler-Heisenberg-Schwinger Lagrangian [28,29,30,31], which to the order $\mathcal{O}\left(s^2{B}_{\mu \nu }^2\right)$ yields for the curly bracket:

$$\left\{\cdots \right\}\simeq \frac{1}{{\left(4\pi s\right)}^2}·\frac{s^2}{3}\sum _{\alpha <\beta }^{}{B}_{\alpha \beta }^2.$$

The $B_{\mu \nu }$-integration can now be performed by using the value of ${\pi }^3$ for the full solid angle in six dimensions. That yields for the heavy-quark condensate ${\langle \overline{\psi }\psi \rangle }_{\mathrm{heavy}}=-\frac{\partial }{\partial m}\langle \mathsf{\Gamma }\left[A_{\mu }^a\right]\rangle$ the following expression:

$${\langle \overline{\psi }\psi \rangle }_{\mathrm{heavy}}=-\frac{4N}{{\left(8\pi G\right)}^3}·\frac{1}{{\left(4\pi \right)}^2}·\frac{{\pi }^3}{3}·m{\int }_0^{\infty }ds{\mathrm{e}}^{-m^{2}s}{\int }_0^{1/s}dB{B}^7{\mathrm{e}}^{-\frac{B^2}{8G}}.$$

Because of the factor ${\mathrm{e}}^{-m^{2}s}$, the essential values of proper time here are $s\lesssim \frac{1}{m^2}$, so that $\frac{1}{s}\gtrsim m^2$. For this reason, in the leading large-m approximation, one can replace the upper limit of $1/s$ in the last integral by $+\infty$, thereby decoupling the s- and the B-integrals from each other. The B-integral then reads

$${\int }_0^{\infty }dB{B}^7{\mathrm{e}}^{-\frac{B^2}{8G}}=12288G^4,$$

which leads to the expression

$${\langle \overline{\psi }\psi \rangle }_{\mathrm{heavy}}=-\frac{2NG}{{\pi }^{2}m}=-\frac{\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle }{48{\pi }^{2}m}.$$

(12)

Thus, we have recovered the heavy-quark condensate

$${\langle \overline{\psi }\psi \rangle }_{\mathrm{heavy}}{|}_{N=3}=-\frac{{\alpha }_s\langle {\left(F_{\mu \nu }^a\right)}^2\rangle }{12\pi m},\mathrm{where}{\alpha }_s=\frac{g^2}{4\pi },$$

which was obtained in Refs. [5,6] through a direct calculation of the quark polarization operator in a constant gauge field. We see that the “area-squared” law (6) for the non-perturbative contribution to the average of a small-sized Wilson loop is consistent with the large-mass limit of the quark condensate. Note also that, for arbitrary number of colors $N>1$, one has $\langle {\left(F_{\mu \nu }^a\right)}^2\rangle =\mathcal{O}\left(N^2\right)$, while $g^2=\mathcal{O}\left(N^{-1}\right)$, so that ${\langle \overline{\psi }\psi \rangle }_{\mathrm{heavy}}=\mathcal{O}\left(N\right)$.

Also in the heavy-quark limit, one can calculate the mixed quark-gluon condensate [32] $\langle g\overline{\psi }{\mathsf{\Sigma }}_{\mu \nu }{F}_{\mu \nu }\psi \rangle$, where ${\mathsf{\Sigma }}_{\mu \nu }=\frac{1}{4i}[{\gamma }_{\mu },{\gamma }_{\nu }]$. As follows from its definition, this condensate measures the average interaction of the color-magnetic moment of a quark with the vacuum gluonic fields. Much as the usual quark condensate, the mixed condensate plays an important role in the QCD sum rules [33,34,35,36,37]. In order to calculate the mixed condensate, one can make use of its world-line representation [38]:

$${\langle g\overline{\psi }{\mathsf{\Sigma }}_{\mu \nu }{F}_{\mu \nu }\psi \rangle }_{\mathrm{heavy}}=\frac{4g^2}{V}\left({\gamma }_{\mu }\frac{\partial }{\partial x_{\mu }\left(0\right)}-m\right)\times$$

$$\times {\int }_0^{\infty }ds{\mathrm{e}}^{-m^{2}s}{\int }_P^{}\mathcal{D}{x}_{\mu }{\int }_A^{}\mathcal{D}{\psi }_{\mu }{\mathrm{e}}^{-{\int }_0^{s}d\tau \left(\frac{1}{4}{\dot{x}}_{\mu }^2+\frac{1}{2}{\psi }_{\mu }{\dot{\psi }}_{\mu }\right)}{\mathsf{\Sigma }}_{\mu \nu }{\mathsf{\Sigma }}_{\lambda \rho }\times$$

$$\times {\int }_0^{s}d\tau \mathrm{tr}\langle F_{\mu \nu }\left(x\left(0\right)\right)F_{\lambda \rho }\left(x\left(\tau \right)\right)\rangle 〈\mathcal{P}exp\left[ig{\int }_0^{s}d\tau T^a\left(A_{\mu }^a{\dot{x}}_{\mu }-{\psi }_{\mu }{\psi }_{\nu }{F}_{\mu \nu }^a\right)\right]〉,$$

(13)

where V denotes the four-volume occupied by the system. Note that Equation (13) differs from Equation (9) by an additional global factor of 2. This factor stems from the expression for the mixed condensate in terms of a variational derivative of the averaged closed-path quark propagator with respect to the coupling g, where the latter is temporarily made x-dependent (cf. Ref. [38]):

$${\langle g\overline{\psi }{\mathsf{\Sigma }}_{\mu \nu }{F}_{\mu \nu }\psi \rangle }_{\mathrm{heavy}}=\mathrm{tr}\langle g\left(x\right){\mathsf{\Sigma }}_{\mu \nu }{F}_{\mu \nu }\left(x\right)S(x,x|A_{\mu }^a)\rangle =2\mathrm{tr}\frac{\delta }{\delta g\left(0\right)}\langle S(x,x|A_{\mu }^a)\rangle ,$$

$$\mathrm{where}S(x,y|A_{\mu }^a)={(m+{\gamma }_{\mu }{D}_{\mu })}_{x,y}^{-1},D_{\mu }\equiv {\partial }_{\mu }-ig{A}_{\mu }^a{T}^a.$$

Furthermore, similarly to Equation (11), the entire dependence of Equation (13) on the gauge field can be reduced to that of the Wilson-loop average. This can be done by virtue of the area-derivative operator as follows:

$$\mathrm{tr}\langle F_{\mu \nu }\left(x\left(0\right)\right)F_{\lambda \rho }\left(x\left(\tau \right)\right)\rangle 〈\mathcal{P}exp\left[ig{\int }_0^{s}d\tau T^a\left(A_{\mu }^a{\dot{x}}_{\mu }-{\psi }_{\mu }{\psi }_{\nu }{F}_{\mu \nu }^a\right)\right]〉=$$

$$=-\frac{4}{g^2}\frac{{\delta }^2}{\delta {\sigma }_{\mu \nu }\left(x\left(0\right)\right)\delta {\sigma }_{\lambda \rho }\left(x\left(\tau \right)\right)}exp\left(-2{\int }_0^{s}d\tau {\psi }_{\mu }{\psi }_{\nu }\frac{\delta }{\delta {\sigma }_{\mu \nu }\left(x\left(\tau \right)\right)}\right)\langle W\left[x_{\mu }\right]\rangle .$$

(14)

Next, one represents the coordinate $x_{\mu }\left(\tau \right)$ as $x_{\mu }\left(\tau \right)=y_{\mu }+z_{\mu }\left(\tau \right)$, where the vector $y_{\mu }=\frac{1}{s}{\int }_0^{s}d\tau x_{\mu }\left(\tau \right)$ describes the position of a trajectory, while the vector-function $z_{\mu }\left(\tau \right)$ describes its shape (cf. Refs. [19,20,21,22,23,24,25]). The $d^{4}y$-integration yields then the factor of V, which cancels with $1/V$ on the right-hand side of Equation (13). Furthermore, the Wilson-loop average $\langle W\left[z_{\mu }\right]\rangle$ can be again represented through Equation (8) as an integral over an auxiliary antisymmetric-tensor field $B_{\mu \nu }$. The variational derivatives in Equation (14) recover then this field, so that the emerging $B_{\mu \nu }$-integral reads

$$\frac{1}{{\left(8\pi G\right)}^3}\left(\prod _{\mu <\nu }^{}{\int }_{-\infty }^{+\infty }d{B}_{\mu \nu }{\mathrm{e}}^{-\frac{B_{\mu \nu }^2}{8G}}\right)B_{\mu \nu }{B}_{\lambda \rho }\sum _{\alpha <\beta }^{}{B}_{\alpha \beta }^2=128G^2({\delta }_{\mu \lambda }{\delta }_{\nu \rho }-{\delta }_{\mu \rho }{\delta }_{\nu \lambda }).$$

The resulting expression for the mixed condensate in terms of the gluon condensate has the form [32]

$${\langle g\overline{\psi }{\mathsf{\Sigma }}_{\mu \nu }{F}_{\mu \nu }\psi \rangle }_{\mathrm{heavy}}=-\frac{1}{144{\pi }^{2}N}\frac{{\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle }^2}{m^3}.$$

As one can readily see, ${\langle g\overline{\psi }{\mathsf{\Sigma }}_{\mu \nu }{F}_{\mu \nu }\psi \rangle }_{\mathrm{heavy}}$ scales with the number of colors as $\mathcal{O}\left(N\right)$, which is the same scaling as in the case of ${\langle \overline{\psi }\psi \rangle }_{\mathrm{heavy}}$. Thus, the ratio of the two condensates,

$$\frac{{\langle g\overline{\psi }{\mathsf{\Sigma }}_{\mu \nu }{F}_{\mu \nu }\psi \rangle }_{\mathrm{heavy}}}{{\langle \overline{\psi }\psi \rangle }_{\mathrm{heavy}}}=\frac{\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle }{3N{m}^2},$$

becomes N-independent in the large-N limit.

With the decrease of the current quark mass m downwards $1/\lambda$, variations of the gauge field inside the quark trajectory produce corrections to Equation (12). We will now follow Ref. [39] to calculate these corrections by using the approach proposed in Ref. [22]. (For the bosonic case, see [40].) An advantageous feature of this approach is that it provides a closed formula for the effective action of a fermion moving in an arbitrary Abelian gauge field. Furthermore, it is known that, in addition to the confining interactions of stochastic gluonic fields, there also exist non-confining non-perturbative interactions of those fields [41]. Below, we briefly address possible influence of such interactions on the heavy-quark condensate, and notice that, for the purely exponential two-point correlation function of gluonic field strengths, the heavy-quark condensate does not depend on the non-confining non-perturbative interactions.

We start our analysis with considering confining interactions of the stochastic gluonic fields. Within the stochastic vacuum model [11,12,13,14,15,16], the corresponding part of the Wilson-loop average is given by the following area-dependent expression:

$$\langle W\left(C\right)\rangle \simeq exp\left[-\frac{\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle }{96N}{\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x\right){\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x^{\prime }\right){\mathrm{e}}^{-|x-x^{\prime }|/\lambda }\right].$$

(15)

This expression generalizes Equation (5) to the case where the mean size of $S_{\min }$ exceeds the vacuum correlation length λ. We choose the surface element $d{\sigma }_{\mu \nu }$ in the form of an oriented, infinitely thin triangle built up of the position vector $z_{\mu }\left(\tau \right)$ and the differential element $d{z}_{\mu }={\dot{z}}_{\mu }d\tau$, namely $d{\sigma }_{\mu \nu }\left(z\right)=\frac{1}{2}(z_{\mu }{\dot{z}}_{\nu }-z_{\nu }{\dot{z}}_{\mu })d\tau$. We further apply an elementary Fourier transform ${\int }_x{\mathrm{e}}^{-\mu \left|x\right|+ipx}=\frac{12{\pi }^2\mu }{{(p^2+{\mu }^2)}^{5/2}}$, where from now on $\mu \equiv \frac{1}{\lambda }$, and ${\int }_x\equiv \int d^{4}x$. This yields for the Wilson-loop average the following expression:

$$\langle W\left(C\right)\rangle =\int \left[\prod _{\mu <\nu }^{}\mathcal{D}{B}_{\mu \nu }{\mathrm{e}}^{-\frac{N}{{\pi }^2\mu \langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle }{\int }_x{B}_{\mu \nu }{(-{\partial }^2+{\mu }^2)}^{5/2}{B}_{\mu \nu }}\right]{\mathrm{e}}^{\frac{i}{2}{\int }_0^{s}d\tau B_{\mu \nu }\left(z\right)z_{\mu }{\dot{z}}_{\nu }}.$$

(16)

One recognizes in the latter exponential factor a Wilson loop corresponding to the Abelian gauge field

$$A_{\nu }\left(z\right)=\frac{1}{2}{z}_{\mu }{B}_{\mu \nu }\left(z\right).$$

(17)

The strength tensor of this field is $F_{\mu \nu }={\partial }_{\mu }{A}_{\nu }-{\partial }_{\nu }{A}_{\mu }=B_{\mu \nu }+C_{\mu \nu }$, where $C_{\mu \nu }\left(z\right)=\frac{1}{2}{z}_{\lambda }({\partial }_{\mu }{B}_{\lambda \nu }-{\partial }_{\nu }{B}_{\lambda \mu })$. Furthermore, owing to the Abelian Stokes’ theorem, it is just the strength tensor $F_{\mu \nu }$ which again automatically appears in the quark spin term of the effective action, being recovered by the operator $\frac{\delta }{\delta {\sigma }_{\mu \nu }}$ in Equation (11). We can now apply the aforementioned closed formula for the effective action of a fermion moving in an arbitrary Abelian field [22]. It yields

$$\langle \mathsf{\Gamma }\left[A_{\mu }^a\right]\rangle =2N{\int }_0^{\infty }\frac{ds}{s}\frac{{\mathrm{e}}^{-m^{2}s}}{{\left(4\pi \right)}^2}{〈{\int }_x{F}_{\mu \nu }\left(x\right)\mathsf{F}\left(\xi \right)F_{\mu \nu }\left(x\right)〉}_B.$$

(18)

Here, the factor of volume V has not yet been divided out, the average ${\langle \cdots \rangle }_B$ is defined by the square bracket in Equation (16), and the Abelian covariant derivative entering the formfactor $\mathsf{F}\left(\xi \right)$ reads $D_{\mu }={\partial }_{\mu }-i{A}_{\mu }$. In the spinor case at issue, this formfactor has the form [22] $\mathsf{F}\left(\xi \right)=\frac{f\left(\xi \right)-1}{2\xi }-\frac{1}{4}f\left(\xi \right)$, where $f\left(\xi \right)={\int }_0^{1}du{\mathrm{e}}^{u(1-u)\xi }$ and $\xi =s{D}_{\mu }^2$. For what follows, we find it convenient to get rid of ξ in the denominator of $\mathsf{F}\left(\xi \right)$, by representing this formfactor as

$$\mathsf{F}\left(\xi \right)=\frac{1}{2}{\int }_0^{1}du\left[u(1-u){\int }_0^{1}d\alpha {\mathrm{e}}^{\alpha u(1-u)\xi }-\frac{1}{2}{\mathrm{e}}^{u(1-u)\xi }\right].$$

To calculate the thus emerging averages of the form ${\langle {\int }_x{F}_{\mu \nu }\left(x\right){\mathrm{e}}^{\alpha u(1-u)\xi }{F}_{\mu \nu }\left(x\right)\rangle }_B$, let us consider the following simplified expression: ${\int }_x{\langle F_{\mu \nu }\left(x\right){\mathrm{e}}^{s{D}_{\mu }^2}{F}_{\mu \nu }\left(x\right)\rangle }_B$. Using the path-integral representation for the operator ${\mathrm{e}}^{s{D}_{\mu }^2}$, we have

$${\int }_x{\langle F_{\mu \nu }\left(x\right){\mathrm{e}}^{s{D}_{\mu }^2}{F}_{\mu \nu }\left(x\right)\rangle }_B={\int }_x{\int }_y\underset{r\left(0\right)=x}{\overset{r\left(s\right)=x+y}{\int }}\mathcal{D}{r}_{\mu }{\mathrm{e}}^{-\frac{1}{4}{\int }_0^{s}d\tau {\dot{r}}_{\mu }^2}{\langle F_{\mu \nu }\left(x\right)exp\left(i{\int }_x^{x+y}d{r}_{\mu }{\mathcal{A}}_{\mu }\right)F_{\mu \nu }(x+y)\rangle }_B.$$

Next, owing to the translation invariance of the B-average, we have

$${\langle F_{\mu \nu }\left(x\right)exp\left(i{\int }_x^{x+y}d{r}_{\mu }{\mathcal{A}}_{\mu }\right)F_{\mu \nu }(x+y)\rangle }_B={\langle F_{\mu \nu }\left(0\right)exp\left(i{\int }_0^{y}d{r}_{\mu }{\mathcal{A}}_{\mu }\right)F_{\mu \nu }\left(y\right)\rangle }_B.$$

Shifting now trajectories $r_{\mu }\left(\tau \right)$ by the vector $-x$, we obtain

$${\int }_x{\langle F_{\mu \nu }\left(x\right){\mathrm{e}}^{s{D}_{\mu }^2}{F}_{\mu \nu }\left(x\right)\rangle }_B=V{\int }_y\underset{r\left(0\right)=0}{\overset{r\left(s\right)=y}{\int }}\mathcal{D}{r}_{\mu }{\mathrm{e}}^{-\frac{1}{4}{\int }_0^{s}d\tau {\dot{r}}_{\mu }^2}{\langle F_{\mu \nu }\left(0\right)exp\left(i{\int }_0^{y}d{r}_{\mu }{\mathcal{A}}_{\mu }\right)F_{\mu \nu }\left(y\right)\rangle }_B,$$

where V again denotes the 4-volume occupied by the system. Furthermore, to stay within the initial two-point approximation, the phase factor $exp\left(i{\int }_0^{y}d{r}_{\mu }{\mathcal{A}}_{\mu }\right)$ should be approximated by unity. Indeed, the use of the formfactor $\mathsf{F}\left(\xi \right)$ corresponds to accounting for only two $F_{\mu \nu }$’s, while the Taylor expansion of the phase factor $exp\left(i{\int }_0^{y}d{r}_{\mu }{\mathcal{A}}_{\mu }\right)$ would yield correlation functions of more than two $F_{\mu \nu }$’s. The path integrals over $r_{\mu }\left(\tau \right)$ get then reduced to the Green function of the heat equation in 4D, and we obtain for the B-average from Equation (18):

$${〈{\int }_x{F}_{\mu \nu }\left(x\right)\mathsf{F}\left(\xi \right)F_{\mu \nu }\left(x\right)〉}_B=\frac{V}{2{\left(4\pi s\right)}^2}{\int }_0^{1}du{\int }_y\left(\frac{4s}{y^2}-\frac{1}{2{\left[u(1-u)\right]}^2}\right){\mathrm{e}}^{-\frac{y^2}{4u(1-u)s}}{\langle F_{\mu \nu }\left(0\right)F_{\mu \nu }\left(y\right)\rangle }_B,$$

(19)

where the elementary α-integration has been performed. A straightforward calculation of the correlation function ${\langle F_{\mu \nu }\left(0\right)F_{\mu \nu }\left(y\right)\rangle }_B$ yields further the condensate $\langle \overline{\psi }\psi \rangle =-\frac{1}{V}\frac{\partial }{\partial m}\langle \mathsf{\Gamma }\left[A_{\mu }^a\right]\rangle$ for the quark masses m down to $m=\mu$. Referring the reader for details to Ref. [39], we present here the final result for the corresponding ratio

$$J\left(r\right)\equiv \frac{\langle \overline{\psi }\psi \rangle }{{\langle \overline{\psi }\psi \rangle }_{\mathrm{heavy}}},$$

(20)

where ${\langle \overline{\psi }\psi \rangle }_{\mathrm{heavy}}$ is given by Equation (12), and $r\equiv \frac{m}{\mu }$. It reads

$$J\left(r\right)=\frac{3r^2}{4}{\int }_0^1\frac{du}{1-a^2}\left\{4+{\left(\frac{a}{r}\right)}^2·\frac{2a^2+1}{1-a^2}+\frac{3}{a^2-1}-{\left(\frac{a}{2r}\right)}^2·\frac{13a^2+2}{{(a^2-1)}^2}+\right.$$

$$\left.+\frac{arccos(1/a)}{{(a^2-1)}^{3/2}}\left[\frac{3a^4}{r^2}-5a^2+2+{\left(\frac{a}{2r}\right)}^2·\frac{3a^2(a^2+4)}{a^2-1}\right]\right\},$$

where $a\equiv \frac{r}{\sqrt{u(1-u)}}$. For $r\gg 1$, the leading large-r terms $4+{\left(\frac{a}{r}\right)}^2·\frac{2a^2+1}{1-a^2}$ yield $J\left(r\right)\to 1$. For arbitrary r’s, the remaining u-integration can be done numerically (cf. Ref. [39]). In particular, we obtain $J\left(1\right)\simeq 0.36$, which means a 64%-decrease in the value of the quark condensate for $m\simeq \mu$.

Physicswise, only the values of $J\left(r\right)$ corresponding to $m=m_c$, $m_b$, and $m_t$ are relevant. We use the standard heavy-quark masses $m_c\simeq 1.3\mathrm{GeV}$, $m_b\simeq 4.2\mathrm{GeV}$, and $m_t\simeq 173\mathrm{GeV}$. The vacuum correlation length in full QCD with light flavors [42], $\lambda \simeq 0.34\mathrm{fm}$, corresponds to $\mu \simeq 580\mathrm{MeV}$. This yields

$$J(m_c/\mu )\simeq 0.60,J(m_b/\mu )\simeq 0.84,J(m_t/\mu )\simeq 0.996\mathrm{in}\mathrm{full}\mathrm{QCD}.$$

(21)

For the alternative case of quenched QCD, i.e., the SU(3) Yang-Mills theory, the vacuum correlation length is [43,44,45] $\lambda \simeq 0.22\mathrm{fm}$, which corresponds to $\mu \simeq 897\mathrm{MeV}$. For this value of μ, we have

$$J(m_c/\mu )\simeq 0.47,J(m_b/\mu )\simeq 0.77,J(m_t/\mu )\simeq 0.993\mathrm{in}\mathrm{quenched}\mathrm{QCD}.$$

(22)

The obtained decrease of $J\left(r\right)$ with the decrease of m shows that, if one keeps calculating the quark condensate via Equation (12) while decreasing m, the value of m to be used in that equation should be larger than just the current quark mass. Furthermore, the sets of values (21) and (22) illustrate the degree of accuracy of Equation (12) for various heavy flavors and various values of the vacuum correlation length λ. Since the case of heavy quarks considered here is intermediate between QCD with light quarks and quenched QCD, the genuine value of $J(m_{\mathrm{f}}/\mu )$, for a given heavy flavor f, lies somewhere in between the two corresponding values of $J(m_{\mathrm{f}}/\mu )$ given by Equations (21) and (22). In any case, we can conclude that Equation (12) is inapplicable to the c-quark, since it can develop up to 53%-corrections [cf. $J(m_c/\mu )$ from Equation (22)].

In addition to the confining interactions of the stochastic background Yang-Mills fields, which lead to the Wilson-loop average in the form of Equation (15), there also exist non-confining non-perturbative interactions of those fields. To account for the non-confining non-perturbative interactions, one represents the full two-point correlation function of gluonic field strengths as [11,12,13,14,15,16,41]

$$\langle g^2{F}_{\mu \nu }^a\left(0\right)F_{\lambda \rho }^b\left(x\right)\rangle =\frac{\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle }{12}·\frac{{\delta }^{ab}}{N^2-1}·\{\kappa ({\delta }_{\mu \lambda }{\delta }_{\nu \rho }-{\delta }_{\mu \rho }{\delta }_{\nu \lambda })+$$

$$+\frac{1-\kappa }{2}\left[{\partial }_{\mu }(x_{\lambda }{\delta }_{\nu \rho }-x_{\rho }{\delta }_{\nu \lambda })+{\partial }_{\nu }(x_{\rho }{\delta }_{\mu \lambda }-x_{\lambda }{\delta }_{\mu \rho })\right]\}{\mathrm{e}}^{-\mu \left|x\right|}.$$

(23)

Here, $\kappa \in [0,1]$ is some parameter, which defines the relative strength of the confining and the non-confining non-perturbative interactions. Lattice simulations in the SU(3) Yang-Mills theory yield the value of $\kappa =0.83\pm 0.03$ (cf. Ref. [41]), which means that the relative contribution of the non-confining non-perturbative interactions amounts to at most 20%. Expressing the Wilson-loop average via the correlation function $\langle g^2{F}_{\mu \nu }^a\left(0\right)F_{\lambda \rho }^b\left(x\right)\rangle$ by means of the non-Abelian Stokes’ theorem and the cumulant expansion, and using the parametrization (23), one obtains the following generalization of Equation (15) (Note that the non-confining non-perturbative interactions produce in the Wilson-loop average a term with the double contour integral, which initially has the form $\frac{1-\kappa }{2}{\oint }_{C}d{x}_{\mu }{\oint }_{C}d{x}_{\mu }^{\prime }\underset{{(x-x^{\prime })}^2}{\overset{\infty }{\int }}d\tau {\mathrm{e}}^{-\mu \sqrt{\tau }}$. The corresponding expression in Equation (24) resulted from the τ-integration in this formula.):

$$\langle W\left(C\right)\rangle =exp\{-\frac{\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle }{96N}\times$$

$$\times \left[\kappa {\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x\right){\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x^{\prime }\right)+\frac{1-\kappa }{{\mu }^2}{\oint }_{C}d{x}_{\mu }{\oint }_{C}d{x}_{\mu }^{\prime }(1+\mu |x-x^{\prime }\left|\right)\right]{\mathrm{e}}^{-\mu |x-x^{\prime }|}\}.$$

(24)

This expression explicitly demonstrates that the term $\propto (1-\kappa )$ in Equation (23) corresponds to the non-confining interactions of the background fields. Furthermore, as has been demonstrated in Ref. [39], for the case of exponential correlations of such interactions, they do not affect the value of the quark condensate. The question of whether this result is specific to the exponential ansatz for the correlation function $\langle g^2{F}_{\mu \nu }^a\left(0\right)F_{\lambda \rho }^b\left(x\right)\rangle$, or it holds for some other ansätze (such as, e.g., the Gaussian one) equally well, requires a separate study.

Altogether, the above-presented analysis shows that Equation (12) applies with a good accuracy only to the t-quark. Instead, the corrections are ∼20% for the b-quark, while for the c-quark they can reach $50\%$, thereby making Equation (12) in this case inapplicable. As we have also seen, once the continuously varied current quark mass m reaches the value of the inverse vacuum correlation length μ, the absolute value of the quark condensate decreases by 64% compared to the value provided by Equation (12). Furthermore, we have used in our calculations the most general ansatz for the Wilson-loop average, which is provided by the stochastic vacuum model, and accounts for the confining, as well as the non-perturbative non-confining, interactions of the stochastic gluonic fields. In particular, for the most simple, exponential, parametrization of the two-point correlation function of gluonic field strengths, the heavy-quark condensate turns out to be independent of the non-confining non-perturbative interactions of the stochastic gluonic fields. Furthermore, we have not considered s-, d-, and u-quarks, whose masses are all smaller than μ. The reason for this restriction is that, for these light quarks, the effect of spontaneous breaking of chiral symmetry starts to play an important role, resulting in the appearance of a significant self-energy contribution to the dynamical constituent quark mass. Thus, since such a self-energy contribution cannot be consistently calculated within the present approach, we had to consider only the heavy quarks, in which case this contribution can be safely disregarded in comparison with the current quark mass. Nevertheless, as we have just seen, even in the case of such heavy quarks as b and c, the corrections to Equation (12) are substantial.

Digression: Casimir Scaling of k-String Tensions

In any SU(N) gauge theory, complex conjugation is known to define a representation conjugated to the fundamental one. For $N\ge 3$, the fundamental representation, under which quarks are transformed, and the representation conjugated to it, under which antiquarks are transformed, are mutually independent. Objects transforming under higher representations of the group SU(N) can be constructed from the direct products of the quark and the antiquark fields. The difference between the number of these quark and antiquark fields modulo N is called the $\mathcal{N}$-ality of a given representation of SU(N), and is characterized by an integer $k=1,\dots ,N-1$. Thus, representations with $\mathcal{N}$-alities $k\le N/2$ and $N-k$ are related to each other via the complex conjugation, which corresponds to the replacement of quarks by antiquarks and vice versa. Consequently, confining strings associated with these representations have equal tensions. Note further that objects transforming under a representation of $\mathcal{N}$-ality k carry charge k with respect to the center group $Z_N$ of SU(N). Therefore, since the $Z_N$ symmetry is unbroken in the hadronic phase, states with different $\mathcal{N}$-alities do not mix with each other. Accordingly, objects transforming trivially under the center group $Z_N$ define a representation of the zero $\mathcal{N}$-ality. Each $\mathcal{N}$-ality-zero representation is contained in the tensor product of some number of adjoint representations. A Wilson loop $W\left(C\right)$ which describes a static object transforming under an $\mathcal{N}$-ality-zero representation possesses the vacuum expectation value $\langle W\left(C\right)\rangle$ which does not exhibit an area law for large contours C. The reason for this is the fact that such an object can be screened by gluons, since those transform under the adjoint representation. Instead, representations which are important for confinement are given by rank-k antisymmetric tensors. The eigenvalue of the quadratic Casimir operator of a rank-k antisymmetric representation reads $C_2^{\left(k\right)}=\frac{k(N-k)(N+1)}{2N}$. In particular, rank-1 antisymmetric representation is just the fundamental representation of the group SU(N). In general, any representation of a non-zero $\mathcal{N}$-ality is contained in a direct product of a certain rank-k antisymmetric representation and some number of adjoint representations. Accordingly, any color source belongs to one of N classes, each of which is characterized by the corresponding rank-k antisymmetric representation. Therefore, for all possible representations of the color source, there exist only N string tensions ${\sigma }_k$’s to characterize confinement. Of those, ${\sigma }_1$ is the string tension corresponding to the fundamental representation, while ${\sigma }_N=0$. The quantity ${\sigma }_k$ can be interpreted as a tension of the so-called k-string, which is a confining string interconnecting k quarks with k antiquarks. As was already discussed, the equality ${\sigma }_k={\sigma }_{N-k}$ takes place, owing to which only [$N/2$] of all string tensions ${\sigma }_k$’s are mutually independent. Thus, for any SU(N) gauge theory, the dynamics of confinement is encoded in these [$N/2$] numbers.

Clearly, a k-string can only be stable provided ${\sigma }_k<k{\sigma }_1$. In the large-N limit, interactions between strings composing a k-string are suppressed, so that ${\sigma }_k\to k{\sigma }_1$ at $N\to \infty$ and a fixed k. In the 2D Yang-Mills theory, one can prove that ${\sigma }_k\propto C_2^{\left(k\right)}$, which is a consequence of the fact that confinement in 1+1 dimensions stems from the one-gluon exchange between the static sources. For this reason, the ratio $\frac{{\sigma }_k}{{\sigma }_1}$ in the 2D Yang-Mills theory obeys exactly the so-called Casimir-scaling formula [46], $\frac{{\sigma }_k}{{\sigma }_1}=\frac{k(N-k)}{N-1}$, i.e., indeed $\frac{{\sigma }_k}{{\sigma }_1}<k$. Nevertheless, in the 4D SU(N) Yang-Mills theory with $N=4$ and $N=6$, lattice data [47,48,49] on the $\frac{{\sigma }_k}{{\sigma }_1}$-ratio show that corrections to the Casimir scaling are of the order of $10\%$, while corrections to the so-called Sine scaling [50,51,52,53], $\frac{{\sigma }_k}{{\sigma }_1}=\frac{sin\frac{\pi k}{N}}{sin\frac{\pi }{N}}$, amount to only a few percent. The Sine-scaling ratio has been found analytically in supersymmetric gauge theories [50], as well as through a possible duality between such gauge theories and string theories [51,52,53,54], but not in the 4D Yang-Mills theory itself. The principal difference of the Sine scaling from the Casimir scaling is reflected in the N-dependence of the leading correction to the large-N limit $\frac{{\sigma }_k}{{\sigma }_1}\to k$. Namely, for the Sine scaling, this correction reads $\frac{{\pi }^2}{6}\frac{k(1-k^2)}{N^2}$, being therefore $\mathcal{O}(\frac{1}{N^2})$, while in the case of the Casimir scaling it reads $\frac{k(1-k)}{N}$, thereby behaving with N as $\mathcal{O}(\frac{1}{N})$. Physicswise, this correction yields the strength of pairwise attractive interactions between the k strings that constitute a k-string. This is the reason why the parametric N-dependence of such a leading correction to the large-N limit of the k-string tension is important. However, the current level of accuracy of lattice simulations does not allow one to unambiguously decide in favor of either the $\mathcal{O}(\frac{1}{N^2})$- or the $\mathcal{O}(\frac{1}{N})$-behavior of this correction. On the theory side too, there is no reason to expect either the Sine or the Casimir scaling to be an exact result in the 4D non-supersymmetric Yang-Mills theory. Nevertheless, as has been shown in Ref. [55], it is the Casimir scaling which takes place as the leading result in the realistic 4D [U(1)]${}^{N-1}$-invariant dual Abelian Higgs model of confinement.

Here, we demonstrate that it is the Casimir scaling which takes place for the k-string tensions in the Gaussian ensemble of the stochastic Yang-Mills fields at $k\sim N$. To this end, we generalize Equation (15) to the case of the k-th power of the Wilson loop in the fundamental representation. This will allow us to obtain the k-string tension ${\sigma }_k$ in the approximation where only the two-point correlation functions of gluonic field strengths are taken into account (i.e., the ensemble of the stochastic Yang-Mills fields can be regarded as Gaussian). Equation (15) can be obtained by expressing the Wilson-loop average via the non-Abelian Stokes’ theorem and the cumulant expansion truncated at its second-order term, where the latter is parametrized as (Cf. the generalization (23) of this parametrization to the case where the non-confining non-perturbative interactions of the stochastic Yang-Mills fields are taken into account. This generalization is provided through the parameter κ, whose value can differ from 1. The corresponding generalization of the Wilson-loop average (15) is given by Equation (24).)

$$\langle F_{\mu \nu }^a\left(x\right)F_{\lambda \rho }^b\left(x^{\prime }\right)\rangle =\frac{{\delta }^{ab}}{N^2-1}·\frac{1}{12}\left({\delta }_{\mu \lambda }{\delta }_{\nu \rho }-{\delta }_{\mu \rho }{\delta }_{\nu \lambda }\right)·\langle {(F_{\alpha \beta }^a)}^2\rangle {\mathrm{e}}^{-|x-x^{\prime }|/\lambda },$$

(25)

and the normalization factors in this equation are obvious. Indeed, substituting Equation (25) into the corresponding expression for the Wilson-loop average,

$$\langle W\left(C\right)\rangle \simeq \frac{1}{N}\mathrm{tr}exp\left[-\frac{g^2}{8}{T}^a{T}^b{\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x\right){\int }_{S_{\min }}d{\sigma }_{\lambda \rho }\left(x^{\prime }\right)\langle F_{\mu \nu }^a\left(x\right)F_{\lambda \rho }^b\left(x^{\prime }\right)\rangle \right],$$

and using the relation (Clearly, for some representation r other than the fundamental one, the factor of $\frac{N^2-1}{2N}$ in this relation gets replaced by the corresponding quadratic Casimir operator $C_r$.) $T^a{T}^a=\frac{N^2-1}{2N}·\widehat{1}$, we recover Equation (15). In order to calculate $\langle {\left(W\left(C\right)\right)}^k\rangle$, we introduce the short-hand notations $\mathrm{Tr}\equiv \frac{1}{N}\mathrm{tr}$ and ${\widehat{f}}_j\equiv \frac{ig}{2}{\int }_{S}d{\sigma }_{\mu \nu }{F}_{\mu \nu }^a{T}_j^a$, where the index j numerates the copy of the fundamental representation. We have

$$\langle {\left(W\left(C\right)\right)}^k\rangle ={\mathrm{Tr}}_1\cdots {\mathrm{Tr}}_k{\mathcal{P}}_S\langle {\mathrm{e}}^{{\widehat{f}}_1}\otimes {\widehat{1}}_2\otimes \cdots \otimes {\widehat{1}}_k·{\widehat{1}}_1\otimes {\mathrm{e}}^{{\widehat{f}}_2}\otimes {\widehat{1}}_3\otimes \cdots \otimes {\widehat{1}}_k\cdots {\widehat{1}}_1\otimes {\widehat{1}}_2\otimes \cdots \otimes {\widehat{1}}_{k-1}\otimes {\mathrm{e}}^{{\widehat{f}}_k}\rangle ,$$

where ⊗ denotes the tensor product of matrices belonging to different copies, and ${\mathcal{P}}_S$ denotes the surface ordering. The unit matrices can further be exponentiated as

$$\langle {\left(W\left(C\right)\right)}^k\rangle ={\mathrm{Tr}}_1\cdots {\mathrm{Tr}}_k{\mathcal{P}}_S\langle {\mathrm{e}}^{{\widehat{f}}_1\otimes {\widehat{1}}_2\otimes \cdots \otimes {\widehat{1}}_k}·{\mathrm{e}}^{{\widehat{1}}_1\otimes {\widehat{f}}_2\otimes \cdots \otimes {\widehat{1}}_k}\cdots {\mathrm{e}}^{{\widehat{1}}_1\otimes {\widehat{1}}_2\otimes \cdots \otimes {\widehat{f}}_k}\rangle .$$

Performing now the average in the Gaussian approximation, we have

$$\langle {\left(W\left(C\right)\right)}^k\rangle ={\mathrm{Tr}}_1\cdots {\mathrm{Tr}}_{k}exp\left\{\frac{1}{2}\left[\langle {\widehat{f}}_1·{\widehat{f}}_1\otimes {\widehat{1}}_2\otimes \cdots \otimes {\widehat{1}}_k\rangle +\langle {\widehat{1}}_1\otimes {\widehat{f}}_2·{\widehat{f}}_2\otimes {\widehat{1}}_3\otimes \cdots \otimes {\widehat{1}}_k\rangle +\cdots +\right.\right.$$

$$+\langle {\widehat{1}}_1\otimes {\widehat{1}}_2\otimes \cdots \otimes {\widehat{f}}_k·{\widehat{f}}_k\rangle +2\left(\langle {\widehat{f}}_1\otimes {\widehat{f}}_2\otimes {\widehat{1}}_3\otimes \cdots \otimes {\widehat{1}}_k\rangle +\langle {\widehat{f}}_1\otimes {\widehat{1}}_2\otimes {\widehat{f}}_3\otimes \cdots \otimes {\widehat{1}}_k\rangle +\cdots +\right.$$

$$\left.\left.\left.+\langle {\widehat{f}}_1\otimes {\widehat{1}}_2\otimes {\widehat{1}}_3\otimes \cdots \otimes {\widehat{f}}_k\rangle +\cdots +\langle {\widehat{1}}_1\otimes \cdots \otimes {\widehat{1}}_{k-2}\otimes {\widehat{f}}_{k-1}\otimes {\widehat{f}}_k\rangle \right)\right]\right\}.$$

(26)

Here, for every j, averages of the type $\langle {\widehat{f}}_j·{\widehat{f}}_j\rangle$ are the same as for a single Wilson loop. Indeed, by virtue of the parametrization (25), we have

$$\langle {\widehat{f}}_j·{\widehat{f}}_j\rangle =-\frac{g^2}{4}{T}^a{T}^b{\int }_{S}d{\sigma }_{\mu \nu }\left(x\right){\int }_{S}d{\sigma }_{\lambda \rho }\left(x^{\prime }\right)\langle F_{\mu \nu }^a\left(x\right)F_{\lambda \rho }^b\left(x^{\prime }\right)\rangle =-\frac{\chi }{N}·\widehat{1},$$

where we have introduced the notation $\chi \equiv \frac{1}{48}\langle {(g{F}_{\alpha \beta }^a)}^2\rangle {\int }_{S}d{\sigma }_{\mu \nu }\left(x\right){\int }_{S}d{\sigma }_{\mu \nu }\left(x^{\prime }\right){\mathrm{e}}^{-|x-x^{\prime }|/\lambda }$. Rather, for averages of the type $\langle {\widehat{f}}_i\otimes {\widehat{f}}_j\rangle$, where $i\ne j$, we have

$$\langle {\widehat{f}}_i\otimes {\widehat{f}}_j\rangle =-\frac{2\chi }{N^2-1}{T}^a\otimes T^a.$$

To handle the exponential of the matrix $T^a\otimes T^a$, we use projection operators $P_s$ and $P_a$, which decompose the direct product of two fundamental representations of SU(N) into the irreducible representations as $N\otimes N=\frac{N(N+1)}{2}\oplus \overline{\frac{N(N-1)}{2}}$. These operators obey the relation $T^a\otimes T^a=\frac{N-1}{2N}{P}_s-\frac{N+1}{2N}{P}_a$, which yields for Equation (26) the following expression:

$$\langle {\left(W\left(C\right)\right)}^k\rangle ={\mathrm{e}}^{-\frac{k}{2N}·\chi }\left[\frac{N+1}{2N}{\mathrm{e}}^{-k(k-1)·\frac{N-1}{2N}·\frac{\chi }{N^2-1}}+\frac{N-1}{2N}{\mathrm{e}}^{k(k-1)·\frac{N+1}{2N}·\frac{\chi }{N^2-1}}\right]=$$

$$={\mathrm{e}}^{-\frac{k}{2N}·\chi }\left[\frac{N+1}{2N}{\mathrm{e}}^{-\frac{k(k-1)}{2N(N+1)}\chi }+\frac{N-1}{2N}{\mathrm{e}}^{\frac{k(k-1)}{2N(N-1)}\chi }\right],$$

(27)

where we have also used the relations ${\mathrm{Tr}}_i{\mathrm{Tr}}_j{P}_s=\frac{N+1}{2N}$ and ${\mathrm{Tr}}_i{\mathrm{Tr}}_j{P}_a=\frac{N-1}{2N}$.

Note that, in a similar way, one can calculate the correlation function of k Wilson loops in the fundamental representation, defined at different contours $C_1,\dots ,C_k$. For concreteness, let us consider the correlation function of two Wilson loops, which are defined at contours $C_1$ and $C_2$ lying in parallel planes. For the likely oriented contours, Equation (27) for $k=2$ goes over to (cf. Refs. [56,57])

$$\langle W\left(C_1\right)W\left(C_2\right)\rangle ={\mathrm{e}}^{-\frac{1}{2N}({\chi }_{11}+{\chi }_{22})}\left[\frac{N+1}{2N}{\mathrm{e}}^{-\frac{{\chi }_{12}}{N(N+1)}}+\frac{N-1}{2N}{\mathrm{e}}^{\frac{{\chi }_{12}}{N(N-1)}}\right],$$

(28)

where ${\chi }_{mn}\equiv \frac{\langle {(g{F}_{\alpha \beta }^a)}^2\rangle }{48}{\int }_{S_m}d{\sigma }_{\mu \nu }\left(x\right){\int }_{S_n}d{\sigma }_{\mu \nu }\left(x^{\prime }\right){\mathrm{e}}^{-|x-x^{\prime }|/\lambda }$, and $S_1$, $S_2$ are the minimal (i.e., plane) surfaces encircled respectively by $C_1$ and $C_2$. With the increase of the separation between the planes, ${\chi }_{12}$ gradually vanishes, and the correlation function (28) gets factorized into the product of two Wilson-loop averages (15). Note that, within the above-described approach, which is based on the non-Abelian Stokes’ theorem and the cumulant expansion, the correlation function (28) is given entirely in terms of the integrals over the surfaces $S_1$ and $S_2$, rather than in terms of an integral over the common minimal surface interconnecting the contours $C_1$ and $C_2$ (e.g., a catenoid in the case of circular coaxial contours $C_1$ and $C_2$ of the same radius). For this reason, unlike the classical-physics example of soap films, in the present case one cannot define any critical separation between the contours, such that for smaller separations a common minimal surface is formed, while for separations larger than the critical one this common surface goes over into $S_1$ and $S_2$.

Rewriting now Equation (27) in an equivalent form

$$\langle {\left(W\left(C\right)\right)}^k\rangle =\frac{N+1}{2N}{\mathrm{e}}^{-\frac{k(N+k)}{2N(N+1)}\chi }+\frac{N-1}{2N}{\mathrm{e}}^{-\frac{k(N-k)}{2N(N-1)}\chi },$$

we notice that, in the limit $k\sim N\gg 1$ of interest, the first of the two exponentials on the right-hand side of this expression is suppressed in comparison with the second one. Thus, in this limit, one obtains Casimir scaling, $\frac{{\sigma }_k}{{\sigma }_1}\simeq \frac{k(N-k)}{N-1}$.

2. Electroweak Phase Transition as a Vacuum Instability

One of the outstanding problems of modern particle physics and cosmology is the description of the electroweak (EW) phase transition [58,59,60,61]. This phase transition has been found to be first-order if the Higgs mass is smaller than 80 GeV, while becoming a smooth crossover for larger values of the Higgs mass [58,59,60,61]. In this Section, we present an analytic calculation of the vacuum instability, which is caused by the “paramagnetic” interaction of the gauge boson with the stochastic background Yang-Mills fields. In this way, we show that this instability may occur at the temperature of 160 GeV, which is the modern lattice value for the critical temperature of the EW phase transition at the Higgs mass of 125 GeV [62,63]. (See also Ref. [64], where the behavior of several thermodynamic quantities across the phase transition has been found by combining 3-loop perturbative calculations with the lattice simulations. An advantageous feature of such an approach is that it allows one to avoid the infrared problems, which are present in the purely perturbative calculations around the critical temperature.) The vacuum instability is associated with the vanishing vacuum-energy density [65], and occurs due to the fact that the “paramagnetic” interaction has the sign opposite to that of the “diamagnetic” interaction (experienced also by the Higgs boson and the ghost), while being larger than the “diamagnetic” interaction in its absolute value. We perform our calculations within the effective 3D theory, which can be obtained from the 4D theory by integrating out the degrees of freedom with momenta $\sim g_{T}T$. Here $g_T$ is the weak coupling constant, as it would have been at temperature T if the 4D theory were not reduced to the 3D one. Thus, the characteristic momenta in the resulting 3D theory can reach the order of $\mathcal{O}\left(g_T^{2}T\right)$, which is however much smaller than $\mathcal{O}\left(g_{T}T\right)$. The Lagrangian of the 3D theory has the form [58,59,60,61]

$$\mathcal{L}=\frac{1}{4}{(F_{ij}^a)}^2+{\left(D_i\mathsf{\Phi }\right)}^{†}\left(D_i\mathsf{\Phi }\right)+m_3^2{\mathsf{\Phi }}^{†}\mathsf{\Phi }+{\lambda }_3{\left({\mathsf{\Phi }}^{†}\mathsf{\Phi }\right)}^2,$$

(29)

where the Yang-Mills field-strength tensor and the covariant derivative in the fundamental representation contain a gauge coupling $g_3^2=g_T^{2}T$. The mass of the gauge boson and the ghost is $m=\frac{g_3{\phi }_T}{2}$, while the mass of the Higgs boson is given by the formula $m_H^2=3{\lambda }_3{\phi }_T^2+m_3^2$. Here, the dimensionality of the parameter ${\phi }_T$ is (mass)${}^{1/2}$, and the dimensionality of the Higgs coupling ${\lambda }_3$ is (mass).

The inverse vacuum correlation length M in the theory (29) is approximately equal to the mass $\mathsf{M}$ of the so-called one-gluon gluelump, which is the bound state of a gluon in the field of a hypothetical static adjoint source. Such a bound state corresponds to a rectangular contour C, whose temporal extension is much larger than its spatial extension. The adjoint string interconnecting the gluon with the static source, breaks upon the creation of a glueball and the subsequent recombination process. This leads to the formation of two one-gluon gluelumps, which correspond to the perimeter-law exponential ${\mathrm{e}}^{-\mathsf{M}L}$ in the adjoint Wilson loop, where L is the length of the contour C. Accordingly, in the limit of large number of colors N, the adjoint Wilson loop is assumed to be of the form [66,67,68,69]

$$\langle W\left(C\right)\rangle =\frac{{\mathrm{e}}^{-\sigma S}+\frac{1}{N^2}{\mathrm{e}}^{-\mathsf{M}L}}{1+\frac{1}{N^2}},$$

(30)

where the normalization condition $\langle W\left(0\right)\rangle =1$ has been imposed. In Equation (30), S denotes the area of the minimal surface ${\mathsf{\Sigma }}_C$ bounded by C, and σ is the string tension in the adjoint representation. When the number of colors N is not large, but is rather equal to 2, one can check the degree to which the mixing between the string and gluelump states leads to the change in the large-N value of $\frac{1}{N^2}$ for the relative coefficient between the area- and the perimeter-law terms in Equation (30). We have performed such a check at zero temperature, by calculating numerically the static potential in the adjoint representation and observing its flattening at a certain critical separation between the static sources, which corresponds to the breaking of the adjoint string [70]. To this end, we have adopted the above-discussed Casimir scaling to express the zero-temperature adjoint string tension ${\sigma }_{\mathrm{a}}$ at $N=2$ via the known fundamental string tension at $N=3$ as ${\sigma }_{\mathrm{a}}=\frac{3}{2}{\sigma }_{\mathrm{f}}{|}_{N=3}$, where ${\sigma }_{\mathrm{f}}{|}_{N=3}\simeq {\left(0.44\mathrm{GeV}\right)}^2$. We have also used for the gluelump mass the inverse to the lattice value of 0.16 fm for the vacuum correlation length in the SU(2) Yang-Mills theory [71]. With these values at hand, we have checked the most sensitive case, where the extension of the flat contour C in the 4-th direction is larger than R only by a factor of 2. As a result, we have observed only few-percent relative corrections to the static potential when varying the said coefficient of $\frac{1}{N^2}=\frac{1}{4}$ by a large factor of 100. For this reason, we use for this coefficient the same value of $\frac{1}{4}$ throughout the subsequent calculations.

Note also that, at $T=0$, not only a one- but also a two-gluon gluelump can be formed, so that the masses of the one- and the two-gluon gluelumps define the correlation lengths of, respectively, non-confining and confining stochastic Yang-Mills fields [72,73,74,75,76,77]. Rather, the vacuum of the high-T dimensionally-reduced 3D Yang-Mills theory at issue is characterized by only one correlation length, which describes the large-distance exponential fall-off of the correlation function $\langle H_i^a\left(\mathbf{x}\right)H_j^b\left(\mathbf{0}\right)\rangle$ of the chromo-magnetic fields. Hence, we set M equal to the inverse of this correlation length, which yields [78,79] $M=\frac{\lambda T}{2\pi }$. Here, $\lambda ={\mathsf{g}}_T^{2}N$ is the so-called ’t Hooft coupling, which stays finite in the large-N limit, and ${\mathsf{g}}_T$ is the weak coupling at arbitrary N.

Similarly to the previous Section, the general strategy of our calculation is based on a reduction of the effective action of a given particle (i.e., a gauge boson, a ghost, or a Higgs boson) propagating in the stochastic Yang-Mills fields, to an equivalent effective action corresponding to some auxiliary Abelian field with a Gaussian action. The possibility for such a reduction is provided by the fact that the whole dependence on the non-Abelian gauge field in the effective action is encoded in the form of the area- and the perimeter-law terms entering Equation (30). The auxiliary Abelian field appears then in the course of the regularization of the area S of the minimal surface ${\mathsf{\Sigma }}_C$ and the perimeter L of the contour C. Furthermore, while for L such a regularization is straightforward (owing to the one-dimensionality of the contour), for S one can rather adopt various parametrizations in terms of C, which all render the path integral in the effective action calculable [7,8,9,39,80]. From all such parametrizations, we find the one used in Ref. [39] mostly suitable for the present analysis, as it explicitly accounts for the finiteness of the vacuum correlation length $1/M$.

We organize this Section as follows. In the next Subsection, we introduce regularized expressions for the area- and the perimeter-law terms that enter Equation (30), along with the corresponding auxiliary Abelian fields. We further proceed to the calculation of the para- and the diamagnetic contributions to the effective action, which are produced by the gauge boson and the ghost. In SubSection 2.2, we calculate additionally the contribution produced to the effective action by the Higgs boson. In the same SubSection 2.2, we use the so-obtained full expression for the effective action in order to find the critical behavior of the vacuum-energy density. In this way, we find this behavior consistent with the assumption about the value of 160 GeV for the critical temperature of the EW phase transition. In Appendix A, we present some technical details of the calculations whose results are used in the main text.

2.1. The Dia- and the Paramagnetic Contributions to the Effective Action

As outlined above, we start our analysis with regularizing the area- and the perimeter laws in Equation (30), so as to be further able to calculate the path integral which enters the effective action and involves the Wilson loop. The regularized expressions have the form

$$\sigma S\simeq c{M}^4{\int }_{{\mathsf{\Sigma }}_C}d{\sigma }_{ij}\left(\mathbf{x}\right){\int }_{{\mathsf{\Sigma }}_C}d{\sigma }_{ij}\left({\mathbf{x}}^{\prime }\right){\mathrm{e}}^{-M|\mathbf{x}-{\mathbf{x}}^{\prime }|}\mathrm{and}ML\simeq \frac{M^2}{2}{\oint }_{C}d{x}_i{\oint }_{C}d{x}_i^{\prime }{\mathrm{e}}^{-M|\mathbf{x}-{\mathbf{x}}^{\prime }|}.$$

(31)

While the proof of the second relation is straightforward, the parameter c from the first relation can be found by using the corresponding expression for the string tension [15] $\sigma =2c{M}^2\int d^{2}z{\mathrm{e}}^{-\left|\mathbf{z}\right|}=4\pi c{M}^2$. To this end, we use the known exact formula for the fundamental string tension [78,79], ${\sigma }_{\mathrm{f}}=\frac{N^2-1}{8\pi N^2}{\left(\lambda T\right)}^2$, and obtain the adjoint string tension via the Casimir scaling: $\sigma =\frac{C_{\mathrm{a}}}{C_{\mathrm{f}}}{\sigma }_{\mathrm{f}}=\frac{{\left(\lambda T\right)}^2}{4\pi }$, where the quadratic Casimir operators of the adjoint and the fundamental representations of SU(N) have the form $C_{\mathrm{a}}=N$ and $C_{\mathrm{f}}=\frac{N^2-1}{2N}$, respectively. (In particular, $\sigma \to 2{\sigma }_{\mathrm{f}}$ for $N\gg 1$, in agreement with the factorization of the adjoint Wilson loop in the large-N limit (cf. Ref. [68]).) Comparison of the two expressions for σ yields then the coefficient c as $c={\left(\frac{\lambda T}{4\pi M}\right)}^2=\frac{1}{4}$.

With this value of c at hand, we follow further the method of the previous Section to represent the so-regularized area- and perimeter laws in terms of the functional integrals over the auxiliary antisymmetric-tensor and vector fields, $B_{ij}$ and $h_i$, as

$$exp\left[-\frac{M^4}{4}{\int }_{{\mathsf{\Sigma }}_C}d{\sigma }_{ij}\left(\mathbf{x}\right){\int }_{{\mathsf{\Sigma }}_C}d{\sigma }_{ij}\left({\mathbf{x}}^{\prime }\right){\mathrm{e}}^{-M|\mathbf{x}-{\mathbf{x}}^{\prime }|}\right]=$$

$$=\int \left[\prod _{i<j}\mathcal{D}{B}_{ij}{\mathrm{e}}^{-\frac{1}{8\pi M^5}{\int }_{\mathbf{x}}{B}_{ij}{(-{\partial }^2+M^2)}^2{B}_{ij}}\right]{\mathrm{e}}^{\frac{i}{2}{\int }_{\mathbf{x}}{B}_{ij}{\mathsf{\Sigma }}_{ij}}\equiv {\langle {\mathrm{e}}^{\frac{i}{2}{\int }_{\mathbf{x}}{B}_{ij}{\mathsf{\Sigma }}_{ij}}\rangle }_B$$

and

$$exp\left[-\frac{M^2}{2}{\oint }_{C}d{x}_i{\oint }_{C}d{x}_i^{\prime }{\mathrm{e}}^{-M|\mathbf{x}-{\mathbf{x}}^{\prime }|}\right]=\int \mathcal{D}{h}_i{\mathrm{e}}^{-\frac{1}{16\pi M^3}{\int }_{\mathbf{x}}{h}_i{(-{\partial }^2+M^2)}^2{h}_i+i{\int }_{\mathbf{x}}{h}_i{j}_i}\equiv {\langle {\mathrm{e}}^{i{\int }_{\mathbf{x}}{h}_i{j}_i}\rangle }_h.$$

(32)

Here $j_i\equiv j_i(\mathbf{x};C)={\oint }_{C}d{z}_i\delta (\mathbf{x}-\mathbf{z})$ is the conserved current associated with the contour C, and ${\int }_{\mathbf{x}}\equiv \int d^{3}x$. Furthermore, similarly to the previous Section, we choose the surface element $d{\sigma }_{ij}$ in the form of an oriented, infinitely thin triangle built up of the position vector $z_i\left(\tau \right)$ and the differential element $d{z}_i={\dot{z}}_{i}d\tau$, namely $d{\sigma }_{ij}\left(\mathbf{z}\right)=\frac{1}{2}(z_i{\dot{z}}_j-z_j{\dot{z}}_i)d\tau$. The surface tensor ${\mathsf{\Sigma }}_{ij}\equiv {\mathsf{\Sigma }}_{ij}(\mathbf{x};C)$, defined as ${\mathsf{\Sigma }}_{ij}(\mathbf{x};C)={\int }_{{\mathsf{\Sigma }}_C}d{\sigma }_{ij}\left(\mathbf{z}\right)\delta (\mathbf{x}-\mathbf{z})$, takes then the form ${\mathsf{\Sigma }}_{ij}=\frac{1}{2}{\int }_0^{s}d\tau (z_i{\dot{z}}_j-z_j{\dot{z}}_i)\delta (\mathbf{x}-\mathbf{z}\left(\tau \right))$. With this expression for ${\mathsf{\Sigma }}_{ij}$, the exponential ${\mathrm{e}}^{\frac{i}{2}{\int }_{\mathbf{x}}{B}_{ij}{\mathsf{\Sigma }}_{ij}}$ can be written as ${\mathrm{e}}^{\frac{i}{2}{\int }_{\mathbf{x}}{B}_{ij}{\mathsf{\Sigma }}_{ij}}={\mathrm{e}}^{i{\int }_{\mathbf{x}}{\mathcal{A}}_i{j}_i}$. Here, we have introduced, instead of $B_{ij}$, an auxiliary vector field ${\mathcal{A}}_i\left(\mathbf{x}\right)=\frac{1}{2}{x}_j{B}_{ji}\left(\mathbf{x}\right)$, whose strength tensor ${\mathcal{F}}_{ij}={\partial }_i{\mathcal{A}}_j-{\partial }_j{\mathcal{A}}_i$ has the form ${\mathcal{F}}_{ij}=B_{ij}+C_{ij}$, where $C_{ij}\left(\mathbf{x}\right)=\frac{1}{2}{x}_k({\partial }_i{B}_{kj}-{\partial }_j{B}_{ki})$. Thus, the Wilson loop (30) can altogether be written as

$$\langle W\left(C\right)\rangle =\frac{{\langle {\mathrm{e}}^{i{\int }_{\mathbf{x}}{\mathcal{A}}_i{j}_i}\rangle }_B+\frac{1}{N^2}{\langle {\mathrm{e}}^{i{\int }_{\mathbf{x}}{h}_i{j}_i}\rangle }_h}{1+\frac{1}{N^2}}.$$

(33)

Let us now consider the one-loop effective action for a spinless adjoint particle of mass m (e.g., a ghost). It has the form

$$\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle ={\int }_0^{\infty }\frac{ds}{s}{\mathrm{e}}^{-m^{2}s}{\int }_P^{}\mathcal{D}\mathbf{z}{\mathrm{e}}^{-\frac{1}{4}{\int }_0^{s}d\tau {\dot{\mathbf{z}}}^2}\langle W\left(C\right)\rangle ,$$

(34)

where the contour C is parameterized by the vector-function $\mathbf{z}\left(\tau \right)$, and P stands for the periodic boundary conditions, i.e., ${\int }_P\equiv {\int }_{\mathbf{z}\left(s\right)=\mathbf{z}\left(0\right)}$. Rather, in the world-line representation of a spinning particle (e.g., a gauge boson), an additional term $\propto F_{ij}^a{T}^a$, where $T^a$ is an SU(N)-generator in the adjoint representation, appears. As was discussed, this term can be recovered by acting on the Wilson loop with the area-derivative operator $\frac{\delta }{\delta {\sigma }_{ij}\left(\mathbf{z}\right)}$, which allows one, also in the spinning case, to reduce the gauge-field dependence of the effective action to that of the Wilson loop only. In what follows, we will again use for the effective action the known closed-form expression (That is, this expression is valid to all orders in s, being therefore suitable for the studies of the infra-red physics. It can be obtained by using either the standard covariant perturbation theory, or the world-line formalism [22].), which corresponds to two $F_{ij}^a$’s standing in the pre-exponent [22,40]:

$$\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle =-\frac{4}{{\left(4\pi \right)}^{3/2}}{\int }_0^{\infty }\frac{ds}{\sqrt{s}}{\mathrm{e}}^{-m^{2}s}{\int }_{\mathbf{x}}{\mathrm{tr}}_c\langle (F_{ij}^a{T}^a)\left[f-\frac{1}{4}·\frac{f-1}{\xi }\right](F_{ij}^a{T}^a)\rangle ,$$

(35)

where $f\equiv f\left(\xi \right)={\int }_0^{1}du{\mathrm{e}}^{u(1-u)\xi }$ and $\xi =s{D}_i^2$. This effective action represents the sum of contributions produced by the ghost and the gauge boson. Hence, for the above-given reason, one can replace each of the two $(F_{ij}^a{T}^a)$’s in Equation (35) either by ${\mathcal{F}}_{ij}$ or by $H_{ij}\equiv {\partial }_i{h}_j-{\partial }_j{h}_i$, depending on whether we consider the area- or the perimeter-law term in $\langle W\left(C\right)\rangle$. Simultaneously, once the parametrization (33) for $\langle W\left(C\right)\rangle$ is used in Equation (35), one should replace $A_i^a{T}^a$ in $D_i$ either by ${\mathcal{A}}_i$ or by $h_i$, and remove “${\mathrm{tr}}_c$”. Starting then with the area-law term, we consider the following path-integral representation for the corresponding contribution to the effective action, which is similar to the one of the previous Section:

$${\int }_{\mathbf{x}}{\langle {\mathcal{F}}_{ij}\left(\mathbf{x}\right)\left[f-\frac{1}{4}·\frac{f-1}{\xi }\right]{\mathcal{F}}_{ij}\left(\mathbf{x}\right)\rangle }_B=$$

$$={\int }_{\mathbf{x}}{\langle {\mathcal{F}}_{ij}\left(\mathbf{x}\right){\int }_0^{1}du\left[{\mathrm{e}}^{u(1-u)\xi }-\frac{1}{4}u(1-u){\int }_0^{1}d\alpha {\mathrm{e}}^{\alpha u(1-u)\xi }\right]{\mathcal{F}}_{ij}\left(\mathbf{x}\right)\rangle }_B=$$

$$=V{\int }_0^{1}du{\int }_{\mathbf{y}}\left[\underset{\mathbf{r}\left(0\right)=\mathbf{0}}{\overset{\mathbf{r}\left(u\right(1-u\left)s\right)=\mathbf{y}}{\int }}\mathcal{D}\mathbf{r}{\mathrm{e}}^{-\frac{1}{4}\underset{0}{\overset{u(1-u)s}{\int }}d\tau {\dot{\mathbf{r}}}^2}-\frac{1}{4}u(1-u){\int }_0^{1}d\alpha \underset{\mathbf{r}\left(0\right)=\mathbf{0}}{\overset{\mathbf{r}\left(\alpha u\right(1-u\left)s\right)=\mathbf{y}}{\int }}\mathcal{D}\mathbf{r}{\mathrm{e}}^{-\frac{1}{4}\underset{0}{\overset{\alpha u(1-u)s}{\int }}d\tau {\dot{\mathbf{r}}}^2}\right]\times$$

$$\times {\langle {\mathcal{F}}_{ij}\left(\mathbf{0}\right)exp\left(i{\int }_{\mathbf{0}}^{\mathbf{y}}d{r}_k{\mathcal{A}}_k\right){\mathcal{F}}_{ij}\left(\mathbf{y}\right)\rangle }_B,$$

(36)

where V is the volume occupied by the system. Approximating the phase factor $exp\left(i{\int }_{\mathbf{0}}^{\mathbf{y}}d{r}_k{\mathcal{A}}_k\right)$ in Equation (36) by unity, so as to respect the initial two-point approximation, we reduce the path integrals to the Green function of the heat equation. Noticing further that ${\langle {\mathcal{F}}_{ij}\left(\mathbf{0}\right){\mathcal{F}}_{ij}\left(\mathbf{y}\right)\rangle }_B={\langle B_{ij}\left(\mathbf{0}\right)B_{ij}\left(\mathbf{y}\right)\rangle }_B+{\langle B_{ij}\left(\mathbf{0}\right)C_{ij}\left(\mathbf{y}\right)\rangle }_B$, we use the following explicit expressions for the latter correlation functions:

$${\langle B_{ij}\left(\mathbf{0}\right)B_{ij}\left(\mathbf{y}\right)\rangle }_B=3M^4{\mathrm{e}}^{-M\left|\mathbf{y}\right|},{\langle B_{ij}\left(\mathbf{0}\right)C_{ij}\left(\mathbf{y}\right)\rangle }_B=-M^5\left|\mathbf{y}\right|{\mathrm{e}}^{-M\left|\mathbf{y}\right|}.$$

(37)

In the same way, we treat the contribution produced to $\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle$ by the perimeter law, which corresponds to the average ${\langle {\mathrm{e}}^{i{\int }_{\mathbf{x}}{h}_i{j}_i}\rangle }_h$ in Equation (33). The corresponding correlation function ${\langle H_{ij}\left(\mathbf{0}\right)H_{ij}\left(\mathbf{y}\right)\rangle }_h$ can be calculated through the average ${\langle h_i\left(\mathbf{0}\right)h_j\left(\mathbf{y}\right)\rangle }_h=M^2{\delta }_{ij}{\mathrm{e}}^{-M\left|\mathbf{y}\right|}$, and reads

$${\langle H_{ij}\left(\mathbf{0}\right)H_{ij}\left(\mathbf{y}\right)\rangle }_h=4M^3\left(\frac{2}{\left|\mathbf{y}\right|}-M\right){\mathrm{e}}^{-M\left|\mathbf{y}\right|}.$$

Furthermore, consistency with the left-hand side of Equation (32) requires to keep here the full exponential ${\mathrm{e}}^{-M\left|\mathbf{y}\right|}$, rather than the first terms of its Taylor series, which means that we effectively consider the infra-red limit of $M\left|\mathbf{y}\right|>1$. For this reason, the ultra-violet contribution $\frac{2}{\left|\mathbf{y}\right|}$ can be disregarded in comparison with M, i.e., we can use the approximate expression

$${\langle H_{ij}\left(\mathbf{0}\right)H_{ij}\left(\mathbf{y}\right)\rangle }_h\simeq -4M^4{\mathrm{e}}^{-M\left|\mathbf{y}\right|},$$

(38)

which looks similar to Equation (37). By representing then the effective action (35) as a sum of the dia- and the paramagnetic contributions (cf. Ref. [65]), $\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle ={\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{dia}}+{\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{para}}$, with

$${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{dia}}=\frac{1}{{\left(4\pi \right)}^{3/2}}{\int }_0^{\infty }\frac{ds}{\sqrt{s}}{\mathrm{e}}^{-m^{2}s}{\int }_{\mathbf{x}}{\mathrm{tr}}_c\langle (F_{ij}^a{T}^a)\frac{f-1}{\xi }(F_{ij}^a{T}^a)\rangle$$

and

$${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{para}}=-\frac{4}{{\left(4\pi \right)}^{3/2}}{\int }_0^{\infty }\frac{ds}{\sqrt{s}}{\mathrm{e}}^{-m^{2}s}{\int }_{\mathbf{x}}{\mathrm{tr}}_c\langle (F_{ij}^a{T}^a)f(F_{ij}^a{T}^a)\rangle ,$$

we rewrite these contributions in terms of the B- and h-averages. For the diamagnetic contribution, we have

$${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{dia}}=\frac{V{M}^3}{1+\frac{1}{N^2}}{\int }_0^{1}du\left(I_1+\frac{1}{N^2}{I}_2\right),$$

(39)

where

$$I_1=\frac{1}{{\left(4\pi \right)}^{3/2}}{\int }_0^{\infty }\frac{ds}{\sqrt{s}}{\mathrm{e}}^{-m^{2}s}{\int }_{\mathbf{x}}{\langle {\mathcal{F}}_{ij}\left(\mathbf{x}\right)\frac{f-1}{\xi }{\mathcal{F}}_{ij}\left(\mathbf{x}\right)\rangle }_B$$

and

$$I_2=\frac{1}{{\left(4\pi \right)}^{3/2}}{\int }_0^{\infty }\frac{ds}{\sqrt{s}}{\mathrm{e}}^{-m^{2}s}{\int }_{\mathbf{x}}{\langle H_{ij}\left(\mathbf{x}\right)\frac{f-1}{\xi }{H}_{ij}\left(\mathbf{x}\right)\rangle }_h.$$

Similarly, for the paramagnetic contribution, we have

$${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{para}}=\frac{V{M}^3}{1+\frac{1}{N^2}}{\int }_0^{1}du\left(J_1+\frac{1}{N^2}{J}_2\right),$$

where

$$J_1=-\frac{4}{{\left(4\pi \right)}^{3/2}}{\int }_0^{\infty }\frac{ds}{\sqrt{s}}{\mathrm{e}}^{-m^{2}s}{\int }_{\mathbf{x}}{\langle {\mathcal{F}}_{ij}\left(\mathbf{x}\right)f{\mathcal{F}}_{ij}\left(\mathbf{x}\right)\rangle }_B$$

and

$$J_2=-\frac{4}{{\left(4\pi \right)}^{3/2}}{\int }_0^{\infty }\frac{ds}{\sqrt{s}}{\mathrm{e}}^{-m^{2}s}{\int }_{\mathbf{x}}{\langle H_{ij}\left(\mathbf{x}\right)f{H}_{ij}\left(\mathbf{x}\right)\rangle }_h.$$

It turns out that the integrals $I_1$, $I_2$, $J_1$, and $J_2$ can be calculated analytically. The details of this calculation are presented in Appendix A, where N is already set equal to 2. Performing the remaining u-integrations in Equations (A1) and (A2) numerically, one can see that ${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{dia}}$ is a positive, monotonically increasing function of the parameter $\gamma =\frac{M}{m}$, while ${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{para}}$ is a negative and a monotonically decreasing function of γ. Both functions vanish at $\gamma \to 0$, i.e., in the limit of $m\to \infty$, as it should be. Remarkably, their ratio in this limit stays finite, being equal to $-24$. With the increase of γ, i.e., with the decrease of m, the absolute value of the ratio decreases, approaching at $\gamma \to \infty$ the expected value of 16 (cf. Ref. [65]). For illustration, we plot in Figure 1 the absolute value of the ratio ${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{para}}/{\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{dia}}$ in the range of $\gamma \in (0,10)$.

2.2. Higgs Contribution to the Effective Action. The Vacuum Instability

As is known [65], the negative sign of ${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{para}}$ can lead to the vacuum instability associated with the vanishing vacuum-energy density. (As can be shown (cf. Ref. [65]), the fermionic contribution to $\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle$ does not affect this instability. Namely, the instability is related to the fact that the gauge bosons have spin 1, while it does not occur for spin-$\frac{1}{2}$ particles. For this reason, the fermionic contribution to the effective action can be disregarded altogether.) The temperature at which this instability occurs, can be identified with the critical temperature of the EW phase transition. The condition of the vacuum instability has the form [65]

$$\frac{\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }{V}+\frac{G}{4}=0,$$

(40)

where G is the gluon condensate in the theory (29). By using the fact that G is equal to the condensate in the high-temperature dimensionally-reduced Yang-Mills theory, we can express it in terms of M. To this end, we recall the following regularization of $\sigma S$, which is suggested by the stochastic vacuum model at $T=0$ and arbitrary N [see Equation (15)]:

$$\sigma S\simeq \frac{\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle }{96N}{\int }_{{\mathsf{\Sigma }}_C}d{\sigma }_{\mu \nu }\left(x\right){\int }_{{\mathsf{\Sigma }}_C}d{\sigma }_{\mu \nu }\left(x^{\prime }\right){\mathrm{e}}^{-M_0|x-x^{\prime }|}.$$

(41)

Here $\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle$ and $M_0$ are respectively the zero-temperature gluon condensate and the zero-temperature inverse vacuum correlation length, where the values of indices μ and ν run from 1 to 4. In the deconfinement phase (but at temperatures smaller than the temperature of dimensional reduction), where only the chromo-magnetic part of the gluon condensate is non-vanishing, the counterpart of this expression for a spatial surface ${\mathsf{\Sigma }}_C$ has the form

$$\sigma S\simeq \frac{\langle {(g{F}_{ij}^a)}^2\rangle }{96N}{\int }_{{\mathsf{\Sigma }}_C}d{\sigma }_{ij}\left(x\right){\int }_{{\mathsf{\Sigma }}_C}d{\sigma }_{ij}\left(x^{\prime }\right){\mathrm{e}}^{-M|x-x^{\prime }|}.$$

Comparing this equation at $N=2$ with Equation (31), we obtain for the chromo-magnetic condensate: $\langle {(g{F}_{ij}^a)}^2\rangle =96·2·c{M}^4=48M^4$. Accordingly, in the dimensionally-reduced theory at issue, the condensate has the form $G=48M^4/T$. Noticing now that $M=g_T^{2}T/\pi$ for $N=2$, and using Equations (A1) and (A2), we can write condition (40) explicitly as

$$1+\frac{4\sqrt{\pi }}{15g_T^2\gamma }{\int }_0^{1}du\left\{\frac{arccosa}{{(1-a^2)}^{3/2}}\left[a^2-2{\gamma }^2·\frac{2+a^2}{1-a^2}\right]+\frac{a}{1-a^2}\left[2{\gamma }^2·\frac{4-5a^2+a^4}{{(1-a^2)}^2}-a^2\right]\right\}=0,$$

(42)

where $a=\gamma \sqrt{u(1-u)}$.

Let us further account for the Higgs contribution to the effective action. Given that the Higgs boson is as scalar particle as the ghost, we can start with the following estimate of this contribution: ${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_H=-\frac{2}{3}{\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{gh}}{|}_{m\to m_H}$. Here the factor $-\frac{2}{3}$ accounts for the difference in the statistics of the Higgs boson and the ghost, as well as for the fact that the ghost transforms under the adjoint representation, while the Higgs boson transforms under the fundamental representation. Noticing further that ${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{gh}}=-2{\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{dia}}$, we can write the above estimate as ${\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_H=\frac{4}{3}{\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle }_{\mathrm{dia}}{|}_{m\to m_H}$. However, this estimate needs to be corrected further. Indeed, since the Higgs boson transforms under the fundamental representation, one should use here Equation (39) with the term $\frac{1}{N^2}{I}_2$ (as well as the term $\frac{1}{N^2}$ in the denominator) subtracted. For the same reason, σ in Equation (31) should be replaced by ${\sigma }_{\mathrm{f}}$, which changes the value of the coefficient c from $\frac{1}{4}$ to $\frac{3}{32}$. As one can see, this leads to the following modification of Equation (37):

$${\langle B_{ij}\left(\mathbf{0}\right)B_{ij}\left(\mathbf{y}\right)\rangle }_B=\frac{9}{8}{M}^4{\mathrm{e}}^{-M\left|\mathbf{y}\right|},{\langle B_{ij}\left(\mathbf{0}\right)C_{ij}\left(\mathbf{y}\right)\rangle }_B=-\frac{3}{8}{M}^5\left|\mathbf{y}\right|{\mathrm{e}}^{-M\left|\mathbf{y}\right|},$$

which results in the overall factor of $\frac{3}{8}$. Once brought together, these corrections lead to the following additional Higgs contribution, which is to be included into the curly brackets in Equation (42) (As discussed in Ref. [65], there also exists an additional contribution to $\langle \mathsf{\Gamma }\left[A_i^a\right]\rangle$, which stems from the higher-order operator $i{g}_T\left({\partial }_{\mu }{A}_{\mu }^a\right)(\phi †T^a\mathsf{\Phi }-{\mathsf{\Phi }}^{†}{T}^a\phi )$, where φ is the quantum fluctuation of the Higgs field. The effect produced by this interaction is that the vacuum instability may take place at somewhat larger values of the vacuum correlation length $1/M$, i.e., at smaller temperatures. Nevertheless, as has been shown in Ref. [65], this contribution is numerically negligible):

$$\frac{5}{8}\frac{\gamma }{{\gamma }_H}\left[\frac{arccos{a}_H}{{(1-a_H^2)}^{3/2}}+2a_H-\frac{a_H}{1-a_H^2}-\frac{\pi }{2}\right],$$

(43)

where ${\gamma }_H=\frac{M}{m_H}$ and $a_H={\gamma }_H\sqrt{u(1-u)}$.

We evaluate now numerically the so-obtained full left-hand side of Equation (42) as a function of T. Owing to the extremely slow logarithmic evolution of ${\alpha }_{\mathrm{w}}=\frac{g_T^2}{4\pi }$, one usually approximates it by a constant [58,59,60,61] ${\alpha }_{\mathrm{w}}\simeq \frac{1}{30}$, which corresponds to the value $g_T\simeq 0.65$ of the weak coupling. We need further a certain realistic parametrization for the T-dependence of ${\phi }_T$ in the formula $m=\frac{g_3{\phi }_T}{2}$. Noticing that $1/m$ should vanish as the inverse correlation length, we choose this dependence in the form ${\phi }_T=\phi {(1-T/T_c)}^{\nu }$. Assuming the crossover behavior to occur as a limiting case of the first-order phase transition, we set $\nu =1/3$, which corresponds to the so-called weak first-order phase transition [81,82]. Furthermore, by using the relation [62,63] $\frac{v^2}{T}=2\langle {\phi }^{†}\phi \rangle$, where $v\simeq 246\mathrm{GeV}$ is the Higgs v.e.v. in the zero-temperature 4D EW theory, we estimate φ as $\phi =\frac{v}{\sqrt{2T}}$. Next, owing to the proportionality of the Higgs mass to the Higgs v.e.v., we can parametrize the critical behavior of the Higgs mass as $m_H=m_H{|}_{T=0}{(1-T/T_c)}^{\nu }$, where $m_H{|}_{T=0}=125\mathrm{GeV}$. Assuming then the value of $T_c=160\mathrm{GeV}$ (cf. Refs. [58,59,60,61,62,63]), we calculate numerically the left-hand side of Equation (42). In Figure 2, we plot this quantity both with and without the Higgs contribution (43). Quite remarkably, we find that it vanishes at the same temperature $T=T_c$, which can be viewed as an indication of self-consistency of the whole calculation.

In conclusion of this Section, our study of the EW phase transition has been using the fact that at temperatures close to the critical temperature of the phase transition, the dynamics of the Standard Model is described by the effective 3D theory (29), where spatial confinement plays an important role. In order to account for this phenomenon, one needs to explore the phase transition by means of the adequate non-perturbative methods. In this way, it has been suggested in Ref. [65] to use for the analysis the one-loop effective action, where the non-perturbative dynamics is encoded in the form of a Wilson loop. We have performed an analytic calculation of the corresponding effective action, which contains the contributions of a gauge boson, a ghost, and a Higgs boson. While the Higgs boson transforms under the fundamental representation, so that the corresponding Wilson loop contains only the area-law term, this is not the case for the gauge boson and the ghost. The latter transform under the adjoint representation, which leads to the appearance in the corresponding Wilson loop of a separate perimeter-law term. By using some known results of lattice phenomenology [71,72,73,74,75,76,77], we have regularized the area- and the perimeter-law terms, and used the auxiliary-field formalism in order to make the corresponding path integral, which enters the effective action, analytically calculable. More specifically, we have calculated that path integral by reducing it to the one of a particle in the Abelian background field, and using for the latter one the known closed-form expression [22,40] which contains two field-strength tensors in the pre-exponent. By separating the thus obtained contributions to the effective action into a “paramagnetic” part, emerging from the interactions of the background gauge fields with the color-spin degrees of freedom of the gauge boson, and the remaining “diamagnetic” part, we have first reproduced the earlier result of Ref. [65]. This result states that, in the limit where the ghost- and the gauge-boson mass m is much smaller than the inverse vacuum correlation length M, the absolute value of the paramagnetic contribution to the effective action exceeds the diamagnetic contribution by a factor of 16. Moreover, within the present formalism, we have obtained the ratio of the para- and the diamagnetic contributions for arbitrary values of the parameter $M/m$ (cf. Figure 1). Furthermore, the dominance of the paramagnetic contribution over the diamagnetic one, which persists also after accounting for the Higgs-boson part of the diamagnetic contribution, leads to the vacuum instability associated with the vanishing full vacuum-energy density. This vacuum instability is self-consistently shown to occur at $T=160\mathrm{GeV}$, which is the commonly accepted present-day value for the critical temperature of the EW phase transition. This observation indicates that the phase transition in the EW theory can be associated with the discussed vacuum instability.

3. A Semi-Classical Analogue of the Relation between the QCD Condensates

Nowadays, it is quite evident that the two main non-perturbative phenomena in QCD, confinement and chiral symmetry breaking, are not independent of each other. A clear indication in favor of this statement stems from the following relation that holds between the order parameters of these two phenomena, which are the chiral and the chromo-electric gluon condensates, $\langle \overline{\psi }\psi \rangle$ and $\langle {(g{E}_i^a)}^2\rangle$ [83,84,85]:

$$\langle \overline{\psi }\psi \rangle \propto -\lambda \langle {\left(g{E}_i^a\right)}^2\rangle .$$

(44)

In this expression, λ is the so-called chromo-electric vacuum correlation length, while the concrete value of the order-$\mathcal{O}\left(1\right)$ non-universal proportionality coefficient is unimportant. Clearly, this relation is specific for a stochastic vacuum, such as that of QCD, while it cannot hold in a vacuum characterized by some constant chromo-electric field $E_i^a$. That is, such a “classical” vacuum cannot support either of the two phenomena. Nevertheless, it is legitimate to pose a question of whether other classical fields can lead to the condensation of quantum fields so as to yield a formula for the corresponding condensates similar to Equation (44). In this Section, we follow Ref. [86] to show that a simple example of such a classical field is provided by the gravitational field of a spherically symmetric object of a constant energy density. For our illustrative purposes, it suffices to calculate the condensate of a scalar field, which is minimally coupled to this gravitational field. In the vanishing-mass limit of the scalar field, the role of the vacuum correlation length is played by the size of the object. We will obtain a formula for the scalar-field condensate through the second order in the curvature $\mathcal{R}$ and the Ricci tensor ${\mathcal{R}}_{\mu \nu }$. Within the present analogy, these two quantities play the role of classical counterparts of $E_i^a$. In order to derive such a formula, we use the known closed-form expression [40] for the one-loop effective action of a scalar field in the gravitational background.

Hence, let us consider a real-valued massive scalar field $\varphi \left(x\right)$ interacting with the gravitational field $g_{\mu \nu }\left(x\right)$. The corresponding Euclidean action has the form

$$S=\frac{1}{2}\int d^{4}x\sqrt{g}\varphi (-\square +m^2)\varphi ,\mathrm{where}\square \varphi =\frac{1}{\sqrt{g}}{\partial }_{\mu }\left(\sqrt{g}{g}^{\mu \nu }{\partial }_{\nu }\varphi \right),\mathrm{and}g\equiv det{g}_{\mu \nu }.$$

Integrating over the field ϕ, one obtains the following effective action:

$$\mathsf{\Gamma }\left[g_{\mu \nu }\right]\equiv -ln\int \mathcal{D}\varphi {\mathrm{e}}^{-S}=\frac{1}{2}lndet(-\square +m^2).$$

In Ref. [40], a closed-form expression for $\mathsf{\Gamma }\left[g_{\mu \nu }\right]$ has been obtained through the second order in curvature:

$$\mathsf{\Gamma }\left[g_{\mu \nu }\right]=\frac{1}{2}{\int }_0^{\infty }\frac{ds}{s}\int \frac{d^{4}x\sqrt{g}}{{\left(4\pi s\right)}^2}\times$$

$$\times \left\{1+s\left(\frac{\mathcal{R}}{6}-m^2\right)+\frac{{\left(ms\right)}^2}{2}\left(m^2-\frac{\mathcal{R}}{6}\right)+s^2\left[m^2{f}_2\mathcal{R}+\mathcal{R}{f}_3\mathcal{R}+{\mathcal{R}}_{\mu \nu }{f}_4{\mathcal{R}}^{\mu \nu }\right]\right\}.$$

(45)

In this formula,

$$f_2\equiv -\frac{f}{6}+\frac{f-1}{2\xi },f_3\equiv \frac{f}{32}-\frac{f-1}{8\xi }-\frac{f_4}{8},f_4\equiv \frac{f-1-\frac{1}{6}\xi }{{\xi }^2},$$

(46)

where $f\equiv {\int }_0^{1}du{\mathrm{e}}^{u(1-u)\xi }$, and $\xi \equiv s\square$. In what follows, we will be interested in the small-m limit, where the s-series in Equation (45) can be recovered from the full factor ${\mathrm{e}}^{-m^{2}s}$. To the order $\mathcal{O}\left(s^2\right)$, this approximation yields

$$\mathsf{\Gamma }\left[g_{\mu \nu }\right]\simeq \frac{1}{2}{\int }_0^{\infty }\frac{ds}{s}{\mathrm{e}}^{-m^{2}s}\int \frac{d^{4}x\sqrt{g}}{{\left(4\pi s\right)}^2}\left\{1+\frac{s\mathcal{R}}{6}+s^2\left[\frac{m^2\mathcal{R}}{12}+m^2{f}_2\mathcal{R}+\mathcal{R}{f}_3\mathcal{R}+{\mathcal{R}}_{\mu \nu }{f}_4{\mathcal{R}}^{\mu \nu }\right]\right\}.$$

(47)

As we will see below, the $\mathcal{O}\left(s^2\right)$-term provides the leading contribution to the scalar-field condensate in the physically interesting case where the size of a spherical object that produces the gravitational field is sufficiently small, being therefore analogous to the chromo-electric vacuum correlation length in QCD. The relation

$$\langle {\varphi }^2\rangle =-2\frac{\partial \mathsf{\Gamma }\left[g_{\mu \nu }\right]/\partial m^2}{\int d^{4}x\sqrt{g}}$$

(48)

allows us further to obtain, from the effective action (47), the scalar-field condensate $\langle {\varphi }^2\rangle$. Clearly, this technique is conceptually similar to the one discussed in Section 1, which allows one to relate the quark and the gluon condensates away from the heavy-quark limit. Nevertheless, as it has already been mentioned, QCD differs from the present model due to the fact that the non-perturbative Yang-Mills fields, being quantum in their nature, should be averaged over, whereas the gravitational field in the present case is classical. (Had QCD been equivalent to the present model, one could make a conjecture that the average over the Yang-Mills fields can play the role of a specific classical background, which, unlike the standard classical chromo-electric field, is capable to produce chiral condensate.) Furthermore, unlike Section 1, we work here in the small-m limit.

In order to calculate the effective action (47), we use the following integral representations of the formfactors $f_2$, $f_3$, and $f_4$:

$$f_2=-\frac{1}{6}{\int }_0^{1}du{\mathrm{e}}^{u(1-u)\xi }+\frac{1}{2}{\int }_0^{1}duu(1-u){\int }_0^{1}d\alpha {\mathrm{e}}^{\alpha u(1-u)\xi },$$

(49)

$$f_3=\frac{1}{32}\left[{\int }_0^{1}du{\mathrm{e}}^{u(1-u)\xi }-4{\int }_0^{1}duu(1-u){\int }_0^{1}d\alpha [1+u(1-u)(1-\alpha )]{\mathrm{e}}^{\alpha u(1-u)\xi }\right],$$

(50)

$$f_4={\int }_0^{1}du{\left[u(1-u)\right]}^2{\int }_0^{1}d\alpha (1-\alpha ){\mathrm{e}}^{\alpha u(1-u)\xi }.$$

(51)

These representations are similar to those which were used in Section 1 and Section 2. Some details of their derivation from Equation (46) can be found in Ref. [86]. Differentiating then the effective action (47) with respect to $m^2$, we have

$$-\frac{\partial \mathsf{\Gamma }\left[g_{\mu \nu }\right]}{\partial m^2}\simeq$$

$$\simeq \frac{1}{2}{\int }_0^{\infty }ds{\mathrm{e}}^{-m^{2}s}\int \frac{d^{4}x\sqrt{g}}{{\left(4\pi s\right)}^2}\left\{1+\frac{s\mathcal{R}}{12}-s{f}_2\mathcal{R}+s^2\left[\frac{m^2\mathcal{R}}{12}+m^2{f}_2\mathcal{R}+\mathcal{R}{f}_3\mathcal{R}+{\mathcal{R}}_{\mu \nu }{f}_4{\mathcal{R}}^{\mu \nu }\right]\right\}.$$

(52)

Furthermore, as is known [87], the ultraviolet-divergent terms $1+\frac{s\mathcal{R}}{12}-s{f}_2\mathcal{R}$ in Equation (52) can be renormalized by adding to $\mathsf{\Gamma }\left[g_{\mu \nu }\right]$ the sum of the bare Einstein-Hilbert and the cosmological-constant actions. As a result, the gravitational constant gets renormalized independently of the cosmological constant. In the present context, the net effect of this renormalization procedure amounts to using the standard value of the gravitational constant, $G=6.7·{10}^{-39}{\mathrm{GeV}}^{-2}$, and subtracting from Equation (52) the three aforementioned ultraviolet-divergent terms. Moreover, since the terms $s^2\left[\frac{m^2\mathcal{R}}{12}+m^2{f}_2\mathcal{R}\right]$ in Equation (52) are suppressed in the small-m limit of interest, they will be henceforth disregarded. Thus, we arrive at the following intermediate expression:

$$-\frac{\partial \mathsf{\Gamma }\left[g_{\mu \nu }\right]}{\partial m^2}\simeq \frac{1}{2{\left(4\pi \right)}^2}{\int }_0^{\infty }ds{\mathrm{e}}^{-m^{2}s}\int d^{4}x\sqrt{g}\left(\mathcal{R}{f}_3\mathcal{R}+{\mathcal{R}}_{\mu \nu }{f}_4{\mathcal{R}}^{\mu \nu }\right).$$

(53)

As outlined above, we consider now the metric corresponding to the inner part of a spherically symmetric object of radius R, which is filled with the matter of a constant energy density ε. In this case, the metric itself, as well as the associated scalar curvature $\mathcal{R}$ and the Ricci tensor ${\mathcal{R}}_{\mu \nu }$, are the functions of the radial coordinate $\mathbf{r}$. In Appendix B, we summarize some known facts about this metric, and calculate the corresponding scalar curvature. As follows from this calculation, the adopted approximation, where the terms $\frac{m^2\mathcal{R}}{12}+m^2{f}_2\mathcal{R}$ in Equation (52) are disregarded in comparison with the terms $\mathcal{R}{f}_3\mathcal{R}+{\mathcal{R}}_{\mu \nu }{f}_4{\mathcal{R}}^{\mu \nu }$, is translated in the inequality

$$m\ll \sqrt{\epsilon G}.$$

(54)

Furthermore, since we are interested in the terms $\mathcal{O}\left({\mathcal{R}}^2\right)$, $\mathcal{O}\left({\mathcal{R}}_{\mu \nu }{\mathcal{R}}^{\mu \nu }\right)$, and not in higher-curvature terms, the formfactors $f_3$ and $f_4$ in Equation (53) can be taken at $\xi =0$. Indeed, in the case of $\xi \ne 0$, one has, for instance in the term $f_3\mathcal{R}$, the following heat-kernel integral:

$${\left.{\mathrm{e}}^{\tau \mathsf{\Delta }}\mathcal{R}\right|}_{\mathrm{x}}=\frac{1}{{\left(4\pi \tau \right)}^{3/2}}\int d^{3}y{\mathrm{e}}^{-\frac{{\mathrm{y}}^2}{4\tau }}\mathcal{R}(\mathbf{x}-\mathbf{y}),$$

where $\tau =u(1-u)$ or $\tau =\alpha u(1-u)$. Once expanded in $\mathbf{y}$, the scalar curvature $\mathcal{R}(\mathbf{x}-\mathbf{y})$ yields $\mathcal{R}\left(\mathbf{x}\right)+y^2\mathcal{O}\left({\mathcal{R}}^2\left(\mathbf{x}\right)\right)$, since the term linear in $\mathbf{y}$ vanishes upon the integration. This observation shows that we should restrict ourselves to the leading term of the heat-kernel expansion, which is indeed equivalent to setting $\xi =0$. Then the formfactors $f_3$ and $f_4$ become just numbers, namely

$$f_3{|}_{\xi =0}=\frac{1}{120},f_4{|}_{\xi =0}=\frac{1}{60}.$$

(55)

Moreover, since we consider the propagation of the field ϕ only to radial distances $r\le R$, the upper limit of the proper-time integration should be restricted by some value $s_{\max }=\gamma R^2$, where γ is a dimensionless constant of the order of unity. The s-integration yields then a factor of $\frac{1}{m^2}\left(1-{\mathrm{e}}^{-\gamma m^2{R}^2}\right)$, which ensures that $\langle {\varphi }^2\rangle \to 0$ in the limit of $R\to 0$.

We can now proceed to the calculation of $\langle {\varphi }^2\rangle$. By making use of Equations (48) and (53), as well as the values (55), it can be represented in the form

$$\langle {\varphi }^2\rangle =\frac{1-{\mathrm{e}}^{-\gamma m^2{R}^2}}{120{\left(4\pi \right)}^2{m}^2}·\frac{I_2+2I_3}{I_1},$$

(56)

where we have denoted by $I_1$, $I_2$, and $I_3$ the following Euclidean integrals:

$$I_1\equiv \underset{r<R}{\int }{d}^{3}x\sqrt{g},I_2\equiv \underset{r<R}{\int }{d}^{3}x\sqrt{g}{\mathcal{R}}^2,I_3\equiv \underset{r<R}{\int }{d}^{3}x\sqrt{g}{\mathcal{R}}_{\mu \nu }{\mathcal{R}}^{\mu \nu }.$$

(57)

Using the inequalities (54) and (B10), we observe that $mR\ll 1$, which, in turn, yields

$$\langle {\varphi }^2\rangle \simeq \frac{\gamma R^2}{120{\left(4\pi \right)}^2}·\frac{I_2+2I_3}{I_1}.$$

Referring the reader for the details to Ref. [86], we use here the result of the calculation of the integrals $I_1$, $I_2$, and $I_3$ in the limit of $y\ll 1$, where $y=c{R}^2$, and the parameter c is defined by Equation (B8). It yields

$$\langle {\varphi }^2\rangle \simeq \frac{\gamma R^2}{120{\left(4\pi \right)}^2}·27c^2=\frac{\gamma }{10}{\left(\epsilon GR\right)}^2.$$

(58)

Thus, in the small-m limit (54), we have obtained an m-independent expression for the condensate $\langle {\varphi }^2\rangle$. Its parametric dependence resembles that of the chiral quark condensate in QCD, given by Equation (44). Indeed, according to Equations (B8) and (B17), $\mathcal{R}\propto \epsilon G$, so that $\langle {\varphi }^2\rangle \propto {\mathcal{R}}^2$, in agreement with the initial Equation (53). Thus, ${\mathcal{R}}^2$ is analogous to $\langle {(g{E}_i^a)}^2\rangle$ in what concerns the role played by the gravitational and the chromo-electric fields as catalysts of, respectively, the ϕ-field and the quark-field condensation. We also notice that, in the limit of $y\ll 1$ at issue, the size R of the object is analogous to the correlation length λ, i.e., to the shortest distance scale of the non-perturbative chromo-electric vacuum. The parameter y resembles the quantity [11,12,13,14,15,16] ${\lambda }^2{\langle {\left(g{E}_i^a\right)}^2\rangle }^{1/2}\simeq 0.1$, whose numerical smallness, provided by the smallness of λ, ensures the convergence of the cumulant expansion in QCD. In our model, the smallness of y ensures that the leading contribution to $\langle {\varphi }^2\rangle$ is given by Equation (53).

Let us now consider the back-reaction produced by the $\langle {\varphi }^2\rangle$-condensate on the energy density ε. To this end, we notice that the contribution of the condensate (58) to the v.e.v. of the trace of the energy-momentum tensor is $\langle T^{\mu }{}_{\mu }\rangle =m^2\langle {\varphi }^2\rangle \propto {\left(m\epsilon GR\right)}^2$. This contribution should be added to the classical expression for the trace of the energy-momentum tensor, $T^{\mu }{}_{\mu }=\epsilon -3p$. As one can see, for the case of sufficiently small y’s at issue, $3p$ entering $T^{\mu }{}_{\mu }$ can be disregarded w.r.t. ε. Indeed, the pressure p, given by (B7), can be approximated by its average value

$$\overline{p}\equiv \frac{1}{R}{\int }_0^{R}drp=\frac{\epsilon }{2}{\int }_0^1\frac{dt}{\sqrt{t}}\frac{\sqrt{1-yt}-\sqrt{1-y}}{3\sqrt{1-y}-\sqrt{1-yt}},$$

where we have changed the integration variable r to $t\equiv \frac{z}{y}$, with z defined through Equation (B9). Calculating the latter integral, we obtain

$$\frac{3\overline{p}}{\epsilon }=-3+6·\frac{3\sqrt{1-y}\left[\mathrm{arccot}\sqrt{\frac{8-9y}{y}}+arctan\left(3\sqrt{\frac{y}{8-9y}}\right)\right]-\sqrt{8-9y}arcsin\sqrt{y}}{\sqrt{9y+1+\frac{1}{y-1}}}.$$

This is a monotonic function of y, which increases linearly at $y\ll 1$, remaining smaller than 1 for $y\lesssim 0.65$. Thus, for $y\ll 1$, one can indeed approximate $T^{\mu }{}_{\mu }$ by ε. Accordingly, it is ε that receives through the quantum correction $\langle T^{\mu }{}_{\mu }\rangle$ a contribution $\propto {\left(m\epsilon GR\right)}^2$.

Recalculating now $\langle {\varphi }^2\rangle$ with the so-corrected ε, and further iterating this procedure, we arrive at the equation

$$\frac{d{\epsilon }_n}{dn}\propto {\left(mGR\right)}^2{\epsilon }_n^2,$$

(59)

where we have denoted by n the cardinal of iterations, and approximately replaced $\mathsf{\Delta }n$ by $dn$. The solution to Equation (59), $-\frac{1}{{\epsilon }_n}+\frac{1}{\epsilon }\propto {\left(mGR\right)}^{2}n$, yields a value of n,

$$n_{*}\propto \frac{1}{\epsilon ·{\left(mGR\right)}^2},$$

(60)

which looks critical in the sense that ${\epsilon }_n\to \infty$ for $n\to n_{*}$. However, this does not happen, i.e., the energy density does not experience an infinite increase for $n\to n_{*}$. Indeed, as soon as ${\epsilon }_n$ starts increasing, the radius of the object also becomes n-dependent, and scales according to Equation (B10) as $R_n\sim \frac{1}{\sqrt{{\epsilon }_{n}G}}$. Therefore, Equation (59) at $n\to n_{*}$ takes the form $\frac{d{\epsilon }_n}{dn}\propto m^{2}G{\epsilon }_n$, with the solution

$${\epsilon }_n\propto \epsilon ·{\mathrm{e}}^{m^{2}Gn}.$$

(61)

The mass m of the ϕ-field is bounded from above according to the inequality (54), where ε should not be replaced by ${\epsilon }_n$, since ε is in any case smaller than ${\epsilon }_n$. Therefore, Equation (60) yields

$$n_{*}\gtrsim \frac{1}{{\epsilon }^2{G}^3{R}^2}\gtrsim \frac{1}{\epsilon G^2}.$$

Substituting this estimate into Equation (61), we obtain

$${\epsilon }_n\propto \epsilon ·{\mathrm{e}}^{\mathrm{const}·\frac{m^2}{\epsilon G}}\sim \epsilon ,$$

where $\mathrm{const}\sim 1$, and at the last step we have used the inequality (54) once again. Thus, at $n\to n_{*}$, ${\epsilon }_n$ remains of the order of ε, i.e., an infinite increase of ${\epsilon }_n$ does not occur.

In conclusion, we have provided an interesting semi-classical analogue of the relation that holds in QCD between the chiral and the chromo-electric condensates. Within this analogue, the role of the chromo-electric condensate is played by the squared curvature of the classical gravitational field produced by a spherically symmetric object of a constant energy density, while the role of the chromo-electric vacuum correlation length is played by the size of that object, which is considered to be sufficiently small. Finally, by estimating the back-reaction of the so-obtained scalar-field condensate on the energy density ε, recalculating the condensate with such a corrected ε, and iterating this procedure, we have shown that the resulting ${\epsilon }_n$ remains of the order of the initial ε, i.e., no instability of the system, which could be associated with an infinite increase of its energy density, occurs.

4. Free Energy of the Gluon Plasma in the High-Temperature Limit

In this Section, we address an important issue regarding the leading correction to the Stefan-Boltzmann law for the free-energy density of the gluon plasma at high temperatures. As we will see, this correction has the order (In this Section, we denote for brevity the finite-temperature Yang-Mills coupling $g_T$ simply as g.) $\mathcal{O}\left(g^2\right)$ for $N\sim 1$, while this order changes to $\mathcal{O}\left({\lambda }^{3/2}\right)$ for $N\gg 1$, where $\lambda =g^{2}N$ is the finite-temperature ’t Hooft coupling. The corrections to the Stefan-Boltzmann law stem from the spatial confinement of gluons constituting the plasma, as well as from the Polyakov loop. For our analysis, we will use the method developed in Refs. [88,89]. We start with representing the partition function of the finite-temperature Euclidean Yang-Mills theory in the form

$$\mathcal{Z}\left(T\right)=〈\int \mathcal{D}{a}_{\mu }^{a}exp\left[-\frac{1}{4g^2}{\int }_0^{\beta }d{x}_4{\int }_V^{}{d}^{3}x{\left(F_{\mu \nu }^a\left[A\right]\right)}^2\right]〉,$$

(62)

where $\beta \equiv 1/T$, and V is the three-dimensional volume occupied by the system. In Equation (62), we have modeled spatial confinement of the $a_{\mu }^a$-gluons by means of the stochastic background fields $B_{\mu }^a$. For this purpose, the full Yang-Mills field $A_{\mu }^a$ has been represented as a sum $A_{\mu }^a=B_{\mu }^a+a_{\mu }^a$, and the stochastic field $B_{\mu }^a$ has been averaged over with some measure $\langle \dots \rangle$. Clearly, at finite temperature T, both the $a_{\mu }^a$-and the $B_{\mu }^a$-fields obey the periodic boundary conditions $a_{\mu }^a(\mathbf{x},\beta )=a_{\mu }^a(\mathbf{x},0)$ and $B_{\mu }^a(\mathbf{x},\beta )=B_{\mu }^a(\mathbf{x},0)$. Integrating over the $a_{\mu }^a$-gluons in the Gaussian approximation, and disregarding for simplicity gluon spin degrees of freedom, one obtains

$$\mathcal{Z}\left(T\right)=\langle {\left\{det\left[-{\left(D_{\mu }^a\left[B\right]\right)}^2\right]\right\}}^{-\frac{1}{2}·2(N^2-1)}\rangle =\langle exp\left\{-(N^2-1)\mathrm{Tr}ln\left[-{\left(D_{\mu }^a\left[B\right]\right)}^2\right]\right\}\rangle ,$$

(63)

with the covariant derivative ${\left(D_{\mu }\left[B\right]f_{\nu }\right)}^a={\partial }_{\mu }{f}_{\nu }^a+f^{abc}{B}_{\mu }^b{f}_{\nu }^c$. Equation (63) includes color degrees of freedom of the $a_{\mu }^a$-gluons, and accounts for their $2(N^2-1)$ physical polarizations. In the one-loop approximation for the $a_{\mu }^a$-field, this equation can be simplified further:

$$\mathcal{Z}\left(T\right)\simeq exp\left\{-(N^2-1)\langle \mathrm{Tr}ln\left[-{\left(D_{\mu }^a\left[B\right]\right)}^2\right]\rangle \right\}.$$

(64)

In Equation (64), “Tr” includes the trace “tr” over color indices and the functional trace over space-time coordinates.

The free-energy density $F\left(T\right)$ is defined through the standard formula

$$\beta VF\left(T\right)=-ln\mathcal{Z}\left(T\right).$$

(65)

Using further for $ln\left[-{\left(D_{\mu }^a\left[B\right]\right)}^2\right]$ the proper-time representation, one has

$$F\left(T\right)=-(N^2-1)·2\sum _{n=1}^{\infty }{\int }_0^{\infty }\frac{ds}{s}\int \mathcal{D}{z}_{\mu }{\mathrm{e}}^{-\frac{1}{4}{\int }_0^{s}d\tau {\dot{z}}_{\mu }^2}\langle W\left[z_{\mu }\right]\rangle .$$

(66)

The integration in Equation (66) is performed over trajectories $z_{\mu }\left(\tau \right)$ which obey the periodic boundary conditions: $z_4\left(s\right)=z_4\left(0\right)+\beta n$ and $\mathbf{z}\left(s\right)=\mathbf{z}\left(0\right)$. The vector-function $z_{\mu }\left(\tau \right)$ describes therefore only the shape of the trajectory, while the factor $\beta V$ on the left-hand side of Equation (65) stems from the integration over positions of the trajectories. Furthermore, the summation over the winding number n yields a factor of 2, which accounts for winding modes with $n<0$. The zero-temperature part of the free-energy density, corresponding to the zeroth winding mode, has been subtracted. Finally, the Wilson loop that enters Equation (66), reads $W\left[z_{\mu }\right]=\frac{1}{N^2-1}\mathrm{tr}\mathcal{P}exp\left(i\oint d{z}_{\mu }{B}_{\mu }\right)$, where $B_{\mu }=B_{\mu }^a{t}^a$, and ${\left(t^a\right)}^{bc}=-i{f}^{abc}$ is a generator of the adjoint representation of the group SU(N).

According to the lattice data [44,45], the correlation function $\langle g^2{H}_i\left(x\right)H_k\left(x^{\prime }\right)\rangle$ exceeds by an order of magnitude the correlation function $\langle g^2{E}_i\left(x\right)H_k\left(x^{\prime }\right)\rangle$. This fact allows one to approximately factorize $\langle W\left[z_{\mu }\right]\rangle$ as $\langle W\left[z_{\mu }\right]\rangle \simeq \langle W\left[\mathbf{z}\right]\rangle {\displaystyle \prod _{n=-\infty }^{+\infty }}\langle P^n\rangle$, where $\langle W\left[\mathbf{z}\right]\rangle =〈\frac{1}{N^2-1}\mathrm{tr}\mathcal{P}exp\left(i\oint d{z}_k{B}_k\right)〉$ is the averaged purely spatial Wilson loop, and $\langle P^n\rangle =〈\frac{1}{N^2-1}\mathrm{tr}\mathcal{T}exp\left(in{\int }_0^{\beta }d{z}_4{B}_4\right)〉$ is a generalization of the Polyakov loop to the case of n windings. Upon this factorization, the world-line integral over $z_4\left(\tau \right)$ in Equation (66) becomes that of a free particle, which yields

$$F\left(T\right)=-2(N^2-1)\sum _{n=1}^{\infty }{\int }_0^{\infty }\frac{ds}{s}\frac{{\mathrm{e}}^{-\frac{{\beta }^2{n}^2}{4s}}}{\sqrt{4\pi s}}\langle P^n\rangle \oint \mathcal{D}\mathbf{z}{\mathrm{e}}^{-\frac{1}{4}{\int }_0^{s}d\tau {\dot{\mathbf{z}}}^2}\langle W\left[\mathbf{z}\right]\rangle .$$

(67)

In order to calculate the world-line integral over $\mathbf{z}\left(\tau \right)$, we use again Equation (30) for the Wilson-loop average in the adjoint representation. This average can be written in the form

$$\langle W\left[\mathbf{z}\right]\rangle =\frac{1}{1+\frac{1}{N^2}}\left({\mathrm{e}}^{-\sigma \mathsf{\Sigma }}+\frac{1}{N^2}{\mathrm{e}}^{-c·g^2\frac{N}{3}T\sqrt{\mathsf{\Sigma }}}\right).$$

(68)

Here, Σ is the area of the minimal surface bounded by the contour $\mathbf{z}\left(\tau \right)$, and c is some positive dimensionless constant, which will be determined below. Furthermore, Equation (68) obeys the normalization condition $\langle W\left[\mathbf{0}\right]\rangle =1$. The second exponential on the right-hand side of Equation (68) represents the perimeter law ${\mathrm{e}}^{-mL}$, where $L={\int }_0^{s}d\tau \left|\dot{\mathbf{z}}\right|$ is the length of the contour $\mathbf{z}\left(\tau \right)$, and the constant m has the dimensionality of mass. Here, we have substituted L by $\sqrt{\mathsf{\Sigma }}$, and parametrized m through the soft scale $g^{2}NT$ as $m=c·g^2\frac{N}{3}T$. The spatial string tension σ in the adjoint representation can be expressed in terms of the spatial string tension ${\sigma }_{\mathrm{f}}$ in the fundamental representation by means of Casimir scaling: $\frac{\sigma }{{\sigma }_{\mathrm{f}}}=\frac{2N^2}{N^2-1}$. This ratio is equal to 9/4 for $N=3$, while going to 2 in the large-N limit. At temperatures $T>T_{*}$ of interest, where $T_{*}$ is the temperature of dimensional reduction, one can express ${\sigma }_{\mathrm{f}}$ in terms of the string tension in the 3D Yang-Mills theory, which was calculated analytically in Refs. [78,79]. The corresponding expression for ${\sigma }_{\mathrm{f}}$ reads (Note that, for $N=3$, the coefficient $\frac{1}{\pi }\simeq 0.32$ in this formula agrees remarkably well with the value of $0.{566}^2$, which was used in Refs. [90,91] for the parametrization of ${\sigma }_{\mathrm{f}}$ at high temperatures.) ${\sigma }_{\mathrm{f}}=\frac{N^2-1}{8\pi }{\left(g^{2}T\right)}^2$, which yields the following spatial string tension in the adjoint representation: $\sigma =\frac{1}{4\pi }{\left(g^{2}NT\right)}^2$.

Hence, the free-energy density (67) can be written in the form $F=F_1+F_2$, where the term $F_1$ corresponds to the exponential ${\mathrm{e}}^{-\sigma \mathsf{\Sigma }}$ from Equation (68), while the term $F_2$ corresponds to the exponential ${\mathrm{e}}^{-c·g^2\frac{N}{3}T\sqrt{\mathsf{\Sigma }}}$ from the same equation. Clearly, in the large-N limit, $F_1\gg F_2$ due to the relative factor of $\frac{1}{N^2}$, so that the thermodynamics of the gluon plasma in that limit is fully determined by spatial confinement. Therefore, let us start with calculating the world-line integral $I\equiv \oint \mathcal{D}\mathbf{z}{\mathrm{e}}^{-\frac{1}{4}{\int }_0^{s}d\tau {\dot{\mathbf{z}}}^2-\sigma \mathsf{\Sigma }}$, which enters the term $F_1$. To this end, we implement for the minimal area Σ the following ansatz: $\mathsf{\Sigma }=\frac{1}{2}{\int }_0^{s}d\tau |\mathbf{z}\times \dot{\mathbf{z}}|$. It corresponds to a parasol-shaped surface made of thin segments. Furthermore, since ${\int }_0^{s}d\tau \mathbf{z}=0$, the point where the segments merge is the origin. Therefore, the chosen ansatz for Σ automatically selects from all cone-shaped surfaces bounded by $\mathbf{z}\left(\tau \right)$ the one of the minimal area. We use further the so-called non-backtracking approximation $\mathsf{\Sigma }\simeq \sqrt{{\mathbf{f}}^2}$, where $\mathbf{f}\equiv \frac{1}{2}{\int }_0^{s}d\tau (\mathbf{z}\times \dot{\mathbf{z}})$. This approximation is widely used in order to simplify the parametrizations of minimal surfaces allowing for an analytic calculation of the corresponding world-line integrals [7,8,9]. (In general, $\frac{1}{2}{\int }_0^{s}d\tau |\mathbf{z}\times \dot{\mathbf{z}}|$ can be larger than $\sqrt{{\mathbf{f}}^2}$. This happens if, in the course of its evolution in spatial directions, the gluon performs backward and/or non-planar motions. Once this happens, the vector product $(\mathbf{z}\times \dot{\mathbf{z}})$ changes its direction, and the integral ${\int }_0^{s}d\tau (\mathbf{z}\times \dot{\mathbf{z}})$ receives mutually cancelling contributions.) One can now calculate the integral I by representing the exponential ${\mathrm{e}}^{-\sigma \mathsf{\Sigma }}$ as ${\mathrm{e}}^{-\sigma \mathsf{\Sigma }}={\int }_0^{\infty }\frac{d\lambda }{\sqrt{\pi \lambda }}{\mathrm{e}}^{-\lambda -\frac{{\sigma }^2{\mathbf{f}}^2}{4\lambda }}$, and introducing further an auxiliary space-independent magnetic field $\mathbf{H}$ according to the formula

$${\mathrm{e}}^{-A{\mathbf{f}}^2}=\frac{1}{{\left(4\pi A\right)}^{3/2}}\int d^{3}H{\mathrm{e}}^{-\frac{{\mathbf{H}}^2}{4A}+i\mathbf{H}\mathbf{f}},\mathrm{where}A>0.$$

(69)

The world-line integral gets then reduced to the one for a spinless particle of an electric charge 1 interacting with the constant magnetic field $\mathbf{H}$, i.e., to the bosonic Euler-Heisenberg-Schwinger Lagrangian, which has the form [28,29,30,31]

$$\oint \mathcal{D}\mathbf{z}{\mathrm{e}}^{-\frac{1}{4}{\int }_0^{s}d\tau {\dot{\mathbf{z}}}^2+i\mathbf{H}\mathbf{f}}=\frac{1}{{\left(4\pi s\right)}^{3/2}}\frac{Hs}{sinh\left(Hs\right)}.$$

(70)

Integrating further over λ, we obtain for the world-line integral at issue:

$$I=\frac{\sigma }{2{\pi }^{5/2}\sqrt{s}}{\int }_0^{\infty }dH\frac{H^3/sinh\left(Hs\right)}{{(H^2+{\sigma }^2)}^2}.$$

In the case of $N=3$, the corresponding free-energy density reads

$$F_1{|}_{N=3}=-\frac{18\sigma }{5{\pi }^3}\sum _{n=1}^{\infty }{\int }_0^{\infty }\frac{ds}{s^2}{\mathrm{e}}^{-\frac{{\beta }^2{n}^2}{4s}}\langle P^n\rangle {\int }_0^{\infty }dH\frac{H^3/sinh\left(Hs\right)}{{(H^2+{\sigma }^2)}^2}.$$

To perform the perturbative expansion of this expression, we introduce a dimensionless integration variable $h=H/\sigma$. Furthermore, we notice that, in the high-temperature limit of interest, $\langle P^n\rangle \simeq \langle P\rangle$, where [92]

$$\langle P\rangle =1+\mathcal{O}\left(g^3\right).$$

(71)

To find the order of the leading g-dependent term of the perturbative expansion, we use the approximation $sinh\left(\sigma hs\right)\simeq \sigma hs·\left(1+\frac{{\left(\sigma hs\right)}^2}{6}\right)$, which yields for $F_1{|}_{N=3}$ the following expression:

$$F_1{|}_{N=3}\simeq -\frac{9\langle P\rangle }{10{\pi }^2}\sum _{n=1}^{\infty }{\int }_0^{\infty }\frac{ds}{s^3}·\frac{{\mathrm{e}}^{-\frac{{\beta }^2{n}^2}{4s}}}{{\left(1+\frac{\sigma s}{\sqrt{6}}\right)}^2}.$$

Approximating further the sum over winding modes by the first two terms, we obtain

$$F_1{|}_{N=3}\simeq -\frac{9\langle P\rangle T^4}{10{\pi }^2}\left[17-\frac{10}{\sqrt{6}}\sigma {\beta }^2+\mathcal{O}\left({\left(\sigma {\beta }^2\right)}^2\right)\right].$$

(72)

Clearly, since $\sigma =\mathcal{O}\left(g^4\right)$, the obtained term $-\frac{10}{\sqrt{6}}\sigma {\beta }^2$ also has the order $\mathcal{O}\left(g^4\right)$. Nevertheless, due to Equation (71), the order of the leading g-dependent term of the perturbative expansion of $F_1{|}_{N=3}$ is 3, rather than 4.

We proceed now to the calculation of the free-energy density $F_2$ for $N=3$, which will allow us to find the value of the constant c in Equation (68). The corresponding world-line integral $\oint \mathcal{D}\mathbf{z}{\mathrm{e}}^{-\frac{1}{4}{\int }_0^{s}d\tau {\dot{\mathbf{z}}}^2-c{g}^{2}T\sqrt{\mathsf{\Sigma }}}$ can be calculated by using again the approximation $\mathsf{\Sigma }\simeq \sqrt{{\mathbf{f}}^2}$. The fourth root in the so-emerging exponential, ${\mathrm{e}}^{-c{g}^{2}T\sqrt[4]{{\mathbf{f}}^2}}$, can be got rid of by using two identical auxiliary integrations as follows:

$${\mathrm{e}}^{-c{g}^{2}T\sqrt[4]{{\mathbf{f}}^2}}=\frac{1}{\pi }{\int }_0^{\infty }\frac{d\lambda }{\sqrt{\lambda }}{\int }_0^{\infty }\frac{d\mu }{\sqrt{\mu }}{\mathrm{e}}^{-\lambda -\mu -\frac{{\left(c{g}^{2}T\right)}^4{\mathbf{f}}^2}{64{\lambda }^2\mu }}.$$

Introducing now once again the auxiliary magnetic field $\mathbf{H}$ according to the formula (69), we obtain for the exponential at issue the following representation:

$${\mathrm{e}}^{-c{g}^{2}T\sqrt[4]{{\mathbf{f}}^2}}=\frac{64}{{\pi }^{5/2}}\frac{1}{{\left(c{g}^{2}T\right)}^6}{\int }_0^{\infty }d\lambda {\lambda }^{5/2}{\mathrm{e}}^{-\lambda }{\int }_0^{\infty }d\mu \mu {\mathrm{e}}^{-\mu }\int d^{3}H{\mathrm{e}}^{-\frac{16{\lambda }^2\mu }{{\left(c{g}^{2}T\right)}^4}{\mathbf{H}}^2+i\mathbf{H}\mathbf{f}}.$$

Performing the functional $\mathbf{z}$-integration as in Equation (70), and integrating further over μ, which can be done analytically, we obtain the following intermediate expression:

$$F_2{|}_{N=3}=-\frac{256T^4}{{\pi }^{7/2}}{\xi }^2\sum _{n=1}^{\infty }{\int }_0^{\infty }\frac{ds}{s^2}{\mathrm{e}}^{-\frac{n^2}{4s}}\langle P^n\rangle {\int }_0^{\infty }dh\frac{h^3}{sinh\left({\xi }^{2}hs\right)}{\int }_0^{\infty }d\lambda \frac{{\lambda }^{5/2}{\mathrm{e}}^{-\lambda }}{{(16{\lambda }^2{h}^2+1)}^2}.$$

(73)

Here, we have denoted $\xi \equiv c{g}^2$, $h\equiv H/{\left(\xi T\right)}^2$, and made s dimensionless by rescaling it as $s_{\mathrm{new}}=T^2{s}_{\mathrm{old}}$. By using the approximation $sinh\left({\xi }^{2}hs\right)\simeq {\xi }^{2}hs[1+{\left({\xi }^{2}hs\right)}^2/6]$, we have

$$F_2{|}_{N=3}\simeq -\frac{256T^4}{{\pi }^{7/2}}\sum _{n=1}^{\infty }{\int }_0^{\infty }\frac{ds}{s^3}{\mathrm{e}}^{-\frac{n^2}{4s}}{\int }_0^{\infty }d\lambda {\lambda }^{5/2}{\mathrm{e}}^{-\lambda }{\int }_0^{\infty }dh\frac{h^2}{{(16{\lambda }^2{h}^2+1)}^2}·\frac{1}{1+{\left({\xi }^{2}hs\right)}^2/6}.$$

The h-integration in this formula can be performed analytically, which yields

$$F_2{|}_{N=3}\simeq -\frac{16T^4}{{\pi }^{5/2}}\sum _{n=1}^{\infty }{\int }_0^{\infty }\frac{ds}{s^3}{\mathrm{e}}^{-\frac{n^2}{4s}}{\int }_0^{\infty }d\lambda \frac{{\lambda }^{3/2}{\mathrm{e}}^{-\lambda }}{{(4\lambda +{\xi }^{2}s/\sqrt{6})}^2}.$$

Approximating again the sum over winding modes by the first two terms, we further have

$${\int }_0^{\infty }\frac{ds}{s^3}\left({\mathrm{e}}^{-\frac{1}{4s}}+{\mathrm{e}}^{-\frac{1}{s}}\right){\int }_0^{\infty }d\lambda \frac{{\lambda }^{3/2}{\mathrm{e}}^{-\lambda }}{{(4\lambda +{\xi }^{2}s/\sqrt{6})}^2}=\frac{17\sqrt{\pi }}{16}-\frac{27{\pi }^{3/2}}{128·6^{1/4}}·\xi +\mathcal{O}\left({\xi }^2\right).$$

This yields the sought free-energy density

$$F_2{|}_{N=3}\simeq -\frac{\langle P\rangle T^4}{10{\pi }^2}\left(17-\frac{27\pi }{8·6^{1/4}}·c{g}^2\right).$$

(74)

Once brought together, Equations (72) and (74) yield

$$F{|}_{N=3}\simeq -\frac{\langle P\rangle T^4}{{\pi }^2}\left[17-\frac{27\pi }{80·6^{1/4}}·c{g}^2-\frac{9}{\sqrt{6}}\sigma {\beta }^2+\mathcal{O}\left({\left(\sigma {\beta }^2\right)}^2\right)\right].$$

(75)

The two leading terms of this expression can be compared with the known perturbative expansion of the free-energy density [93,94],

$$F_{\mathrm{pert}}{|}_{N=3}=-\frac{8{\pi }^2{T}^4}{45}\left[1-\frac{15g^2}{16{\pi }^2}+\mathcal{O}\left(g^3\right)\right].$$

(76)

Comparing the leading term of Equation (75), $-\frac{17T^4}{{\pi }^2}\simeq -1.72T^4$, with the Stefan-Boltzmann expression represented by the leading term of Equation (76), $-\frac{8{\pi }^2{T}^4}{45}\simeq -1.75T^4$, we conclude that the above-used approximation of the full sum over winding modes by the $(n=1)$- and the $(n=2)$-terms is very good. Comparing further with each other the $\mathcal{O}\left(g^2\right)$-terms of Equations (75) and (76), we obtain:

$$c=\frac{80\pi }{27·6^{3/4}}\simeq 2.4.$$

(77)

By using Equation (68), we proceed now to arbitrary N, which yields

$$F_1=-\frac{T^4}{8{\pi }^2}·\frac{N^2(N^2-1)}{N^2+1}\left(1+\frac{{\lambda }^{3/2}}{8\pi \sqrt{3}}\right)\left(17-\frac{5{\lambda }^2}{2\pi \sqrt{6}}+\mathcal{O}\left({\lambda }^4\right)\right)$$

and

$$F_2=-\frac{T^4}{8{\pi }^2}·\frac{N^2-1}{N^2+1}\left(1+\frac{{\lambda }^{3/2}}{8\pi \sqrt{3}}\right)\left(17-\frac{9\pi c}{8·6^{1/4}}·\lambda +\mathcal{O}\left({\lambda }^2\right)\right).$$

Here, $\lambda =g^{2}N$ is the so-called ’t Hooft coupling, which stays finite in the large-N limit, and we have approximated the Polyakov loop as (cf. Equation (71)) [95] $\langle P\rangle \simeq 1+\frac{{\lambda }^{3/2}}{8\pi \sqrt{3}}$. According to the lattice data, this appproximation holds at temperatures as high as $T\gtrsim 10T_c$, while $\langle P\rangle$ stays $\lesssim 0.9$ even at $T=4T_c$ (cf. the lattice data quoted in Ref. [95]). Hence, for the rest of this Section, we imply the said limit of very high temperatures, where the analytic coupling-dependent expression for $\langle P\rangle$ can only be provided.

Accordingly, we obtain for the full free-energy density $F=F_1+F_2$:

$$F=$$

$$=-\frac{T^4}{8{\pi }^2}·\frac{N^2(N^2-1)}{N^2+1}\left(1+\frac{{\lambda }^{3/2}}{8\pi \sqrt{3}}\right)\left[17\left(1+\frac{1}{N^2}\right)-\frac{9\pi c}{8·6^{1/4}}·\frac{\lambda }{N^2}-\frac{5{\lambda }^2}{2\pi \sqrt{6}}+\mathcal{O}\left({\lambda }^4\right)+\mathcal{O}\left(\frac{{\lambda }^2}{N^2}\right)\right].$$

In the large-N limit of this expression, the c-dependent term, which corresponds to the leading perturbative correction from Equation (76), gets $\frac{1}{N^2}$-suppressed in comparison with the $\mathcal{O}\left({\lambda }^2\right)$-term, which corresponds to the leading σ-dependent correction from Equation (72). This result follows, of course, from the relative factor of $\frac{1}{N^2}$ between the perimeter- and the area-law exponentials in Equation (68). The large-N limit of the free-energy density thus reads

$$F=-\frac{T^4{N}^2}{8{\pi }^2}\left(1+\frac{{\lambda }^{3/2}}{8\pi \sqrt{3}}\right)\left[17-\frac{5{\lambda }^2}{2\pi \sqrt{6}}+\mathcal{O}\left(\frac{10}{N^2}\right)+\mathcal{O}\left(\frac{10\lambda }{N^2}\right)+\mathcal{O}\left({\lambda }^4\right)+\mathcal{O}\left(\frac{{\lambda }^2}{N^2}\right)\right].$$

Accordingly, $-\frac{F}{T^4}$ could have a maximum corresponding to the most probable configuration of the system, once the relation $\frac{17{\lambda }^{3/2}}{8\pi \sqrt{3}}=\frac{5{\lambda }^2}{2\pi \sqrt{6}}$ would hold, i.e., at $\lambda =\frac{289}{8}$. However, since this value of λ is much larger than unity, it lies outside the range of applicability of the λ-expansion, so that such a maximum of $-\frac{F}{T^4}$ is not realized. Thus, the main qualitative result of our study is that the leading coupling-dependent correction to the Stefan-Boltzmann expression, while being $\mathcal{O}\left(\lambda \right)$ for $N\sim 1$, goes over to a $\mathcal{O}\left({\lambda }^{3/2}\right)$-term of the opposite sign for $N\gg 1$ and $T\gtrsim 10T_c$.

5. Yang-Mills Running Coupling with an IR-Stable Fixed Point

In this Section, we will follow Ref. [96] to provide one more illustration of the importance of confinement effects in the one-loop vacuum amplitudes. Namely, by using a parametrization of the Wilson-loop average with the minimal-area law, we will calculate the polarization operator of a gluon which propagates in the confining background. This calculation yields an infra-red finiteness (also called infra-red freezing) of the running strong coupling. It turns out that this phenomenon takes place not only in the hadronic phase, but also in the deconfinement phase, where it is caused by the above-discussed spatial confinement. In particular, we will demonstrate how the momentum scale defining the onset of freezing can be obtained both analytically and numerically. As a phenomenological application of these results, we will evaluate the non-perturbative contribution to the so-called thrust variable, which makes the value of this variable closer to the experimental one.

To perform the mentioned calculation of the polarization operator, we parametrize the minimal area entering the area law of the Wilson-loop average, in terms of the Wilson-loop contour. The main idea of such a parametrization is to convert the proper time in the path integral into a length coordinate, which we denote as τ. After that, one can naturally parametrize the minimal area ${\mathsf{\Sigma }}_{\min }$ as a transverse distance $\left|\mathbf{r}\right(\tau \left)\right|$ integrated along this coordinate:

$${\mathsf{\Sigma }}_{\min }={\int }_0^{R}d\tau \left|\mathbf{r}\left(\tau \right)\right|.$$

(78)

With such a parametrization of ${\mathsf{\Sigma }}_{\min }$, we calculate the polarization operator of a gluon as a function of the distance R. In the absence of the confining background, i.e., for a free gluon, the polarization operator yields the standard Yang-Mills one-loop running coupling [97]. In the presence of the background, the running coupling ${\alpha }_s\left(p\right)=\frac{g^2\left(p\right)}{4\pi }$ goes to a constant in the infra-red region [96]:

$${\alpha }_s\left(p\right)=\frac{4\pi }{bln\frac{p^2}{{\mathsf{\Lambda }}^2}}\to \frac{4\pi }{\tilde{b}ln\frac{p^2+m^2}{{\mathsf{\Lambda }}^2}}.$$

(79)

Here,

$$b=\frac{11}{3}N$$

is the absolute value of the first coefficient of the Yang-Mills β-function, and $m\propto {\left(\mathrm{string}\mathrm{tension}\right)}^{1/2}$ is a non-perturbative mass parameter. Below, we will calculate the values of m and $\tilde{b}$ both analytically and numerically. Note also that, in the absence of the confining background, one can use various regularization schemes (such as the $MS$ or the $\overline{MS}$ ones) [98], which is no longer possible once the confining background is taken into account. (Indeed, within these schemes, the gluon propagator at a certain value of $p^2$, goes over to the propagator of a free gluon, which is not possible in the presence of the background.) To circumvent this problem, one uses the definition of renormalized ${\alpha }_s\left(p\right)$ through the static potential [99,100], where its Fourier image ${\alpha }_s\left(R\right)$ defines the Coulomb interaction $-\frac{4}{3}\frac{{\alpha }_s\left(R\right)}{R}$ at the distance R between the static quarks. Since the static potential can be extracted from the Wilson-loop average, it looks natural to adopt this definition within the above-mentioned approach, where the Wilson-loop average appears in the course of the calculation of the polarization operator [99]. Hence, it is this definition of ${\alpha }_s\left(p\right)$ which is assumed in the below calculation.

Let us first recollect some steps of the derivation of ${\alpha }_s\left(p\right)$, which is based on the integration over free quantum fluctuations of the Yang-Mills field [97,98]. The renormalization-group procedure starts with splitting the total Yang-Mills field $A_{\mu }^a$ into a low-momentum part ${\overline{A}}_{\mu }^a$ and a quantum fluctuation $a_{\mu }^a$, whose momentum is thus larger than that of the field ${\overline{A}}_{\mu }^a$. Hence, one can substitute the ansatz $A_{\mu }^a={\overline{A}}_{\mu }^a+a_{\mu }^a$ into the bare Yang-Mills action,

$$S_0\left[A\right]=\frac{1}{4g_0^2}\int d^{4}x{\left(F_{\mu \nu }^a\left[A\right]\right)}^2,$$

where $F_{\mu \nu }^a\left[A\right]$ is the Yang-Mills field-strength tensor, and $g_0$ is the bare coupling. Fixing the so-called background Feynman gauge ${\left(D_{\mu }{a}_{\mu }\right)}^a=0$, one adds to the action the term $S_{\mathrm{g}.\mathrm{f}.}=\frac{1}{2g_0^2}\int d^{4}x{\left[{\left(D_{\mu }{a}_{\mu }\right)}^a\right]}^2$, where ${\left(D_{\mu }{a}_{\nu }\right)}^a={\partial }_{\mu }{a}_{\nu }^a+f^{abc}{\overline{A}}_{\mu }^b{a}_{\nu }^c$. The full action then reads

$$S_0+S_{\mathrm{g}.\mathrm{f}.}=\frac{1}{4g_0^2}\int d^{4}x\left\{{\left(F_{\mu \nu }^a\left[\overline{A}\right]\right)}^2-4a_{\nu }^a{\left(D_{\mu }{F}_{\mu \nu }\left[\overline{A}\right]\right)}^a-2a_{\mu }^a\left[{\delta }_{\mu \nu }{\left(D^2\right)}^{ac}+2f^{abc}{F}_{\mu \nu }^b\left[\overline{A}\right]\right]a_{\nu }^c+\mathcal{O}\left(a^3\right)\right\}.$$

To perform the one-loop renormalization of the running coupling $g\left(p\right)$, one should integrate out the $a_{\mu }^a$-gluons which propagate in the ${\overline{A}}_{\mu }^a$-background. The full renormalized effective action can be written as

$$S=S_0+S^{\mathrm{dia}}+S^{\mathrm{para}}=\int \frac{d^{4}p}{{\left(2\pi \right)}^4}\frac{1}{4g^2\left(p\right)}{F}_{\mu \nu }^a\left(p\right)F_{\mu \nu }^a(-p).$$

(80)

In Equation (80), the diamagnetic part of the effective action for the ${\overline{A}}_{\mu }^a$-fields has the form

$$S^{\mathrm{dia}}=\mathrm{tr}ln(-D^2)=\mathrm{tr}\left[{(-{\partial }^2)}_{xx}^{-1}{\mathsf{\Delta }}^{\left(2\right)}\left(x\right)-\frac{1}{2}{(-{\partial }^2)}_{xy}^{-1}{\mathsf{\Delta }}^{\left(1\right)}\left(y\right){(-{\partial }^2)}_{yx}^{-1}{\mathsf{\Delta }}^{\left(1\right)}\left(x\right)\right],$$

(81)

where we have used the decomposition

$$-D^2=-{\partial }^2+{\mathsf{\Delta }}^{\left(1\right)}+{\mathsf{\Delta }}^{\left(2\right)}\mathrm{with}{\mathsf{\Delta }}^{\left(1\right)}\left(x\right)\equiv i{t}^a\left({\partial }_{\mu }{\overline{A}}_{\mu }^a+2{\overline{A}}_{\mu }^a{\partial }_{\mu }\right),{\mathsf{\Delta }}^{\left(2\right)}\left(x\right)\equiv {\left({\overline{A}}_{\mu }^a{t}^a\right)}^2,$$

and ${\left(t^a\right)}^{bc}=-i{f}^{abc}$ being a generator of the group SU(N) in the adjoint representation. The paramagnetic part of the effective action in Equation (80) reads

$$S^{\mathrm{para}}=\frac{1}{2}\mathrm{tr}ln\left[1+{(-{\partial }^2)}^{-1}\left(2f^{abc}{F}_{\mu \nu }^b\right)\right]=-\frac{1}{4}{(-{\partial }^2)}_{xy}^{-1}\left(2f^{abc}{F}_{\mu \nu }^b\left(y\right)\right){(-{\partial }^2)}_{yx}^{-1}\left(2f^{adc}{F}_{\mu \nu }^d\left(x\right)\right).$$

(82)

The action $S^{\mathrm{dia}}$ is associated with the orbital motion of the $a_{\mu }^a$-gluons in the chromo-magnetic field of the low-energy ${\overline{A}}_{\mu }^a$-gluons. The diamagnetic interaction of the $a_{\mu }^a$-and the ${\overline{A}}_{\mu }^a$-gluons leads to the screening of a test color charge, and exists in the Abelian case as well. Rather, the action $S^{\mathrm{para}}$ is associated with the alignment of “color spins” of the $a_{\mu }^a$-gluons along the chromo-magnetic field of the ${\overline{A}}_{\mu }^a$-gluons. The corresponding paramagnetic interaction of the ${\overline{A}}_{\mu }^a$-field with the color degrees of freedom of the $a_{\mu }^a$-gluons leads to the antiscreening of a test color charge, being a specific property of the non-Abelian gauge theories (cf. Section 2). The dia- and the paramagnetic parts of the one-loop effective action can be written as

$$S^{\left\{\genfrac{}{}{0pt}{}{\mathrm{dia}}{\mathrm{para}}\right\}}=N·\left\{\genfrac{}{}{0pt}{}{\frac{1}{12}}{(-1)}\right\}·\int \frac{d^{4}p}{{\left(2\pi \right)}^4}{F}_{\mu \nu }^a\left(p\right)F_{\mu \nu }^a(-p){\mathsf{\Pi }}_{\mathrm{free}}\left(p^2\right),$$

(83)

where

$${\mathsf{\Pi }}_{\mathrm{free}}\left(p^2\right)\equiv \frac{1}{16{\pi }^2}ln\frac{{\mathsf{\Lambda }}_0^2}{p^2}=\int \frac{d^{4}q}{{\left(2\pi \right)}^4}\frac{1}{q^2{(q+p)}^2}$$

(84)

is the free scalar polarization operator. Thus, the absolute value of the paramagnetic contribution to the one-loop vacuum polarization in 4D is 12 times larger than the value of the diamagnetic contribution. That is, the Yang-Mills vacuum is a color-paramagnetic medium (The color-paramagnetism of the Yang-Mills vacuum should not be confused with the usual diamagnetism of the QCD vacuum in the hadronic phase. That is, an external magnetic field decreases in the QCD vacuum due to the fact that the latter is populated by the electrically charged ${\pi }^{\pm }$-mesons. Such mesons thus form a diamagnetic medium. Instead, the quark-gluon plasma acts as an (ordinary) paramagnetic medium, which enhances an external magnetic field [101].), which antiscreens a test color charge, i.e., enhances its color field. For this reason, the Yang-Mills theory possesses the property of asymptotic freedom. In fact, Equation (80) yields

$$\frac{1}{g^2\left(p\right)}=\frac{1}{g_0^2}-\frac{b}{16{\pi }^2}ln\frac{{\mathsf{\Lambda }}_0^2}{p^2},$$

(85)

where $g_0\equiv g\left({\mathsf{\Lambda }}_0\right)$, and ${\mathsf{\Lambda }}_0$ is a bare cut-off. Introducing a renormalized cut-off $\mathsf{\Lambda }={\mathsf{\Lambda }}_{0}exp\left(-\frac{8{\pi }^2}{b{g}_0^2}\right)$, one arrives at the standard result [102,103] ${\alpha }_s\left(p\right)=\frac{4\pi }{bln\frac{p^2}{{\mathsf{\Lambda }}^2}}$.

In reality, however, gluon fluctuations $a_{\mu }^a$’s propagating along the loop, are self-confined rather than free, which leads to the formation of a colored bound state of two $a_{\mu }^a$-gluons. This confining interaction is associated with the formation of the string in the octet channel. Using Casimir scaling, we can relate the corresponding string tension to that of the fundamental representation as (see e.g., Ref. [104]) $\sigma =\frac{9}{8}{\sigma }_{\mathrm{f}}=0.225{\mathrm{GeV}}^2$ with ${\sigma }_{\mathrm{f}}={\left(440\mathrm{MeV}\right)}^2$. Confinement of gluon fluctuations $a_{\mu }^a$’s is likely to be caused by the interaction of these fluctuations with the stochastic background fields $B_{\mu }^a$’s, whose momenta are even smaller than the momenta of the ${\overline{A}}_{\mu }^a$-gluons [99]. These additional background fields can be accounted for by substituting into the Yang-Mills action a modified ansatz $A_{\mu }^a=B_{\mu }^a+{\overline{A}}_{\mu }^a+a_{\mu }^a$. Accordingly, the definition of the one-loop effective action, $S=-ln{\langle \int \mathcal{D}{a}_{\mu }^a{\mathrm{e}}^{-S\left[A\right]}\rangle }_B$, includes now some gauge- and $O\left(4\right)$-invariant average ${\langle \dots \rangle }_B$ over the background fields. In the course of the calculation of $\int \mathcal{D}{a}_{\mu }^a{\mathrm{e}}^{-S\left[A\right]}$, the presence of the background leads to the substitution ${\partial }^2\to D^2\left[B\right]$ in Equations (81) and (82), where ${\left(D_{\mu }\left[B\right]a_{\nu }\right)}^a={\partial }_{\mu }{a}_{\nu }^a+f^{abc}{B}_{\mu }^b{a}_{\nu }^c$. In the diagrammatic language, this means that we are still considering the same one-loop diagrams with only two external lines of the ${\overline{A}}_{\mu }^a$-field, but the internal $a_{\mu }^a$-loop of those diagrams receives now infinitely many contributions of the $B_{\mu }^a$-field. The latter appear through the world-line representation of the operator ${\left(D^2\left[B\right]\right)}^{-1}$, which contains the phase factor of the $a_{\mu }^a$-gluon.

Neglecting for simplicity the spin degrees of freedom of the $a_{\mu }^a$-gluon, we notice that the scalar polarization operator in the coordinate representation has the form

$${\mathsf{\Pi }}_{\mathrm{free}}\left(x\right)=\int \frac{d^{4}p}{{\left(2\pi \right)}^4}{\mathrm{e}}^{ipx}{\mathsf{\Pi }}_{\mathrm{free}}\left(p^2\right)=D_0^2\left(x\right),\mathrm{where}{D}_0\left(x\right)=\frac{1}{4{\pi }^2{x}^2}.$$

(86)

In the presence of the background field $B_{\mu }^a$, this expression goes over to

$${\langle \mathsf{\Pi }(x,y|B)\rangle }_B\equiv \mathsf{\Pi }(x,y)={\langle \mathrm{tr}{\left(D^2\left[B\right]\right)}_{xy}^{-1}{\left(D^2\left[B\right]\right)}_{yx}^{-1}\rangle }_B.$$

The world-line representation of this average reads

$$\mathsf{\Pi }(x,y)={\int }_0^{\infty }ds{\int }_0^{\infty }d\overline{s}\int {\left(\mathcal{D}{z}_{\mu }\right)}_{xy}{\left(\mathcal{D}{\overline{z}}_{\mu }\right)}_{yx}{\mathrm{e}}^{-{\int }_0^{s}d\lambda \frac{{\dot{z}}_{\mu }^2}{4}-{\int }_0^{\overline{s}}d\overline{\lambda }\frac{{\dot{\overline{z}}}_{\mu }^2}{4}}\times$$

$$\times {〈\mathrm{tr}\mathcal{P}exp\left[i\left({\int }_0^{s}d\lambda {\dot{z}}_{\mu }{B}_{\mu }^a\left(z\right)t^a+{\int }_0^{\overline{s}}d\lambda {\dot{\overline{z}}}_{\mu }{B}_{\mu }^a\left(\overline{z}\right)t^a\right)\right]〉}_B,$$

(87)

where we have introduced the short-hand notation $\int {\left(\mathcal{D}{z}_{\mu }\right)}_{xy}\equiv \underset{z\left(0\right)=y}{\overset{z\left(s\right)=x}{\int }}\mathcal{D}{z}_{\mu }$. Due to confinement of the $a_{\mu }^a$-gluon, the Wilson-loop average entering Equation (87) should respect the minimal-area law: ${\langle \dots \rangle }_B\simeq {\mathrm{e}}^{-\sigma {\mathsf{\Sigma }}_{\min }}$, where ${\mathsf{\Sigma }}_{\min }={\mathsf{\Sigma }}_{\min }[z,\overline{z}]$ is the area of the minimal surface encircled by the paths $z_{\mu }$ and ${\overline{z}}_{\mu }$. Furthermore, it turns out that the parametrization of ${\mathsf{\Sigma }}_{\min }$ through Equation (78) leads to the restoration of the translation invariance of $\mathsf{\Pi }(x,y)$.

To calculate $\mathsf{\Pi }(x,y)$, we introduce the “center-of-mass” and the relative coordinate of the gluons (cf. Ref. [99]):

$$u_{\mu }=\frac{\overline{s}{z}_{\mu }+s{\overline{z}}_{\mu }}{s+\overline{s}},r_{\mu }=z_{\mu }-{\overline{z}}_{\mu }.$$

We define further the new integration variables, which have the dimensionality of mass:

$$\mu =\frac{|x-y|}{2s},\overline{\mu }=\frac{|x-y|}{2\overline{s}}.$$

(88)

Then the kinetic terms of the gluons take the form

$${\int }_0^s\frac{{\dot{z}}_{\mu }^2}{4}d\lambda +{\int }_0^{\overline{s}}\frac{{\dot{\overline{z}}}_{\mu }^2}{4}d\overline{\lambda }=\frac{1}{2}{\int }_0^{R}d\tau \left[(\mu +\overline{\mu }){\dot{u}}_{\mu }^2+{\mu }_r{\dot{r}}_{\mu }^2\right].$$

Here $R\equiv |x-y|$, the variable τ denotes the distance to the point x along the line which passes through x and y, and ${\mu }_r\equiv \frac{\mu \overline{\mu }}{\mu +\overline{\mu }}$ is the “reduced mass”. In order to calculate the world-line integral entering $\mathsf{\Pi }(x,y)$, we apply to Equation (78) the Cauchy-Schwarz inequality in the form

$${\mathsf{\Sigma }}_{\min }={\int }_0^{R}d\tau \left|\mathbf{r}\left(\tau \right)\right|\le {\left(R{\int }_0^{R}d\tau {\mathbf{r}}^2\right)}^{1/2}.$$

(89)

Getting rid of the power 1/2 in ${\mathrm{e}}^{-\sigma {\mathsf{\Sigma }}_{\min }}$ through the integration over an auxiliary parameter λ, we obtain

$$\mathsf{\Pi }(x,y)\simeq \frac{R^2}{4}{\int }_0^{\infty }\frac{d\mu }{{\mu }^2}{\int }_0^{\infty }\frac{d\overline{\mu }}{{\overline{\mu }}^2}\int {\left(\mathcal{D}{u}_{\mu }\right)}_{xy}{\left(\mathcal{D}{r}_{\mu }\right)}_{00}\times$$

$$\times {\int }_0^{\infty }\frac{d\lambda }{\sqrt{\pi \lambda }}exp\left[-\lambda -\frac{\mu +\overline{\mu }}{2}{\int }_0^{R}d\tau {\dot{u}}_{\mu }^2-\frac{{\mu }_r}{2}{\int }_0^{R}d\tau {\dot{r}}_{\mu }^2-\frac{{\sigma }^{2}R}{4\lambda }{\int }_0^{R}d\tau {\mathbf{r}}^2\right].$$

(90)

The integrals over $u_{\mu }\left(\tau \right)$ and $r_4\left(\tau \right)$ in Equation (90) can be recognized as those of a free particle, while the integral over $\mathbf{r}\left(\tau \right)$ is the one of the harmonic oscillator. Using the known expressions for these quantum-mechanical integrals, and introducing the variables

$$\xi \equiv \sigma R^{3/2}/\sqrt{2{\mu }_r\lambda }\mathrm{instead}\mathrm{of}\lambda ,a=\mu R/2\mathrm{instead}\mathrm{of}\mu ,b=\overline{\mu }R/2\mathrm{instead}\mathrm{of}\overline{\mu },$$

we can write down the polarization operator in the form

$$\mathsf{\Pi }(x,y)\equiv \mathsf{\Pi }\left(R\right)=\frac{\sigma }{16{\pi }^{9/2}{R}^2}f\left(\sigma R^2\right),$$

where

$$f\left(\sigma R^2\right)={\int }_0^{\infty }\frac{d\xi }{\sqrt{\xi }{sinh}^{3/2}\xi }{\int }_0^{\infty }dadb\sqrt{\frac{a+b}{ab}}exp\left[-a-b-{\left(\frac{\sigma R^2}{2\xi }\right)}^2\frac{a+b}{ab}\right].$$

(91)

When integrating analytically over a and b, it has been found (cf. Ref. [96]) that, owing to the $a\leftrightarrow b$ symmetry, the corresponding saddle-point equations can be solved even if the pre-exponential is lifted to the exponential. However, the leading large-distance asymptotic behavior of $\mathsf{\Pi }\left(R\right)$ stems from the mere substitution of the saddle-point values appearing without this lifting, $a=b\simeq \frac{\sigma R^2}{2\xi }$, into the pre-exponential. Such a substitution yields

$$\mathsf{\Pi }\left(R\right)\simeq \frac{{\sigma }^{3/2}}{16{\pi }^{7/2}R}{\int }_0^{\infty }\frac{d\xi }{\xi {sinh}^{3/2}\xi }{\mathrm{e}}^{-2\sigma R^2/\xi }.$$

(92)

The latter integral can be evaluated by splitting the integration region into two intervals:

$$\mathsf{\Pi }\left(R\right)\equiv \frac{{\sigma }^{3/2}}{16{\pi }^{7/2}R}(I_1+I_2),$$

(93)

where

$$I_1\simeq {\int }_0^1\frac{d\xi }{{\xi }^{5/2}}{\mathrm{e}}^{-2\sigma R^2/\xi },I_2\simeq 2^{3/2}{\int }_1^{\infty }\frac{d\xi }{\xi }{\mathrm{e}}^{-\frac{3\xi }{2}-\frac{2\sigma R^2}{\xi }}.$$

Then, at $\sigma R^2\gg 1$, one obtains that $I_1=\mathcal{O}\left({\mathrm{e}}^{-2\sigma R^2}\right)$ is a subleading term. Rather, the integral $I_2$ is saturated by its saddle point, ${\xi }_{\mathrm{s}.\mathrm{p}.}=2R\sqrt{\sigma /3}$, which lies inside the integration region. This yields

$$\mathsf{\Pi }\left(R\right)\simeq \frac{{\sigma }^{5/4}}{2^{5/2}·3^{1/4}·{\pi }^3·R^{3/2}}{\mathrm{e}}^{-2\sqrt{3\sigma }R}\mathrm{at}\sigma R^2\gg 1.$$

(94)

Let us now consider the opposite limit of small distances. In this limit, one obtains

$$\int {\left(\mathcal{D}{u}_{\mu }\right)}_{xy}{\left(\mathcal{D}{r}_{\mu }\right)}_{00}exp\left[-\frac{\mu +\overline{\mu }}{2}{\int }_0^{R}d\tau {\dot{u}}_{\mu }^2-\frac{{\mu }_r}{2}{\int }_0^{R}d\tau {\dot{r}}_{\mu }^2-\frac{{\sigma }^{2}R}{4\lambda }{\int }_0^{R}d\tau {\mathbf{r}}^2\right]\to$$

$$\to {\left[\frac{\mu \overline{\mu }}{{\left(2\pi R\right)}^2}\right]}^{2}exp\left[-\frac{(\mu +\overline{\mu })R}{2}\right]\mathrm{at}\sigma R^2\ll 1.$$

Since this expression does not depend on λ anymore, the remaining λ-integration in Equation (90) results in a factor of 1. Equation (90) then goes over to the polarization operator of a free scalar particle, Equation (86):

$$\mathsf{\Pi }(x,y)\to \frac{1}{64{\pi }^4{R}^4}{\int }_0^{\infty }d\mu {\int }_0^{\infty }d\overline{\mu }exp\left[-\frac{(\mu +\overline{\mu })R}{2}\right]={\mathsf{\Pi }}_{\mathrm{free}}(x-y)\mathrm{at}\sigma R^2\ll 1.$$

(95)

Thus, we obtain a formula for the polarization operator, which interpolates between the limits (94) and (95):

$$\mathsf{\Pi }\left(R\right)=\frac{1}{16{\pi }^4{R}^4}{\mathrm{e}}^{-A\sqrt{\sigma R^2}}\left(1+B{\left(\sigma R^2\right)}^{5/4}\right),$$

(96)

where the analytic values of the coefficients read $A=2\sqrt{3}\simeq 3.46$, $B=\frac{2^{3/2}\pi }{3^{1/4}}\simeq 6.75$.

The integral in Equation (91) can also be calculated numerically, which yields the fitted value of $A=3.38$. As we see, it is very close to the above-obtained analytic value of 3.46. The corresponding analytic and numerical values of the mass parameter m in Equation (79) read

$$m_{\mathrm{an}}=2\sqrt{3\sigma }=1.64\mathrm{GeV},m_{\mathrm{num}}=3.38\sqrt{\sigma }=1.60\mathrm{GeV}.$$

(97)

Note that the analytically calculated value depends on the dimensionality D of space-time as

$$m_{\mathrm{an}}=2\sqrt{(D-1)\sigma }.$$

(98)

This is readily seen from the saddle point of the integral $I_2$ in Equation (93) by noticing that, in D dimensions, ${sinh}^{3/2}\xi \to {sinh}^{(D-1)/2}\xi$ in Equation (92). This result is a direct consequence of the minimal-area ansatz (89) that we have used. (However, in case of $D\ne 4$ dimensions, there is no asymptotic freedom, and in a space with $D>4$ dimensions, point particles are unlikely to be confined.)

The infra-red behavior of the running coupling, defined by Equation (79), stems from the replacement of the free-gluon polarization operator, Equation (84), by that of a confined gluon,

$${\mathsf{\Pi }}_{\mathrm{conf}}\left(p^2\right)\equiv \frac{1}{16{\pi }^2}ln\frac{{\mathsf{\Lambda }}_0^2}{p^2+m^2}.$$

(99)

Therefore, it looks instructive to compare the inverse Fourier image of this desired exact expression with our approximate result, Equation (96). We have

$$\int d^{4}p{\mathrm{e}}^{-ipx}ln\frac{{\mathsf{\Lambda }}_0^2}{p^2+m^2}={\int }_0^{\infty }\frac{ds}{s}\int d^{4}p{\mathrm{e}}^{-ipx}\left[{\mathrm{e}}^{-s(p^2+m^2)}-{\mathrm{e}}^{-s{\mathsf{\Lambda }}_0^2}\right]=$$

$$={\pi }^2{\int }_0^{\infty }\frac{ds}{s^3}{\mathrm{e}}^{-\frac{x^2}{4s}-m^{2}s}=\frac{8{\pi }^2{m}^2}{x^2}{K}_2\left(m\right|x\left|\right),$$

where $K_2$ is a Macdonald function. (When deriving the third formula in this chain, we have used the obvious fact that $\int d^{4}p{\mathrm{e}}^{-ipx}=0$ for $x\ne 0$.) Therefore, Equation (99) in the coordinate representation reads

$${\mathsf{\Pi }}_{\mathrm{conf}}\left(R\right)=\frac{m^2}{32{\pi }^4{R}^2}{K}_2\left(mR\right).$$

(100)

The short-distance asymptotic limit of this formula coincides with Equation (95). As for the large-distance limit, we see that Equation (100) has the same exponential fall-off as the obtained result (96). The only difference between the large-distance behavior of Equations (100) and (96) stems from the pre-exponential factor, which plays however only a secondary role.

It is also instructive to compare Equation (99) with the Fourier image of Equation (96). Because of the logarithmic divergency, we compare the derivatives of the two expressions, which are UV-finite. That is, we compare

$$\frac{d{\mathsf{\Pi }}_{\mathrm{conf}}\left(p^2\right)}{d{p}^2}=-\frac{1}{16{\pi }^2}\frac{1}{p^2+m^2}$$

with the derivative $\frac{d\mathsf{\Pi }\left(p^2\right)}{d{p}^2}$, where

$$\mathsf{\Pi }\left(p^2\right)=\frac{4{\pi }^2}{\left|p\right|}{\int }_0^{\infty }dR{R}^2{J}_1\left(\right|p\left|R\right)\mathsf{\Pi }\left(R\right)$$

is the Fourier image of $\mathsf{\Pi }\left(R\right)$, and $J_{\nu }$’s denote the Bessel functions. Introducing the dimensionless variables $x=\sqrt{\sigma }R$ and $q=\left|p\right|/\sqrt{\sigma }$, one has

$$\frac{d\mathsf{\Pi }\left(p^2\right)}{d{p}^2}\equiv \frac{F\left(q\right)}{\sigma },\mathrm{where}$$

$$F\left(q\right)=\frac{1}{8{\pi }^2{q}^2}{\int }_0^{\infty }\frac{dx}{x}\left[\frac{1}{2}\left(J_0\left(qx\right)-J_2\left(qx\right)\right)-\frac{J_1\left(qx\right)}{qx}\right]\left(1+B{x}^{5/2}\right){\mathrm{e}}^{-Ax}.$$

(101)

In Ref. [96], this function has been calculated numerically in the momentum range $0.02\mathrm{GeV}\le \left|p\right|\le 5\mathrm{GeV}$, which corresponds to $q\in [0.04,10.54]$. With a high accuracy, the result can be fitted by the function

$$F_{\mathrm{fit}}\left(q\right)=-\frac{1}{16{\pi }^2}\frac{d_1}{q^2+d_2^2},$$

(102)

which yields the values $d_1=0.585$ and $d_2=2.491$. Therefore, the value of m corresponding to the coefficient $d_2$, reads

$$m=d_2\sqrt{\sigma }=1.18\mathrm{GeV}.$$

(103)

It is closer to the phenomenological estimates ($\simeq 1\mathrm{GeV}$) [105,106] than the values (97). Furthermore, the coefficient $d_1$ yields a numerical prediction for the parameter $\tilde{b}$ from Equation (79):

$$\tilde{b}\simeq 0.585b=6.435.$$

(104)

Thus, the numerical analysis in the momentum representation yields an effective decrease of b in the infra-red region.

Let us now extend our analysis to the deconfinement phase, i.e., to the temperatures $T>T_c$. As has already been discussed, large spatial Wilson loops in this phase still exhibit the area law, which is the essence of the so-called magnetic- or spatial-confinement phenomenon. Accordingly, if the points x and y are separated by a time-like interval, the two $a_{\mu }^a$-gluons that propagate from one of these points to the other, are not confined. Consequently, ${\alpha }_s\left(p^2\right)$ at $p^2>0$ is given by the perturbative formula without freezing. Rather, if the points x and y are separated by a space-like interval, magnetic confinement of $a_{\mu }^a$-gluons takes place, so that one can expect to find the freezing of ${\alpha }_s\left(p^2\right)$ at $p^2<0$. We will start our analysis with the calculation of the polarization operator of a gluon in the SU(3) Yang-Mills theory at temperatures higher than the temperature of dimensional reduction, $T>T_{*}\simeq 2T_c$. As has been discussed, the Yang-Mills theory becomes a superrenormalizable theory in three spatial dimensions, where the renormalization of the dimensionful coupling $g_T\sqrt{T}$ is exact in one loop. Clearly, this coupling is not asymptotically free, being therefore unrelated to the study of freezing. Accordingly, the purpose of this calculation is to prepare the subsequent analysis of ${\alpha }_s$ in the temperature range $T_c<T<T_{*}$.

Let us first consider the propagation of a free particle at temperature T from the origin to the point $R_{\mu }=(\mathbf{R},R_4)$. The corresponding propagator reads

$${(-{\partial }^2)}_{R,0}^{-1}={\int }_0^{\infty }ds\sum _n^{}\frac{1}{{\left(4\pi s\right)}^2}exp\left[-\frac{{\mathbf{R}}^2+{(R_4+\beta n)}^2}{4s}\right]=$$

$$={\int }_0^{\infty }ds\sum _n^{}\frac{1}{2T\sqrt{\pi s}}exp\left[-\frac{{(R_4+\beta n)}^2}{4s}\right]·T\int {\left(\mathcal{D}\mathbf{z}\right)}_{\mathbf{R}\mathbf{0}}exp\left(-{\int }_0^s\frac{{\dot{\mathbf{z}}}^2}{4}d\lambda \right),$$

(105)

where ${\displaystyle \sum _n}\equiv {\displaystyle \sum _{n=-\infty }^{+\infty }}$. Performing the Poisson resummation, one has

$$\frac{1}{2T\sqrt{\pi s}}\sum _n^{}exp\left[-\frac{{(R_4+\beta n)}^2}{4s}\right]=\sum _n^{}exp\left(-{\omega }_n^{2}s+i{\omega }_n{R}_4\right),\mathrm{where}{\omega }_n=2\pi nT.$$

(106)

In the limit of $T\to \infty$, only the zeroth term on the right-hand side of Equation (106) survives, which means the dimensional reduction. The sum goes to 1, and

$${(-{\partial }^2)}_{R,0}^{-1}\to T{\int }_0^{\infty }ds\int {\left(\mathcal{D}\mathbf{z}\right)}_{\mathbf{R}\mathbf{0}}exp\left(-{\int }_0^s\frac{{\dot{\mathbf{z}}}^2}{4}d\lambda \right)=\frac{T}{4\pi L},$$

(107)

where $L\equiv \left|\mathbf{R}\right|$. With the effect of magnetic confinement included, the polarization operator at $T>T_{*}$ reads [cf. Equations (87) and (89)]:

$$\mathsf{\Pi }(\mathbf{x},\mathbf{y}|T)=T^2{\int }_0^{\infty }ds{\int }_0^{\infty }d\overline{s}I(s,\overline{s}),\mathrm{where}$$

$$I(s,\overline{s})\simeq \int {\left(\mathcal{D}\mathbf{z}\right)}_{\mathbf{x}\mathbf{y}}{\left(\mathcal{D}\overline{\mathbf{z}}\right)}_{\mathbf{y}\mathbf{x}}exp\left[-{\int }_0^s\frac{{\dot{\mathbf{z}}}^2}{4}d\lambda -{\int }_0^{\overline{s}}\frac{{\dot{\overline{\mathbf{z}}}}^2}{4}d\overline{\lambda }-{\sigma }_s{\left(L{\int }_0^{L}d\tau {\mathbf{w}}^2\right)}^{1/2}\right].$$

(108)

Here ${\sigma }_s$ is the spatial string tension, $\mathbf{R}=\mathbf{y}-\mathbf{x}$, and $\mathbf{w}$ is a two-dimensional vector orthogonal to $\mathbf{R}$. This polarization operator can be calculated in a way which is similar to the calculation of the polarization operator at zero temperature. Referring the reader for details to Ref. [96], we present here the final result for the polarization operator at $T>T_{*}$:

$$\mathsf{\Pi }(\mathbf{x},\mathbf{y}|T)\simeq \frac{\sqrt{{\sigma }_s}{T}^2}{4{\pi }^{2}L}{\mathrm{e}}^{-mL}\mathrm{at}{\sigma }_s{L}^2\gg 1,$$

(109)

$$\mathsf{\Pi }(\mathbf{x},\mathbf{y}|T)\simeq {\left(\frac{T}{4\pi L}\right)}^2\mathrm{at}{\sigma }_s{L}^2\ll 1.$$

(110)

The mass parameter m entering Equation (109) reads $m=2\sqrt{2{\sigma }_s}$. It reproduces Equation (98) at $D=3$, as expected. Furthermore, Equation (110) is nothing but the polarization operator of a free scalar particle, which is given by the square of the propagator (107).

As has been outlined above, we proceed now to the physically more interesting range of temperatures $T_c<T<T_{*}$, and calculate there the polarization operator $\mathsf{\Pi }(x,y|T)$, where x and y are four-vectors. Similarly to the case of $T=0$, the exponential fall-off of $\mathsf{\Pi }(x,y|T)$ at such temperatures, which stems from magnetic confinement, is relevant to the freezing of ${\alpha }_s$. The world-line representation of $\mathsf{\Pi }(x,y|T)$ at $T_c<T<T_{*}$ can be constructed by using Equations (105) and (108) as

$$\mathsf{\Pi }(x,y|T)=\frac{1}{4\pi }{\int }_0^{\infty }\frac{ds}{\sqrt{s}}{\int }_0^{\infty }\frac{d\overline{s}}{\sqrt{\overline{s}}}I(s,\overline{s})\sum _{n,k}^{}exp\left\{-\frac{1}{4}\left[\frac{{(R_4+\beta n)}^2}{s}+\frac{{(R_4+\beta k)}^2}{\overline{s}}\right]\right\}.$$

(111)

At large distances, characterized by the inequality ${\sigma }_s{L}^2\gg 1$, the result for this quantity reads (cf. Ref. [96]):

$$\mathsf{\Pi }(x,y|T)\simeq \frac{{\sigma }_{s}L}{8\sqrt{2}{\pi }^3{R}^3}{\mathrm{e}}^{-2\sqrt{2{\sigma }_{s}LR}},$$

(112)

where $R\equiv \sqrt{L^2+R_4^2}$. Recalling that magnetic confinement holds only if $R_{\mu }$ is a space-like vector, one can always place the points x and y along some spatial axis, which makes L and R equal. The resulting formula for the polarization operator at large distances takes the form

$$\mathsf{\Pi }(x,y|T)\simeq \frac{{\sigma }_s}{8\sqrt{2}{\pi }^3{R}^2}{\mathrm{e}}^{-2\sqrt{2{\sigma }_s}R}\mathrm{at}{\sigma }_s{R}^2\gg 1.$$

(113)

Thus, at temperatures $T_c<T<T_{*}$, freezing of ${\alpha }_s\left(p^2\right)$ takes place at the temperature-dependent momentum scale, which is analytically defined as $m=2\sqrt{2{\sigma }_s}$. The factor of 2 under the square root in this formula is the number of spatial dimensions minus one, in accordance with Equation (98).

In the opposite limit of small distances, the polarization operator at $L=R$ reads (In particular, at $T\to \infty$, this result goes over to Equation (110), as expected.)

$$\mathsf{\Pi }(x,y|T)\to {\left[\frac{T}{4\pi R}coth\left(\pi TR\right)\right]}^2\mathrm{at}{\sigma }_s{R}^2\ll 1.$$

(114)

Using the value of $T_c=270\mathrm{MeV}$, one can prove that, at $T_c<T<T_{*}$, the condition ${\sigma }_s{R}^2\ll 1$ implies also the inequality $TR\ll 1$. For this reason, at such temperatures, Equation (114) can be approximated by its zero-temperature counterpart, $\frac{1}{{\left(4{\pi }^2{R}^2\right)}^2}$. Therefore, the formula for the polarization operator, which interpolates between this short-distance limit and the large-distance limit (113), reads

$$\mathsf{\Pi }(x,y|T)=\frac{{\mathrm{e}}^{-2\sqrt{2{\sigma }_s}R}}{{\left(4{\pi }^2{R}^2\right)}^2}\left(1+\pi \sqrt{2}{\sigma }_s{R}^2\right)\mathrm{at}{T}_c<T<T_{*}.$$

Note that this expression depends on temperature only implicitly, namely through ${\sigma }_s\left(T\right)$. Furthermore, the analysis of the polarization operator $\mathsf{\Pi }(x,y|T)$ in the momentum representation can be performed similarly to the zero-temperature case. Specifically, the factor $\left(1+B{x}^{5/2}\right){\mathrm{e}}^{-Ax}$ in Equation (101) should be replaced by $\left(1+\pi \sqrt{2}{x}^2\right){\mathrm{e}}^{-2\sqrt{2}x}$. Using further the fitting function (102), we get the values $d_1=0.711$ and $d_2=2.289$. The mass parameter, which defines the onset of freezing, is given by the formula $m=d_2\sqrt{{\sigma }_s}$. For T varying from $T_c$ to $T_{*}\simeq 2T_c$, m varies from 1.09 GeV to 1.56 GeV. The parameter $\tilde{b}\simeq 0.711b=7.821$ turns out to be larger than its zero-temperature counterpart (104).

Both at $T=0$ and at $T_c<T<T_{*}$, freezing modifies Equation (85) as

$$\frac{1}{g^2\left(p\right)}=\frac{1}{g_0^2}-\frac{\tilde{b}}{16{\pi }^2}ln\frac{{\mathsf{\Lambda }}_0^2+m^2}{p^2+m^2}.$$

Defining the renormalized cut-off Λ through the bare one ${\mathsf{\Lambda }}_0$ as

$$\mathsf{\Lambda }=\sqrt{{\mathsf{\Lambda }}_0^2+m^2}exp\left(-\frac{2\pi }{{\alpha }_s\left({\mathsf{\Lambda }}_0\right)\tilde{b}}\right),$$

one ensures that ${\alpha }_s\left(p\right)=\frac{4\pi }{\tilde{b}ln\frac{p^2+m^2}{{\mathsf{\Lambda }}^2}}$ takes the value ${\alpha }_s\left({\mathsf{\Lambda }}_0\right)$ at $p={\mathsf{\Lambda }}_0$. It is then natural to choose a sufficiently large momentum scale ${\mathsf{\Lambda }}_0$, so that ${\alpha }_s\left({\mathsf{\Lambda }}_0\right)$ is essentially unaffected by freezing and by the finite-temperature effects, and to match ${\alpha }_s\left({\mathsf{\Lambda }}_0\right)$ to the experimental value. Choosing ${\mathsf{\Lambda }}_0=3\mathrm{GeV}$, where [107] ${\alpha }_s\simeq 0.2524$, one obtains for ${\alpha }_s\left(p\right)$ with freezing the following values of ${\alpha }_s\left(0\right)$ at $T=0$, $T_c$, $1.5T_c$, and $2T_c$: 0.341, 0.377, 0.357, and 0.333, respectively. For a given p, one observes an increase of ${\alpha }_s\left(p\right)$ at $T=T_c$ with respect to ${\alpha }_s\left(p\right)$ at $T=0$, and a subsequent decrease of ${\alpha }_s\left(p\right)$ with the increase of T. Furthermore, the obtained values of ${\alpha }_s\left(0\right)$ are very close to those that follow from other approaches where the infra-red finiteness of ${\alpha }_s\left(p\right)$ has been studied. For instance, Ref. [108] yields ${\alpha }_s\left(0\right)\simeq 0.35$, Ref. [109] predicts the values of ${\alpha }_s\left(0\right)$ in the range $0.3<{\alpha }_s\left(0\right)<0.35$, Ref. [110] suggests the value ${\alpha }_s\left(0\right)=0.4$.

As has been pointed out in Refs. [111,112], the infra-red behavior of ${\alpha }_s\left(p\right)$ can be used to estimate non-perturbative corrections to various QCD observables. The parameters essential for this estimate are the first few integral moments of ${\alpha }_s\left(p\right)$ in the infra-red region. One class of such observables consists of the so-called event-shape variables, such as the thrust $\mathcal{T}$:

$$\mathcal{T}=\underset{\mathbf{n}}{max}\frac{{\sum }_i\left|{\mathbf{p}}_i\mathbf{n}\right|}{{\sum }_i\left|{\mathbf{p}}_i\right|}.$$

In this formula, the vectors ${\mathbf{p}}_i$’s are the momenta of final-state hadrons, and $\mathbf{n}$ is a unit vector which maximizes $\mathcal{T}$. If the momenta of hadrons form an almost collinear jet, then the maximum is achieved for the vector $\mathbf{n}$ which is directed along the jet axis. The value of $\mathcal{T}$ becomes equal to 1 for an ideal pencil-like jet. However, due to the radiation of gluons, the observed $\mathcal{T}$ is different from 1. In perturbative QCD, such corrections stem from the hard-gluon radiation [113]. Specifically, the observable $1-\mathcal{T}$ has the following expansion in powers of the perturbative strong coupling:

$$1-\mathcal{T}{|}_{\mathrm{pert}}=0.334{\alpha }_s^{\mathrm{pert}}\left(p\right)+1.02{\left({\alpha }_s^{\mathrm{pert}}\left(p\right)\right)}^2+\mathcal{O}\left({\left({\alpha }_s^{\mathrm{pert}}\left(p\right)\right)}^3\right).$$

However, in addition to the gluon radiation in the perturbative region, the soft-gluon radiation below the scale ${\mu }_{\mathrm{IR}}\simeq 3\mathrm{GeV}$ also takes place, resulting in the so-called hadronization correction:

$$1-\mathcal{T}=1-\mathcal{T}{|}_{\mathrm{pert}}+\frac{2\lambda }{p}.$$

(115)

The non-perturbative higher-twist contribution λ can be related to the lowest integral moment of ${\alpha }_s\left(p\right)$ with freezing as [112]

$$\lambda =\frac{C_{\mathrm{f}}}{\pi }{\int }_0^{3\mathrm{GeV}}dp{\alpha }_s\left(p\right),$$

where $C_{\mathrm{f}}=4/3$ is the quadratic Casimir operator of the fundamental representation of the group SU(3). Using for ${\alpha }_s\left(p\right)=\frac{4\pi }{\tilde{b}ln\frac{p^2+m^2}{{\mathsf{\Lambda }}^2}}$ the above-obtained parameters m and $\tilde{b}$ at $T=0$ and $T=T_c$, We have

$$\lambda {|}_{T=0}=0.376\mathrm{GeV},\lambda {|}_{T=0.27\mathrm{GeV}}=0.395\mathrm{GeV}.$$

(116)

As a possible phenomenological application of these results, let us evaluate the thrust of a two-jet event at the scale $p=M_Z$. Using the value ${\alpha }_s^{\mathrm{pert}}\left(M_Z\right)=0.12$ (cf. Ref. [107]), one obtains first the purely perturbative contribution:

$$1-\mathcal{T}{|}_{\mathrm{pert}}=0.334{\alpha }_s^{\mathrm{pert}}\left(M_Z\right)+1.02{\left({\alpha }_s^{\mathrm{pert}}\left(M_Z\right)\right)}^2\simeq 0.055,$$

(117)

which underestimates the experimental value measured at LEP [114], $1-\mathcal{T}\simeq 0.068$. Rather, inclusion of the higher-twist contribution according to Equation (115), with the parameter λ given by Equation (116), yields for the full quantity

$$1-\mathcal{T}{|}_{T=0}=0.063,$$

(118)

which is closer to the aforementioned experimental value. Furthermore, for the case of hadronization taking place in the quark-gluon plasma, the value of $\lambda {|}_{T=0.27\mathrm{GeV}}$ from Equation (116) yields at the $M_Z$-scale

$$1-\mathcal{T}{|}_{T=0.27\mathrm{GeV}}=0.064.$$

(119)

Thus, although ${\alpha }_s\left(p\right)$ with freezing enhances the perturbative contribution (117), the corresponding full values (118) and (119) are nevertheless still somewhat smaller than the experimental one.

In conclusion of this Section, we have analyzed the infra-red freezing of the running strong coupling in the confinement and the deconfinement phases. The direct calculation of the polarization operator of a gluon confined by the stochastic background fields has been done by using the parametrization (89) for the minimal area. This parametrization reduces the corresponding world-line integral to that of the three-dimensional harmonic oscillator. In the deconfinement phase, at temperatures $T_c<T<T_{*}\simeq 2T_c$, freezing of ${\alpha }_s\left(p\right)$ also takes place due to the so-called magnetic confinement. Upon the calculation of the world-line integral, we have found the momentum scales at which the freezing phenomenon occurs in the hadronic and the deconfinement phases. Analytically, these scales are given by one and the same formula (98), where $D=4$ at $T=0$ and $D=3$ at $T>T_c$. Numerically, the values of the freezing scales following from the fits in the momentum representation are smaller than the corresponding analytic ones and closer to the phenomenological value of 1 GeV [cf. Equation (103)]. Finally, we have estimated the physical effect produced by freezing on the thrust variable. The corresponding non-perturbative contribution, which arises due to the soft-gluon radiation, brings the purely perturbative value of this quantity closer to the experimental one.

6. Correlation Lengths of the Stochastic Yang-Mills Fields

As has been discussed in Section 1, stochastic vacuum model [11,12,13,14,15,16] puts in practice the idea that it is the stochasticity of non-perturbative background Yang-Mills fields that provides confinement. This idea is implemented through the assumption that these stochastic fields form a Gaussian ensemble, where the two-point correlation function $\langle \mathrm{tr}{F}_{\mu \nu }^a\left(x\right)T^a{\mathsf{\Phi }}_{x{x}^{\prime }}{F}_{\lambda \rho }^b\left(x^{\prime }\right)T^b{\mathsf{\Phi }}_{x^{\prime }x}\rangle$ contains the tensor structure essential for confinement. The corresponding part of this correlation function is parametrized through the gluon condensate and the vacuum correlation length, which yields the Wilson-loop average in the form of Equation (15). The thus emerging confinement is quantified by the string tension, which for a given representation r of the group SU(N), reads ${\sigma }_r\propto C_r{\lambda }^2\langle {\left(g{F}_{\mu \nu }^a\right)}^2\rangle$. Here, $C_r$ is the quadratic Casimir operator of the representation r, and the proportionality coefficient slightly depends on a particular parametrization of the two-point correlation function.

In spite of the numerical support provided by the dedicated lattice simulations [41,42,43,44,45,71,115], the stochastic vacuum model needs also a theoretical support in the form of analytic studies of the two-point correlation function. Such calculations have been initially performed in various Abelian models, where confinement is provided by the monopole condensation. Those include the 3D Georgi-Glashow model and the 4D compact QED [116,117], as well as the 4D dual Abelian Higgs model in the London limit and its [U(1)]${}^{N-1}$-symmetric generalization [118,119,120,121,122]. In this way, it has been explicitly shown that the vacuum correlation lengths in all these models are equal to the inverse masses of the corresponding dual vector bosons. From this point of view, the only difference between the 3D Georgi-Glashow model (or the 4D compact QED) and the 4D dual Abelian-Higgs-type models in the London limit is a particular mechanism of generation of the dual-boson mass. Specifically, in the case of the monopole-antimonopole plasma of the 3D Georgi-Glashow model, this mechanism is provided by the Debye screening, while being just the Higgs mechanism in the case of the 4D dual Abelian-Higgs-type models. Moreover, as will be shown at the end of this Section, the latter models can be modified in such a way as to allow for the existence of two different vacuum correlation lengths, which can make them even closer to QCD.

A lattice simulation indicating that the nonperturbative-nonconfining and confining stochastic background Yang-Mills fields can have different correlation lengths has been done in Ref. [123]. Analytic studies of the vacuum correlation length in QCD have been initiated in Ref. [124], suggesting that a direct relation between this length and the gluon condensate can exist. Next, Refs. [74,75,76,77] analytically explored the possibility for the ensemble of stochastic background Yang-Mills fields to have two different correlation lengths. That has been done within the theory of so-called gluelumps. As has already been mentioned in Section 2, gluelumps can be thought of as bound states of gluons in the field of a hypothetical infinitely heavy adjoint source [72,73]. In the Yang-Mills theory, they define the correlation lengths of the two-point function of gluonic field strengths in the same way as, in full QCD, physically existing heavy-light mesons define the correlation length of a non-local quark condensate [125] $\langle \overline{\psi }\left(x\right){\mathsf{\Phi }}_{x{x}^{\prime }}\psi \left(x^{\prime }\right)\rangle$. (Here, ${\mathsf{\Phi }}_{x{x}^{\prime }}$ denotes a phase factor along the straight-line path, which is provided by the heavy-quark propagator.) Unlike the fundamental representation of the group SU(N), the adjoint representation allows for two different types of heavy-light objects: those with a single gluon, called one-gluon gluelumps, and those with two gluons, called two-gluon gluelumps. The first case is similar to the aforementioned non-local quark condensate, whereas the second case is conceptually different, as it corresponds to two gluons connected together with the heavy source by three fundamental strings. (For brevity, we call here strings that interconnect charges transforming under the fundamental representation of the group SU(N) simply as fundamental strings. This term should not be confused with its counterpart used in String Theory. In the same manner, we use the term “adjoint string” for strings that interconnect charges transforming under the adjoint representation of SU(N).) As has been shown in Refs. [74,75,76,77], the Green functions of one- and two-gluon gluelumps define respectively the correlation lengths of nonperturbative-nonconfining and confining self-interactions of the stochastic Yang-Mills fields. Here, we follow Ref. [77] to provide the full calculation of the world-line integrals, which represent those Green functions. Such a calculation turns out to be possible by virtue of an effective parametrization of minimal areas swept out by the strings in the gluelumps, analogous to Equation (89). Furthermore, in Appendix C, we calculate the Green function of the one-gluon gluelump for the case where the heavy adjoint source evolves along a trajectory which deviates from the straight-line one. This calculation allows us to show that the correlation length of nonperturbative-nonconfining self-interactions of the background fields decreases with the deformation of that trajectory, in agreement with the lattice results of Ref. [115].

The correlation lengths of nonperturbative-nonconfining and confining self-interactions of the background fields describe the fall-off of, respectively, the functions $D_1\left(x\right)$ and $D\left(x\right)$, which appear in the following generalization of the parametrization (23) (see Refs. [11,12,13,14,15,16]):

$$\frac{g^2}{N}\langle \mathrm{tr}{F}_{\mu \nu }^a\left(x\right)T^a{\mathsf{\Phi }}_{x{x}^{\prime }}{F}_{\lambda \rho }^b\left(x^{\prime }\right)T^b{\mathsf{\Phi }}_{x^{\prime }x}\rangle =({\delta }_{\mu \lambda }{\delta }_{\nu \rho }-{\delta }_{\mu \rho }{\delta }_{\nu \lambda })D(x-x^{\prime })+$$

$$+\frac{1}{2}\left[{\partial }_{\mu }^x\left({(x-x^{\prime })}_{\lambda }{\delta }_{\nu \rho }-{(x-x^{\prime })}_{\rho }{\delta }_{\nu \lambda }\right)+{\partial }_{\nu }^x\left({(x-x^{\prime })}_{\rho }{\delta }_{\mu \lambda }-{(x-x^{\prime })}_{\lambda }{\delta }_{\mu \rho }\right)\right]D_1(x-x^{\prime }).$$

(120)

We calculate the functions $D_1\left(x\right)$ and $D\left(x\right)$ analytically, by using their relations to the Green functions of the one- and the two-gluon gluelumps. These relations read [74,75,76,77]

$$D_1\left(x\right)=-4g^2{C}_{\mathrm{f}}\frac{dG\left(x\right)}{d{x}^2},D\left(x\right)=\frac{g^4(N-1)}{2}S\left(x\right),$$

(121)

where $C_{\mathrm{f}}=(N^2-1)/\left(2N\right)$ is the quadratic Casimir operator of the fundamental representation. Equations (121) hold at distances $\left|x\right|\ge \mathcal{O}\left({\sigma }_{\mathrm{f}}^{-1/2}\right)$, where ${\sigma }_{\mathrm{f}}$ is the string tension in the fundamental representation. At such distances, perturbative contributions are negligible, and the Wilson-loop averages appearing in the Green functions $G\left(x\right)$ and $S\left(x\right)$ exhibit the area law. (In the case of a one-gluon gluelump, we omit the perimeter-law contribution associated with the breaking of a sufficiently long adjoint string through the creation of a glueball. That is, the string in the one-gluon gluelump is supposed to be short enough to make the subleading processes of this kind impossible.) We start our analysis with the Green function of a one-gluon gluelump,

$$G\left(x\right)={\int }_0^{\infty }ds\int {\left(\mathcal{D}{z}_{\mu }\right)}_{x0}exp\left(-{\int }_0^{s}d\lambda \frac{{\dot{z}}_{\mu }^2}{4}-\sigma {\mathsf{\Sigma }}_{\min }\right),$$

(122)

where the minimal surface of area ${\mathsf{\Sigma }}_{\min }$ is swept out by the adjoint string of tension σ. This string connects the gluon to a heavy source which transforms under the adjoint representation of the group SU(N). While the effects of non-staticity of the source are considered in Appendix C, here we assume the source to be static, i.e., evolving entirely along the $x_4$-axis. This means that only the $x_4$-coordinate of the point x in Equation (122) is non-vanishing, i.e., $x=(\mathbf{0},L)$. Hence, we need to find $G\left(x\right)$ and $S\left(x\right)$ for such an x, i.e., to express these functions in terms of $L\equiv \left|x\right|$. We calculate the world-line integral from Equation (122) by approximating ${\mathsf{\Sigma }}_{\min }$ via the Cauchy-Schwarz inequality as [cf. Equation (89)]

$${\mathsf{\Sigma }}_{\min }={\int }_0^{L}d\tau \left|\mathbf{z}\left(\tau \right)\right|\le \sqrt{L{\int }_0^{L}d\tau {\mathbf{z}}^2}.$$

(123)

This approximation reduces the world-line integral to that of the harmonic oscillator of a variable frequency. Indeed, using the formula ${\mathrm{e}}^{-\sqrt{A}}={\int }_0^{\infty }\frac{d\lambda }{\sqrt{\pi \lambda }}exp\left(-\lambda -\frac{A}{4\lambda }\right)$ with $A>0$, and changing the proper-time variable as $s\to \mu =\frac{L}{2s}$, We have

$$G\left(x\right)\simeq \frac{L}{2}{\int }_0^{\infty }\frac{d\mu }{{\mu }^2}\int {\left(\mathcal{D}{z}_4\right)}_{L0}\oint \mathcal{D}\mathbf{z}{\int }_0^{\infty }\frac{d\lambda }{\sqrt{\pi \lambda }}exp\left(-\lambda -\frac{\mu }{2}{\int }_0^{L}d\tau {\dot{z}}_{\mu }^2-\frac{{\sigma }^{2}L}{4\lambda }{\int }_0^{L}d\tau {\mathbf{z}}^2\right),$$

(124)

where $\oint \mathcal{D}\mathbf{z}$ denotes the integration over trajectories which start and end up at $\mathbf{z}=\mathbf{0}$. The kinetic term $\frac{\mu {\dot{z}}_{\mu }^2}{2}$ means, of course, that the auxiliary parameter μ can be viewed as an effective gluon mass. Furthermore, the approximate equality “≃” in Equation (124) is understood in the sense of approaching the upper limit for ${\mathsf{\Sigma }}_{\min }$ in Equation (123).

As one can now see, the resulting world-line integral over $z_4\left(\tau \right)$ is that of a free particle, namely

$$\int {\left(\mathcal{D}{z}_4\right)}_{L0}exp\left(-\frac{\mu }{2}{\int }_0^{L}d\tau {\dot{z}}_4^2\right)=\sqrt{\frac{\mu }{2\pi L}}{\mathrm{e}}^{-\mu L/2}.$$

Rather, the harmonic-oscillator-type integral over $\mathbf{z}\left(\tau \right)$ reads

$$\oint \mathcal{D}\mathbf{z}exp\left(-\frac{\mu }{2}{\int }_0^{L}d\tau {\dot{\mathbf{z}}}^2-\frac{{\sigma }^{2}L}{4\lambda }{\int }_0^{L}d\tau {\mathbf{z}}^2\right)={\left[\frac{\omega }{4\pi sinh\left(\frac{\omega L}{2\mu }\right)}\right]}^{3/2},$$

(125)

where $\omega \equiv \sigma \sqrt{2\mu L/\lambda }$ is the frequency of the oscillator. Changing further the integration variable λ to $\xi \equiv \sigma L^{3/2}/\sqrt{2\mu \lambda }$, we can perform the μ-integration as

$$G\left(x\right)=\sigma \sqrt{\frac{L}{32{\pi }^5}}{\int }_0^{\infty }\frac{d\mu }{\sqrt{\mu }}{\int }_0^{\infty }\frac{d\xi }{\sqrt{\xi }{sinh}^{3/2}\xi }exp\left(-\frac{\mu L}{2}-\frac{{\sigma }^2{L}^3}{2\mu {\xi }^2}\right)=\frac{\sigma }{4{\pi }^2}{\int }_0^{\infty }\frac{d\xi }{\sqrt{\xi }{sinh}^{3/2}\xi }{\mathrm{e}}^{-l/\xi },$$

(126)

where $l\equiv \sigma L^2$. In the infra-red limit $l\gg 1$ of interest, the latter integral can be evaluated analytically as follows:

$${\int }_0^1\frac{d\xi }{{\xi }^2}{\mathrm{e}}^{-l/\xi }+2^{3/2}{\int }_1^{\infty }\frac{d\xi }{\sqrt{\xi }}{\mathrm{e}}^{-\frac{3\xi }{2}-\frac{l}{\xi }}\simeq \frac{{\mathrm{e}}^{-l}}{l}+2^{3/2}{\int }_0^{\infty }\frac{d\xi }{\sqrt{\xi }}{\mathrm{e}}^{-\frac{3\xi }{2}-\frac{l}{\xi }}=\frac{{\mathrm{e}}^{-l}}{l}+4\sqrt{\frac{\pi }{3}}{\mathrm{e}}^{-\sqrt{6\sigma }L}.$$

(127)

Here, the replacement of the lower limit of integration in the second integral by 0 was legitimate, since the saddle-point value ${\xi }_{*}=\sqrt{2l/3}$, in the limit $l\gg 1$ at issue, is larger than 1, so that the contribution of the integration region $0<\xi <1$ to the whole integral is exponentially suppressed. Thus, the leading result in the large-$\left|x\right|$ limit reads

$$G\left(x\right)\simeq \frac{\sigma }{\sqrt{3{\pi }^3}}{\mathrm{e}}^{-\sqrt{6\sigma }\left|x\right|}.$$

(128)

By means of Equation (121), it yields the following function $D_1\left(x\right)$:

$$D_1\left(x\right)=g^2{C}_{\mathrm{f}}{\left(\frac{2\sigma }{\pi }\right)}^{3/2}\frac{{\mathrm{e}}^{-\sqrt{6\sigma }\left|x\right|}}{\left|x\right|}.$$

(129)

This expression can be compared with the one from Ref. [76],

$${\tilde{D}}_1\left(x\right)=g^2{C}_{\mathrm{f}}\frac{M_0\sigma }{2\pi }·\frac{{\mathrm{e}}^{-M_0\left|x\right|}}{\left|x\right|},$$

(130)

where the value $M_0\simeq 1.5\mathrm{GeV}$ was obtained from the Schrödinger equation with the linear potential. To make such a comparison, we use Casimir scaling to evaluate the string tension in the adjoint representation of the group SU(3) as $\sigma =\frac{9}{4}{\sigma }_{\mathrm{f}}$, where the value of the string tension in the fundamental representation reads ${\sigma }_{\mathrm{f}}={\left(440\mathrm{MeV}\right)}^2$. This yields

$$\mathrm{mass}\mathrm{of}\mathrm{the}1g\mathrm{gluelump}=\sqrt{6\sigma }\simeq 1.6\mathrm{GeV},$$

(131)

which turns out to be very close to the above-quoted value of $M_0$. In particular, in the large-N limit, one has $\sigma \to 2{\sigma }_{\mathrm{f}}$, which leads to an extremely good coincidence between the exponentials in Equations (129) and (130), since in that case $\sqrt{6\sigma }\simeq 1.5\mathrm{GeV}$.

We proceed now to the Green function of the two-gluon gluelump:

$$S\left(x\right)=\frac{L^2}{4}{\int }_0^{\infty }\frac{d\mu }{{\mu }^2}{\int }_0^{\infty }\frac{d\overline{\mu }}{{\overline{\mu }}^2}\int {\left(\mathcal{D}{z}_{\mu }\right)}_{x0}\int {\left(\mathcal{D}{\overline{z}}_{\mu }\right)}_{x0}exp\{-\frac{\mu }{2}{\int }_0^{L}d\tau {\dot{z}}_{\mu }^2-\frac{\overline{\mu }}{2}{\int }_0^{L}d\tau {\dot{\overline{z}}}_{\mu }^2-$$

$$-{\sigma }_{\mathrm{f}}{\int }_0^{L}d\tau \left(\left|\mathbf{z}\right|+|\overline{\mathbf{z}}|+|\mathbf{z}-\overline{\mathbf{z}}|\right)\}.$$

The use of the Cauchy-Schwarz inequality (123) separately for each of the three distances in the last integral would involve integrations over three auxiliary parameters. In order to reduce this number to just one, it is more useful to apply first the Cauchy-Schwarz inequality in the form

$$\frac{1}{n}\sum _{i=1}^n{A}_i\le \sqrt{\frac{1}{n}\sum _{i=1}^n{A}_i^2},\mathrm{where}{A}_i>0.$$

For $n=3$, it yields

$$\left|\mathbf{z}\right|+|\overline{\mathbf{z}}|+|\mathbf{z}-\overline{\mathbf{z}}|\le \sqrt{3}·\sqrt{{\mathbf{z}}^2+{\overline{\mathbf{z}}}^2+{(\mathbf{z}-\overline{\mathbf{z}})}^2}.$$

Since the common square root for the three distances is now assembled, we can again apply the Cauchy-Schwarz inequality in the form of Equation (123) with only one auxiliary integration:

$$S\left(x\right)\simeq \frac{L^2}{4}{\int }_0^{\infty }\frac{d\mu }{{\mu }^2}{\int }_0^{\infty }\frac{d\overline{\mu }}{{\overline{\mu }}^2}\int {\left(\mathcal{D}{z}_{\mu }\right)}_{x0}\int {\left(\mathcal{D}{\overline{z}}_{\mu }\right)}_{x0}{\int }_0^{\infty }\frac{d\lambda }{\sqrt{\pi \lambda }}exp\{-\lambda -\frac{\mu }{2}{\int }_0^{L}d\tau {\dot{z}}_{\mu }^2-$$

$$-\frac{\overline{\mu }}{2}{\int }_0^{L}d\tau {\dot{\overline{z}}}_{\mu }^2-\frac{3{\sigma }_{\mathrm{f}}^{2}L}{4\lambda }{\int }_0^{L}d\tau \left[{\mathbf{z}}^2+{\overline{\mathbf{z}}}^2+{(\mathbf{z}-\overline{\mathbf{z}})}^2\right]\}.$$

(132)

Similarly to the case of the one-gluon gluelump, integrations over $z_4\left(\tau \right)$ and ${\overline{z}}_4\left(\tau \right)$ in this formula yield

$$\int {\left(\mathcal{D}{z}_4\right)}_{L0}\int {\left(\mathcal{D}{\overline{z}}_4\right)}_{L0}exp\left[-\frac{1}{2}{\int }_0^{L}d\tau (\mu {\dot{z}}_4^2+\overline{\mu }{\dot{\overline{z}}}_4^2)\right]=\frac{\sqrt{\mu \overline{\mu }}}{2\pi L}exp\left[-\frac{(\mu +\overline{\mu })L}{2}\right].$$

We observe now that, if the terms ${\mathbf{z}}^2$ and ${\overline{\mathbf{z}}}^2$ in Equation (132) were absent, the world-line integral would be that of two mutually interacting but otherwise free particles, which was calculated in the previous Section. Here, however, we have to deal with two (also mutually interacting) harmonic oscillators. Nevertheless, such a world-line integral $\oint \mathcal{D}\mathbf{z}\oint \mathcal{D}\overline{\mathbf{z}}$ can still be calculated. This turns out to be possible by virtue of the known fact that two positively definite quadratic forms (which are the kinetic and the potential energies) can be diagonalized simultaneously. Introducing new trajectories $\mathbf{u}\left(\tau \right)$ and $\mathbf{v}\left(\tau \right)$ according to the formulae

$$\mathbf{z}=\mathbf{u}+\alpha \mathbf{v},\overline{\mathbf{z}}=\beta \mathbf{u}+\mathbf{v},$$

(133)

we obtain the diagonalization conditions

$$\mu \alpha +\overline{\mu }\beta =0,\alpha +\beta =(\alpha -1)(\beta -1).$$

The solution to these equations is straightforward, and reads

$$\beta (\mu ,\overline{\mu })=-\frac{\mu }{\overline{\mu }}·\alpha (\mu ,\overline{\mu }),\mathrm{where}\alpha (\mu ,\overline{\mu })=1-\frac{\overline{\mu }}{\mu }+\sqrt{\frac{\overline{\mu }}{\mu }+{\left(1-\frac{\overline{\mu }}{\mu }\right)}^2}.$$

(134)

The kinetic- and the potential-energy terms take then the form

$$\mu {\dot{\mathbf{z}}}^2+\overline{\mu }{\dot{\overline{\mathbf{z}}}}^2=(\mu +\overline{\mu }{\beta }^2){\dot{\mathbf{u}}}^2+(\mu {\alpha }^2+\overline{\mu }){\dot{\mathbf{v}}}^2,$$

$${\mathbf{z}}^2+{\overline{\mathbf{z}}}^2+{(\mathbf{z}-\overline{\mathbf{z}})}^2=2\left[({\beta }^2-\beta +1){\mathbf{u}}^2+({\alpha }^2-\alpha +1){\mathbf{v}}^2\right].$$

The Jacobian stemming from the change of trajectories in Equation (133) is, of course, also μ- and $\overline{\mu }$-dependent, namely one has

$$\mathcal{D}\mathbf{z}\mathcal{D}\overline{\mathbf{z}}=|1-\alpha \beta |\mathcal{D}\mathbf{u}\mathcal{D}\mathbf{v}=\left(1+{\alpha }^2\frac{\mu }{\overline{\mu }}\right)\mathcal{D}\mathbf{u}\mathcal{D}\mathbf{v}.$$

Thus, we have reduced the world-line integral $\oint \mathcal{D}\mathbf{z}\oint \mathcal{D}\overline{\mathbf{z}}$ to the product of two integrals for non-interacting harmonic oscillators. Introducing finally the dimensionless variables

$$\nu \equiv \frac{\mu }{\sqrt{{\sigma }_{\mathrm{f}}}},\overline{\nu }\equiv \frac{\overline{\mu }}{\sqrt{{\sigma }_{\mathrm{f}}}},d\equiv \sqrt{{\sigma }_{\mathrm{f}}}L,$$

we can write down the result in the following form:

$$S\left(x\right)=\frac{3^{3/2}{\sigma }_{\mathrm{f}}^2}{64{\pi }^{9/2}}·d^{5/2}{\int }_0^{\infty }\frac{d\nu }{{\nu }^{3/2}}{\int }_0^{\infty }\frac{d\overline{\nu }}{{\overline{\nu }}^{3/2}}\left(1+{\alpha }^2\frac{\nu }{\overline{\nu }}\right){\int }_0^{\infty }\frac{d\lambda }{{\lambda }^2}exp\left[-\lambda -\frac{(\nu +\overline{\nu })d}{2}\right]\times$$

$$\times {\left[({\beta }^2-\beta +1)({\alpha }^2-\alpha +1)(\nu +\overline{\nu }{\beta }^2)(\overline{\nu }+\nu {\alpha }^2)\right]}^{3/4}\times$$

$$\times {sinh}^{-3/2}\left(d^{3/2}\sqrt{\frac{3}{\lambda }·\frac{{\beta }^2-\beta +1}{\nu +\overline{\nu }{\beta }^2}}\right)·{sinh}^{-3/2}\left(d^{3/2}\sqrt{\frac{3}{\lambda }·\frac{{\alpha }^2-\alpha +1}{\overline{\nu }+\nu {\alpha }^2}}\right),$$

where $\alpha \equiv \alpha (\nu ,\overline{\nu })$ and $\beta \equiv \beta (\nu ,\overline{\nu })$ are given by the Equation (134). The remaining ordinary integrations in this formula have been done numerically (cf. Ref. [77]). The so-calculated quantity $-\frac{ln(S\left(x\right)/{\sigma }_{\mathrm{f}}^2)}{d}$ asymptotically approaches the value of 6.0 at large d’s. Therefore, for ${\sigma }_{\mathrm{f}}={\left(440\mathrm{MeV}\right)}^2$, we obtain

$$\mathrm{mass}\mathrm{of}\mathrm{the}2g\mathrm{gluelump}=6\sqrt{{\sigma }_{\mathrm{f}}}\simeq 2.6\mathrm{GeV}.$$

(135)

This result is remarkably close to 2.56 GeV, which is the value obtained in Ref. [76] within the Hamiltonian approach. Thus, we have analytically shown that nonperturbative-nonconfining and confining self-interactions of the stochastic background Yang-Mills fields can have two different correlation lengths. Namely, we have obtained these lengths from the Green functions of the one- and the two-gluon gluelumps. This has been done by means of an analytic calculation of the world-line integrals describing one and two gluons (inter)connected by strings with the static source of the gauge field. The resulting inverse correlation lengths are given by Equations (131) and (135). Their numerical values agree remarkably well with those obtained by means of a different, Hamiltonian, approach. In Appendix C, we have also demonstrated that, in accordance with the lattice results of Ref. [115], the correlation length of nonperturbative-nonconfining interactions of the background fields decreases with the deformation of the path interconnecting the two points in the correlation function (120).

Finally, as has been announced at the beginning of this Section, let us discuss a possible modification of the 4D dual Abelian Higgs model, which can allow for the existence of different correlation lengths of the functions $D\left(x\right)$ and $D_1\left(x\right)$. For concreteness, we consider the London limit of this model. In accordance with the above-discussed results on the vacuum correlation lengths in QCD, our aim is to modify the dual Abelian Higgs model in such a way as to make the correlation length of the function $D_1\left(x\right)$ larger than the correlation length of the function $D\left(x\right)$. This aim can be achieved by assuming that, in addition to the condensate of magnetic monopoles, the vacuum of the model contains also the condensate of electrically charged particles. While the monopole condensate is described by the magnetically charged dual Higgs field, the condensate of electrically charged particles can be described by some other, electrically charged, Higgs field. The Higgs mechanism yields then the masses of the vector boson, as well as of the dual vector boson, which we denote respectively as $M_A$ and $M_B$. For a test particle, which is immersed into such a medium and carries an electric charge e, the vacuum-to-vacuum transition amplitude is given by the product of the ’t Hooft- and the Wilson-loop averages, $\langle H\left(C\right)\rangle \langle W\left(C\right)\rangle$. The ’t Hooft-loop average reads (For simplicity, we omit here the contribution of closed dual strings. In the leading approximation, such strings can be neglected due to their small sizes (cf. Ref. [118,119,120,121,122]).)

$$\langle H\left(C\right)\rangle =\int \mathcal{D}{B}_{\mu }exp\left\{-\int d^{4}x\left[\frac{1}{4}{(F_{\mu \nu }\left[B\right]+e{\tilde{\mathsf{\Sigma }}}_{\mu \nu })}^2+\frac{M_B^2}{2}{B}_{\mu }^2\right]\right\},$$

(136)

where $F_{\mu \nu }\left[B\right]={\partial }_{\mu }{B}_{\nu }-{\partial }_{\nu }{B}_{\mu }$ is the strength tensor of the dual field $B_{\mu }$, and ${\mathsf{\Sigma }}_{\mu \nu }(x,C)={\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x\left(\mathbf{w}\right)\right)\delta (x-x\left(\mathbf{w}\right))$ is the surface tensor defined at the minimal surface $S_{\min }$ bounded by the contour C. (Here, $\mathbf{w}$ is a 2D vector, and $x_{\mu }\left(\mathbf{w}\right)$ is a vector-function that parametrizes the surface $S_{\min }$.) Performing in Equation (136) the $B_{\mu }$-integration, and using the relation ${\partial }_{\nu }{\mathsf{\Sigma }}_{\mu \nu }=j_{\mu }$, where $j_{\mu }(x,C)\equiv {\oint }_{C}d{x}_{\mu }\left(\tau \right)\delta (x-x\left(\tau \right))$ is the current associated with the contour C, We have

$$\langle H\left(C\right)\rangle =exp\left\{-\frac{e^2}{2}\left[\frac{M_B^2}{2}{\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(x\right){\int }_{S_{\min }}d{\sigma }_{\mu \nu }\left(y\right)+{\oint }_{C}d{x}_{\mu }{\oint }_{C}d{y}_{\mu }\right]D_{M_B}(x-y)\right\},$$

(137)

where $D_m\left(x\right)=m{K}_1\left(m\right|x\left|\right)/(4{\pi }^2\left|x\right|)$ is the Yukawa propagator in 4D. The Wilson-loop average,

$$\langle W\left(C\right)\rangle =\int \mathcal{D}{A}_{\mu }exp\left\{-\int d^{4}x\left[\frac{1}{4}{F}_{\mu \nu }^2\left[A\right]+\frac{M_A^2}{2}{A}_{\mu }^2+ie{A}_{\mu }{j}_{\mu }\right]\right\},$$

can be calculated through the straightforward $A_{\mu }$-integration, which yields

$$\langle W\left(C\right)\rangle =exp\left[-\frac{e^2}{2}{\oint }_{C}d{x}_{\mu }{\oint }_{C}d{y}_{\mu }{D}_{M_A}(x-y)\right].$$

(138)

As follows from Equations (137) and (138), the function $D\left(x\right)$ describing confining interactions, falls off at large distances as $D\left(x\right)\sim {\mathrm{e}}^{-M_B\left|x\right|}$, while the function $D_1\left(x\right)$ describing nonconfining interactions, falls off as $D_1\left(x\right)\sim {\mathrm{e}}^{-M_B\left|x\right|}+{\mathrm{e}}^{-M_A\left|x\right|}$. Therefore, for $M_A<M_B$, the correlation length $1/M_A$ of the nonconfining interactions is larger than the correlation length $1/M_B$ of the confining interactions. Thus, the difference in the vacuum correlation lengths of the confining and the nonconfining interactions of the stochastic background Yang-Mills fields can be accommodated into the dual-superconductor model of confinement through the extension of this model by the condensate of electrically charged particles.