# Energy-driven disorder in the mean field QCD

Sergei N. Nedelko<EMAIL_ADDRESS>Vladimir E. Voronin
<EMAIL_ADDRESS>Bogoliubov Laboratory of Theoretical Physics, JINR,
141980 Dubna, Russia

(2020)

###### Abstract

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 $\langle g^{2}F^{2}\rangle$,
$\langle(g^{2}\tilde{F}F)^{2}\rangle$ and $\langle\bar{\psi}\psi\rangle$ are
relevant to the anomalous breakdown of scale and $U_{\mathrm{A}}(1)$
symmetries, spontaneous breaking of chiral $SU_{\mathrm{L}}(N_{f})\times
SU_{\mathrm{R}}(N_{f})$ symmetry.

In principle, the condensates could be described in terms of the background
(vacuum) gluon fields within the self-consistent 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}={\rm 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 field

$\displaystyle\hat{B}_{\mu}=-\frac{1}{2}\hat{n}B_{\mu\nu}x_{\nu},\
\tilde{B}_{\mu\nu}=\pm B_{\mu\nu}\ B_{\mu\rho}=B^{2}\delta_{\nu\rho},\
\hat{n}=t_{3}\cos\xi+t_{8}\sin\xi,$ (1)

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 $\langle g^{2}F^{2}\rangle$ 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

$\displaystyle\lim_{V\to\infty}V^{-1}\int\limits_{V}d^{4}xg^{2}F^{a}_{\mu\nu}F^{a}_{\mu\nu}=B_{\mathrm{vac}}^{2}\not=0,$
$\displaystyle\lim_{V\to\infty}V^{-1}\int\limits_{V}d^{4}x|g^{2}\tilde{F}^{a}_{\mu\nu}F^{a}_{\mu\nu}|\simeq
B_{\mathrm{vac}}^{2}.$

For domain wall network with sufficiently thin boundaries between domains,
both gauge field invariants $F^{2}$ and $\tilde{F}F=\pm 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 [17, 25, 26, 13, 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 space-time 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_{\rm 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 [21, 13]. It is important that interaction of gluon and
quark quasi-zero modes ( zero modes in the limit $R\to\infty$) 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\to\infty$, fixed by the asymptotic freedom and correct strong
field limit, to zero at $R\to 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

$\displaystyle Z[B_{\rm vac}]=N\int\limits_{{\cal
F}}DA\int\limits_{\Psi}D\bar{\psi}D\psi\exp\\{-S[A,\bar{\psi},\psi]\\},$

assumes certain choice of the functional spaces of integration over gluon and
quark fields. If one allows nonzero gluon condensates then the functional
space ${\cal F}$ has to be subjected to an appropriate condition, for instance

$\displaystyle{\cal
F}=\left\\{A:\lim_{V\to\infty}\frac{1}{V}\int\limits_{V}d^{4}xg^{2}F^{a}_{\mu\nu}(x)F^{a}_{\mu\nu}(x)=B_{\rm
vac}^{2}\right\\}.$ (2)

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_{\mu}^{a}$ into the
background fields $B_{\mu}^{a}$ with extensive classical action specified by
Eq. (2) and the fluctuations $Q_{\mu}^{a}$ in the background $B_{\mu}^{a}$
supplemented by the background gauge condition $D(B)Q=0$ leads to the
representation [32]

$\displaystyle Z[B_{\rm vac}]$ $\displaystyle=$ $\displaystyle
N^{\prime}\int\limits_{{\cal B}}DB\exp\left\\{-S^{V}_{\rm eff}[B]\right\\}$
$\displaystyle=$ $\displaystyle N^{\prime\prime}\int\limits_{{\cal
B}}DB\int\limits_{{\cal
Q}}DQ\int\limits_{\Psi}D\bar{\psi}D\psi\det[D(B)D(B+Q)]\delta[D(B)Q]\exp\\{-S_{V}[B+Q,\bar{\psi},\psi]\\}.$

At this step, functional spaces ${\cal Q}$ and $\Psi$ are restricted by the
condition

$\displaystyle\lim_{V\to\infty}\left(S_{V}[B+Q,\bar{\psi},\psi]-S_{V}[B]\right)<\infty,$
(4)

which excludes long-range fields from the set of fluctuations. Integral over
the quark $\psi$ and gluon $Q$ fluctuations defines effective action
$S^{V}_{\rm eff}[B]$ for a given background field $B$. Whether the constant
$B_{\rm vac}$ is nonzero has to be determined by the minima of the quantum
effective action. In the infinite volume limit $V\to\infty$ the global minima
of $S^{V}_{\rm 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_{\mu}^{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 ${\cal B}$ to configurations with constant strength.
Representation (II) 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
[33, 34, 35, 36, 37, 28]). 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

$\displaystyle\frac{\delta S^{V}_{\rm eff}[B]}{\delta B}=0$

”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 subspace $\tilde{\cal
B}\subset{\cal B}$. Given this, in the infinite volume limit one may reduce
the integration over background fields $B$ in (II) to the subspace of dominant
(vacuum) fields and arrive at the mean field representation of the QCD
partition function

$\displaystyle Z[B_{\rm
vac}]=N^{\prime\prime}\lim_{V\to\infty}\int\limits_{\tilde{\cal
B}}D\sigma_{B}\int\limits_{{\cal
Q}}DQ\int\limits_{\Psi}D\bar{\psi}D\psi\det[D(B)D(B+Q)]\delta[D(B)Q]$
$\displaystyle\times\exp\\{-S_{V}[B+Q,\bar{\psi},\psi]\\},$ (5)

where $D\sigma_{B}$ is a measure of integration over the space of vacuum gauge
fields. A treatment of these possibly soliton-like vacuum fields
$B\in\bar{\cal B}$ in the functional integral (II) 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 (II) normalization constant $N^{\prime\prime}$ 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 $\Lambda_{\rm
QCD}$, that is $B_{\rm vac}\propto\Lambda^{2}_{\rm 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 subspace $\tilde{\cal B}$.

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 [21, 17, 13, 18]

$\displaystyle\mathcal{L}_{\mathrm{eff}}=-\frac{1}{4\Lambda^{2}}\left(D^{ab}_{\nu}F^{b}_{\rho\mu}D^{ac}_{\nu}F^{c}_{\rho\mu}+D^{ab}_{\mu}F^{b}_{\mu\nu}D^{ac}_{\rho}F^{c}_{\rho\nu}\right)-U_{\mathrm{eff}},$
$\displaystyle
U_{\mathrm{eff}}=\frac{\Lambda^{4}}{12}\mathrm{Tr}\left(C_{1}\breve{f}^{2}+\frac{4}{3}C_{2}\breve{f}^{4}-\frac{16}{9}C_{3}\breve{f}^{6}\right)$
(6)

indicated an intrinsic possibility for disordered ground state of QCD. Here
$\Lambda$ is a scale, $F^{a}_{\mu\nu}$ is the standard strength tensor for
$SU_{\rm c}(3)$ color gauge field,
$\breve{f}_{\mu\nu}=T^{a}F^{a}_{\mu\nu}/\Lambda^{2}$,
$D^{ab}_{\mu}=\delta^{ab}\partial_{\mu}-i\breve{A}^{ab}_{\mu}$. 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

$\breve{A}^{k}_{\mu}=-\frac{1}{2}\breve{n}_{k}F_{\mu\nu}x_{\nu},\,\tilde{F}_{\mu\nu}=\pm
F_{\mu\nu},\quad F_{\mu\nu}F_{\mu\nu}=b^{2}\Lambda^{4},\quad
b_{\mathrm{vac}}^{2}=\frac{-C_{2}+\sqrt{C_{2}^{2}+3C_{1}C_{3}}}{3C_{3}},$

where the matrix $\breve{n}_{k}$ belongs to the Cartan subalgebra of $su(3)$

$\breve{n}_{k}=T^{3}\ \cos\left(\xi_{k}\right)+T^{8}\
\sin\left(\xi_{k}\right),\quad\xi_{k}=\frac{2k+1}{6}\pi,\,k=0,1,\dots,5.$ (7)

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 $\langle g^{2}F^{2}\rangle$. For
instance, if all parameters of the field besides angle $\omega$ between
chromoelectric and chromomagnetic fields are put to the vacuum values then
initial GL Lagrangian describes a sine-Gordon field $\omega$,

$\mathcal{L}_{\mathrm{eff}}=-\frac{1}{2}\Lambda^{2}b_{\mathrm{vac}}^{2}\partial_{\mu}\omega\partial_{\mu}\omega-
b_{\mathrm{vac}}^{4}\Lambda^{4}\left(C_{2}+3C_{3}b_{\mathrm{vac}}^{2}\right){\sin^{2}\omega}.$

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.

Figure 1: Domain wall solution of Eqs. (8): gauge field interpolates between
two vacuum states, namely self-dual and anti-self-dual configurations and
$\xi$ going from one boundary of Weyl chamber to another, field strength $b$
has a dip inside the wall.

More general domain wall correspond to $\omega$, $\xi$ and $b$ varying
simultaneously. In this case equations of motion read:

$\displaystyle-6b^{\prime}\omega^{\prime}+b^{3}\sin
2\omega\left(C_{3}b^{2}(\cos
6\xi+10)+3C_{2}\right)-3b\omega^{\prime\prime}=0,$
$\displaystyle-15b^{\prime\prime}+3b\left(-4C_{1}+2\omega^{\prime
2}+5\xi^{\prime 2}\right)-12C_{2}b^{3}(\cos 2\omega-3)-2C_{3}b^{5}(3\cos
2\omega-5)(\cos 6\xi+10)=0,$ $\displaystyle 2C_{3}b^{6}(3\cos 2\omega-5)\sin
6\xi-15b\left(2b^{\prime}\xi^{\prime}+b\xi^{\prime\prime}\right)=0.$ (8)

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}$ [40,
13], like samples shown in Fig. 2.

Figure 2: Two-dimensional slices of topological charge density for various
superpositions of domain wall solutions for Eqs.(8). Red and blue colors
correspond to Abelian self-dual and anti-self-dual field.

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 [41, 28].

The functional space $\tilde{\mathcal{B}}$ 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.

##  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

$\displaystyle\exp(-V_{R}F(B,R))={\rm Ren}\ N\int\limits_{{\cal
Q}}DQ\int\limits_{\Psi}D\bar{\psi}D\psi\int\limits_{\cal C}D\bar{c}Dc\
\exp\left\\{-S_{V_{R}}[Q;\bar{\psi},\psi;\bar{c},c;B]\right\\},$ (9)

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

$\displaystyle F(0,R)=0.$ (10)

Renormalization prescription is specified below. Functional spaces $\cal{Q}$,
$\Psi$ and $\cal{C}$ contain the quark, gauge and ghost fields subject to the
bag-like boundary conditions

$\displaystyle\left.\left(i\not{\eta}e^{i\alpha\gamma_{5}}-1\right)\psi(x)\right|_{x\in\partial
V_{R}}=0,$ (11) $\displaystyle\left.\breve{n}Q_{\mu}(x)\right|_{x\in\partial
V_{R}}=0,$ (12) $\displaystyle\left.\breve{n}c(x)\right|_{x\in\partial
V_{R}}=0.$ (13)

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,

$Q(x)=\sum_{n}q_{n}Q^{(n)}(x),\quad
c(x)=\sum_{n}c_{n}C_{n}(x),\quad\psi(x)=\sum_{n}\theta_{n}\Psi_{n}(x),$

over orthogonal normalized complete set functions in $\Psi$, $\cal{Q}$ and
$\cal{C}$. It is convenient to diagonalize quadratic part of the action, using
the eigenfunctions of the corresponding differential operators,

$\displaystyle\hat{\not{D}}\Psi_{n}=\lambda_{n}^{\text{q}}\Psi_{n},$ (14)
$\displaystyle-\breve{D}^{2}C_{n}=\lambda_{n}^{\text{gh}}C_{n},$ (15)
$\displaystyle\left[-\breve{D}^{2}\delta_{\mu\nu}+2i\breve{B}_{\mu\nu}\right]Q^{(n)}_{\nu}=\lambda_{n}^{\text{gl}}Q^{(n)}_{\mu},$
(16)

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 $\delta U$ to the classical
action

$U_{\textrm{cl}}=\int_{V_{R}}d^{4}x\frac{1}{4g^{2}}\left(B^{a}_{\mu\nu}\right)^{2}=\frac{\pi^{2}B^{2}R^{4}}{2g^{2}}.$

is given by the determinants:

$\exp\left(-\delta
U\right)=N\left[\det\left(-\breve{D}^{2}\delta_{\mu\nu}+2i\breve{B}_{\mu\nu}\right)\right]^{-\frac{1}{2}}\
\det\left(-\breve{D}^{2}\right)\ \left[\det i\not{D}\right]^{N_{f}},$ (17)

where $N_{f}$ is the number of massless quark flavors. Renormalized functional
determinants are calculated below by means of analytical regularization

$\mathrm{Tr}\log\Delta=-\frac{d}{ds}\sum_{n}\left.\lambda_{n}^{-s}\vphantom{\frac{1}{1}}\right|_{s=0}=-\left.\frac{d}{ds}\zeta(s)\right|_{s=0},$

where $\lambda_{n}$ are eigenvalues of operator $\Delta$. Computation of
$\zeta(s)$ is based on the method summarized in papers [42, 43].

###  III.1 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

$M\left(\frac{k}{2}+1-m-\frac{\lambda^{2}}{2v_{a}B},k+2,\frac{v_{a}BR^{2}}{2}\right)=0,k=0,1,2,\dots,\
m=-\frac{k}{2},-\frac{k}{2}+1,\dots,\frac{k}{2},$ (18)

where $v_{a}$ is the absolute value of the $a$-th nonzero eigenvalue of the
matrix $\breve{n}=n^{a}T^{a}$, and $T^{a}$ are generators in the adjoint
representation. The eigenmodes corresponding to the zero eigenvalues of
$\breve{n}$ 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

$\delta
U^{\mathrm{gh}}=-\sum_{kmr}\mathrm{Tr}\ln\frac{\lambda_{kmr}^{2}(vB,R)}{\lambda_{kmr}^{2}(0,R)}=\left.\frac{d}{ds}\zeta_{\mathrm{gh}}(s)\right|_{s=0},$
(19)

where $\mathrm{Tr}$ denotes summation over $v^{a}$. Dimensionless quantities

$\lambda=\lambda/\mu,\quad B=B/\mu^{2},\quad R=R\mu$ (20)

have been introduced using auxiliary renormalization scale $\mu$. In the limit
$B\to 0$ Eq. (18) transforms to (see Appendix D)

$(k+1)!\left(\frac{\lambda R}{2}\right)^{-k-1}J_{k+1}(\lambda R)=0,$

which defines eigenvalues $\lambda(0,R)$ in Eq. (19). Zeta function can be
written in the form [42, 43]

$\displaystyle\zeta^{\text{gh}}(s)=\mathrm{Tr}\frac{\sin\pi
s}{\pi}\sum_{k=0}^{\infty}\sum_{m=-\frac{k}{2}}^{\frac{k}{2}}(k+1)\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\Psi_{\mathrm{gh}}(k+1,m,t,\breve{B},R),$
$\displaystyle\Psi_{\mathrm{gh}}(k,m,t,\breve{B},R)=\log\frac{\exp\left(-\frac{\breve{B}R^{2}}{4}\right)M\left(\frac{k+1}{2}-m+\frac{t^{2}}{2\breve{B}},k+1,\frac{\breve{B}R^{2}}{2}\right)}{k!\left(\frac{tR}{2}\right)^{-k}I_{k}(tR)}.$

This expression is still not suitable for analytical continuation to $s\to 0$,
since intervals of convergence of the integral and sum do not overlap. To make
it ready for continuation to $s\to 0$, several terms of asymptotic series in
$k$ of the integrand are added and subtracted:

$\displaystyle\zeta^{\text{gh}}(s)$ $\displaystyle=$
$\displaystyle\mathrm{Tr}\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\left[\sum_{m=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi_{\mathrm{gh}}(k,m,kt,\breve{B},R)-\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{gh}}(t,\breve{B},R)}{k^{i}}\right]$
(21) $\displaystyle+\mathrm{Tr}\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{gh}}(t,\breve{B},R)}{k^{i}}.$

The first term is analytical at $s\to 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\to 0$ in the complex plane (see
Appendix A for details). The final expression for $\delta U^{\mathrm{gh}}$
looks as

$\displaystyle\delta U^{\mathrm{gh}}(B,R)$ $\displaystyle=$
$\displaystyle-4\sum_{k=1}^{\infty}k\left[\sum_{m=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi_{\text{gh}}\left(k,m,0,\frac{\sqrt{3}B}{2},R\right)-\frac{3}{4}B^{2}R^{4}\frac{1}{48}\left(1-\frac{1}{k}+\frac{1}{k^{2}}\right)\right]$
(22) $\displaystyle+\frac{B^{2}R^{4}}{48}(2-3\gamma+3\log 2-3\log
R)-\frac{3B^{4}R^{8}}{10240}.$

The first term is a convergent sum, ready for numerical computation.

The ghost contribution $F_{\mathrm{gh}}(B,R)=\delta
U_{\mathrm{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.

Figure 3: Contribution of the Faddeev-Popov determinant to the free energy
density. Domain size $R$ and field strength $B$ are given in units of
renormalization scale $\mu$, Eq. (20). If $R$ is sufficiently large, then a
minimum exists at nonzero field strength $B$. For any $B$ the ghost
contribution is minimal in the infinite volume limit $R\to\infty$.

###  III.2 Gauge field fluctuations

Eigenvalues of relevant gluon operator

$-\breve{D}^{2}\delta_{\mu\nu}+2i\breve{B}_{\mu\nu}$ (23)

differ from the the ghost eigenvalues just by an overall shifts due to the
term $2i\breve{B}_{\mu\nu}$. Though in the infinite space-time this shift
leads to the presence of the infinitely many exact gluon zero modes, for
$R<\infty$ 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<\infty$, 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 $\delta U^{\mathrm{gl}}$:

$\displaystyle\delta U^{\mathrm{gl}}(B,R)$ $\displaystyle=$ $\displaystyle
4\sum_{k=1}^{\infty}k\left[\sum_{m=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi_{\text{gl}}\left(k,m,0,\frac{\sqrt{3}B}{2},R\right)-\frac{3}{4}B^{2}R^{4}\frac{1}{24}\left(1-\frac{1}{k}-\frac{5}{k^{2}}\right)\right]$
(24) $\displaystyle-$ $\displaystyle\frac{B^{2}R^{4}}{48}(31+30\gamma-30\log
2+30\log R)+\frac{3B^{4}R^{8}}{5120}.$

Figure 4: Dependence of one-loop gluon contribution to the free energy density
on $R$ and $B$ (in units of renormalization scale $\mu$, Eq. (20)). The strong
field and large size volume limits are incorrect due to the presence of
infinitely many gluon quasi-zero modes.

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}\to\infty$ (see Fig. 5), which occurs both in the strong field and
thermodynamic limits.

Figure 5: Dependence of gluon eigenvalues (see Eq. (38)) on $BR^{2}$. In the
limit $BR^{2}\to\infty$, infinite number of zero modes emerges.

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” $\mu_{0}^{2}=\bar{\varkappa}B$, which
provides one with an appropriate Gaussian measure. Here $\bar{\varkappa}$ 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]

$\displaystyle\mu_{0}^{2}=\bar{\varkappa}B,\
\bar{\varkappa}=\frac{g^{2}}{24\pi^{2}}.$

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

$\displaystyle\bar{\lambda}_{\mathrm{eff}}^{2}(B,R)=\bar{\lambda}^{2}(B,R)+\varkappa(BR^{2})B,$
$\displaystyle\lim_{R\to\infty}\bar{\lambda}^{2}(B,R)=0,\ \
\lim_{z\to\infty}\varkappa(z)=\bar{\varkappa}>0,$ (25)

where $\bar{\lambda}^{2}(B,R)$ are quasi-zero eigenvalues of operator (23),
and the function $\varkappa(BR^{2})$ should approach $\bar{\varkappa}$ 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 $\varkappa$ 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
$\varkappa(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 $\varkappa(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,

$\displaystyle\lim_{z\to 0}\varkappa(z)=0,$ (26)

follows from the scaling of all eigenvalues at small $R$

$\displaystyle\lim_{R\to 0}\lambda^{2}(z)\propto\lim_{R\to 0}1/R^{2}=0,$

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.

Figure 6:  Wavy line corresponds to gluon quasi-zero modes, dashed line
corresponds to quark quasi-zero modes. The filled circles denote all possible
diagrams with normal quark and gluon modes.

A trial function $\varkappa(z)$ can be taken in the form

$\varkappa(z)=\frac{2}{\pi}\left[\arctan\exp\left(\frac{z-z_{0}}{a}\right)+\arctan\exp\left(\frac{z_{0}-z}{a}\right)-2\arctan\exp\left(-\frac{z_{0}}{a}\right)\right],$
(27)

which additionally to above restrictions at $z\to 0$ and $z\to\infty$ 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
$\varkappa$ is plotted in Fig. 7.

Incorporation of effective “mass” $\varkappa(BR^{2})$ leads to the following
gluon contribution to the effective action (see Appendix A for details of
calculation):

$\displaystyle\delta U_{\varkappa}^{\mathrm{gl}}(B,R)$ $\displaystyle=$
$\displaystyle
4\sum_{k=1}^{\infty}k\left[\sum_{m=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\text{gl}}(k,m,0,\frac{\sqrt{3}B}{2},R)-\frac{3}{4}B^{2}R^{4}\frac{1}{24}\left(1-\frac{1}{k}-\frac{5}{k^{2}}\right)\right]$
(28) $\displaystyle-\frac{B^{2}R^{4}}{48}(31+30\gamma-30\log 2+30\log
R)+\frac{3B^{4}R^{8}}{5120}$
$\displaystyle+4\sum_{k=0}^{\infty}\left[(k+1)\log\frac{\lambda_{\uparrow
k\frac{k}{2}0}^{2}(B,R)+\varkappa(BR^{2})B}{\lambda_{\uparrow
k\frac{k}{2}0}^{2}(B,R)}+(k+1)\log\frac{\lambda_{\downarrow
k\frac{-k}{2}0}^{2}(B,R)-\varkappa(BR^{2})B}{\lambda_{\downarrow
k\frac{-k}{2}0}^{2}(B,R)}\right].$

Figure 7: The profile of function $\varkappa(z)$ used in Fig. 8.

Figure 8: Improved gluon contribution to the free energy with effective ”mass”
$\varkappa$ given by Eq. (27). Parameters of $\varkappa$ are taken such that
the minimum would be well-pronounced,$z_{0}=100$ and $a=10$. All quantities
are measured in units of renormalization scale $\mu$, Eq. (20).

The corresponding free energy density is plotted in Fig. 8. The correct
behavior of the free energy density for $BR^{2}\to\infty$, 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.

###  III.3 Quark contribution

To complete the calculation for full QCD with massless quarks, we have to
study the quark contribution to the effective potential

$\delta U^{\mathrm{q}}(B,R)=-N_{f}\mathrm{Tr}\ln\frac{i\not D}{\left.i\not
D\right|_{B=0}}=-N_{f}\sum_{k,j,n}\mathrm{Tr}\ln\frac{i\lambda_{kjn}(B,R)}{i\lambda_{kjn}(0,R)}=N_{f}\left.\frac{d}{ds}\zeta^{\mathrm{q}}(s)\right|_{s=0},$

where $N_{f}$ is number of quark flavors, $\lambda(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).

Figure 9: Dependence of quark eigenvalues on $BR^{2}$. At $BR^{2}\to\infty$,
infinite number of zero modes emerges. Normal modes come in positive-negative
pairs $\lambda^{\pm}$. Note that $\lambda^{+}\not{=}-\lambda^{-}$ for finite
$BR^{2}$, but the equality is restored in the limit $BR^{2}\to\infty$. An
infinite number of quasi-zero modes are not chiral for finite $R$ and become
chiral zero-modes for asymptotically large $BR^{2}$.

Zeta function $\zeta^{\mathrm{q}}(s)$ can be split into two parts [44]

$\displaystyle\zeta^{\mathrm{q}}(s)=\cos(\pi
s)\zeta_{\not{D}^{2}}\left(\frac{s}{2}\right)-i\sin(\pi s)\eta(s),$ (29)
$\displaystyle\left.\frac{d}{ds}\zeta^{\mathrm{q}}(s)\right|_{s=0}=\frac{1}{2}\zeta_{\not{D}^{2}}^{\prime}(0)-i\pi\eta(0).$

For the purpose of the present study we need only $\zeta_{\not{D}^{2}}$.
Parity-odd term $\eta(s)$ contributes to the imaginary part of the effective
potential and can be related to the $U_{\rm A}(1)$ anomalous breakdown [22].
Parity-even part $\zeta_{\not{D}^{2}}$ can be written in the following form:

$\displaystyle\zeta_{\not{D}^{2}}(s)=\mathrm{Tr}\frac{\sin\pi
s}{\pi}\sum_{k=0}^{\infty}\sum_{j_{3}=-\frac{k+1}{2}}^{\frac{k+1}{2}}(k+1)\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\Psi^{\mathrm{q}}(k+1,j_{3},t,\hat{B},R),$
(30)
$\displaystyle\Psi^{\mathrm{q}}(k+1,j,t,\hat{B},R)=\log\frac{A(-it,k,j_{3},\hat{B},R)A(it,k,j_{3},\hat{B},R)}{A(-it,k,j_{3},0,R)A(it,k,j_{3},0,R)},$

where $A(\lambda,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

$\displaystyle\delta U^{\mathrm{q}}(B,R)$ $\displaystyle=$ $\displaystyle
N_{f}\left\\{-\frac{1}{2}\mathrm{Tr}\sum_{k=1}^{\infty}(k+1)\left[\sum_{j=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\text{q}}(k,j,0,\hat{B},R)-\hat{B}^{2}R^{4}\frac{1}{12}\left(1-\frac{2}{k^{2}}\right)\right]\right.$
(31) $\displaystyle+\left.\frac{B^{2}R^{4}}{144}(5+6\gamma-6\log
2+\pi^{2}+6\log
R)-B^{4}R^{8}\frac{1}{30720}\vphantom{\sum_{\frac{1}{1}}^{\frac{1}{1}}}\right\\}.$

Figure 10: Dependence of one-loop quark contribution to the free energy
density on $R$ and $B$, $N_{f}=1$. All quantities are given in units of
renormalization scale $\mu$, Eq. (20). The strong-field and large-size regimes
are incorrect.

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

$\displaystyle\delta U_{\varkappa}^{\mathrm{q}}(B,R)$ $\displaystyle=$
$\displaystyle
N_{f}\left\\{-\frac{1}{2}\mathrm{Tr}\sum_{k=1}^{\infty}(k+1)\left[\sum_{j_{3}=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\text{q}}(k,j,0,\hat{B},R)-\hat{B}^{2}R^{4}\frac{1}{12}\left(1-\frac{2}{k^{2}}\right)\right]\right.$
(32) $\displaystyle\left.+\frac{B^{2}R^{4}}{144}(5+6\gamma-6\log
2+\pi^{2}+6\log
R)-B^{4}R^{8}\frac{1}{30720}\vphantom{\sum_{\frac{1}{1}}^{\frac{1}{1}}}\right\\}$
$\displaystyle-\frac{N_{f}}{2}\mathrm{Tr}\sum_{k=0}^{\infty}\left[(k+1)\log\frac{\lambda_{k\frac{k}{2}0}^{2}(B,R)+\varkappa(BR^{2})B}{\lambda_{k\frac{k}{2}0}^{2}(B,R)}+(k+1)\log\frac{\lambda_{k\frac{-k}{2}0}^{2}(B,R)-\varkappa(BR^{2})B}{\lambda_{k\frac{-k}{2}0}^{2}(B,R)}\right],$

similarly to Eq. (28). The effective potential is shown in Fig. 11.

Figure 11: Improved quark contribution to the free energy density $R$ and $B$,
$N_{f}=1$, in units of renormalization scale $\mu$, Eq. (20). Function
$\varkappa$ is given by Eq. (27) ($z_{0}=100,a=10$). correct strong-field and
large-size limits are restored.

Combining ghost, improved gluon and quark contributions, one finds

$\displaystyle\delta U_{\varkappa}(B,R)=\delta
U_{\varkappa}^{\mathrm{gl}}(B,R)+\delta U^{\mathrm{gh}}(B,R)+\delta
U_{\varkappa}^{\mathrm{q}}(B,R)=$
$\displaystyle=4\sum_{k=1}^{\infty}k\left[\sum_{m=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\text{gl}}(k,m,0,\frac{\sqrt{3}B}{2},R)-\frac{3}{4}B^{2}R^{4}\frac{1}{24}\left(1-\frac{1}{k}-\frac{5}{k^{2}}\right)\right]$
$\displaystyle+4\sum_{k=0}^{\infty}\left[(k+1)\log\frac{\lambda_{\mathrm{gl},k\frac{k}{2}0}^{2}(B,R)+\varkappa(BR^{2})B}{\lambda_{\mathrm{gl},k\frac{k}{2}0}^{2}(B,R)}+(k+1)\log\frac{\lambda_{\mathrm{gl},k\frac{-k}{2}0}^{2}(B,R)-\varkappa(BR^{2})B}{\lambda_{\mathrm{gl},k\frac{-k}{2}0}^{2}(B,R)}\right]$
$\displaystyle-4\sum_{k=1}^{\infty}k\left[\sum_{m=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\text{gh}}(k,m,0,\frac{\sqrt{3}B}{2},R)-\frac{3}{4}B^{2}R^{4}\frac{1}{48}\left(1-\frac{1}{k}+\frac{1}{k^{2}}\right)\right]$
$\displaystyle-\frac{N_{f}}{2}\mathrm{Tr}\sum_{k=1}^{\infty}(k+1)\left[\sum_{j_{3}=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\text{q}}(k,j,0,\hat{B},R)-\hat{B}^{2}R^{4}\frac{1}{12}\left(1-\frac{2}{k^{2}}\right)\right]$
$\displaystyle-\frac{N_{f}}{2}\mathrm{Tr}\sum_{k=0}^{\infty}\left[(k+1)\log\frac{\lambda_{\mathrm{q},k\frac{k}{2}0}^{2}(B,R)+\varkappa(BR^{2})B}{\lambda_{\mathrm{q},k\frac{k}{2}0}^{2}(B,R)}+(k+1)\log\frac{\lambda_{\mathrm{q},k\frac{-k}{2}0}^{2}(B,R)-\varkappa(BR^{2})B}{\lambda_{\mathrm{q},k\frac{-k}{2}0}^{2}(B,R)}\right]$
$\displaystyle-\frac{B^{2}R^{4}}{48}(29+33\gamma-33\log 2+33\log
R)+\frac{3B^{4}R^{8}}{10120}+N_{f}\left\\{\frac{B^{2}R^{4}}{144}(5+6\gamma-6\log
2+\pi^{2}+6\log
R)-\frac{B^{4}R^{8}}{30720}\vphantom{\sum_{\frac{1}{1}}^{\frac{1}{1}}}\right\\}.$
(33)

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.

Figure 12: Total quantum correction to the free energy density
$F_{\varkappa}(B,R)$ given by Eq. (33) with $\varkappa$ given by Eq. (27)
($z_{0}=30,a=5$). Left-hand side plot is for pure gluodynamics ($N_{f}=0$),
right-hand side plot is for full QCD with two massless quark flavors
($N_{f}=2$), in units of renormalization scale $\mu$.

###  III.4 One-loop beta function

The effective action $U=U^{\rm cl}+\delta U_{\varkappa}$ should not depend on
renormalization scale $\mu$ [45, 46], that is

$\mu\frac{d}{d\mu}U=0,$ (34)

with classical action ($V_{R}$ is four-dimensional ball of radius $R$)

$U^{\textrm{cl}}=\int_{V_{R}}d^{4}x\frac{1}{4g^{2}}\left(B^{a}_{\mu\nu}\right)^{2}=\frac{\pi^{2}B^{2}R^{4}}{2g^{2}},$

and $\delta U_{\varkappa}$ given by Eq.(33). Since only terms containing $\log
R$ contribute to the Eq.(34), one obtains

$\mu\frac{d}{d\mu}\left[\frac{\pi^{2}B^{2}R^{4}}{2g^{2}}+\frac{B^{2}R^{4}}{16}\left(-11+\frac{2}{3}N_{f}\right)\log
R\right]=0,$

that is equivalent to the equation

$\mu\frac{d}{d\mu}g=-\frac{g^{3}}{16\pi^{2}}\left(11-\frac{2}{3}N_{f}\right),$

exposing the correct one-loop $\beta$-function of QCD with $N_{\rm 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_{\rm 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].

###### Acknowledgements.

We acknowledge useful discussions with Michael Bordag, Irina Pirozhenko,
Victor Braguta and Artem Roenko.

##  Appendix A Zeta function for ghost and gluon fields

###  A.1 Ghosts

We start with the expression for $\zeta^{\text{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

$\displaystyle
M\left(\frac{m+1}{2}+\frac{k}{2}-j+n+\frac{k^{2}t^{2}}{4},1+k+m,z\right)\sim$
$\displaystyle\frac{2^{k+m}\Gamma(k+m+1)}{k^{k+m}\sqrt{2\pi
k}(1+t^{2}z)^{\frac{1}{4}}}\exp\left(k\sqrt{1+t^{2}z}-(k+m)\log\left[\sqrt{1+t^{2}z}+1\right]-\frac{2j}{k}\frac{z}{\sqrt{1+t^{2}z}+1}+\frac{z}{2}\right)\sum_{i=0}^{\infty}\frac{A_{i}(z)}{k^{i}}$
(35)

for fixed $z\geqslant 0$ and arbitrary constants $m,n$. The coefficients
$A_{i}(z)$ are found with the help of recursion relation ($i\geqslant
0,A_{0}=1$)

$A_{i+1}(z)=-\frac{1}{2}\phi-\frac{1}{2}\frac{dA_{i}}{dz}\frac{dz}{d\xi}+\int
dz\frac{d\xi}{dz}\left(\psi+\frac{1}{2}\frac{d\phi}{dz}\frac{dz}{d\xi}-\frac{1}{4}\phi^{2}\right)A_{i},$

where

$\frac{d\xi}{dz}=\frac{1}{2z}\sqrt{1+t^{2}z},\quad\phi=\frac{2}{1+t^{2}z}\left(m-\frac{2j}{k}\right),\quad\psi=\frac{m^{2}+4nz}{1+t^{2}z}+\frac{4t^{2}z+4z^{2}-t^{4}z^{2}+8t^{2}z^{3}+4t^{4}z^{4}}{(1+t^{2}z)^{3}}.$

The constants of integration are fixed by the requirement

$\lim_{t\to\infty}A_{i}(z)=0,\quad i\geqslant 1.$

Asymptotic expansion for modified Bessel function $I$ is given by (see [50],
10.41.3, 10.41.7, 10.41.9)

$I_{k}(kt)=\frac{1}{\sqrt{2\pi
k}(1+t^{2})^{\frac{1}{4}}}\exp\left(k\sqrt{1+t^{2}}+k\log\frac{t}{1+\sqrt{1+t^{2}}}\right)\sum_{i=0}^{\infty}\frac{U_{i}(p)}{k^{i}},$
(36)

where

$U_{k+1}(p)=\frac{1}{2}p^{2}(1-p^{2})U_{k}^{\prime}(p)+\frac{1}{8}\int_{0}^{p}dp^{\prime}(1-5p^{\prime
2})U_{k}(p^{\prime}),\quad p=\sqrt{1+t^{2}}.$

Now, we substitute formulas (35) and (36) into $\Psi^{\mathrm{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_{i}^{\mathrm{gh}}$

$\displaystyle u_{0}^{\mathrm{gh}}(t,B,R)$ $\displaystyle=$ $\displaystyle
B^{2}R^{4}\frac{R^{4}t^{4}-2R^{2}t^{2}-4+4\left(1+R^{2}t^{2}\right)^{1/2}}{24R^{4}t^{4}\left(1+R^{2}t^{2}\right)^{1/2}},$
$\displaystyle u_{1}^{\mathrm{gh}}(t,B,R)$ $\displaystyle=$
$\displaystyle-B^{2}R^{4}\frac{1+3R^{2}t^{2}}{48\left(1+R^{2}t^{2}\right)},$
$\displaystyle u_{2}^{\mathrm{gh}}(t,B,R)$ $\displaystyle=$
$\displaystyle\frac{B^{2}R^{4}}{192\left(1+R^{2}t^{2}\right)^{1/2}}\left[\frac{16}{R^{2}t^{2}}+\frac{R^{2}t^{2}(32+7R^{2}t^{2})}{(1+R^{2}t^{2})^{3}}+32\frac{1-\left(1+R^{2}t^{2}\right)^{1/2}}{R^{4}t^{4}}\right]$
$\displaystyle+\frac{B^{4}R^{8}}{1920\left(1+R^{2}t^{2}\right)^{1/2}}\left[-\frac{128}{R^{6}t^{6}}+\frac{32}{R^{4}t^{4}}-\frac{16}{R^{2}t^{2}}+\frac{10+13R^{2}t^{2}}{(1+R^{2}t^{2})^{2}}-256\frac{1-\left(1+R^{2}t^{2}\right)^{1/2}}{R^{8}t^{8}}\right].$

The sums over $k$ are calculated for $\mathfrak{R}s>1$ and analytically
continued in terms of Riemann $\zeta$ function

$\sum_{k=1}^{\infty}\frac{k^{1-2s}}{k^{i}}=\zeta(i-1+2s).$

to the strip $0<\mathfrak{R}s<1$ where integrals over $t$ converge at $t\to
0$. One obtains

$\displaystyle\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{gh}}(t,B,R)}{k^{i}}$
$\displaystyle=$ $\displaystyle
B^{2}R^{4+2s}\left\\{\vphantom{\frac{1}{1}}\right.\zeta(-1+2s)\frac{-\Gamma(2-s)\Gamma(1/2+s)}{24\sqrt{\pi}(2+s)}+\zeta(2s)\frac{\pi
s(1-2s)}{48\sin\pi s}$
$\displaystyle\left.+\zeta(1+2s)\frac{\left(6-s(2+s)(11+5s)\right)\Gamma(1-s)\Gamma(3/2+s)}{72\sqrt{\pi}(2+s)}\right\\}$
$\displaystyle-B^{4}R^{8+2s}\zeta(1+2s)\frac{s\Gamma(2-s)\Gamma(3/2+s)}{480\sqrt{\pi}(4+s)}.$

The expansion of counterterms in powers of $s$ around $s=0$ is

$\displaystyle\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{gh}}(t,B,R)}{k^{i}}$
$\displaystyle=$ $\displaystyle
B^{2}R^{4}\left[-\frac{1}{96}+\frac{1}{288}(4-6\gamma+\log 64-6\log
R)s\right]+B^{4}R^{8}\left[-\frac{1}{7680}s\right]$ $\displaystyle+O(s^{2}),$

where $\gamma=0.5772156649\dots$ is Euler’s constant. For $su(3)$ and
$\breve{n}$ given by (7)

$\mathrm{Tr}\breve{n}^{2}=3,\quad\mathrm{Tr}\breve{n}^{4}=4\frac{9}{16}=\frac{9}{4}.$

Since the whole spectrum is invariant with respect to $B\to-B$, trace over
color leads to factor 4. Evaluating the derivative of $\zeta^{\mathrm{gh}}(s)$
with respect to $s$ at $s=0$ one arrives at the expression (22).

###  A.2 Gluons

With Dirichlet boundary condition for color-charged modes the whole set of
eigenvalues is determined by the equations [21]

$M\left(\frac{k}{2}+1-m\pm
1-\frac{\lambda^{2}}{2n_{a}B},k+2,\frac{n_{a}BR^{2}}{2}\right)=0,\quad
k=0,1,2,\dots,\quad m=-\frac{k}{2},-\frac{k}{2}+1,\dots,\frac{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

$\displaystyle\zeta^{\text{gl}}(s)$ $\displaystyle=$
$\displaystyle\mathrm{Tr}\frac{\sin\pi
s}{\pi}\sum_{k=0}^{\infty}\sum_{m=-\frac{k}{2}}^{\frac{k}{2}}2(k+1)\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\Psi^{\mathrm{gl}}(k+1,m,t,\breve{B},R)$
$\displaystyle=$ $\displaystyle 2\mathrm{Tr}\left\\{\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\left[\sum_{m=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\mathrm{gl}}(k,m,kt,\breve{B},R)-\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{gl}}(t,\breve{B},R)}{k^{i}}\right]\right.$
$\displaystyle\left.+\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{gl}}(t,\breve{B},R)}{k^{i}}\right\\},$
$\displaystyle\Psi^{\mathrm{gl}}(k,m,t,B,R)$ $\displaystyle=$
$\displaystyle\log\frac{\exp\left(-\frac{BR^{2}}{2}\right)M\left(\frac{k+1}{2}-m+1+\frac{t^{2}}{2B},k+1,\frac{BR^{2}}{2}\right)M\left(\frac{k+1}{2}-m-1+\frac{t^{2}}{2B},k+1,\frac{BR^{2}}{2}\right)}{\left(k!\left(\frac{tR}{2}\right)^{-k}I_{k}(tR)\right)^{2}}.$

We substitute formulas (35) and (36) into $\Psi^{\mathrm{gl}}$, expand it in
powers of $k^{-1}$, sum over $m$ and compute derivative with respect to $t$.
Coefficients of asymptotic expansion $u_{i}^{\mathrm{gl}}$ are given by

$\displaystyle u_{0}^{\mathrm{gl}}(t,B,R)$ $\displaystyle=$ $\displaystyle
B^{2}R^{4}\frac{R^{4}t^{4}-2R^{2}t^{2}-4+4\left(1+R^{2}t^{2}\right)^{1/2}}{12R^{4}t^{4}\left(1+R^{2}t^{2}\right)^{1/2}},$
$\displaystyle u_{1}^{\mathrm{gl}}(t,B,R)$ $\displaystyle=$
$\displaystyle-B^{2}R^{4}\frac{1+3R^{2}t^{2}}{24\left(1+R^{2}t^{2}\right)},$
$\displaystyle u_{2}^{\mathrm{gl}}(t,B,R)$ $\displaystyle=$
$\displaystyle\frac{B^{2}R^{4}}{96\left(1+R^{2}t^{2}\right)^{1/2}}\left[-\frac{80}{R^{2}t^{2}}+\frac{R^{2}t^{2}(32+7R^{2}t^{2})}{(1+R^{2}t^{2})^{3}}-160\frac{1-\left(1+R^{2}t^{2}\right)^{1/2}}{R^{4}t^{4}}\right]$
$\displaystyle+\frac{B^{4}R^{8}}{960\left(1+R^{2}t^{2}\right)^{1/2}}\left[-\frac{128}{R^{6}t^{6}}+\frac{32}{R^{4}t^{4}}-\frac{16}{R^{2}t^{2}}+\frac{10+13R^{2}t^{2}}{(1+R^{2}t^{2})^{2}}-256\frac{1-\left(1+R^{2}t^{2}\right)^{1/2}}{R^{8}t^{8}}\right].$

One obtains

$\displaystyle\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{gl}}(t,B,R)}{k^{i}}$
$\displaystyle=$ $\displaystyle
B^{2}R^{4+2s}\left\\{\vphantom{\frac{1}{1}}\right.\zeta(-1+2s)\frac{-\Gamma(2-s)\Gamma(1/2+s)}{12\sqrt{\pi}(2+s)}+\zeta(2s)\frac{\pi
s(1-2s)}{24\sin\pi s}$
$\displaystyle+\left.\zeta(1+2s)\frac{\left(30+s(2+s)(11+5s)\right)\Gamma(1-s)\Gamma(3/2+s)}{36\sqrt{\pi}(2+s)}\right\\}$
$\displaystyle-B^{4}R^{8+2s}\zeta(1+2s)\frac{s\Gamma(2-s)\Gamma(3/2+s)}{240\sqrt{\pi}(4+s)}.$

Now, counterterms can be expanded in powers of $s$

$\displaystyle\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{gl}}(t,B,R)}{k^{i}}$
$\displaystyle=$ $\displaystyle
B^{2}R^{4}\left[\frac{5}{48}+\frac{1}{144}(31+30\gamma-30\log 2+30\log
R)s\right]+B^{4}R^{8}\left[-\frac{1}{3840}s\right]$ $\displaystyle+O(s^{2}).$

Finally,

$\displaystyle\delta
U^{\mathrm{gl}}(B,R)=-\frac{1}{2}\left.\frac{d}{ds}\zeta^{\text{gl}}(s)\right|_{s=0}$
$\displaystyle=$ $\displaystyle
4\sum_{k=1}^{\infty}k\left[\sum_{m=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\text{gl}}(k,m,0,\frac{\sqrt{3}B}{2},R)-\frac{3}{4}B^{2}R^{4}\frac{1}{24}\left(1-\frac{1}{k}-\frac{5}{k^{2}}\right)\right]$
$\displaystyle-\frac{B^{2}R^{4}}{48}(31+30\gamma-30\log 2+30\log
R)+\frac{3B^{4}R^{8}}{5120}.$

###  A.3 Contribution of quasi-zero modes

In this section, we calculate the contribution of gluon quasi-zero modes
$\lambda_{\uparrow k\frac{k}{2}0}$ to the effective potential. These modes
correspond to the smallest-magnitude solutions to the equations

$M\left(-\frac{\lambda^{2}}{2|\breve{B}|},k+2,\frac{|\breve{B}|R^{2}}{2}\right)=0,\quad
k=0,1,2,\dots$ (38)

that reduce to

$\left(\frac{\lambda R}{2}\right)^{-k-1}J_{k+1}(\lambda R)=0,\quad
k=0,1,\dots$

at $B\to 0$. The contribution of gluon quasi-zero modes to the effective
potential can be expressed as

$\displaystyle\delta
U^{\textrm{gl}(0)}(B,R)=\frac{1}{2}\textrm{Tr}\log\frac{\lambda_{\uparrow
k\frac{k}{2}0}^{2}(B,R)}{\lambda_{k0}^{2}(0,R)}=\left.-\frac{1}{2}\frac{d}{ds}\zeta^{\textrm{gl}(0)}(s)\right|_{s=0},$
$\displaystyle\zeta^{\textrm{gl}(0)}(s)=2\mathrm{Tr}\sum_{k=0}^{\infty}(k+1)\left(\lambda_{\uparrow
k\frac{k}{2}0}^{-2s}(B,R)-\lambda_{k0}^{-2s}(0,R)\right),$ (39)

where factor $2$ in the definition of $\zeta^{\textrm{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\to\frac{\sqrt{3}}{2}B$ for adjoint representation of
$su(3)$)

To continue $\zeta^{\textrm{gl}(0)}(s)$ to $s\to 0$, we add and subtract
several terms of asymptotic expansion in $k$ found with the help of formula
(53)

$\displaystyle\zeta^{\textrm{gl}(0)}(s)$ $\displaystyle=$ $\displaystyle
8\sum_{k=0}^{\infty}\left[(k+1)\left(\lambda_{\uparrow
k\frac{k}{2}0}^{-2s}(B,R)-\lambda_{k0}^{-2s}(0,R)\right)\right.$
$\displaystyle\left.-\left(k+1\right)^{-2s}R^{2s}\left(BR^{2}s-2\alpha
BR^{2}s(1+s)(k+1)^{-2/3}+BR^{2}\frac{s}{4}\left(4+BR^{2}(1+2s)\right)(k+1)^{-1}\right)\right]$
$\displaystyle+8\sum_{k=0}^{\infty}\left(k+1\right)^{-2s}R^{2s}\left(BR^{2}s-2\alpha
BR^{2}s(1+s)(k+1)^{-2/3}+BR^{2}\frac{s}{4}\left(4+BR^{2}(1+2s)\right)(k+1)^{-1}\right).$

The first sum is an analytic function for $\mathfrak{R}s>0$. The second sum is
evaluated for $\mathfrak{R}s>1/2$ and analytically continued to $s\to 0$:

$\displaystyle\sum_{k=0}^{\infty}\left(k+1\right)^{-2s}R^{2s}\left(BR^{2}s-2\alpha
BR^{2}s(1+s)(k+1)^{-2/3}+BR^{2}\frac{s}{4}\left(4+BR^{2}(1+2s)\right)(k+1)^{-1}\right)$
$\displaystyle=BR^{2+2s}s\left(\zeta(2s)-2\alpha(1+s)\zeta(2/3+2s)+\zeta(1+2s)\right)+B^{2}R^{4+2s}\frac{s}{4}(1+2s)\zeta(1+2s).$

Finally, we obtain

$\displaystyle\delta U^{\textrm{gl}(0)}(B,R)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left.\frac{d}{ds}\zeta^{\textrm{gl}(0)}(s)\right|_{s=0}$
$\displaystyle=$
$\displaystyle-4\sum_{k=0}^{\infty}\left[-(k+1)\log\frac{\lambda_{\uparrow
k\frac{k}{2}0}^{2}(B,R)}{\lambda_{k0}^{2}(0,R)}-BR^{2}\left(1-2\alpha(k+1)^{-2/3}+(k+1)^{-1}\right)-\frac{B^{2}R^{4}}{4(k+1)}\right]$
$\displaystyle-2BR^{2}\left[-1+2\gamma+2\log
R-4\alpha\zeta\left(\frac{2}{3}\right)\right]-B^{2}R^{4}(1+\gamma+\log R).$

###  A.4 Contribution of quasi-zero modes with the effective “mass”

If one includes the effective “mass” $\varkappa$ for quasi-zero modes in the
considerations of the previous section, the formulas become

$\displaystyle\delta
U_{\varkappa}^{\textrm{gl}(0)}(B,R)=\frac{1}{2}\textrm{Tr}\log\frac{\lambda_{\uparrow
k\frac{k}{2}0}^{2}(B,R)+\varkappa(BR^{2})B}{\lambda_{k0}^{2}(0,R)}=\left.-\frac{1}{2}\frac{d}{ds}\zeta^{\textrm{gl}(0)}(s)\right|_{s=0},$
$\displaystyle\zeta_{\varkappa}^{\textrm{gl}(0)}(s)=2\mathrm{Tr}\sum_{k=0}^{\infty}(k+1)\left(\left(\lambda_{\uparrow
k\frac{k}{2}0}^{2}(B,R)+\varkappa(BR^{2})B\right)^{-s}-\lambda_{k0}^{-2s}(0,R)\right).$
(40)

In analogy to the previous section,

$\displaystyle\zeta_{\varkappa}^{\textrm{gl}(0)}(s)=$
$\displaystyle=8\sum_{k=0}^{\infty}\left[(k+1)\left(\left(\lambda_{\uparrow
k\frac{k}{2}0}^{2}(B,R)+\varkappa(BR^{2})B\right)^{-s}-\lambda_{k0}^{-2s}(0,R)\right)\right.$
$\displaystyle\left.-\left(k+1\right)^{-2s}R^{2s}\left(BR^{2}s-2\alpha
BR^{2}s(1+s)(k+1)^{-2/3}+BR^{2}\frac{s}{4}\left(4+BR^{2}(1+2s)-4\varkappa\left(BR^{2}\right)\right)(k+1)^{-1}\right)\right]$
$\displaystyle+8\sum_{k=0}^{\infty}\left(k+1\right)^{-2s}R^{2s}\left(BR^{2}s-2\alpha
BR^{2}s(1+s)(k+1)^{-2/3}+BR^{2}\frac{s}{4}\left(4+BR^{2}(1+2s)-4\varkappa\left(BR^{2}\right)\right)(k+1)^{-1}\right)$

The first sum is an analytic function for $\mathfrak{R}s>0$. The second sum is
evaluated for $\mathfrak{R}s>1/2$ and analytically continued to $s\to 0$:

$\displaystyle\sum_{k=0}^{\infty}\left(k+1\right)^{-2s}R^{2s}\left(BR^{2}s-2\alpha
BR^{2}s(1+s)(k+1)^{-2/3}+BR^{2}\frac{s}{4}\left(4+BR^{2}(1+2s)-4\varkappa\left(BR^{2}\right)\right)(k+1)^{-1}\right)$
$\displaystyle=BR^{2+2s}s\left(\zeta(2s)-2\alpha(1+s)\zeta(2/3+2s)+\left(1-\varkappa\left(BR^{2}\right)\right)\zeta(1+2s)\right)+B^{2}R^{4+2s}\frac{s}{4}(1+2s)\zeta(1+2s).$

Finally,

$\displaystyle\delta U_{\varkappa}^{\textrm{gl}(0)}(B,R)$ $\displaystyle=$
$\displaystyle\left.-\frac{1}{2}\frac{d}{ds}\zeta_{\varkappa}^{\textrm{gl}(0)}(s)\right|_{s=0}$
$\displaystyle=$
$\displaystyle-4\sum_{k=0}^{\infty}\left[-(k+1)\log\frac{\lambda_{\uparrow
k\frac{k}{2}0}^{2}(B,R)+\varkappa(BR^{2})B}{\lambda_{k0}^{2}(0,R)}\right.$
$\displaystyle\left.-BR^{2}\left(1-2\alpha(k+1)^{-2/3}+\left(1-\varkappa\left(BR^{2}\right)\right)(k+1)^{-1}\right)-\frac{B^{2}R^{4}}{4(k+1)}\right]$
$\displaystyle-2BR^{2}\left[-1+2\left(\gamma+\log
R\right)\left(1-\varkappa\left(BR^{2}\right)\right)-4\alpha\zeta\left(\frac{2}{3}\right)\right]-B^{2}R^{4}(1+\gamma+\log
R).$

###  A.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

$\zeta_{\varkappa}^{\mathrm{gl}}(s)=\zeta^{\mathrm{gl}}(s)-\zeta^{\textrm{gl}(0)}(s)+\zeta_{\varkappa}^{\textrm{gl}(0)}(s),$

where zeta functions in right-hand side are given by Eqs. (A.2),(39) and (40).
The corresponding effective potential is given by

$\displaystyle\delta U_{\varkappa}^{\mathrm{gl}}(B,R)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left.\frac{d}{ds}\zeta_{\varkappa}^{\text{gl}}(s)\right|_{s=0}$
$\displaystyle=$ $\displaystyle
4\sum_{k=1}^{\infty}k\left[\sum_{m=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\text{gl}}(k,m,0,\frac{\sqrt{3}B}{2},R)-\frac{3}{4}B^{2}R^{4}\frac{1}{24}\left(1-\frac{1}{k}-\frac{5}{k^{2}}\right)\right]$
$\displaystyle-\frac{B^{2}R^{4}}{48}(31+30\gamma-30\log 2+30\log
R)+\frac{3B^{4}R^{8}}{5120}$
$\displaystyle+4\sum_{k=0}^{\infty}\left[(k+1)\log\frac{\lambda_{\uparrow
k\frac{k}{2}0}^{2}(B,R)+\varkappa(BR^{2})B}{\lambda_{k\frac{k}{2}0}^{2}(B,R)}+BR^{2}\varkappa\left(BR^{2}\right)(k+1)^{-1}\right]+4BR^{2}\varkappa(BR^{2})(\gamma+\log
R).$

Thus obtained effective potential is not an even function of $B$. To restore
invariance under $B\to-B$, one adds term $\varkappa$ to contribution of modes
$\lambda_{\downarrow k,-\frac{k}{2},0}^{2}(B,R)$ (these modes become quasi-
zero when field $B$ changes direction to the opposite). The function
$\varkappa(z)$ should also be even function (possibly constant). After these
steps one obtains formula (28) for the effective potential.

##  Appendix B Quark eigenmodes

The equation for eigenvalues of quarks field in the presence of homogeneous
(anti-)self-dual gluon field reads

$\not{D}\psi=\lambda\psi.$ (41)

Here

$\displaystyle
D_{\mu}=\partial_{\mu}-i\hat{B}_{\mu},\quad\hat{B}_{\mu}=\hat{n}B_{\mu},\quad
B_{\mu}=-\frac{1}{2}B_{\mu\nu}x_{\nu},\quad
B_{\mu\nu}=\pm\frac{1}{2}\varepsilon_{\mu\nu\alpha\beta}B_{\alpha\beta},$
$\displaystyle\hat{n}=\cos\xi t^{3}+\sin\xi
t^{8},\quad\xi=\frac{\pi}{6}+\frac{\pi}{3}k,k=0,1,\dots,5,$

$t_{i}$ are generators of $SU(3)$ in fundamental representation. The field
strength tensor $B_{\mu\nu}$ may be parametrized as

$B_{ij}=\varepsilon_{ijk}B_{k},\quad B_{i4}=\pm B_{i},\quad
B=\sqrt{B_{1}^{2}+B_{2}^{2}+B_{3}^{2}},$

where “$+$” stands for self-dual field and “$-$” for anti-self-dual field. We
choose anti-Hermitian chiral representation of gamma matrices

$\left\\{\gamma_{\mu},\gamma_{\nu}\right\\}=-2\delta_{\mu\nu},\quad\gamma_{i}=\left(\begin{matrix}0&\sigma_{i}\\\
-\sigma_{i}&0\end{matrix}\right),\quad\gamma_{4}=i\left(\begin{matrix}0&\mathbf{1}\\\
\mathbf{1}&0\end{matrix}\right),\quad\gamma_{5}=\gamma_{1}\gamma_{2}\gamma_{3}\gamma_{4}=\left(\begin{matrix}\mathbf{1}&0\\\
0&\mathbf{-1}\end{matrix}\right).$

It is convenient to introduce projectors

$\Sigma_{\pm}=\frac{1}{2}\left(1\pm\frac{\Sigma_{i}B_{i}}{B}\right),\quad
P_{\pm}=\frac{1}{2}(1\pm\gamma_{5}),$

where

$\Sigma_{i}=\frac{1}{2}\varepsilon_{ijk}\sigma_{jk},\quad\sigma_{\mu\nu}=-\frac{1}{2}\varepsilon_{\mu\nu\alpha\beta}\sigma_{\alpha\beta}\gamma_{5},\quad\sigma_{ij}=\varepsilon_{ijk}\Sigma_{k},\quad\sigma_{i4}=-\gamma_{5}\Sigma_{i}.$

For self-dual or anti-self-dual field $B_{\mu\nu}$ one obtains identities

$\sigma_{\mu\nu}B_{\mu\nu}=4BP_{\mp}\left(\Sigma_{+}-\Sigma_{-}\right),\quad\gamma_{\mu}B_{\mu\nu}x_{\nu}P_{\mp}=-i\not{x}BP_{\mp}\left(\Sigma_{+}-\Sigma_{-}\right),$

where the upper sign corresponds to self-dual field, the lower sign
corresponds to anti-self-dual-field. Color matrix $\hat{n}$ is diagonal, and
for the sake of simplicity we keep the notation $\hat{n}$ for its diagonal
elements.

Acting on equation (41) with the projectors $P_{\pm}=(1\pm\gamma_{5})/2$, one
rewrites it as

$\displaystyle\not{D}P_{+}\psi$ $\displaystyle=$ $\displaystyle\lambda
P_{-}\psi,$ $\displaystyle\not{D}P_{-}\psi$ $\displaystyle=$
$\displaystyle\lambda P_{+}\psi.$

or

$\displaystyle\not{D}\psi_{+}$ $\displaystyle=$
$\displaystyle\lambda\psi_{-},$ $\displaystyle\not{D}\psi_{-}$
$\displaystyle=$ $\displaystyle\lambda\psi_{+},$

where $\psi_{\pm}=P_{\pm}\psi$. Substituting one equation into the other (this
is valid if $\lambda\neq 0$), one finds

$\displaystyle\not{D}^{2}\psi_{\mp}=\lambda^{2}\psi_{\mp},$ (42)
$\displaystyle\lambda\psi_{\pm}=\not{D}\psi_{\mp},$ (43)
$\displaystyle\psi=\psi_{+}+\psi_{-}=\left(\frac{\not{D}}{\lambda}+1\right)\psi_{\mp}.$
(44)

Only one chiral component is independent, the other one is found via Eq. (43).
The expression for $\gamma_{\mu}B_{\mu\nu}x_{\nu}P_{\pm}$
($\gamma_{\mu}B_{\mu\nu}x_{\nu}$ originates from $\not{D}$) is more complex
than $\gamma_{\mu}B_{\mu\nu}x_{\nu}P_{\mp}$, so we find $\psi_{-}$ from
equation (42) in the case of self-dual field and $\psi_{+}$ in the case of
anti-self-dual field. Here and below, if $\pm$ or $\mp$ appears alongside with
Dirac operator $\not{D}$ or its eigenmode $\psi$, the upper sign should be
taken for self-dual field and the lower sign for anti-self-dual field. $\eta$
is the normal to the surface of a sphere, $\eta^{2}=1$.

The analogue of total angular momentum in Euclidean space is [51]

$\displaystyle J_{\mu\nu}$ $\displaystyle=$ $\displaystyle
K_{\mu\nu}+S_{\mu\nu},$ $\displaystyle K_{\mu\nu}$ $\displaystyle=$
$\displaystyle-i\left(x_{\mu}\partial_{\nu}-x_{\nu}\partial_{\mu}\right),\quad
S_{\mu\nu}=\frac{i}{4}\left[\gamma_{\mu},\gamma_{\nu}\right],$ $\displaystyle
J_{i}^{\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{1}{2}\varepsilon_{ijk}J_{jk}\pm
J_{i4}\right)=\frac{1}{2}\left(\frac{1}{2}\varepsilon_{ijk}K_{jk}\pm
K_{i4}\right)+\frac{1}{2}\left(\frac{1}{2}\varepsilon_{ijk}S_{jk}\pm
S_{i4}\right)=K_{i}^{\pm}+\frac{1}{2}P_{\mp}\Sigma_{i}.$

Algebra of operators $K_{\mu\nu},S_{\mu\nu},J_{\mu\nu}$ splits into two
$so(3)$ algebras

$[K_{i}^{+},K_{j}^{+}]=i\varepsilon_{ijk}K_{k}^{+},\quad[K_{i}^{-},K_{j}^{-}]=i\varepsilon_{ijk}K_{k}^{-},\quad[K_{i}^{+},K_{j}^{-}]=0,$

and analogous relations for $J^{\pm}$. One introduces lowering and raising
operators

$\displaystyle\Sigma^{(\pm)}=\frac{1}{2}\left(\Sigma_{1}\pm
i\Sigma_{2}\right),\quad K^{(\pm)}=K_{1}\pm iK_{2},$
$\displaystyle[K_{3},K^{(\pm)}]=\pm K^{(\pm)},\quad[K^{(+)},K^{(-)}]=2K_{3}.$

We choose reference frame such that direction of field coincides with $z$
axis.

$B_{i}=\left\\{0,0,B\right\\}.$

Thus, equation (42) can be cast into the form

$\displaystyle\left\\{-\left[\frac{1}{r^{3}}\partial_{r}r^{3}\partial_{r}-\frac{4}{r^{2}}\mathbf{K}^{2}_{\mp}+2\hat{B}K_{2,1}^{z}-\frac{1}{4}\hat{B}^{2}r^{2}\right]-2\hat{B}\right\\}\psi_{\mp+}$
$\displaystyle=$ $\displaystyle\lambda^{2}\psi_{\mp+},$
$\displaystyle\left\\{-\left[\frac{1}{r^{3}}\partial_{r}r^{3}\partial_{r}-\frac{4}{r^{2}}\mathbf{K}^{2}_{\mp}+2\hat{B}K_{2,1}^{z}-\frac{1}{4}\hat{B}^{2}r^{2}\right]+2\hat{B}\right\\}\psi_{\mp-}$
$\displaystyle=$ $\displaystyle\lambda^{2}\psi_{\mp-},$

where $c$ in $\psi_{cs}$ stands for chirality, and $s$ for spin. It is
convenient to introduce basis spinors $u_{cs}$

$\displaystyle P_{\pm}u_{\pm s}=u_{\pm s},\quad P_{\pm}u_{\mp
s}=0,\quad\Sigma_{\pm}u_{c\pm}=u_{c\pm},\quad\Sigma_{\pm}u_{c\mp}=0,\quad
u^{\dagger}_{cs}u_{c^{\prime}s^{\prime}}=\delta_{cc^{\prime}}\delta_{ss^{\prime}},$
$\displaystyle\Sigma^{(\mp)}u_{c\pm}=u_{c\mp},\quad\Sigma^{(\pm)}u_{c\pm}=0.$

The solutions are

$\displaystyle\psi_{\mp+}^{km_{1}m_{2}}$ $\displaystyle=$
$\displaystyle\left(\frac{\hat{B}r^{2}}{2}\right)^{k/2}e^{-\hat{B}r^{2}/4}M\left(\frac{k}{2}+1-m_{2,1}+\frac{-2\hat{B}-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}r^{2}}{2}\right)Y_{\frac{k}{2}m_{1}m_{2}}(\varphi,\chi,\eta)u_{\mp+},$
(45) $\displaystyle\psi_{\mp-}^{km_{1}m_{2}}$ $\displaystyle=$
$\displaystyle\left(\frac{\hat{B}r^{2}}{2}\right)^{k/2}e^{-\hat{B}r^{2}/4}M\left(\frac{k}{2}+1-m_{2,1}+\frac{2\hat{B}-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}r^{2}}{2}\right)Y_{\frac{k}{2}m_{1}m_{2}}(\varphi,\chi,\eta)u_{\mp-}.$
(46)

Where $Y_{\frac{k}{2}m_{1}m_{2}}$ are spherical harmonics in four-dimensional
Euclidean space [51].

The solutions of (41) may be characterized by eigenvalues of independent
operators $\mathbf{J}^{\mp 2},\ J_{3}^{\mp},\ J_{3}^{\pm},\ L$ (see [52, 53,
54]). Here

$L=i\gamma_{\mu}\gamma_{5}B_{\mu\nu}D_{\nu},$

and upper or lower signs should be taken for self-dual field and anti-self-
dual field, correspondingly. The solutions to the eigenvalue problem (41)

$\left(\frac{\not D}{\lambda}+1\right)\psi_{cs}^{km_{1}m_{2}},$

where $\psi_{cs}^{km_{1}m_{2}}$ ($c=\mp,\ s=\pm$) are solutions of Eq. (42)
given by Eqs. (45) and (46), diagonalize all these operators.

Now, we impose bag boundary condition (11). Substituting Eq. (44) in Eq. (11)
and acting on it with operators $P_{\pm}$, we find that boundary condition
reduces to

$ie^{-i\alpha\gamma_{5}}\not{\eta}\frac{\not
D}{\lambda}\psi_{\mp}(x)=\psi_{\mp}(x),\quad x^{2}=R^{2},$

where $\psi_{\mp}$ are solutions of Eq. (42). 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}\to m_{1},\alpha\to-\alpha$). 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}=\frac{k}{2}+\frac{1}{2}$, one solution for minimal
value $j_{3}=-\frac{k}{2}-\frac{1}{2}$ and two solutions for each intermediate
$j_{3}=-\frac{k}{2}+\frac{1}{2},\dots,\frac{k}{2}-\frac{1}{2}$ (for $k>0$).
Below we consider these three cases.

###  B.1 Intermediate values of $j_{3}$

As described above, we consider linear combination

$\psi_{-}=C_{1}f(k,m_{2}+1,r)Y_{\frac{k}{2},m_{1},m_{2}}u_{-+}+C_{2}f(k,m_{2},r)Y_{\frac{k}{2},m_{1},m_{2}+1}(\varphi,\chi,\eta)u_{--}$

for

$m_{2}=j_{3}-\frac{1}{2}=-\frac{k}{2},-\frac{k}{2}+1,\dots,\frac{k}{2}-1,\quad
k=1,2,\dots,$

where we introduced notation

$f(k,m,r)=\left(\frac{\hat{B}r^{2}}{2}\right)^{k/2}e^{-\hat{B}r^{2}/4}M\left(\frac{k}{2}+1-m+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}r^{2}}{2}\right).$

The boundary condition reads

$\displaystyle\frac{1}{\lambda}\left\\{C_{1}\left[-\partial_{r}-\frac{\hat{B}R}{2}+\frac{2}{R}m_{2}\right]f(k,m_{2}+1,R)+C_{2}\sqrt{\frac{k}{2}\left(\frac{k}{2}+1\right)-m_{2}(m_{2}+1)}\frac{2}{R}f(k,m_{2},R)\right\\}=$
$\displaystyle=-ie^{-i\alpha}C_{1}f(k,m_{2}+1,R),$
$\displaystyle\frac{1}{\lambda}\left\\{C_{2}\left[-\partial_{r}+\frac{\hat{B}R}{2}-\frac{2}{R}(m_{2}+1)\right]f(k,m_{2},R)+C_{1}\sqrt{\frac{k}{2}\left(\frac{k}{2}+1\right)-m_{2}(m_{2}+1)}\frac{2}{R}f(k,m_{2}+1,R)\right\\}=$
$\displaystyle=-ie^{-i\alpha}C_{2}f(k,m_{2},R).$

Using identities (see [55], 13.4.10,13.4.11)

$\displaystyle aM(a+1,b,z)=aM(a,b,z)+zM^{\prime}(a,b,z),$
$\displaystyle(b-a)M(a-1,b,z)=(b-a-z)M(a,b,z)+zM^{\prime}(a,b,z),$

we bring these equations to the form

$\displaystyle\frac{C_{1}}{\lambda}\left[\left(-\frac{\lambda^{2}}{2\hat{B}}+ie^{-i\alpha}\frac{\lambda
R}{2}\right)M\left(\frac{k}{2}-m_{2}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)-\right.$
$\displaystyle\left.-\left(\frac{k}{2}-m_{2}+\frac{-\lambda^{2}}{2\hat{B}}\right)M\left(\frac{k}{2}+1-m_{2}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)\right]+$
$\displaystyle+\frac{C_{2}}{\lambda}\sqrt{\frac{k}{2}\left(\frac{k}{2}+1\right)-m_{2}(m_{2}+1)}M\left(\frac{k}{2}+1-m_{2}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)=0,$
(47)
$\displaystyle\frac{C_{2}}{\lambda}\left[\left(\frac{\lambda^{2}}{2\hat{B}}+ie^{-i\alpha}\frac{\lambda
R}{2}\right)M\left(\frac{k}{2}+1-m_{2}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)-\right.$
$\displaystyle\left.-\left(\frac{k}{2}+1+m_{2}+\frac{\lambda^{2}}{2\hat{B}}\right)M\left(\frac{k}{2}-m_{2}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)\right]+$
$\displaystyle+\frac{C_{1}}{\lambda}\sqrt{\frac{k}{2}\left(\frac{k}{2}+1\right)-m_{2}(m_{2}+1)}M\left(\frac{k}{2}-m_{2}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)=0.$
(48)

This system is uniform with respect to $C_{1},C_{2}$ and has nontrivial
solution if the determinant of the system is zero. The latter requirement
leads to equation for eigenvalues:

$\displaystyle\frac{1}{\lambda^{2}}\left[\left(-\frac{\lambda^{2}}{2\hat{B}}+ie^{-i\alpha}\frac{\lambda
R}{2}\right)M\left(\frac{k}{2}+\frac{1}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)-\right.$
$\displaystyle\left.-\left(\frac{k}{2}+\frac{1}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}}\right)M\left(\frac{k}{2}+\frac{3}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)\right]\times$
$\displaystyle\times\left[\left(\frac{\lambda^{2}}{2\hat{B}}+ie^{-i\alpha}\frac{\lambda
R}{2}\right)M\left(\frac{k}{2}+\frac{3}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)-\right.$
$\displaystyle\left.-\left(\frac{k}{2}+\frac{1}{2}+j_{3}+\frac{\lambda^{2}}{2\hat{B}}\right)M\left(\frac{k}{2}+\frac{1}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)\right]-$
$\displaystyle-\frac{1}{\lambda^{2}}\left[\frac{k}{2}\left(\frac{k}{2}+1\right)-j_{3}^{2}+\frac{1}{4}\right]M\left(\frac{k}{2}+\frac{1}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)M\left(\frac{k}{2}+\frac{3}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)=0.$

where substitution $m_{2}\to j_{3}-\frac{1}{2}$ is made.

###  B.2 Maximal value of $j_{3}$

In this case there is only one function with $j_{3}=\frac{k}{2}+\frac{1}{2}$
at given $k=0,1,2,\dots$

$\psi_{-}=C_{1}f(k,m_{2}+1,r)Y_{\frac{k}{2},m_{1},m_{2}}u_{-+},\quad
m_{2}=j_{3}-\frac{1}{2}=\frac{k}{2}.$

Boundary condition in this case can be obtained from Eq. (B.1) with $C_{2}=0$.
After simplification we find

$M\left(-\frac{\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)-ie^{i\alpha}\frac{\lambda
R}{2(k+2)}M\left(1-\frac{\lambda^{2}}{2\hat{B}},k+3,\frac{\hat{B}R^{2}}{2}\right)=0.$

###  B.3 Minimal value of $j_{3}$

There is also only one function with $j_{3}=-\frac{k}{2}-\frac{1}{2}$ at given
$k=0,1,2,\dots$

$\psi_{-}=C_{2}f(k,m_{2},r)Y_{\frac{k}{2},m_{1},m_{2}+1}(\varphi,\chi,\eta)u_{--},\quad
m_{2}=j_{3}-\frac{1}{2}=-\frac{k}{2}-1.$

Boundary condition can be obtained from Eq. (B.1) 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\left(\frac{\lambda^{2}}{2\hat{B}},k+2,-\frac{\hat{B}R^{2}}{2}\right)-ie^{i\alpha}\frac{\lambda
R}{2(k+2)}M\left(1+\frac{\lambda^{2}}{2\hat{B}},k+3,-\frac{\hat{B}R^{2}}{2}\right)=0.$

##  Appendix C Zeta function for quark fields

We put chiral angle $\alpha=\pi/2$, so the Dirac operator is Hermitian and
eigenvalues are real. Zeta function $\zeta(s)$ of fermions is split into two
parts

$\displaystyle\zeta^{\mathrm{q}}(s)=\cos(\pi
s)\zeta_{\not{D}^{2}}\left(\frac{s}{2}\right)-i\sin(\pi s)\eta(s),$
$\displaystyle\zeta_{\not{D}^{2}}(s)=\mathrm{Tr}\sum_{k,j_{3},n}(k+1)\left(\frac{1}{\left|\lambda_{kj_{3}n}(\hat{B},R)\right|^{2s}}-\frac{1}{\left|\lambda_{kj_{3}n}(0,R)\right|^{2s}}\right),$
$\displaystyle\eta(s)=\mathrm{Tr}\sum_{k,j_{3},n}(k+1)\left(\frac{\mathrm{sign}{\lambda_{kj_{3}n}(\hat{B},R)}}{\left|\lambda_{kj_{3}n}(\hat{B},R)\right|^{s}}-\frac{\mathrm{sign}{\lambda_{kj_{3}n}(0,R)}}{\left|\lambda_{kj_{3}n}(0,R)\right|^{s}}\right).$

Eigenvalues $\lambda_{kjn}$ are found from equations
$A(\lambda,k,j_{3},\hat{B},R|\alpha)=0$. The set of eigenvalues for
$j_{3}=\pm\frac{k+1}{2}$ is determined by the equations

$\displaystyle
A(\lambda,k,\frac{k+1}{2},\hat{B},R)=M\left(-\frac{\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)-\frac{iR\lambda}{2(k+2)}M\left(1-\frac{\lambda^{2}}{2\hat{B}},k+3,\frac{\hat{B}R^{2}}{2}\right)=0,$
$\displaystyle
A(\lambda,k,-\frac{k+1}{2},\hat{B},R)=M\left(\frac{\lambda^{2}}{2\hat{B}},k+2,-\frac{\hat{B}R^{2}}{2}\right)-\frac{iR\lambda}{2(k+2)}M\left(1+\frac{\lambda^{2}}{2\hat{B}},k+3,-\frac{\hat{B}R^{2}}{2}\right)=0.$

At $B=0$ the equations take the form

$(k+1)!\left(\frac{\lambda R}{2}\right)^{-k-1}\left[J_{k+1}(\lambda
R)-iJ_{k+2}(\lambda R)\right]=0.$

The eigenvalues for $j_{3}=-\frac{k-1}{2},\dots,\frac{k-1}{2},k\geqslant 1$
are found from equation

$\displaystyle
A(\lambda,k,j_{3},\hat{B},R)=\exp\left(-\frac{\hat{B}R^{2}}{2}\right)\frac{1}{\lambda^{2}}\left[i\frac{\lambda
R}{2}M\left(\frac{k}{2}+\frac{1}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)-\right.$
$\displaystyle\left.-\left(\frac{k}{2}+\frac{1}{2}-j_{3}\right)M\left(\frac{k}{2}+\frac{3}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)+\frac{\lambda^{2}R^{2}}{4(k+2)}M\left(\frac{k}{2}+\frac{3}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+3,\frac{\hat{B}R^{2}}{2}\right)\right]\times$
$\displaystyle\times\left[i\frac{\lambda
R}{2}M\left(\frac{k}{2}+\frac{3}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)\right.$
$\displaystyle\left.-\left(\frac{k}{2}+\frac{1}{2}+j_{3}\right)M\left(\frac{k}{2}+\frac{1}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)+\frac{\lambda^{2}R^{2}}{4(k+2)}M\left(\frac{k}{2}+\frac{3}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+3,\frac{\hat{B}R^{2}}{2}\right)\right]$
$\displaystyle-\exp\left(-\frac{\hat{B}R^{2}}{2}\right)\frac{1}{\lambda^{2}}\left[\frac{k}{2}\left(\frac{k}{2}+1\right)-j_{3}^{2}+\frac{1}{4}\right]\times$
$\displaystyle\times
M\left(\frac{k}{2}+\frac{1}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)M\left(\frac{k}{2}+\frac{3}{2}-j_{3}+\frac{-\lambda^{2}}{2\hat{B}},k+2,\frac{\hat{B}R^{2}}{2}\right)=0.$

Note the factor $\exp\left(-\frac{\hat{B}R^{2}}{2}\right)$ which does not
affect the solutions of the equation, but makes the equation invariant with
respect to $j_{3}\to-j_{3},\ B\to-B$ (via Kummer transformation). So, the
decomposition of $\zeta(s)$ and the free energy will contain only even powers
of $B$. At $B\to 0$ this equation transforms to

$\displaystyle
A(\lambda,k,j_{3},0,R|\alpha)=\frac{1}{\lambda^{2}}(k+1)!\left(\frac{\lambda
R}{2}\right)^{-k-1}\left[i\frac{\lambda R}{2}J_{k+1}(\lambda
R)-\left(\frac{k}{2}+\frac{1}{2}-j\right)J_{k+1}(\lambda R)+\frac{\lambda
R}{2}J_{k+2}(\lambda R)\right]\times$
$\displaystyle\times(k+1)!\left(\frac{\lambda
R}{2}\right)^{-k-1}\left[i\frac{\lambda R}{2}J_{k+1}(\lambda
R)-\left(\frac{k}{2}+\frac{1}{2}+j\right)J_{k+1}(\lambda R)+\frac{\lambda
R}{2}J_{k+2}(\lambda R)\right]$
$\displaystyle-\frac{1}{\lambda^{2}}\left[\frac{k}{2}\left(\frac{k}{2}+1\right)-j^{2}+\frac{1}{4}\right]\left[(k+1)!\left(\frac{\lambda
R}{2}\right)^{-k-1}J_{k+1}(\lambda R)\right]^{2}=0.$

It follows that

$\zeta_{\not{D}^{2}}(s)=\zeta_{\not{D}^{2}}^{(1)}(s)+\zeta_{\not{D}^{2}}^{(2)}(s),$

where the first term summarizes contributions of the eigenvalues with maximal
and minimal projections of total angular momentum $j_{3}=\pm\frac{k+1}{2}$,
and the second term with intermediate values
$j_{3}=-\frac{k-1}{2},\dots,\frac{k-1}{2}$. It is easily seen that

$A(\lambda,k,j_{3},\hat{B},R)=A(\lambda,k,-j_{3},-\hat{B},R),$

so the invariance of the spectrum under transformation $B\to-B$ is manifest.

###  C.1 Contribution of maximal and minimal $j_{3}$

The following considerations are analogous to Appendix A. We start with the
representation

$\displaystyle\zeta_{\not{D}^{2}}^{(1)}(s)$ $\displaystyle=$
$\displaystyle\mathrm{Tr}\left\\{\frac{\sin\pi
s}{\pi}\sum_{k=0}^{\infty}\sum_{j_{3}=\pm\frac{k+1}{2}}(k+1)\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\Psi^{\mathrm{q}}(k+1,j_{3},t,\hat{B},R)\right\\}$
$\displaystyle=$ $\displaystyle\mathrm{Tr}\left\\{\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\left[\sum_{j_{3}=\pm\frac{k}{2}}\Psi^{\mathrm{q}}(k,j_{3},kt,\hat{B},R)-\sum_{i=1}^{2}\frac{u_{i}^{\mathrm{q}}(t,\hat{B},R)}{k^{i}}\right]\right.$
$\displaystyle+\left.\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}+\sum_{i=1}^{2}\frac{u_{i}^{\mathrm{q}}(t,\hat{B},R)}{k^{i}}\right\\}.$

The functions $u_{i}^{\mathrm{q}}$ are given by

$\displaystyle u_{1}^{\mathrm{q}}(t,\hat{B},R)$ $\displaystyle=$
$\displaystyle\hat{B}^{2}R^{4}\frac{-8-4R^{2}t^{2}+R^{4}t^{4}+8\sqrt{1+R^{2}t^{2}}}{6R^{4}t^{4}\sqrt{1+R^{2}t^{2}}},$
$\displaystyle u_{2}^{\mathrm{q}}(t,\hat{B},R)$ $\displaystyle=$
$\displaystyle\hat{B}^{2}R^{4}\frac{-8-12R^{2}t^{2}-3R^{4}t^{4}+8\sqrt{1+R^{2}t^{2}}+8R^{2}t^{2}\sqrt{1+R^{2}t^{2}}}{4R^{4}t^{4}\left(1+R^{2}t^{2}\right)^{\frac{3}{2}}}.$

Evaluating the sums and integrals, one finds

$\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=1}^{2}\frac{u_{i}^{\mathrm{q}}(t,\hat{B},R)}{k^{i}}=\hat{B}^{2}R^{4+2s}\left\\{\zeta(2s)\frac{s\Gamma(1-s)\Gamma\left(s+\frac{1}{2}\right)}{2\sqrt{\pi}(s+2)}-\zeta(1+2s)\frac{s\Gamma(1-s)\Gamma\left(s+\frac{3}{2}\right)}{2\sqrt{\pi}(s+2)}\right\\}.$

Expansion of the counterterm in powers of $s$ is

$\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}k^{1-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{q}}(t,\hat{B},R)}{k^{i}}=-\frac{1}{16}\hat{B}^{2}R^{4}s+O(s^{2}),$

and the corresponding contribution to the effective potential reads

$\frac{1}{2}\left.\frac{d}{ds}\zeta^{(1)}_{\not{D}^{2}}\left(s\right)\right|_{s=0}=-\frac{1}{32}\mathrm{Tr}\hat{B}^{2}R^{4}=-\frac{1}{64}B^{2}R^{4}.$

###  C.2 Contribution of intermediate $j_{3}$

The starting expression in this case is

$\displaystyle\zeta_{\not{D}^{2}}^{(2)}(s)$ $\displaystyle=$
$\displaystyle\mathrm{Tr}\left\\{\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}\sum_{j_{3}=-\frac{k-1}{2}}^{\frac{k-1}{2}}(k+1)\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\Psi^{\mathrm{q}}(k+1,j_{3},t,\hat{B},R)\right\\}$
$\displaystyle=$ $\displaystyle\mathrm{Tr}\left\\{\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}(k+1)k^{-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\left[\sum_{j_{3}=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\mathrm{q}}(k+1,j_{3},kt,\hat{B},R)-\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{q}}(t,\hat{B},R)}{k^{i}}\right]\right.$
$\displaystyle\left.+\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}(k+1)k^{-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{q}}(t,\hat{B},R)}{k^{i}}\right\\},$

where

$\displaystyle u_{0}^{\mathrm{q}}(t,B,R)$ $\displaystyle=$
$\displaystyle\hat{B}^{2}R^{4}\frac{R^{4}t^{4}-2R^{2}t^{2}-4+4\left(1+R^{2}t^{2}\right)^{1/2}}{6R^{4}t^{4}\left(1+R^{2}t^{2}\right)^{1/2}},$
$\displaystyle u_{1}^{\mathrm{q}}(t,\hat{B},R)$ $\displaystyle=$
$\displaystyle\hat{B}^{2}R^{4}\frac{-3R^{4}t^{4}-8+8\left(1+R^{2}t^{2}\right)^{1/2}-12R^{2}t^{2}+8R^{2}t^{2}\left(1+R^{2}t^{2}\right)^{1/2}}{6t^{4}\left(1+R^{2}t^{2}\right)^{\frac{3}{2}}},$
$\displaystyle u_{2}^{\mathrm{q}}(t,\hat{B},R)$ $\displaystyle=$
$\displaystyle\frac{\hat{B}^{2}R^{4}}{48\left(1+R^{2}t^{2}\right)^{1/2}}\left[-\frac{32}{R^{2}t^{2}}+\frac{R^{2}t^{2}(8+13R^{2}t^{2})}{(1+R^{2}t^{2})^{3}}-64\frac{1-\left(1+R^{2}t^{2}\right)^{1/2}}{R^{4}t^{4}}\right]$
$\displaystyle+\frac{\hat{B}^{4}R^{8}}{480\left(1+R^{2}t^{2}\right)^{1/2}}\left[-\frac{128}{R^{6}t^{6}}+\frac{32}{R^{4}t^{4}}-\frac{16}{R^{2}t^{2}}+\frac{10+13R^{2}t^{2}}{(1+R^{2}t^{2})^{2}}-256\frac{1-\left(1+R^{2}t^{2}\right)^{1/2}}{R^{8}t^{8}}\right].$

We find

$\displaystyle\sum_{k=1}^{\infty}(k+1)k^{-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{q}}(t,\hat{B},R)}{k^{i}}=$
$\displaystyle\hat{B}^{2}R^{4+2s}\left\\{\vphantom{\frac{1}{1}}\right.\zeta(-1+2s)\frac{-\Gamma(2-s)\Gamma(1/2+s)}{6\sqrt{\pi}(2+s)}+\zeta(2s)\frac{-(1+2s^{2})\Gamma(1-s)\Gamma(1/2+s)}{6\sqrt{\pi}(2+s)}+$
$\displaystyle\left.\zeta(1+2s)\frac{(12+s(4-s)(1+s))\Gamma(1-s)\Gamma(3/2+s)}{18\sqrt{\pi}(2+s)}\right\\}-\hat{B}^{4}R^{8+2s}\zeta(1+2s)\frac{s\Gamma(2-s)\Gamma(3/2+s)}{240\sqrt{\pi}(4+s)}.$

The small-$s$ expansion is

$\displaystyle\frac{\sin\pi
s}{\pi}\sum_{k=1}^{\infty}(k+1)k^{-2s}\int_{0}^{\infty}\frac{dt}{t^{2s}}\frac{d}{dt}\sum_{i=0}^{2}\frac{u_{i}^{\mathrm{q}}(t,\hat{B},R)}{k^{i}}$
$\displaystyle=\hat{B}^{2}R^{4}\left[\frac{1}{12}+\frac{1}{144}(29+24\gamma-24\log
2+4\pi^{2}+24\log
R)s\right]+\hat{B}^{4}R^{8}\left[-\frac{1}{1920}s\right]+O(s^{2}).$

And we find second part of the effective potential

$\displaystyle\frac{1}{2}\left.\frac{d}{ds}\zeta_{\not{D}^{2}}^{(2)}\left(s\right)\right|_{s=0}=-\frac{1}{2}\mathrm{Tr}\sum_{k=1}^{\infty}(k+1)\left[\sum_{j_{3}=-\frac{k-1}{2}}^{\frac{k-1}{2}}\Psi^{\text{q}}(k,j_{3},0,\hat{B},R)-\hat{B}^{2}R^{4}\frac{1}{12}\left(1-\frac{2}{k^{2}}\right)\right]$
$\displaystyle+\frac{B^{2}R^{4}}{576}(29+24\gamma-24\log 2+4\pi^{2}+24\log
R)-\frac{B^{4}R^{8}}{30720}.$

Combining together two contributions $\zeta_{\not{D}^{2}}^{(1)}$ and
$\zeta_{\not{D}^{2}}^{(2)}$ considered above, one obtains Eq. (31).

##  Appendix D Connection between Kummer and Bessel functions

We find the desired limit from definition of Kummer function via series (see
[55], 13.1.2)

$\lim_{z\to 0}M\left(a+\frac{b}{z},c,dz\right)=\lim_{z\to
0}\sum_{n=0}^{\infty}\frac{\left(a+\frac{b}{z}\right)_{n}(dz)^{n}}{(c)_{n}n!}=\sum_{n=0}^{\infty}\frac{(bd)^{n}}{(c)_{n}n!}={{}_{0}F_{1}}(c,bd).$

Here $(c)_{n}$ is Pochhammer symbol:

$(c)_{0}=1,\quad(c)_{n}=c(c+1)\cdots(c+n-1).$

Comparing series expansions of Bessel functions ([55], 9.1.10, 9.6.10) and
hypergeometric function ${}_{0}F_{1}$, we find

$\displaystyle{{}_{0}F_{1}}(a,z)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\frac{z^{k}\Gamma(a)}{k!\Gamma(a+k)}=(a-1)!\left(\sqrt{z}\right)^{1-a}\left(\sqrt{z}\right)^{a-1}\sum_{k=0}^{\infty}\frac{(\sqrt{z})^{2k}}{k!\Gamma(a-1+k+1)}=(a-1)!\left(\sqrt{z}\right)^{1-a}I_{a-1}(2\sqrt{z}),$
$\displaystyle{{}_{0}F_{1}}(a,-z)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\frac{(-z)^{k}\Gamma(a)}{k!\Gamma(a+k)}=(a-1)!\left(\sqrt{z}\right)^{1-a}\left(\sqrt{z}\right)^{a-1}\sum_{k=0}^{\infty}\frac{(-1)^{k}(\sqrt{z})^{2k}}{k!\Gamma(a-1+k+1)}=(a-1)!\left(\sqrt{z}\right)^{1-a}J_{a-1}(2\sqrt{z}).$

for $z\geqslant 0$.

##  Appendix E Asymptotic expansion for quasi-zero eigenvalues

Let $x$ stands for the smallest solution of the equation. The equation

$M\left(-\frac{x(z)^{2}}{2z},k+2,\frac{z}{2}\right)=0,\quad k=0,1,2,\dots$
(49)

reduces to

$\left(\frac{x(0)}{2}\right)^{-k-1}J_{k+1}(x(0))=0,\quad k=0,1,\dots$

at $z\to 0$. The solution can be sought in the form of a power series in $z$

$x(z)=\sum_{i=0}^{\infty}\frac{z^{i}}{i!}\frac{\partial^{i}}{\partial
z^{i}}x^{(i)}(0).$ (50)

The derivatives $x^{(i)}(0)$ can be found successively from the equations

$\frac{\partial^{i}}{\partial
z^{i}}\left.M\left(-\frac{x(z)^{2}}{2z},k+2,\frac{z}{2}\right)\right|_{z=0}=0.$

We find

$\displaystyle\left.\frac{dx(z)}{dz}\right|_{z=0}$ $\displaystyle=$
$\displaystyle-\frac{k+2}{2x(0)},$
$\displaystyle\left.\frac{d^{2}x(z)}{dz^{2}}\right|_{z=0}$ $\displaystyle=$
$\displaystyle-\frac{k^{2}+8k+12-x^{2}(0)}{12x^{3}(0)},$
$\displaystyle\left.\frac{d^{3}x(z)}{dz^{3}}\right|_{z=0}$ $\displaystyle=$
$\displaystyle-\frac{(k+2)(k^{2}+8k+12-x^{2}(0))}{8x^{5}(0)},$

where we used identities (see [55], 9.1.27)

$\displaystyle J_{\nu-1}(z)+J_{\nu+1}=\frac{2\nu}{z}J_{\nu}(z),$
$\displaystyle J^{\prime}_{\nu}(z)=-J_{\nu+1}(z)+\frac{\nu}{z}J_{\nu}(z)$

and the fact that

$J_{k+1}(x(0))=0.$

Next, we use uniform asymptotic expansion of zeros of Bessel functions
$J_{\nu}(z)$ (see [50], Eq. 10.21.vii)

$\rho_{\nu}(t)=\nu\sum_{k=0}^{\infty}\frac{\alpha_{k}}{\nu^{2k/3}},\quad\theta\left(-2^{\frac{1}{3}}\alpha\right)=\pi
t,\quad\alpha_{0}=1,\quad\alpha_{1}=\alpha,\quad\alpha_{2}=\frac{3}{10}\alpha^{2},\dots$
(51)

where $\theta(x)$ is the phase of Airy functions

$\theta(z)=\arctan\frac{\mathrm{Ai}(z)}{\mathrm{Bi}(z)}.$

Quasi-zero solutions of (49) at $z\to 0$ become smallest zeroes of Bessel
functions that correspond to $t=1$. With $t=1$ formula (51) gives the desired
expansion for the first zero of $J_{k+1}(z)$

$x(0)=(k+1)\sum_{i=0}^{\infty}\frac{\alpha_{k}}{(k+1)^{2i/3}},$ (52)

where $\alpha\approx 1.855757$. This series can be once again expanded in
powers of $k$.

According to Eq. (52), $x(0)\sim k$ for $k\gg 1$. One notices that power of
leading-asymptotics term of derivative $x^{(i)}(0)$ is smaller for larger $i$

$x(0)\sim k,\left.\frac{dx(a)}{da}\right|_{a=0}\sim k^{0},\dots$

so we need only several terms in series (50) to find asymptotic expansion of
$x(z)$ in $k$ up to given order:

$x(z)=k+\alpha
k^{1/3}+1+\frac{3\alpha^{2}}{10}k^{-1/3}+\frac{\alpha}{3}k^{-2/3}+\left(\frac{1}{70}-\frac{\alpha^{3}}{350}\right)k^{-1}+z\left[-\frac{1}{2}+\frac{\alpha}{2}k^{-2/3}-\frac{1}{2}k^{-1}\right]+O(k^{-4/3}).$
(53)

## References

  * Callan  _et al._ [1978] C. G. Callan, Jr., R. F. Dashen, and D. J. Gross, Toward a Theory of the Strong Interactions, Phys. Rev. D17, 2717 (1978), [,36(1977)].
  * Diakonov and Petrov [1984] D. Diakonov and V. Yu. Petrov, Instanton Based Vacuum from Feynman Variational Principle, Nucl. Phys. B245, 259 (1984).
  * Eichhorn  _et al._ [2011] A. Eichhorn, H. Gies, and J. M. Pawlowski, Gluon condensation and scaling exponents for the propagators in Yang-Mills theory, Phys. Rev. D83, 045014 (2011), [Erratum: Phys. Rev.D83,069903(2011)], arXiv:1010.2153 [hep-ph] .
  * Savvidy [1977] G. K. Savvidy, Infrared Instability of the Vacuum State of Gauge Theories and Asymptotic Freedom, Phys. Lett. 71B, 133 (1977).
  * Pagels and Tomboulis [1978] H. Pagels and E. Tomboulis, Vacuum of the Quantum Yang-Mills Theory and Magnetostatics, Nucl. Phys. B143, 485 (1978).
  * Minkowski [1978] P. Minkowski, Comment on an Incorrect Interdependence of the PCAC - Quark Masses, the Potential CP Violating Phase in QCD and the Mass of the Alleged Axion, Phys. Lett. 76B, 439 (1978).
  * Leutwyler [1980] H. Leutwyler, Vacuum Fluctuations Surrounding Soft Gluon Fields, Phys. Lett. 96B, 154 (1980).
  * Leutwyler [1981] H. Leutwyler, Constant Gauge Fields and their Quantum Fluctuations, Nucl. Phys. B179, 129 (1981).
  * Trottier and Woloshyn [1993] H. D. Trottier and R. M. Woloshyn, The Savvidy ’ferromagnetic vacuum’ in three-dimensional lattice gauge theory, Phys. Rev. Lett. 70, 2053 (1993), arXiv:hep-lat/9210028 [hep-lat] .
  * Flory [1983a] C. A. Flory, COVARIANT CONSTANT CHROMOMAGNETIC FIELDS AND ELIMINATION OF THE ONE LOOP INSTABILITIES, Submitted to: Phys. Rev. D (1983a).
  * Flory [1983b] C. A. Flory, A Selfdual Gauge Field, Its Quantum Fluctuations, and Interacting Fermions, Phys. Rev. D28, 1425 (1983b).
  * Elizalde and Soto [1985] E. Elizalde and J. Soto, Zeta Regularized Lagrangians for Massive Quarks in Constant Background Mean Fields, Annals Phys. 162, 192 (1985).
  * Nedelko and Voronin [2015] S. N. Nedelko and V. E. Voronin, Domain wall network as QCD vacuum and the chromomagnetic trap formation under extreme conditions, Eur. Phys. J. A51, 45 (2015), arXiv:1403.0415 [hep-ph] .
  * de Forcrand [2007] P. de Forcrand, Localization properties of fermions and bosons, AIP Conf. Proc. 892, 29 (2007), arXiv:hep-lat/0611034 .
  * Moran and Leinweber [2008] P. Moran and D. Leinweber, Buried treasure in the sand of the QCD vacuum, in  _QCD Downunder II_ (2008) arXiv:0805.4246 [hep-lat] .
  * Bruckmann  _et al._ [2012] F. Bruckmann, F. Gruber, N. Cundy, A. Schafer, and T. Lippert, Topology of dynamical lattice configurations including results from dynamical overlap fermions, Phys. Lett. B 707, 278 (2012), arXiv:1107.0897 [hep-lat] .
  * Galilo and Nedelko [2011a] B. V. Galilo and S. N. Nedelko, Weyl group, CP and the kink-like field configurations in the effective SU(3) gauge theory, Phys. Part. Nucl. Lett. 8, 67 (2011a), arXiv:1006.0248 [hep-ph] .
  * George  _et al._ [2013] D. P. George, A. Ram, J. E. Thompson, and R. R. Volkas, Symmetry breaking, subgroup embeddings and the Weyl group, Phys. Rev. D87, 105009 (2013), arXiv:1203.1048 [hep-th] .
  * Efimov and Nedelko [1995] G. V. Efimov and S. N. Nedelko, Nambu-Jona-Lasinio model with the homogeneous background gluon field, Phys. Rev. D51, 176 (1995).
  * Burdanov  _et al._ [1996] Ya. V. Burdanov, G. V. Efimov, S. N. Nedelko, and S. A. Solunin, Meson masses within the model of induced nonlocal quark currents, Phys. Rev. D54, 4483 (1996), arXiv:hep-ph/9601344 [hep-ph] .
  * Kalloniatis and Nedelko [2001] A. C. Kalloniatis and S. N. Nedelko, Confinement and chiral symmetry breaking via domain-like structures in the QCD vacuum, Phys. Rev. D64, 114025 (2001), arXiv:hep-ph/0108010 [hep-ph] .
  * Kalloniatis and Nedelko [2004] A. C. Kalloniatis and S. N. Nedelko, Realization of chiral symmetry in the domain model of QCD, Phys. Rev. D69, 074029 (2004), [Erratum: Phys. Rev.D70,119903(2004)], arXiv:hep-ph/0311357 [hep-ph] .
  * Nedelko and Voronin [2016] S. N. Nedelko and V. E. Voronin, Regge spectra of excited mesons, harmonic confinement and QCD vacuum structure, Phys. Rev. D93, 094010 (2016), arXiv:1603.01447 [hep-ph] .
  * Nedelko and Voronin [2017] S. N. Nedelko and V. E. Voronin, Influence of confining gluon configurations on the $P\to\gamma^{*}\gamma$ transition form factors, Phys. Rev. D95, 074038 (2017), arXiv:1612.02621 [hep-ph] .
  * Bali  _et al._ [2012] G. Bali, F. Bruckmann, G. Endrodi, Z. Fodor, S. Katz, S. Krieg, A. Schafer, and K. Szabo, The QCD phase diagram for external magnetic fields, JHEP 02, 044, arXiv:1111.4956 [hep-lat] .
  * Bali  _et al._ [2014] G. Bali, F. Bruckmann, G. Endrödi, S. Katz, and A. Schäfer, The QCD equation of state in background magnetic fields, JHEP 08, 177, arXiv:1406.0269 [hep-lat] .
  * Bonati  _et al._ [2016] C. Bonati, M. D’Elia, M. Mariti, M. Mesiti, F. Negro, A. Rucci, and F. Sanfilippo, Magnetic field effects on the static quark potential at zero and finite temperature, Phys. Rev. D 94, 094007 (2016), arXiv:1607.08160 [hep-lat] .
  * Kondo  _et al._ [2015] K.-I. Kondo, S. Kato, A. Shibata, and T. Shinohara, Quark confinement: Dual superconductor picture based on a non-Abelian Stokes theorem and reformulations of Yang–Mills theory, Phys. Rept. 579, 1 (2015), arXiv:1409.1599 [hep-th] .
  * Nedelko and Voronin [2019] S. N. Nedelko and V. E. Voronin, Finite size effects in the free energy density for Abelian (anti-)self-dual gluon field in $SU(3)$ gluodynamics, Phys. Part. Nucl. Lett. 16, 553 (2019), arXiv:1906.00432 [hep-ph] .
  * Faddeev [2009] L. D. Faddeev, Mass in Quantum Yang-Mills Theory: Comment on a Clay Millenium problem, (2009), arXiv:0911.1013 [math-ph] .
  * Faddeev [2014] L. D. Faddeev, A couple of methodological comments on the quantum Yang-Mills theory, Theor. Math. Phys. 181, 1638 (2014).
  * Galilo and Nedelko [2011b] B. V. Galilo and S. N. Nedelko, Impact of the strong electromagnetic field on the QCD effective potential for homogeneous Abelian gluon field configurations, Phys. Rev. D84, 094017 (2011b), arXiv:1107.4737 [hep-ph] .
  * Shabanov and Klauder [1999] S. V. Shabanov and J. R. Klauder, Towards a nonperturbative path integral in gauge theories, Phys. Lett. B456, 38 (1999), arXiv:hep-th/9902017 [hep-th] .
  * Shabanov [2000] S. V. Shabanov, Geometry of the physical phase space in quantum gauge systems, Phys. Rept. 326, 1 (2000), arXiv:hep-th/0002043 [hep-th] .
  * Cho [1980] Y. Cho, COLORED MONOPOLES, Phys. Rev. Lett. 44, 1115 (1980), [Erratum: Phys.Rev.Lett. 44, 1566 (1980)].
  * Faddeev and Niemi [2007] L. D. Faddeev and A. J. Niemi, Spin-Charge Separation, Conformal Covariance and the SU(2) Yang-Mills Theory, Nucl. Phys. B776, 38 (2007), arXiv:hep-th/0608111 [hep-th] .
  * Kondo  _et al._ [2006] K. I. Kondo, T. Murakami, and T. Shinohara, Yang-Mills theory constructed from Cho-Faddeev-Niemi decomposition, Prog. Theor. Phys. 115, 201 (2006), arXiv:hep-th/0504107 [hep-th] .
  * Faddeev [2010] L. D. Faddeev, Separation of scattering and selfaction revisited, in  _Subtleties in quantum field theory: Lev Lipatov Festschrift_ , edited by D. Diakonov (2010) pp. 1–6, arXiv:1003.4854 [hep-th] .
  * Faddeev [2006] L. D. Faddeev, Notes on divergences and dimensional transmutation in Yang-Mills theory, Theor. Math. Phys. 148, 986 (2006), [Teor. Mat. Fiz.148,133(2006)].
  * Vachaspati [2006] T. Vachaspati,  _Kinks and Domain walls_ (Cambridge University Press, 2006).
  * Cho  _et al._ [2012] Y. M. Cho, D. G. Pak, P. M. Zhang, and L. P. Zou, Weyl symmetric structure of qcd vacuum, Phys. Rev. D 86, 045025 (2012).
  * Bordag  _et al._ [1996] M. Bordag, E. Elizalde, and K. Kirsten, Heat kernel coefficients of the Laplace operator on the D-dimensional ball, J. Math. Phys. 37, 895 (1996), arXiv:hep-th/9503023 [hep-th] .
  * Kirsten [2001] K. Kirsten,  _Spectral Functions in Mathematics and Physics_ (Chapman and Hall/CRC, 2001).
  * Deser  _et al._ [1997] S. Deser, L. Griguolo, and D. Seminara, Gauge invariance, finite temperature and parity anomaly in D = 3, Phys. Rev. Lett. 79, 1976 (1997), arXiv:hep-th/9705052 [hep-th] .
  * Wiesendanger and Wipf [1994] C. Wiesendanger and A. Wipf, Running coupling constants from finite size effects, Annals Phys. 233, 125 (1994).
  * Cognola  _et al._ [1993] G. Cognola, K. Kirsten, and S. Zerbini, One loop effective potential on hyperbolic manifolds, Phys. Rev. D48, 790 (1993), arXiv:hep-th/9302051 [hep-th] .
  * Braguta  _et al._ [2020] V. Braguta, A. Kotov, D. Kuznedelev, and A. Roenko, Study of the Confinement/Deconfinement Phase Transition in Rotating Lattice SU(3) Gluodynamics, JETP Lett. 112, 6 (2020).
  * Yamamoto and Hirono [2013] A. Yamamoto and Y. Hirono, Lattice QCD in rotating frames, Phys. Rev. Lett. 111, 081601 (2013), arXiv:1303.6292 [hep-lat] .
  * Olver [1997] F. W. J. Olver, Asymptotics and special functions (A K Peters/CRC Press, New York, 1997) Chap. Differential Equations with a Parameter: Expansions in Elementary Functions, pp. 382–386, 1st ed.
  * [50] DLMF, NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/, Release 1.0.23 of 2019-06-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller and B. V. Saunders, eds.
  * Pais [1954] A. Pais, Spherical Spinors in a Euclidean 4-Space, Proceedings of the National Academy of Sciences of the United States of America 40, 835 (1954).
  * Carter and Mclenaghan [1979] B. Carter and R. G. Mclenaghan, Generalized Total Angular Momentum Operator For The Dirac Equation In Curved Space-Time, Phys. Rev. D19, 1093 (1979).
  * Shapovalov [2018] V. N. Shapovalov, Symmetry and Classification of the Dirac–Fock Equation, Theor. Math. Phys. 197, 1572 (2018), [Teor. Mat. Fiz.197,no.2,208(2018)].
  * Breev and Shapovalov [2016] A. I. Breev and A. V. Shapovalov, The Dirac equation in an external electromagnetic field: symmetry algebra and exact integration,  _Proceedings, 23rd International Conference on Integrable Systems and Quantum Symmetries (ISQS-23): Prague, Czech Republic, June 23-27, 2015_ , J. Phys. Conf. Ser. 670, 012015 (2016), arXiv:1509.08612 [math-ph] .
  * Abramowitz and Stegun [1972] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ (U.S. Dept. of Commerce, National Bureau of Standards, Washington, D.C., 1972).