Energy-driven disorder in the mean field QCD

An impact of the finite size effects on the vacuum free energy density of full QCD with $N_{\rm f}$ massless flavors in the presence of homogeneous (anti-)self-dual Abelian background gluon field is studied. The zero temperature free energy density of the four-dimensional spherical domain is computed as a function of the background field strength $B$ and domain radius $R$. Calculation is performed in the one-loop approximation improved by accounting for mixing of the quark and gluon quasi-zero modes with normal modes, with the use of the $\zeta$-function regularization. It is indicated that, under plausible assumption on the character of the mixing, the quantum correction to the free energy density has a minimum as a function of $B$ and $R$. Within the mean field approach to QCD vacuum based on domain wall network representation of the mean field, an existence of the minimum may prevent infinite growth of individual domain, thus protecting the vacuum from the long-range ordering, and, hence, serving as the origin of disorder in the statistical ensemble of domain wall networks, driven by the minimization of the overall free energy of the dominant gauge field configurations.


I. INTRODUCTION
It is generally accepted that physical QCD vacuum can be characterized by various gluon, quark and mixed condensates. The condensates have played an important role in understanding the basic features of hadron physics. In particular, the lowest dimension condensates g 2 F 2 , (g 2F F ) 2 and ψ ψ are relevant to the anomalous breakdown of scale and U A (1) symmetries, spontaneous breaking of chiral SU L (N f ) × SU R (N f ) symmetry.
In principle, the condensates could be described in terms of the background (vacuum) gluon fields within the selfconsistent mean field approach to QCD vacuum. However, necessity to express the condensates, which are vacuum expectation values of the color neutral composite fields, in terms of the vector potential of the background gauge field in pair with the strong coupling regime complicates the actualization of the mean field approach to QCD vacuum.
Depending on base standpoint, the gauge mean field configurations underlying the condensates have been taken in various forms ranging from superposition of quasi-classical gluon configurations like the instanton gas or liquid [1,2] to the fields with constant field strength squared, g 2 F 2 = const, representing the global minimum of the quantum effective action of QCD. Properties of the quantum effective action for homogeneous gluon fields were studied in various approaches [3][4][5][6][7][8][9][10][11][12]. In particular, Leutwyler has demonstrated that the covariantly constant Abelian (anti-)self-dual fieldB B µν x ν ,B µν = ±B µν B µρ = B 2 δ νρ ,n = t 3 cos ξ + t 8 sin ξ, is singled out from other gluon fields as the only gauge field configuration with constant strength which is stable against small gluon and quark fluctuations and, hence, could be considered as likely contender for the global minimum of the effective action [8]. The stability was understood as the absence of the tachyonic modes in the spectrum of small quantum fluctuations in a given background gluon field. An obvious shortcoming of the purely homogeneous gauge field as a candidate for the mean vacuum field is that it would describe a globally ordered vacuum state and thus break all the symmetries of QCD. However, this argument does not apply to the more complicated relevant case when the gluon fields minimizing quantum effective action, the vacuum fields, belong to a set of lumpy configurations corresponding to distributed in R 4 domains of homogeneous Abelian (anti-)self-dual gluon field with size and shape randomly varying around certain mean values [13]. As a whole, such a set can be characterized as the statistical ensemble of gauge fields being homogeneous Abelian (anti-)self-dual almost everywhere in R 4 besides the boundaries between domains where the field appears to be neither homogeneous nor Abelian (anti-)self-dual. It has to be noted that topological charge density distribution in the typical gauge field configurations was studied within Lattice QCD with dynamical quarks with the results supporting the picture of entangled space-time regions of sign-alternating topological charge [14][15][16]. An instance of disordered lumpy configurations is well seen within the Ginzburg-Landau (GL) modeling of quantum effective action of QCD [13,17,18]. In this approach almost everywhere homogeneous Abelian (anti-)self-dual gluon field can be represented as the domain wall networks, arising straightforwardly as soon as the existence of the nonzero scalar gluon condensate g 2 F 2 is assumed. Domain wall networks come out of the structure of the degenerate discrete global minima of the GL effective potential, related to each other via discrete symmetry transformations -CP and Weyl reflections in the root space of color su(3) algebra.
Domain-structured configurations and homogeneous field are almost indistinguishable with respect to important bulk properties like For domain wall network with sufficiently thin boundaries between domains, both gauge field invariants F 2 and F F = ±F 2 are nonzero and constant almost everywhere in R 4 . Almost everywhere homogeneous Abelian (anti-)self-dual gluon fields have been incorporated into hadronization scheme within the domain model of QCD vacuum [19][20][21][22][23][24]. The Abelian (anti-)self-dual nature ensures confinement of both dynamical and static quarks (absence of poles in the propagators of color charged fields as well as fulfillment of the area law for Wilson loop), resolution of the U A (1) problem, and spontaneous breakdown of chiral symmetry. With minimal set of parameters, the domain model provided universal and rather accurate description of the masses and various decay constants of light, heavy-light mesons and heavy quarkonia, including their excited states, as well as some form factors. The picture of QCD vacuum based on Abelian (anti-)self-dual mean field turned out to be suggestive for exposing a catalyzing impact of strong electromagnetic field on quark deconfinement [13,17,[25][26][27].
On the whole, these rather satisfactory phenomenological applications put forward a task to clear up the mechanism behind the balance between the competitive tendencies for a long-range order in the ground state and disorder that may originate from two complementary origins -the topologically stable defects in the background field and existence of a minimum of the effective action with respect to the size of the regions of homogeneity (domain size). Topologically nontrivial gluon field configurations are expected to emerge through the division of the arbitrary gauge field into Abelian and non-Abelian parts, and may be seen in the domain wall network as frustrations of the color and spacetime orientation of the background field at the domain wall junctions. We shortly comment on this mechanism below but do not discuss it in detail here. Discussion of relevant physics can be found in review [28] and references therein.
Existence of a minimum of the effective action with respect to the domain size can be called the energy-driven origin of disorder, which is the main subject of the present paper. The quark and gluon quasi-zero modes characteristic for the covariantly constant Abelian (anti-)self-dual background field in a finite region may play peculiar role in formation of domains. Infinite number of quasi-zero modes become degenerate zero modes in the limit of the infinite domain size and thus potentially lead to an infra-red singular behaviour of the effective potential. In the plain one-loop approximation, contribution of the quark and gluon quasi-zero modes to the effective potential have opposite signs, and their strong concurrent impact on the potentially IR singular behavior of the effective potential becomes manifest in the infinite volume limit.
In the present paper the effective potential for the Abelian (anti-)self-dual field (1) in a four-dimensional spherical domain with radius R is calculated for SU (3) chromodynamics with N f massless quarks using the zeta function regularization, which completes the previously reported results for pure gluodynamics [29]. The quark and gluon fields are subject to the bag-like boundary conditions [13,21]. It is important that interaction of gluon and quark quasi-zero modes ( zero modes in the limit R → ∞) with the normal modes has to be treated beyond one loop as it has been put forward by Leutwyler [8]. This interaction leads to the contribution of the quasi-zero modes to the effective potential being regular in the limit of infinite domain size. This is crucial for overall existence of the thermodynamic limit and, as underlined in [8], consistency of the strong field limit of the effective action with the asymptotic freedom. In fact, a mixing of quasi-zero and normal modes has to be taken into account in all-loop orders as there is no good small parameter, which presents a hardly solvable task especially in the finite region. However, as it was demonstrated in [8], the most significant consequence of the mixing is the emergence of the effective "mass" for the quasi-zero modes. This observation enables sensible modeling of the plausible form of dependence of the "effective mass" on the domain size R, and identification of conditions for the existence of the minimum of the effective potential with respect to the radius of the domain R and field strength B inside domain.
It is shown that formation of the domains with finite size R can be energetically preferable if the effective "mass" falls from asymptotic nonzero value at R → ∞, fixed by the asymptotic freedom and correct strong field limit, to zero at R → 0 as it follows from the dependence of normal modes on the domain size. The final result for the free energy density of full QCD is shown in Fig. 12. The minimum in field strength and domain size is clearly present for a wide class of functional dependence of the zero mode effective "mass" on the domain size.
Finite mean size of the domains in the network is determined by the minimum of the effective action density inside individual domains. Existence of the mean domain size minimizing action density can be seen as a condition for sustainability of a domain wall network. In general, there is an infinite set of networks with degenerate values of quantum effective action that constitutes a statistical ensemble of dominant vacuum gluon field configurations. In accordance with the character of free energy dependence on domain size, the degree of disorder in the ensemble may vary from highly disordered distribution of entangled domains with variable size and shape to the almost periodic distribution of identically sized and shaped domains, reminiscent of spin liquid and (anti-) ferromagnetic state respectively.
In the next section we update the formal framework underlying the domain wall network representation of the background gluon fields in order to adjust it to the formulation of the particular problem studied in this paper. The ghost, gluon and and quark contributions to the free energy density of a spherical domain with Abelian self-dual gauge field are discussed in detail in the third section. Appendices contain quite involved technical details of calculations, which are the result of the present paper in themselves.

II. EFFECTIVE ACTION OF QCD AND THE DOMAIN WALL NETWORKS
The initial, gauge unfixed, Euclidean functional integral representation for QCD partition function assumes certain choice of the functional spaces of integration over gluon and quark fields. If one allows nonzero gluon condensates then the functional space F has to be subjected to an appropriate condition, for instance This definition is reminiscent of the Schrödinger functional representation (see for instance [8,30,31]. The difference is that condition (2) is imposed onto the gauge invariant combination of the gauge fields and has an integral (functional) form. Division of the general gauge fields A a µ into the background fields B a µ with extensive classical action specified by Eq. At this step, functional spaces Q and Ψ are restricted by the condition which excludes long-range fields from the set of fluctuations. Integral over the quark ψ and gluon Q fluctuations defines effective action S V eff [B] for a given background field B. Whether the constant B vac is nonzero has to be determined by the minima of the quantum effective action. In the infinite volume limit V → ∞ the global minima of S V eff [B] dominate the integral over background fields B and thus determine the specific class of gluon field configurations relevant to the self-consistent mean field description of QCD vacuum.
Fields with a constant field strength have to be verified for the role of the vacuum mean field B a µ first. Such a verification has been going since late seventies when chromomagnetic covariantly constant Abelian gauge field has been suggested for the role of QCD vacuum [4]. As has already been mentioned, covariantly constant Abelian (anti-)self-dual gauge field has appeared to be more preferable sample in many aspects, in particular due to the direct relation to confinement of dynamical color charged fields, chiral symmetry breaking, seen also in terms of meson properties through hadronization.
It should be stressed that condition (2) in no way restricts the background field functional space B to configurations with constant strength. Representation (4) assumes division of arbitrary gauge field A into two parts, the background B and fluctuation Q fields with simultaneous gauge fixing for the fluctuation part. If one is going to study Abelian background B then a specific parametrization of the gauge field has to be implemented [28,[33][34][35][36][37]). Both color and space orientations of the background field may appear to be frustrated at some space-time locations thus making manifest topological singularities in the vector potential, which, in general, cover the whole range of defects of various dimensions -domain wall, vortex, monopole and zero-dimensional instanton-like defects.
More importantly, as it has been stressed, in particular, by Faddeev, quantum equations "could have soliton solutions, which are absent in the classical limit. In particular, it is not completely crazy idea that quantum Yang-Mills equations have soliton-like solutions due to the dimensional transmutation" [38].
We have to conclude that representation of the dominant vacuum gauge field configurations by the plain covariantly constant field, corresponding to a globally ordered ground state, is an extreme case. Non-homogeneous at certain space-time locations field configurations with extensive classical action may have lower effective action than constant fields or can be topologically protected and should be taken into consideration anyway.
In general, fields which dominate the functional integral in the infinite volume (thermodynamic) limit belong to functional subspaceB ⊂ B. Given this, in the infinite volume limit one may reduce the integration over background fields B in (4) to the subspace of dominant (vacuum) fields and arrive at the mean field representation of the QCD partition function where Dσ B is a measure of integration over the space of vacuum gauge fields. A treatment of these possibly solitonlike vacuum fields B ∈B in the functional integral (4) must be nonperturbative, while perturbative expansion over fluctuations is likely to be applicable for calculation of various physical quantities [23,24]. It has to be noted that according to condition (4) normalization constant N should contain a factor which cancels a trivial extensive contribution of the background field to the classical action S[B] in the infinite volume limit. Certain prescriptions for regularization and renormalization of UV divergences are assumed. Due to the dimensional transmutation in gauge theories the condensates have to be expressed in terms of internal scale Λ QCD , that is B vac ∝ Λ 2 QCD . Interrelation between dimensional transmutation in gauge theories and gauge field condensates was discussed in [39].
The properties of the effective action are crucial for practical implementation of the described scheme. Details of discussion of the very existence of nontrivial minimum of the quantum effective action for homogeneous fields can be found in papers [3][4][5][6][7][8][9][10][11][12]. Some basic properties have been estimated and can be used for identification of the likely features of the mean field subspaceB.
An important observation was that in the infinite volume limit calculation of the effective potential for covariantly constant gauge fields could not be performed within the plain one-loop approximation because of the presence of infinitely degenerate tachyonic (for chromomagnetic field) or zero modes (for (anti-)self-dual field). One has to improve the one-loop calculation by accounting for interaction of the zero modes with normal modes [8] in order to generate a well-defined Gaussian measure for zero modes. The issues of tachyonic and zero modes were addressed in papers [10][11][12] on the basis of specific scaling properties of the massless QCD, concluding that tachyonic and zero modes seem to be the artifacts of the one-loop approximation. Functional renormalization group calculation of the effective potential for Abelian (anti-)self-dual covariantly constant field [3] does not encounter any problem with gluon zero modes. It has confirmed earlier estimates based on the improved one-loop approximation for the effective potential [8] and existence of a minimum at nonzero homogeneous Abelian (anti-)self-dual gluon field.
Ginzburg-Landau approach to the quantum effective action of QCD with effective Lagrangian [13,17,18,21] indicated an intrinsic possibility for disordered ground state of QCD. Here Λ is a scale, F a µν is the standard strength tensor for SU c (3) color gauge field,f µν = T a F a µν /Λ 2 , D ab µ = δ ab ∂ µ − iȂ ab µ . The effective Lagrangian respects all symmetries of QCD besides scale invariance, and the real constants C i have to be positive to provide a minimum of the effective potential at nonzero gauge field. Given this, one can check that there is a discrete set of global minima corresponding to the covariantly constant Abelian (anti-)self-dual fields where the matrixn k belongs to the Cartan subalgebra of su(3) These minima are connected with each other by discrete parity and Weyl symmetry transformations. One concludes that the domain wall solutions of the quantum equations of motion emerge as soon as the effective action has a global minimum corresponding to nonzero gluon condensate g 2 F 2 . For instance, if all parameters of the field besides angle ω between chromoelectric and chromomagnetic fields are put to the vacuum values then initial GL Lagrangian describes a sine-Gordon field ω, The available standard kink solution describes a planar domain wall between the regions with homogeneous Abelian self-dual and anti-self-dual gluon fields. Topological charge density vanishes on the wall where the chromomagnetic and chromoelectric fields are orthogonal to each other. More general domain wall correspond to ω, ξ and b varying simultaneously. In this case equations of motion read: The plain domain wall solution of these equations is shown in Fig. 1. It corresponds to gauge field which interpolates between two different vacuum configurations. In general, initial GL Lagrangian may produce more nontrivial soliton-like solutions both in Euclidean and Minkowski metrics. A combination of the additive and multiplicative superpositions of the domain walls allow one to generate various domain walls and domain wall networks in R 4 [13,40], like samples shown in Fig. 2. Lagrangian (6) has the simplest form, but its symmetry properties and emergence of the periodic discrete minima as a consequence of the scale invariance breakdown seem to be a general property, qualitatively insensitive to the detailed form of the effective potential. Another form of the strong field behavior of the GL effective Lagrangian can affect the particulars of the kink solution but can hardly influence its very existence and general properties. It has to be noted that the role of Weyl reflections in topology of QCD ground state has been intensively discussed in recent years in the context of dual superconductor picture of confinement [28,41].
The functional spaceB in the integral (5) is assumed to include infinitely many networks with equal values of the free energy. Though implicit implementation of this prescription within the domain model of QCD vacuum demonstrated high phenomenological performance [23,24], the conceptual problem of the network stability remains: since domain wall configurations in four-dimensional Euclidean space are not topologically protected and the presence of kink configuration usually increases action then any network should evolve to a single infinitely large domain. In the next section we study a particular effect which may prevent an infinite growth of a single domain. Red and blue colors correspond to Abelian self-dual and anti-self-dual field.

III. FREE ENERGY DENSITY
Lagrangian (6) does not account for possible finite size effects. Available evaluations of the effective potential for homogeneous gauge field which were performed in the infinite space-time [3][4][5][6][7][8][9][10][11][12], thus, as a matter of fact, implicitly assuming that the free energy density in the region with homogeneous field does not depend on the size of the region. Meanwhile, finite size effects are not excluded and may prevent infinite growth of an individual domain, thus protecting overall stability of the domain wall network configurations.
The present section is devoted to the study of the dependence of the renormalized free energy density F (B, R) of a finite spherical domain of the Abelian (anti-)self-dual homogeneous field on its radius R and the strength B of the field in full QCD with massless quarks, defined by the finite volume partition function where V R is the volume of four-dimensional ball with radius R, which is an idealization of the domain shown in the leftmost picture in Fig.2, and S V R is the gauge-fixed action of QCD in the presence of the background gluon field defined by Eq.(1). The Feynman background gauge is used below. Normalization constant N is fixed by the condition Renormalization prescription is specified below. Functional spaces Q, Ψ and C contain the quark, gauge and ghost fields subject to the bag-like boundary conditions These boundary conditions were discussed in papers [13,21]. The choice assumes that there is a physical boundary of the spherical domain, given by the domain wall illustrated in Fig.1. Bag-like boundary conditions are required by the qualitatively different character of field fluctuations in the bulk of domain (confining self-dual background field) and on the boundary (chromomagnetic field) [13]. The functional integral (9) is defined through decomposition of the quark, gauge and ghost fields, over orthogonal normalized complete set functions in Ψ, Q and C. It is convenient to diagonalize quadratic part of the action, using the eigenfunctions of the corresponding differential operators, −D 2 C n = λ gh n C n , subject to boundary conditions (11), (12) and (13). Indices n denote all relevant quantum numbers as described below. The quark, gluon and ghost spectra are purely discrete for any R if field strength B is nonzero. At finite R, all eigenvalues are nonzero for quark fields, and positive for gauge and ghost fields. The one-loop correction δU to the classical action is given by the determinants: where N f is the number of massless quark flavors. Renormalized functional determinants are calculated below by means of analytical regularization where λ n are eigenvalues of operator ∆. Computation of ζ(s) is based on the method summarized in papers [42,43].
A. Ghost contribution to the free energy density Evaluation of the ghost contribution to the free energy density is rather straightforward. Details of solution of the eigenvalue problem, Eqs. (15) and (13), are given in paper [21]. Ghost eigenfunctions are expressed in terms of confluent hypergeometric function M (a, b, z). The eigenvalues are defined by equation (15), which can be written in the form where v a is the absolute value of the a-th nonzero eigenvalue of the matrixn = n a T a , and T a are generators in the adjoint representation. The eigenmodes corresponding to the zero eigenvalues ofn do not contribute to the free energy density F (B, R) due to the normalization condition (10). Eigenvalues with given k, m, radial number r and color index a are (k + 1)-degenerate. Ghost contribution takes the form where Tr denotes summation over v a . Dimensionless quantities have been introduced using auxiliary renormalization scale µ. In the limit B → 0 Eq. (18) transforms to (see Appendix D) which defines eigenvalues λ(0, R) in Eq. (19). Zeta function can be written in the form [42,43] ζ gh (s) = Tr sin πs π This expression is still not suitable for analytical continuation to s → 0, since intervals of convergence of the integral and sum do not overlap. To make it ready for continuation to s → 0, several terms of asymptotic series in k of the integrand are added and subtracted: The first term is analytical at s → 0. The sums and integrals in the second term are expressed via analytical functions in the regions of s where they converge, and analytically continued to s → 0 in the complex plane (see Appendix A for details). The final expression for δU gh looks as The first term is a convergent sum, ready for numerical computation. The ghost contribution F gh (B, R) = δU gh (B, R)/V R to the free energy density defined by Eq. (22) is shown in Fig. 3. It demonstrates expected behaviour both in the field strength B and domain size R. In the infinite volume limit it approaches expected one-loop ghost contribution.

B. Gauge field fluctuations
Eigenvalues of relevant gluon operator −D 2 δ µν + 2iB µν (23) differ from the the ghost eigenvalues just by an overall shifts due to the term 2iB µν . Though in the infinite space-time this shift leads to the presence of the infinitely many exact gluon zero modes, for R < ∞ the degeneracy of these modes disappears and all corresponding eigenvalues become positive. This subset of modes will be referred below as quasi-zero modes. Since exact zero modes are absent for R < ∞, then, at first glance, one may hope to compute gluon contribution straightaway with the result that is usable at least for not very large values of the domain size R. The calculation goes completely analogously to the computation performed in the previous subsection, and leads to δU gl : Free energy density, as it comes out of Eq. (24), is shown in Fig. 4. One would expect that the character of the strong field limit should be insensitive to the presence of the boundary, since it corresponds to the short distances. Meanwhile, free energy density decreases without bound at large B and fixed R, which does not comply with known results [5,8], and, more generally, with asymptotic freedom. It also does not approach a constant at large R and fixed B, that means the absence of sensible thermodynamic limit in the system. This behavior is due to the manifestation of infinitely many quasi-zero eigenvalues that tend to zero as dimensionless quantity BR 2 → ∞ (see Fig. 5), which occurs both in the strong field and thermodynamic limits.
If all eigenvalues are sufficiently large to provide Gaussian damping in the functional integral (at small BR 2 ), then formula (17) can be considered as justified. The smaller the eigenvalues, the worse the one-loop approximation becomes, and finally it results in a strong field limit that is upside down . It becomes clear that one has to take into account the mixing between normal and quasi-zero modes, that means going beyond one-loop approximation of the free energy at large BR 2 . Calculation in the effective potential for Abelian self-dual field in the infinite space-time gives a guiding prescription [8], based on the observation that if one evaluates functional integral over normal (nonzero) modes first and accounts for their interaction with zero modes, than the obtained effective action has a finite quadratic in zero modes part. In other words, due to the interactions zero modes gain an effective "mass" µ 2 0 =κB, which provides one with an appropriate Gaussian measure. Hereκ is a constant. Schematically, contribution to effective "mass" is shown in Fig. 6. The filled circles denote all possible diagrams which include the propagators of normal gluon and quark modes. For SU (2) gluodynamics the lowest order value of the zero mode "effective mass" takes the value [8] Final result for the free energy density in the infinite volume agrees completely with renormalization group estimate [5] and, as it should be, with asymptotic freedom [8].
For finite volume case, this observation suggests that, in particular, interaction of normal modes with quasi-zero modes generates a shiftλ whereλ 2 (B, R) are quasi-zero eigenvalues of operator (23), and the function κ(BR 2 ) should approachκ both for infinite volume and in the strong field limit. The normal mode propagators, involved into diagrams in Fig.6, can be represented at most as an infinite series over quantum numbers of the modes. Truncation of these series is unreliable since in general diagrams are UV-divergent. Given that, calculation of the dependence of κ on dimensionless quantity z = BR 2 appears to be an extremely complicated task, even in the lowest perturbation order. Moreover, in the absence of small expansion parameter, perturbative expansion can not lead to a decisive result anyway.
Meanwhile, it seems to be possible to identify the general form of κ(z) suitable for qualitative estimate of the available fundamentally different dependencies of the free energy density on the domain size. Two restrictions for the function κ(z) can be identified. The first one is given by Eq. (25). As has already been noted, it follows from the existence of thermodynamic limit and agreement with the asymptotic freedom and the strong field limit. Another restriction, follows from the scaling of all eigenvalues at small R which means that corrections to the effective action of quasi-zero modes coming from the diagrams (Fig. 6) are expected to vanish, in analogy with decoupling of the infinitely heavy particles. A trial function κ(z) can be taken in the form which additionally to above restrictions at z → 0 and z → ∞ also reflects change in the behavior of eigenvalues at certain value of z, when all normal eigenvalues start rising as R decreases, see Fig. 5. Function κ is plotted in Fig. 7. Incorporation of effective "mass" κ(BR 2 ) leads to the following gluon contribution to the effective action (see Appendix A for details of calculation): The corresponding free energy density is plotted in Fig. 8. The correct behavior of the free energy density for BR 2 → ∞, consistent with asymptotic freedom and existence of thermodynamic limit, is restored. It is seen that free energy density acquires a minimum at intermediate values of field strength and domain size.
Free energy density for pure SU (3) gluodynamics for finite domain of Abelian (anti-)self-dual gluon fields is given by sum of Eqs. 22 and 28. The result is illustrated in left-hand side of Fig. 12. Existence of the minimum of the free energy density as a function of two variables is clearly seen.  9. Dependence of quark eigenvalues on BR 2 . At BR 2 → ∞, infinite number of zero modes emerges. Normal modes come in positive-negative pairs λ ± . Note that λ + / = − λ − for finite BR 2 , but the equality is restored in the limit BR 2 → ∞. An infinite number of quasi-zero modes are not chiral for finite R and become chiral zero-modes for asymptotically large BR 2 .

C. Quark contribution
To complete the calculation for full QCD with massless quarks, we have to study the quark contribution to the effective potential where N f is number of quark flavors, λ(B, R) are eigenvalues of Dirac operator in a spherical domain of radius R with homogeneous Abelian (anti-)self-dual field with bag boundary condition (see Appendix B). Zeta function ζ q (s) can be split into two parts [44] ζ q (s) = cos(πs)ζ / D 2 For the purpose of the present study we need only ζ / D 2 . Parity-odd term η(s) contributes to the imaginary part of the effective potential and can be related to the U A (1) anomalous breakdown [22]. Parity-even part ζ / D 2 can be written in the following form: where A(λ, k, j, B, R) = 0 is the equation for eigenvalues. Proceeding in the same manner as in the previous section (see Appendix C for details), we arrive at the expression The quark contribution to the free energy density given by Eq. (31) and shown in Fig. 10 exhibits, due to the quark quasi-zero modes, the inconsistency with the correct strong field and thermodynamic limits similar to inconsistencies in the gluon effective potential given by (24) and illustrated in Fig. 4. Following the reasoning analogous to the case of gluons in the previous subsection, the improved calculation of the quark contribution has to take into account generation of the effective "mass" for quasi-zero eigenmodes of quarks due to the interaction of quark quasi-zero modes with normal gluon and normal quark modes (see right-hand side diagram in Fig. 6). The improved quark contribution reads , (32) similarly to Eq. (28). The effective potential is shown in Fig. 11. Combining ghost, improved gluon and quark contributions, one finds − B 2 R 4 48 (29 + 33γ − 33 log 2 + 33 log R) + 3B 4 R 8 10120 Corresponding free energy density is shown in right-hand side of Fig. 12. The total free energy density demonstrates a well-pronounced minimum as a function of field strength and domain size. Quark contribution does not change the result of pure gluodynamics qualitatively, though the field strength at the minimum is considerably reduced.

D. One-loop beta function
The effective action U = U cl + δU κ should not depend on renormalization scale µ [45,46], that is with classical action (V R is four-dimensional ball of radius R) and δU κ given by Eq. (33). Since only terms containing log R contribute to the Eq.(34), one obtains µ d dµ that is equivalent to the equation exposing the correct one-loop β-function of QCD with N f quark flavors.

IV. DISCUSSION
We have studied, as far as it has been possible with analytical methods, an influence of the finite size effects on the vacuum free energy density of full QCD with N f massless flavors in the presence of homogeneous (anti-)self-dual Abelian background gluon field. The most essential result is illustrated in the right-hand side Fig.12, where the zero temperature free energy density of the four-dimensional spherical domain is plotted as a function of the background field strength B and domain radius R. It indicates that the quantum correction to the free energy density may have a minimum at finite values of B and R. In the domain wall network representation of the vacuum mean field, existence of this minimum means that in the statistically dominant networks the individual domains should have finite size varying near the mean value, and infinite growth of individual domains is prohibited by the minimization of the overall free energy of the configuration. This suggests that the domain wall ensemble should include an infinite number of configurations with degenerate energy.
The character of this energy-driven disorder strongly depends on the details of behavior of the free energy density in the vicinity of the minimum. One may expect that shallow and flat profile near the minimum may lead to strong variations of the geometrical shape as well as deviations of the field strength from the mean value, which will allow the presence of highly irregular networks among the dominant configurations, characterized by the strong entanglement of domains, etc. The deep and steep profile would assume that spatially periodic networks should dominate, bringing the long-range (periodic) order into the mean field configurations.
The result of the present paper, Fig.12, is certainly has a status of preliminary rough estimate, since the validity of one-loop approximation is indeterminate, especially due to the indefiniteness of the treatment of interaction between quasi-zero and normal modes. A straightforward numerical calculation, within the lattice approximation for instance, could be useful. A difficulty for lattice calculation can be caused by the non-standard boundary conditions, which have to represent the physical boundary of a domain. However, this possibility does not look like hopeless in view of the recent lattice QCD calculations for rotating strongly interacting matter, where Dirichlet and Neumann boundary conditions have been implemented [47,48].

ACKNOWLEDGMENTS
We acknowledge useful discussions with Michael Bordag, Irina Pirozhenko, Victor Braguta and Artem Roenko. We start with the expression for ζ gh (s) given by (21). Suitable asymptotic expansion for hypergeometric function M at large k and fixed j/k can be found with the help of the method described in Ref. [49], Chapter 10, §9. We obtain for fixed z 0 and arbitrary constants m, n. The coefficients A i (z) are found with the help of recursion relation (i 0, A 0 = 1) The constants of integration are fixed by the requirement lim t→∞ A i (z) = 0, i 1.
Asymptotic expansion for modified Bessel function I is given by (see [50], 10.41.3, 10.41.7, 10.41.9) where Now, we substitute formulas (A1) and (A2) into Ψ gh in formula (21), expand it in powers of k −1 , sum over m and compute derivative with respect to t. We arrive at the following expressions for u gh The sums over k are calculated for Rs > 1 and analytically continued in terms of Riemann ζ function ∞ k=1 k 1−2s k i = ζ(i − 1 + 2s).
Since the whole spectrum is invariant with respect to B → −B, trace over color leads to factor 4. Evaluating the derivative of ζ gh (s) with respect to s at s = 0 one arrives at the expression (22).

Gluons
With Dirichlet boundary condition for color-charged modes the whole set of eigenvalues is determined by the equations [21] M k 2 and every solution of these equations with given k, m and a is 2(k + 1)-degenerate. Repeating the procedure carried out for ghosts, one finds We substitute formulas (A1) and (A2) into Ψ gl , expand it in powers of k −1 , sum over m and compute derivative with respect to t. Coefficients of asymptotic expansion u gl i are given by Finally,

Contribution of quasi-zero modes
In this section, we calculate the contribution of gluon quasi-zero modes λ ↑k k 2 0 to the effective potential. These modes correspond to the smallest-magnitude solutions to the equations at B → 0. The contribution of gluon quasi-zero modes to the effective potential can be expressed as where factor 2 in the definition of ζ gl(0) (s) originates from two polarizations of quasi-zero gluon modes, and color trace yields factor four. For the sake of brevity we omit color eigenvalue and restore it in the final answer (B → √ 3 2 B for adjoint representation of su (3)) To continue ζ gl(0) (s) to s → 0, we add and subtract several terms of asymptotic expansion in k found with the help of formula (E5) The first sum is an analytic function for Rs > 0. The second sum is evaluated for Rs > 1/2 and analytically continued to s → 0: Finally, we obtain 4. Contribution of quasi-zero modes with the effective "mass" If one includes the effective "mass" κ for quasi-zero modes in the considerations of the previous section, the formulas become In analogy to the previous section, The first sum is an analytic function for Rs > 0. The second sum is evaluated for Rs > 1/2 and analytically continued to s → 0: ∞ k=0 (k + 1) −2s R 2s BR 2 s − 2αBR 2 s(1 + s)(k + 1) −2/3 + BR 2 s 4 4 + BR 2 (1 + 2s) − 4κ BR 2 (k + 1) −1 = BR 2+2s s ζ(2s) − 2α(1 + s)ζ(2/3 + 2s) + 1 − κ BR 2 ζ(1 + 2s) + B 2 R 4+2s s 4 (1 + 2s)ζ(1 + 2s). Finally, 5. Contribution of all eigenmodes with the effective "mass" for quasi-zero modes The desired zeta function corresponding to one-loop correction with the effective "mass" for quasi-zero modes is written as where zeta functions in right-hand side are given by Eqs. (A3),(A5) and (A6). The corresponding effective potential is given by Thus obtained effective potential is not an even function of B. To restore invariance under B → −B, one adds term κ to contribution of modes λ 2 ↓k,− k Acting on equation (B1) with the projectors P ± = (1 ± γ 5 )/2, one rewrites it as where ψ ± = P ± ψ. Substituting one equation into the other (this is valid if λ = 0), one finds Only one chiral component is independent, the other one is found via Eq. (B3). The expression for γ µ B µν x ν P ± (γ µ B µν x ν originates from / D) is more complex than γ µ B µν x ν P ∓ , so we find ψ − from equation (B2) in the case of self-dual field and ψ + in the case of anti-self-dual field. Here and below, if ± or ∓ appears alongside with Dirac operator / D or its eigenmode ψ, the upper sign should be taken for self-dual field and the lower sign for anti-self-dual field. η is the normal to the surface of a sphere, η 2 = 1.
The analogue of total angular momentum in Euclidean space is [51] J µν = K µν + S µν , Algebra of operators K µν , S µν , J µν splits into two so (3) algebras and analogous relations for J ± . One introduces lowering and raising operators We choose reference frame such that direction of field coincides with z axis.
The solutions are ψ km1m2 Where Y k 2 m1m2 are spherical harmonics in four-dimensional Euclidean space [51]. The solutions of (B1) may be characterized by eigenvalues of independent operators J ∓2 , J ∓ 3 , J ± 3 , L (see [52][53][54]). Here and upper or lower signs should be taken for self-dual field and anti-self-dual field, correspondingly. The solutions to the eigenvalue problem (B1) where ψ km1m2 cs (c = ∓, s = ±) are solutions of Eq. (B2) given by Eqs. (B5) and (B6), diagonalize all these operators. Now, we impose bag boundary condition (11). Substituting Eq. (B4) in Eq. (11) and acting on it with operators P ± , we find that boundary condition reduces to where ψ ∓ are solutions of Eq. (B2). The boundary condition breaks symmetry associated with generator L. In order to satisfy boundary condition, we have to mix eigenstates of operator L.
Let us consider the case of self-dual field (the case of anti-self-dual field is obtained via m 2 → m 1 , α → −α). To satisfy boundary condition, linear combination of solutions with equal projection of total angular momenta j 3 is constructed. At given k, there is one solution for maximal value of j 3 = k 2 + 1 2 , one solution for minimal value j 3 = − k 2 − 1 2 and two solutions for each intermediate j 3 = − k 2 + 1 2 , . . . , k 2 − 1 2 (for k > 0). Below we consider these three cases.
The boundary condition reads Using identities (see [55], 13 we bring these equations to the form Boundary condition can be obtained from Eq. (B8) with C 1 = 0. After simplification and Kummer transformation (see [55], 13.1.27) M (a, b, z) = e z M (b − a, b, −z).
we obtain M λ 2
Quasi-zero solutions of (E1) at z → 0 become smallest zeroes of Bessel functions that correspond to t = 1. With t = 1 formula (E3) gives the desired expansion for the first zero of J k+1 (z) where α ≈ 1.855757. This series can be once again expanded in powers of k. According to Eq. (E4), x(0) ∼ k for k 1. One notices that power of leading-asymptotics term of derivative x (i) (0) is smaller for larger i x(0) ∼ k, dx(a) da a=0 ∼ k 0 , . . . so we need only several terms in series (E2) to find asymptotic expansion of x(z) in k up to given order: x(z) = k + αk 1/3 + 1 + 3α