#  Chromoelectric flux tubes within non-Abelian Proca theory

Vladimir Dzhunushaliev<EMAIL_ADDRESS>Department of Theoretical
and Nuclear Physics, Al-Farabi Kazakh National University, Almaty 050040,
Kazakhstan  Academician J. Jeenbaev Institute of Physics of the NAS of the
Kyrgyz Republic, 265 a, Chui Street, Bishkek 720071, Kyrgyzstan  Vladimir
Folomeev<EMAIL_ADDRESS>Academician J. Jeenbaev Institute of Physics of
the NAS of the Kyrgyz Republic, 265 a, Chui Street, Bishkek 720071, Kyrgyzstan
International Laboratory for Theoretical Cosmology, Tomsk State University of
Control Systems and Radioelectronics (TUSUR), Tomsk 634050, Russia

###### Abstract

Flux tube solutions within non-Abelian SU(3) Proca theory with external
sources are obtained. It is shown that such tubes have a longitudinal
chromoelectric field possessing two components (nonlinear and gradient), as
well as a transverse chromomagnetic field whose force lines create concentric
circles with the center on the axis of the tube. The scenario of a possible
relationship between non-Abelian Proca theory and quantum chromodynamics is
considered. In such scenario: (a) the components of color fields have
different behavior: those which are almost classical, and those which are
purely quantum; (b) the second components create a gluon condensate that is a
source of the field for the almost classical components of the Proca field;
(c) Proca masses may appear as a result of an approximate description of the
gluon condensate; (d) the question of gauge invariance is considered. It is
shown that the results obtained are in good agreement with the results of
lattice calculations. We make an assumption that an approximate description of
a flux tube in quantum chromodynamics can be carried out using classical Proca
equations but with a mandatory account of a gluon condensate.

non-Abelian Proca theory, external sources, flux tube, total energy, linear
dependence

###### pacs:

11.90.+t

##  I Introduction

One of manifestations of the confinement phenomenon in quantum chromodynamics
(QCD) is the presence of a tube filled with a longitudinal chromoelectric
field created by quarks located at the ends of the tube. This phenomenon has
been supported by lattice calculations (see, e.g., Refs. [1, 2, 3, 4, 5, 6, 7]
and references therein). In particular, in those works, it was shown that
between static quarks there appears a flux tube threaded by chromoelectric and
chromomagnetic fields, whose structure and a relationhip between the fields
and currents were studied. These calculations are based on the direct
application of numerical methods for nonperturbative quantization, i.e., for
the calculation of the corresponding path integral. However, in order to
clearly appreciate physical processes, one would like to have some approximate
model describing the appearance of the flux tube between quarks filled with a
longitudinal chromoelectric field.

In the present paper we show that in non-Abelian Yang-Mills-Proca theory there
exist such solutions describing flux tubes supported by external sources. In
recent years, there is a growing interest in Proca fields inspired by the
possibility of obtaining within such theories new solutions that can be used
for a description of various physical objects and processes. For example, the
literature in the field suggests gravitating compact configurations [8, 9, 10,
11, 12, 13, 14], various aspects of black hole physics [15, 16], dark matter
models [17, 18], a study of processes at cosmological scales [19, 20], a
consideration of effects arising due to the possible presence of the rest mass
of a photon [21], generalized Proca theories with second-order equations of
motion in a curved spacetime [22]. Also, as applied to modelling flux tubes,
Proca fields have been used by us earlier in constructing various tubelike
configurations filled with electric and magnetic fields [23, 24, 25, 26, 27].

A distinctive feature of the configurations considered in the present work is
that they contain a flux of chromoelectric field directed from one static
source to another. This field consists of two components, one of which is
nonlinear, and therefore it may be called nonperturbative, while another one
is gradient, and it is similar to an electric field in Maxwell’s
electrodynamics. Another feature of the systems obtained here is that the
total energy of such field configurations increases linearly with increasing
distance between the static sources.

These properties of the solutions enable us to assume that there can be some
relationship between non-Abelian Proca theory under consideration and QCD. In
particular, one can assume that non-Abelian Proca theory is some approximation
for nonperturbative QCD. Such an assumption leads to interesting consequences
(for details see the discussion in Sec. V): (a) the components of color fields
have different behavior: those which are almost classical (with a nonzero
quantum average), and those which are purely quantum (with a zero quantum
average); (b) the second components create a gluon condensate that is a source
of the field for the almost classical components of the Proca field; (c) Proca
masses may appear as a result of an approximate description of the gluon
condensate.

The paper is organized as follows. In Secs. II and III, we write down the
general field equations for non-Abelian Proca theory employed. In Sec. IV, we
solve these equations numerically and obtain axially symmetric regular
solutions describing localized flux tube configurations consisting of
chromoelectric and chromomagnetic Proca fields sourced by external charges and
currents. In Sec. VI, we compare our computations with the results of lattice
calculations. In Sec. V, we discuss a possible scenario describing the
relationship between non-Abelian Proca theory under consideration and QCD.
Finally, Sec. VII summarizes the results obtained in the paper.

##  II Non-Abelian SU(3) Yang-Mills-Proca theory

The Lagrangian describing a system consisting of a non-Abelian SU(3) Proca
field $A^{a}_{\mu}$ sourced by color charges and currents can be taken in the
form (hereafter, we work in natural units $c=\hbar=1$)

$\mathcal{L}=-\frac{1}{4}F^{a}_{\mu\nu}F^{a\mu\nu}+\frac{m^{2}}{2}A^{a}_{\mu}A^{a\mu}-A^{a}_{\mu}j^{a\mu}.$
(1)

Here
$F^{a}_{\mu\nu}=\partial_{\mu}A^{a}_{\nu}-\partial_{\nu}A^{a}_{\mu}+gf_{abc}A^{b}_{\mu}A^{c}_{\nu}$
is the field strength tensor for the Proca field of mass $m$, where $f_{abc}$
are the SU(3) structure constants, $g$ is the coupling constant, $j^{a\mu}$
are color currents, $a,b,c=1,2,\dots,8$ are color indices, $\mu,\nu=0,1,2,3$
are spacetime indices.

Using (1), the corresponding field equations can be written in the form:

$D_{\nu}F^{a\mu\nu}-m^{2}A^{a\mu}\equiv\frac{1}{\sqrt{-\mathcal{G}}}\frac{\partial}{\partial
x^{\nu}}\left(\sqrt{-\mathcal{G}}F^{a\mu\nu}\right)+gf_{abc}A^{b}_{\nu}F^{c\mu\nu}-m^{2}A^{a\mu}=-j^{a\mu},$
(2)

where $\mathcal{G}$ is the determinant of the spacetime metric, and the energy
density is

$\begin{split}\varepsilon=&\frac{1}{2}\left(E^{a}_{i}\right)^{2}+\frac{1}{2}\left(H^{a}_{i}\right)^{2}+m^{2}\left(A^{a0}A^{a}_{0}-\frac{1}{2}A^{a}_{\alpha}A^{a\alpha}\right),\end{split}$
(3)

where $i=1,2,3$ and $E^{a}_{i}$ and $H^{a}_{i}$ are the components of the
color electric and magnetic field strengths, respectively.

##  III General equations and Ansatz

In this section we write out the field equations within the theory (1) in some
generalized form in order to analyze the possibility of obtaining flux tube
solutions. Namely, to find a distribution of the non-Abelian Proca field
created by static charges, one can choose the following Ansatz in cylindrical
coordinates $(t,\rho,\varphi,z)$:

$A^{2}_{t}=\frac{f(\rho,z)}{g},\;A^{2}_{\varphi}=\frac{\rho
k(\rho,z)}{g},\;A^{5}_{\rho}=\frac{u(\rho,z)}{g},A^{5}_{z}=\frac{v(\rho,z)}{g},A^{7}_{t}=\frac{h(\rho,z)}{g},\;A^{7}_{\varphi}=\frac{\rho
w(\rho,z)}{g}.$ (4)

Henceforth the upper indices $2,5,7$ are color indices, and one can introduce
the vector $\vec{A}^{5}=\frac{1}{g}\left\\{u(\rho,z),0,v(\rho,z)\right\\}$.
The potentials (4) give the following nonzero components of color electric and
magnetic field strengths:

$\displaystyle E^{2}_{z}=$ $\displaystyle-\frac{hv}{2g}-\frac{f_{,z}}{g},\quad
E^{7}_{z}=\frac{fv}{2g}-\frac{h_{,z}}{g},\quad
E^{2}_{\rho}=-\frac{hu}{2g}-\frac{f_{,\rho}}{g},\quad
E^{7}_{\rho}=\frac{fu}{2g}-\frac{h_{,\rho}}{g},\quad
E^{5}_{\varphi}=\frac{\rho}{2g}\left(hk-fw\right),$ (5) $\displaystyle
H^{2}_{z}=$
$\displaystyle-\frac{2\left(k_{,\rho}+\frac{k}{\rho}\right)+uw}{2g},\quad
H^{7}_{z}=-\frac{2\left(w_{,\rho}+\frac{w}{\rho}\right)-uk}{2g},\quad
H^{2}_{\rho}=\frac{2k_{,z}+vw}{2g},\quad$ $\displaystyle H^{7}_{\rho}=$
$\displaystyle\frac{\left(2w_{,z}-kv\right)}{2g},\quad
H^{5}_{\varphi}=\rho\frac{\left(v_{,\rho}-u_{,z}\right)}{g},$ (6)

where a comma in lower indices denotes differentiation with respect to the
corresponding coordinate.

The Poynting vector created by crossed color electric and magnetic fields is

$S^{i}=\frac{\epsilon^{ijk}}{\sqrt{\gamma}}E^{a}_{j}H^{a}_{k},$

where $\gamma$ is the determinant of the space metric and $\epsilon^{ijk}$ is
the completely antisymmetric Levi-Civita symbol. For the field strengths (5)
and (6), this expression yields the following nonvanishing component:

$\begin{split}g^{2}\rho
S_{\varphi}=&k\left[-\frac{f_{,\rho}}{\rho}-\frac{f}{4}\left(u^{2}+v^{2}\right)+\frac{u}{2}\left(h_{,\rho}-\frac{h}{\rho}\right)+\frac{vh_{,z}}{2}\right]-k_{,\rho}\left(f_{,\rho}+\frac{hu}{2}\right)-k_{,z}\left(f_{,z}+\frac{hv}{2}\right)\\\
+&w\left[-\frac{h_{,\rho}}{\rho}-\frac{h}{4}\left(u^{2}+v^{2}\right)-\frac{u}{2}\left(f_{,\rho}-\frac{f}{\rho}\right)-\frac{vf_{,z}}{2}\right]-w_{,\rho}\left(h_{,\rho}-\frac{fu}{2}\right)-w_{,z}\left(h_{,z}-\frac{fv}{2}\right).\end{split}$
(7)

The Proca equations (2) for the Ansatz (4) have the form:

$\displaystyle
f_{,zz}+f_{,\rho\rho}+\frac{f_{,\rho}}{\rho}-\frac{f}{4}\left(u^{2}+v^{2}+w^{2}\right)+\frac{hkw}{4}+uh_{,\rho}+\frac{hu_{,\rho}}{2}+vh_{,z}+\frac{hv_{,z}}{2}+\frac{hu}{2\rho}-m^{2}f=$
$\displaystyle-gj^{2}_{t},$ (8) $\displaystyle
k_{,zz}+k_{,\rho\rho}+\frac{k_{,\rho}}{\rho}-\frac{k}{\rho^{2}}+\frac{k}{4}\left(h^{2}-u^{2}-v^{2}\right)-\frac{1}{4}fhw+\frac{wu_{,\rho}}{2}+uw_{,\rho}+\frac{wv_{,z}}{2}+vw_{,z}+\frac{uw}{2\rho}-m^{2}k=$
$\displaystyle g\rho j^{2}_{\varphi},$ (9) $\displaystyle
v_{,\rho\rho}+\frac{v_{,\rho}}{\rho}-u_{,z\rho}-\frac{u_{,z}}{\rho}+\frac{v}{4}\left(f^{2}+h^{2}-k^{2}-w^{2}\right)+\frac{hf_{,z}}{2}-\frac{fh_{,z}}{2}+\frac{kw_{,z}}{2}-\frac{wk_{,z}}{2}-m^{2}v=$
$\displaystyle gj^{5}_{v},$ (10) $\displaystyle
u_{,zz}-v_{,z\rho}+\frac{u}{4}\left(f^{2}+h^{2}-k^{2}-w^{2}\right)+\frac{hf_{,\rho}}{2}-\frac{fh_{,\rho}}{2}+\frac{kw_{,\rho}}{2}-\frac{wk_{,\rho}}{2}-m^{2}u=$
$\displaystyle gj^{5}_{u},$ (11) $\displaystyle
h_{,zz}+h_{,\rho\rho}+\frac{h_{,\rho}}{\rho}-\frac{h}{4}\left(u^{2}+v^{2}+k^{2}\right)+\frac{fkw}{4}-uf_{,\rho}-\frac{fu_{,\rho}}{2}-vf_{,z}-\frac{fv_{,z}}{2}-\frac{fu}{2\rho}-m^{2}h=$
$\displaystyle-gj^{7}_{t},$ (12) $\displaystyle
w_{,zz}+w_{,\rho\rho}+\frac{w_{,\rho}}{\rho}-\frac{w}{\rho^{2}}+\frac{w}{4}\left(f^{2}-u^{2}-v^{2}\right)-\frac{1}{4}fhk-\frac{ku_{,\rho}}{2}-uk_{,\rho}-\frac{kv_{,z}}{2}-vk_{,z}-\frac{ku}{2\rho}-m^{2}w=$
$\displaystyle g\rho j^{7}_{\varphi}.$ (13)

Let us calculate the covariant divergence of the Proca equations (2),

$D_{\mu}\left(D_{\nu}F^{a\mu\nu}-m^{2}A^{a\mu}\right)=-D_{\mu}j^{a\mu}.$ (14)

Taking into account Eqs. (8)-(13), we then have from Eq. (14):

$m^{2}\divergence\vec{A}^{5}\equiv
m^{2}\left(v_{,z}+u_{,\rho}+\frac{u}{\rho}\right)=\frac{1}{2}\left(fj^{7}_{t}-hj^{2}_{t}\right)+\frac{1}{2}\rho\left(kj^{7}_{\varphi}-wj^{2}_{\varphi}\right)-\left(j^{5}_{u}\right)_{,\rho}-\frac{j^{5}_{u}}{\rho}-\left(j^{5}_{v}\right)_{,z}.$
(15)

The conservation of the four-current $j^{a\mu}$ implies that the components of
the current must be chosen so as to make the right-hand side of Eq. (15) be
identically equal to zero. Then, using Eq. (15), one can get rid of mixed
derivatives in Eqs. (10) and (11). As a result, the full system of equations
(8)-(13) for the functions $f,h,u,v,k,$ and $w$ [on account of the divergent
equation (15)] takes the form:

$\displaystyle\laplacian{f}-\frac{f}{4}\left(u^{2}+v^{2}+w^{2}\right)+\frac{hkw}{4}+\vec{A}^{5}\cdot\gradient{h}-m^{2}f=$
$\displaystyle-gj^{2}_{t},$ (16)
$\displaystyle\laplacian{k}-\frac{k}{\rho^{2}}+\frac{k}{4}\left(h^{2}-u^{2}-v^{2}\right)-\frac{1}{4}fhw+\vec{A}^{5}\cdot\gradient{w}-m^{2}k=$
$\displaystyle g\rho j^{2}_{\varphi},$ (17)
$\displaystyle\laplacian{v}+\frac{v}{4}\left(f^{2}+h^{2}-k^{2}-w^{2}\right)+\frac{hf_{,z}}{2}-\frac{fh_{,z}}{2}+\frac{kw_{,z}}{2}-\frac{wk_{,z}}{2}-m^{2}v=$
$\displaystyle gj^{5}_{v},$ (18)
$\displaystyle\laplacian{u}-\frac{u}{\rho^{2}}+\frac{u}{4}\left(f^{2}+h^{2}-k^{2}-w^{2}\right)+\frac{hf_{,\rho}}{2}-\frac{fh_{,\rho}}{2}+\frac{kw_{,\rho}}{2}-\frac{wk_{,\rho}}{2}-m^{2}u=$
$\displaystyle gj^{5}_{u},$ (19)
$\displaystyle\laplacian{h}-\frac{h}{4}\left(u^{2}+v^{2}+k^{2}\right)+\frac{fkw}{4}-\vec{A}^{5}\cdot\gradient{f}-m^{2}h=$
$\displaystyle-gj^{7}_{t},$ (20)
$\displaystyle\laplacian{w}-\frac{w}{\rho^{2}}+\frac{w}{4}\left(f^{2}-u^{2}-v^{2}\right)-\frac{1}{4}fhk-\vec{A}^{5}\cdot\gradient{k}-m^{2}w=$
$\displaystyle g\rho j^{7}_{\varphi},$ (21)

where $\laplacian=\partial_{zz}+\partial_{\rho\rho}+\partial_{\rho}/\rho$ is
the Laplacian in cylindrical coordinates.

##  IV Flux tube solutions

To find a flux tube connecting static charges, consider a simplified system
where in Eqs. (16)-(21) the vector potentials $w=k=0$. As a result, there will
be the following reduced set of equations for the functions $f,h,v,$ and $u$:

$\displaystyle\laplacian{f}-\frac{f}{4}\left(u^{2}+v^{2}\right)+\vec{A}^{5}\cdot\gradient{h}-f=$
$\displaystyle-j^{2}_{t},$ (22)
$\displaystyle\laplacian{v}+\frac{v}{4}\left(f^{2}+h^{2}\right)+\frac{hf_{,z}}{2}-\frac{fh_{,z}}{2}-v=$
$\displaystyle j^{5}_{v},$ (23)
$\displaystyle\laplacian{u}-\frac{u}{\rho^{2}}+\frac{u}{4}\left(f^{2}+h^{2}\right)+\frac{hf_{,\rho}}{2}-\frac{fh_{,\rho}}{2}-u=$
$\displaystyle j^{5}_{u},$ (24)
$\displaystyle\laplacian{h}-\frac{h}{4}\left(u^{2}+v^{2}\right)-\vec{A}^{5}\cdot\gradient{f}-h=$
$\displaystyle-j^{7}_{t},$ (25)

and in this case the Poynting vector (7) is equal to zero.

The equations (22)-(25) are written already in terms of the dimensionless
variables

$\left(\rho,z\right)\to\frac{1}{m}\left(\rho,z\right),\quad\left(f,h,u,v\right)\to
m\left(f,h,u,v\right),\quad\left(j^{2}_{t},j^{7}_{t},j^{5}_{u},j^{5}_{v}\right)\to\frac{m^{3}}{g}\left(j^{2}_{t},j^{7}_{t},j^{5}_{u},j^{5}_{v}\right).$

For our purpose, nonzero components of the four-current are chosen by giving
the following ad hoc expressions:

$\displaystyle j^{2}_{t}=$
$\displaystyle-\exp\left(-\frac{\rho^{2}}{\rho_{1}^{2}}\right)\left\\{j_{1}\exp\left[-\frac{(z-z_{0})^{2}}{z_{1}^{2}}\right]+j_{2}\exp\left[-\frac{(z+z_{0})^{2}}{z_{1}^{2}}\right]\right.$
$\displaystyle\left.+j_{\text{glc}}\left[\frac{1}{1+\exp\left[-10(z+z_{0})\right]}+\frac{1}{1+\exp\left[-10(z-z_{0})\right]}-1\right]\right\\},$
(26) $\displaystyle j^{5}_{u}=$
$\displaystyle\rho\left\\{j_{3}\exp\left[-\frac{\rho^{2}}{\rho_{1}^{2}}-\frac{(z-z_{0})^{2}}{z_{1}^{2}}\right]+j_{4}\exp\left[-\frac{\rho^{2}}{\rho_{1}^{2}}-\frac{(z+z_{0})^{2}}{z_{1}^{2}}\right]\right\\}.$
(27)

Here $j_{1,2,3,4}$, $\rho_{1}$, and $z_{1}$ are some parameters characterizing
linear sizes of color charge and current densities; $z_{0}$ characterizes the
distance between static charges; $j_{\text{glc}}$ is a constant determining
the magnitude of a “gluon condensate” (cf. Sec. V). The latter is modeled in
the form of some distribution of constant charge density which fills the space
between the static charges.

The law of conservation of the current $D_{\mu}j^{a\mu}=0$ following from Eq.
(14) implies that in our case Eq. (15) must hold. In our model problem we
solve this equation as follows:

$\displaystyle v_{,z}+u_{,\rho}+\frac{u}{\rho}=$ $\displaystyle 0,$ (28)
$\displaystyle fj^{7}_{t}-hj^{2}_{t}=$ $\displaystyle 0,$ (29)
$\displaystyle\left(j^{5}_{u}\right)_{,\rho}+\frac{j^{5}_{u}}{\rho}+\left(j^{5}_{v}\right)_{,z}=$
$\displaystyle 0.$ (30)

That is, we take $\divergence\vec{A}^{5}=0$, since for the currents $j^{a\mu}$
created by fields (spinor or scalar) the equation $D_{\mu}j^{a\mu}=0$ must
hold. Then, taking into account the expressions (26) and (27), we have from
Eqs. (29) and (30) the following expressions for the components of the four-
current $j^{7}_{t}$ and $j^{5}_{v}$:

$\displaystyle j^{7}_{t}=$ $\displaystyle\frac{h}{f}j^{2}_{t},$ (31)
$\displaystyle j^{5}_{v}=$
$\displaystyle\frac{\sqrt{\pi}z_{1}}{\rho_{1}^{2}}\exp(-\frac{\rho^{2}}{\rho_{1}^{2}})\left(\rho^{2}-\rho_{1}^{2}\right)\left[j_{3}\erf\left(\frac{z-z_{0}}{z_{1}}\right)+j_{4}\erf\left(\frac{z+z_{0}}{z_{1}}\right)\right].$
(32)

###  IV.1 Boundary conditions and numerical approach

The set of nonlinear partial differential equations (22)-(25) is solved
numerically by specifying appropriate boundary conditions on the region of
integration $0\leq\rho<\infty$, $-\infty<z<\infty$. In order to have a color
longitudinal electric field directed along the axis of the tube, we choose the
boundary conditions as follows:

$\begin{split}&\left.\frac{\partial
f}{\partial\rho}\right|_{\rho=0}=\left.\frac{\partial
h}{\partial\rho}\right|_{\rho=0}=\left.\frac{\partial
v}{\partial\rho}\right|_{\rho=0}=0,\quad\left.u\right|_{\rho=0}=0;\\\
&f=h=u=v=0\quad\text{ as }\quad\rho^{2}+z^{2}\to\infty.\end{split}$ (33)

These boundary conditions enable us to obtain solutions regular over the whole
space.

For numerical calculations, it is convenient to introduce new compactified
coordinates

$\bar{\rho}=\frac{2}{\pi}\arctan\rho,\quad\bar{z}=\frac{2}{\pi}\arctan z,$
(34)

the use of which permits one to map the infinite region $(-\infty,\infty)$ to
the finite interval $[-1,1]$.

All results of numerical calculations given below have been obtained using the
Intel MKL PARDISO sparse direct solver [28, 29] and the CESDSOL
library111Complex Equations-Simple Domain partial differential equations
SOLver, a C++ package developed by I. Perapechka, see, e.g., Ref. [30]..
Solutions are sought on a grid of $200\times 400$ points, covering the region
of integration $0\leq\bar{\rho}\leq 1$ and $-1\leq\bar{z}\leq 1$ given by the
compactified coordinates from Eq. (34). This enables us to get solutions with
typical relative errors on the order of $10^{-4}$. Also, in performing
calculations, we have been checking the fulfilment of the divergent equation
(28).

Figure 1: Functions $f,h,u,v$, force lines for the electric $\vec{E}^{2}$,
$\vec{E}^{7}$ and magnetic $\vec{H}^{5}$ fields [in Cartesian coordinates
$(x,y)$], and the energy density $\varepsilon$ from Eq. (3) for the magnetic-
field-dominated system (two top rows) and for the electric-field-dominated
system (two bottom rows) with the choice $z_{0}=5$ (for the values of other
parameters, see the caption of Fig. 2).

###  IV.2 Asymptotic behavior

Even before obtaining numerical solutions, one can estimate their asymptotic
behavior, bearing in mind that all the functions $f,h,u,$ and $v$ are assumed
to be exponentially decaying with distance. For this purpose, it is convenient
to use a spherical coordinate system $\\{r,\theta,\varphi\\}$. In this case,
as $r\rightarrow\infty$, the functions $f,h,u,v\to 0$ exponentially fast. As a
result, from Eqs. (22)-(25), one can find the following asymptotic equations:

$\displaystyle\bigtriangleup_{r,\theta}(f,h,v)-(f,h,v)=$ $\displaystyle 0,$
(35)
$\displaystyle\bigtriangleup_{r,\theta}u-\frac{u}{r^{2}\sin^{2}{\theta}}-u=$
$\displaystyle 0,$ (36)

where $\bigtriangleup_{r,\theta}$ is the Laplacian in the coordinates
$r,\theta$. Eqs. (35) for the functions $f,h,$ and $v$ have obvious solutions
in the form

$(f,h,v)\approx C_{(f,h,v)}\left(Y\right)^{0}_{l_{(f,h,v)}}\frac{e^{-r}}{r},$
(37)

where $\left(Y\right)^{0}_{l_{(f,h,v)}}$ are spherical functions and
$C_{(f,h,v)}$ are constants. In turn, Eq. (36) has a solution similar to (37),
but only with the angular part expressed in terms of special functions (we do
not show this expression here to avoid overburdening the text). The general
solution of Eqs. (35) and (36) represents a superposition of the above
solutions, and it is obtained by summing over $l_{(f,h,v)}$. Note that the
results of numerical calculations obtained below indicate that the angular
part $\left(Y\right)^{0}_{l_{(f,h,v)}}=\text{const}$.

###  IV.3 Numerical results

Depending on the choice of free parameters determining characteristics of the
currents (26) and (27), one can obtain systems with different ratio of
magnitudes of the electric and magnetic fields with different color indices.
As an illustration, we consider two systems where either the magnetic or
electric field dominates. The corresponding numerical results of solving Eqs.
(22)-(25) are given in Fig. 1 for the magnetic-field-dominated system (two top
rows) and for the electric-field-dominated system (two bottom rows).

For us, it is of especial interest to study the dependence of the total
dimensionless energy of the systems under consideration on the distance
between charges,

$W\equiv\frac{W_{\text{NU}}}{m}=2\pi\int\limits_{-\infty}^{\infty}dz\int\limits_{0}^{\infty}\rho\varepsilon
d\rho,$ (38)

[this expression is derived using Eq. (3), and $W_{\text{NU}}$ is a
dimensional energy in natural units] and the flux $\Phi_{z}$ of the color
electric Proca field $\vec{E}^{7}$ through the plane $z=0$,

$\Phi_{z}\equiv g\Phi_{z\text{NU}}=2\pi\int\limits_{0}^{\infty}\rho
E^{7}_{z}d\rho$ (39)

[this expression is derived using Eq. (5) on account of the transition to the
dimensionless electric field strength
$E^{7}_{z}\to\left(m^{2}/g\right)E^{7}_{z}$ (for other components of the
electric and magnetic field strengths the transition to the dimensionless
variables is similar)]. The corresponding results of computations are given in
Fig. 2.

Figure 2: Total energy $W$ from Eq. (38) and the flux of the longitudinal
chromoelectric field $E^{7}_{z}$ from Eq. (39) as functions of $L=2z_{0}$. The
top plots correspond to the magnetic-field-dominated systems with the choice
$j_{1}=j_{2}=j_{3}=0.2,\rho_{1}=2,z_{1}=1$, $j_{4}=-j_{3}$, and
$j_{\text{glc}}=0.1$; the bottom plots correspond to the electric-field-
dominated systems with the choice
$j_{1}=j_{2}=1,j_{3}=0.1,\rho_{1}=3,z_{1}=1$, $j_{4}=-j_{3}$, and
$j_{\text{glc}}=1$.

It is seen from the results obtained that

  * •

The total energy $W$ increases linearly with increasing distance between the
static charges; this is analogous to the behavior of the energy of interaction
between two quarks in QCD when the distance between them increases. At the
same time, the flux of the non-Abelian electric field $\vec{E}^{7}$ through
the plane $z=0$ first increases (modulus) with increasing distance between the
charges, and for large $L$ changes only slightly.

  * •

The behavior of the chromoelectric fields $\vec{E}^{7}$ and $\vec{E}^{2}$ is
fundamentally different. The field $\vec{E}^{7}$ creates a flux directed from
one charge to another, and its $z-$component has both the nonlinear, $fv/2$,
and gradient, $-h_{,z}$, terms directed in different directions and creating a
total nonvanishing flux of this field through the plane $z=0$. For the field
$\vec{E}^{2}$, the situation differs in principle: the flux of this field
through the plane $z=0$ is equal to zero and both components (nonlinear,
$-hv/2$, and gradient, $-f_{,z}$) are equal to zero on the plane $z=0$. Such a
behavior is analogous to the behavior of force lines between charges with
opposite signs in the first case and between charges with the same sign in
Maxwell’s electrodynamics.

  * •

The electric fields $\vec{E}^{2,7}$ and the magnetic field $\vec{H}^{5}$ have
different color indices; this can be written as $E^{a}_{i}H^{a}_{j}=0$ (no
summation over $a$).

  * •

The profile of the longitudinal electric field $E^{7}_{z}$ obtained here has
one node (see Fig. 3) and it differs from the nodeless profiles of the
longitudinal electric field obtained using lattice calculations (see, e.g.,
Ref. [7] where the color electric field is also sourced by static charges –
quarks). One can naively expect that this difference may disappear for another
choice of charge/current densities (26), (27), (31), and (32).

  * •

The force lines of the color magnetic field create circles with a center
coinciding with that of the flux tube.

  * •

The profile of force lines for the electric field $\vec{E}^{7}$ (see Fig. 1)
shows that there are closed force lines, i.e. $\curl{\vec{E}^{7}}\neq 0$. This
is related to the fact that according to (5)
$\vec{E}^{7}=f\vec{A}^{5}/2-\gradient{h}$ and consequently
$\curl{\vec{E}^{7}}=f\vec{H}^{5}/2+\gradient{f}\times\vec{A}^{5}$. Thus the
presence of nonzero magnetic field $\vec{H}^{5}$ results in the appearance of
nonzero $\curl{\vec{E}^{7}}\neq 0$ (the effect related to the nonlinearity of
the non-Abelian field).

  * •

It is seen from Fig. 2 that the flux of the electric field through the plane
$z=0$ changes slowly for quite large distances between the charges.
Apparently, this leads to the appearance of the practically linear dependence
of the total energy on the distance.

  * •

The spatial distribution of the energy density of the systems with the
electric field dominance differs considerably from that of the configurations
with the magnetic field dominance by the presence of two distinct maxima
caused by a dominating contribution coming from the charge density compared
with the current density.

Figure 3: The typical profiles of the longitudinal component $E^{7}_{z}$ of
the electric field in the plane $z=0$ for the configurations with the magnetic
field dominance (shown in two top rows of Fig. 1) and with the electric field
dominance (shown in two bottom rows of Fig. 1).

It is interesting to note that, due to the nonlinearity of the system, even
when the current $j^{5}$ is completely “turned off”, there are configurations
with the flux of the longitudinal electric field which have the magnetic field
that is much smaller than the electric field.

To conclude this section, let us estimate the magnitude of the Proca field
mass $m$ that is needed to obtain a typical value of the tension coefficient
of the string between quarks obtained in lattice calculations. Since we
consider static charges, one can regard the total energy of the system $W$
calculated above as a corresponding interaction potential of the quarks that
is calculated in lattice calculations in QCD. Then, for example, it is seen
from the top left panel of Fig. 2 that for the magnetic-field-dominated
systems our potential increases linearly with distance between charges in the
region $L\gtrsim 2$: $W\sim KL$. On changing to dimensional units, this
expression takes the form [see Eq. (38)]:

$W_{\text{NU}}=K_{\text{NU}}l,$

where $l=L/m$ and $K_{\text{NU}}$ is a string tension coefficient (measured in
natural units). It is known from lattice calculations that
$K_{\text{NU}}\approx 0.2\,\text{GeV}^{2}$ [3]. Then, to obtain this value of
$K_{\text{NU}}$ for the systems shown in the top left panel of Fig. 2 (where
$K\approx 0.38$), it is necessary to choose $m\approx 0.73\,\text{GeV}$. In
turn, the dimensional distance between charges will then be determined by the
expression

$l\approx 0.27\times 10^{-13}\cdot L\,\text{cm}.$

This quantity is in good agreement with a typical distance between quarks
obtained in lattice calculations [3].

##  V Discussion

The results obtained tell us that in non-Abelian Yang-Mills-Proca theory with
an external field source in the form of ad hoc currents there are flux tubes
filled with a longitudinal color electric field and having a total energy that
increases linearly with increasing distance between static charges.

Since similar behavior is inherent in a tube connecting quarks and threaded by
a longitudinal electric field, it is of great interest to see the reasons of
such a behavior of non-Abelian electric Proca fields. This would enable one to
suggest some ideas concerning a possible relationship between these two
theories (for instance, Proca theory could serve as an approximate description
of QCD in the situation under consideration). These fields are sourced by the
charge and current densities (26), (27), (31), and (32). The charge density
$j^{2}_{t}$ is concentrated at the points $z=\pm z_{0}$; this corresponds to
the fact that at these points two charges are located. The static charge
density $j^{7}_{t}$ is distributed over the whole space. The color current
density $j^{5}_{u}$, as well as the charge density $j^{2}_{t}$, is
concentrated at the points $z=\pm z_{0}$. In contrast, the current density
$j^{5}_{z}$ is distributed over the space.

In such a situation, let us emphasize that there are charge/current densities
concentrated near some center that corresponds to a location of a color
charge/current; there are charge/current densities distributed over the space
(this is especially evident for the current density $j^{5}_{z}$). This can be
written in the form

$j^{a}_{\mu}=\left(j^{a}_{\mu}\right)_{1}+\left(j^{a}_{\mu}\right)_{2}.$

Here the four-current $\left(j^{a}_{\mu}\right)_{1}$ is concentrated near the
center of location of static charge/current (“quarks”), while the four-current
$\left(j^{a}_{\mu}\right)_{2}$ is an extended object. Physically, this means
that the four-current $\left(j^{a}_{\mu}\right)_{1}$ is created by a
charge/current localized near some center, while the current
$\left(j^{a}_{\mu}\right)_{2}$ is created by some extended object which, as
will be shown below, is related to a gluon condensate.

Our calculations indicate that in non-Abelian Proca theory there are flux
tubes threaded by a longitudinal electric field and having a total energy that
increases with increasing distance between static charges creating such a
tube. Such result enables one to suppose that non-Abelian Proca theory under
consideration with sources in the form of static charges is some approximation
for QCD in describing a flux tube. In this case the four-current
$\left(j^{a}_{\mu}\right)_{1}$ is created by quarks, while the four-current
$\left(j^{a}_{\mu}\right)_{2}$ is an approximate description of a current
created by a gluon condensate.

The appearance of such a condensate may be approximately imagined as follows.
In a quantum-averaged Lagrangian of SU(3) Yang-Mills theory, there is the
following term:

$\left\langle
f^{abc}\hat{A}^{b}_{\mu}\hat{A}^{c}_{\nu}f^{amn}\hat{A}^{m\mu}\hat{A}^{n\nu}\right\rangle.$
(40)

Here $\left\langle\cdots\right\rangle$ denotes the quantum averaging;
$\hat{A}^{a}_{\mu}(a,b,c=2,5,7)$ belongs to the sub-group SU(2):
$\hat{A}^{a}_{\mu}\in\text{SU(2)}\subset\text{SU(3)}$ and
$\hat{A}^{m}_{\mu}\in\text{SU(3)/SU(2)}$ with $m=1,3,4,6,8$. The quantized
gauge potentials $\hat{A}^{a}_{\mu}$ and $\hat{A}^{m}_{\mu}$ behave
differently:

$\left\langle\hat{A}^{a}_{\mu}\right\rangle\approx
A^{a}_{\mu},\left\langle\hat{A}^{m}_{\mu}\right\rangle=0\text{ but
}\left\langle\hat{A}^{m}_{\mu}\hat{A}^{m\mu}\right\rangle\neq 0.$ (41)

This enables one to say that the four-potential $\hat{A}^{a}_{\mu}$ is almost
classical, while $\hat{A}^{m}_{\mu}$ is purely quantum. In this case the
expression (40) will approximately look like

$\left\langle
f^{abc}\hat{A}^{b}_{\mu}\hat{A}^{c}_{\nu}f^{amn}\hat{A}^{m\mu}\hat{A}^{n\nu}\right\rangle\approx
f^{abc}A^{b}_{\mu}A^{c}_{\nu}\left\langle
f^{amn}\hat{A}^{m\mu}\hat{A}^{n\nu}\right\rangle.$

Notice the following fact. On account of the antisymmetry of the structure
constants $f^{abc}$, the expression $f^{amn}\hat{A}^{m\mu}\hat{A}^{n\nu}$ can
be written in the form:

$f^{amn}\hat{A}^{m\mu}\hat{A}^{n\nu}=\frac{1}{2}f^{amn}\left(\hat{A}^{m\mu}\hat{A}^{n\nu}-\hat{A}^{n\mu}\hat{A}^{m\nu}\right)=\frac{1}{2}f^{amn}\left[\hat{A}^{m\mu},\hat{A}^{n\nu}\right]_{m,n},$
(42)

where $[\cdots]_{m,n}$ is a commutator for indices $m,n$. For free
noninteracting fields similar commutator

$\left[\hat{A}^{m\mu}(t,x^{i}),\hat{A}^{n\nu}(t,y^{j})\right]=0$

for the points $x^{i},y^{j}$ separated by a spacelike interval. For strongly
interacting fields the commutator (42) can be nonzero, since for them there
can exist a static regular distribution of the field in space. This means that

$\left[\hat{A}^{m}_{\mu}(t,x^{i}),\hat{A}^{n}_{\nu}(t,y^{j})\right]_{m,n}=D^{mn}_{\mu\nu}\left(x^{i},y^{j}\right)\neq
0.$

Thus in a quantum-averaged Lagrangian of SU(3) Yang-Mills theory there appears
the following term:

$f^{abc}f^{amn}A^{b}_{\mu}A^{c}_{\nu}D^{mn}_{\mu\nu}\left(x^{i},x^{j}\right).$

In our opinion this term may describe a gluon condensate created by the
quantum components $\hat{A}^{m}_{\mu}$ of a non-Abelian gauge field. Note also
that if the asymptotic behavior of the commutator is

$D^{mn}_{\mu\nu}\left(x^{i}\rightarrow\infty,x^{j}\rightarrow\infty\right)\rightarrow
m^{2},$ (43)

then the Proca mass $m$ is a purely quantum phenomenon appearing in a strongly
nonlinear regime of QCD in nonperturbative quantization (the nonperturbative
generation of the Proca mass). In this connection, note that Ref. [31]
discusses the question of the massless limit in Proca theory, and it is shown
that the apparent discontinuity in the massless limit is only an artefact of
the perturbation theory.

The above discussion justifies the assumption that the four-current
$\left(j^{a}_{\mu}\right)_{1}$ is created by static quarks, while the current
$\left(j^{a}_{\mu}\right)_{2}$ is created by a gluon condesate, i.e., it is a
purely quantum object created by the quantum components $\hat{A}^{m}_{\mu}$.
Note that this assumption, as well as the assumption (41) concerning the
different behavior of different components of the gauge potential
$A^{B}_{\mu}$ (with $B=1,\ldots,8$) and the assumption (43) concerning the
nature of the Proca mass $m$, can be verified using lattice calculations.

A very important and interesting question is that of gauge invariance in the
case under consideration. The results concerning the flux tube were obtained
in non-Abelian SU(2) Yang-Mills-Proca theory that is clearly not a gauge-
invariant theory. In this section we consider an assumption that the present
Proca theory is some approximation for a real quantum SU(3) theory. In such a
case the divergent equation (15) can be thought of as a choice of the gauge in
initial SU(3) theory where a contribution coming from the gluon condensate
created by the fields $\hat{A}^{m}_{\mu}\in\text{SU(3)/SU(2)}$ is taken into
account. This enables one to assume that for such a nonpertubative
quantization the gauge invariance is conserved, but, of course, this
assumption requires a more careful analysis.

The discussion presented in this section enables one to assume that an
approximate description of a flux tube can be made using classical equations
within Yang-Mills-Proca theory where one takes into account the fact that
these fields are sourced not only by color charges but also by a gluon
condensate of quantum fields with zero average value and with nonzero
dispersion.

##  VI Comparison with the results of lattice calculations

Using lattice calculations, it was shown in Ref. [32] that the simulated
chromoelectric field $\vec{E}$ consists of the “nonperturbative”,
$\vec{E}^{\text{NP}}$, and “Coulomb”, $\vec{E}^{\text{C}}$, fields:

$\vec{E}=\vec{E}^{\text{NP}}+\vec{E}^{\text{C}}.$

In order to compare this expression with the expressions for electric and
magnetic fields, let us rewrite the expressions (5) and (6) for these fields
on account of a possible relationship between quantum SU(3) gauge theory and
its approximate description in the form of SU(2) Proca theory discussed in
Sec.V:

$\displaystyle\vec{E}^{2}=$
$\displaystyle-\frac{h\vec{A}^{5}}{2}-\frac{\gradient{f}}{g}+gf_{2mn}\left\langle\hat{A}^{m}_{t}\hat{\vec{A}}^{n}\right\rangle=\left(\vec{E}^{2}\right)_{\text{nl}}+\left(\vec{E}^{2}\right)_{\gradient}+\left(\vec{E}^{2}\right)_{\text{gc}}\quad\text{with}\quad
m,n=1,3,4,6;$ (44) $\displaystyle\vec{E}^{7}=$
$\displaystyle\frac{f\vec{A}^{5}}{2}-\frac{\gradient{h}}{g}+gf_{7mn}\left\langle\hat{A}^{m}_{t}\hat{A}^{n}_{z}\right\rangle=\left(\vec{E}^{7}\right)_{\text{nl}}+\left(\vec{E}^{7}\right)_{\gradient}+\left(\vec{E}^{7}\right)_{\text{gc}}\quad\text{with}\quad
m,n=1,3,4,6,8;$ (45) $\displaystyle H^{5}_{\varphi}=$
$\displaystyle\rho\frac{\left(v_{,\rho}-u_{,z}\right)}{g}+g\rho
f_{5mn}\left\langle\hat{A}^{m}_{\rho}\hat{A}^{n}_{z}\right\rangle=\left(H^{5}_{\varphi}\right)_{\gradient}+\left(H^{5}_{\varphi}\right)_{\text{gc}}\quad\text{with}\quad
m,n=1,3,4,6,8.$ (46)

Here the first terms on the right-hand sides of Eqs. (44) and (45) are
nonlinear terms, the second terms are gradient terms, and the third terms are
associated with the presence of a gluon condensate [recall that the indices
$2,5,7\in\text{SU(2)}\subset\text{SU(3)}$ and $m,n\in\text{SU(3)/SU(2)}$.
Similarly, for the magnetic field (46)].

It is seen that the “nonperturbative” component $\vec{E}^{\text{NP}}$ from
lattice calculations coincides with the nonlinear plus gluon terms of the
total vector of the chromoelectric field $\vec{E}^{7}$,
$\left(\vec{E}^{7}\right)_{\text{nl}}+\left(\vec{E}^{7}\right)_{\text{gc}}$,
and correspondingly the “Coulomb” field $\vec{E}^{\text{C}}$ coincides with
the gradient term $\left(\vec{E}^{7}\right)_{\gradient}$.

In Ref. [32], it was also shown that the magnetic field obtained using lattice
calculations is negligibly small compared with the electric field. It must be
mentioned that in those calculations the comparison was made for magnetic and
electric fields possessing the same color index. In our calculations we deal
with the similar situation: there are nonvanishing color electric fields
$\vec{E}^{2,7}$ and color magnetic fields having the same color index but
equal to zero, $\vec{H}^{2,7}=0$. On the other hand, in our calculations there
is nonzero magnetic field $\vec{H}^{5}\neq 0$ and there is no a corresponding
electric field: $\vec{E}^{5}=0$. This situation can be described as
$E^{a}_{i}H^{a}_{j}=0$ (no summation over $a$). Notice that this result may be
verified using lattice calculations.

As noted in the previous section, in lattice calculations the profile of the
longitudinal chromoelectric field differs from that of obtained by us in that
in our case there is one node. This is related to the fact that there is a
nonzero curl of the field
$\vec{E}^{7}=\left(\vec{E}^{7}\right)_{\text{nl}}+\left(\vec{E}^{7}\right)_{\gradient}$
used in our calculations:
$\curl{\vec{E}^{7}}=f\vec{H}^{5}/2+\gradient{f}\times\vec{A}^{5}$. Due to the
presence of the quite large chromomagnetic field $\vec{H}^{5}$ there appear
closed force lines of the field $\vec{E}^{7}$; this leads to the appearance
both of positive and of negative values of $E^{7}_{z}$. In turn, in lattice
calculations this magnetic field is small, and therefore there are no closed
force lines of the field $\vec{E}^{7}$. Naively we may assume that this
situation with one node can be changed by choosing profiles of charge/current
densities different from those of given by Eqs. (26), (27), (31), and (32).

##  VII Conclusions

Summarizing the results obtained,

  * •

Within non-Abelian Yang-Mills-Proca theory, we have shown that there are
regular solutions supported by static charges and currents. These solutions
describe a flux tube threaded by a color longitudinal electric field with a
color index $a=2,7$, and force lines of a color magnetic field with another
color index $a=5$ are wound on the tube.

  * •

These flux tubes have a linear dependence of their total energy on the
distance between static charges; this enables us to say about the effect
similar to the confinement phenomenon in QCD.

  * •

The above dependence of the total energy on the distance between charges is a
consequence of the presence of the Proca mass which ensures an exponential
decrease of the fields at infinity, as well as because of the nonlinear
interaction between fields inherent in Yang-Mills theory.

  * •

The flux of the electric field through the plane $z=0$ depends strongly on the
length of the tube: when the length increases, the flux becomes a slowly
varying quantity.

  * •

Our calculations indicate that there is a longitudinal chromoelectric field
that ensures the confinement of the charges.

  * •

It is shown that the results obtained are in good agreement with the results
of lattice calculations.

  * •

The question of a possible relationship between non-Abelian Yang-Mills-Proca
theory employed here and quantum chromodynamics is considered.

As a result, we have obtained the distribution of chromoelectric and
chromomagnetic fields bearing a fairly close resemblance to the properties of
the same fields that are needed for the confinement of quarks in QCD. We have
discussed this interesting fact in Sec. V, where a possible scenario for a
situation when non-Abelian Proca theory may serve as an approximate
description of QCD has been considered.

In future studies it would be interesting, by comparing with lattice
calculations, to obtain more realistic profiles of the currents which could
provide better agreement of our model with lattice calculations. For example,
one might verify the existence of a chromomagnetic field $\vec{H}^{5}$ which
is perpendicular to the chromoelectric field $\vec{E}^{2,7}$
($\vec{H}^{5}\perp\vec{E}^{2,7}$) in color space in the sense that
$H^{a}_{i}E^{a}_{j}=0$ (no summation over $a$).

## References

  * [1] A. Di Giacomo, M. Maggiore, and S. Olejnik, “Confinement and Chromoelectric Flux Tubes in Lattice QCD,” Nucl. Phys. B 347, 441 (1990).
  * [2] G. S. Bali, K. Schilling, and C. Schlichter, “Observing long color flux tubes in SU(2) lattice gauge theory,” Phys. Rev. D 51, 5165 (1995).
  * [3] Y. A. Simonov, “Confinement,” Phys. Usp. 39, 313 (1996).
  * [4] P. Cea, L. Cosmai, and A. Papa, “Chromoelectric flux tubes and coherence length in QCD,” Phys. Rev. D 86, 054501 (2012).
  * [5] N. Cardoso, M. Cardoso, and P. Bicudo, “Inside the SU(3) quark-antiquark QCD flux tube: screening versus quantum widening,” Phys. Rev. D 88, 054504 (2013).
  * [6] M. Baker, P. Cea, V. Chelnokov, L. Cosmai, F. Cuteri, and A. Papa, “Isolating the confining color field in the SU(3) flux tube,” Eur. Phys. J. C 79, 478 (2019).
  * [7] M. Baker, V. Chelnokov, L. Cosmai, F. Cuteri, and A. Papa, “The flux tube profile in full QCD,” PoS LATTICE2021, 355 (2022).
  * [8] R. Brito, V. Cardoso, C. A. R. Herdeiro, and E. Radu, “Proca stars: Gravitating Bose-Einstein condensates of massive spin 1 particles,” Phys. Lett. B 752, 291 (2016).
  * [9] C. A. R. Herdeiro, A. M. Pombo, and E. Radu, “Asymptotically flat scalar, Dirac and Proca stars: discrete vs. continuous families of solutions,” Phys. Lett. B 773, 654 (2017).
  * [10] V. Dzhunushaliev and V. Folomeev, “Dirac star in the presence of Maxwell and Proca fields,” Phys. Rev. D 99, 104066 (2019).
  * [11] C. Herdeiro, I. Perapechka, E. Radu, and Y. Shnir, “Asymptotically flat spinning scalar, Dirac and Proca stars,” Phys. Lett. B 797, 134845 (2019).
  * [12] V. Dzhunushaliev and V. Folomeev, “Dirac star with SU(2) Yang-Mills and Proca fields,” Phys. Rev. D 101, 024023 (2020).
  * [13] J. C. Bustillo, N. Sanchis-Gual, A. Torres-Forné, J. A. Font, A. Vajpeyi, R. Smith, C. Herdeiro, E. Radu, and S. H. W. Leong, “GW190521 as a Merger of Proca Stars: A Potential New Vector Boson of $8.7\times 10^{-13}$ eV,” Phys. Rev. Lett. 126, 081101 (2021).
  * [14] V. Dzhunushaliev and V. Folomeev, “Axially symmetric Proca-Higgs boson stars,” Phys. Rev. D 104, 104024 (2021).
  * [15] L. Heisenberg, R. Kase, M. Minamitsuji, and S. Tsujikawa, “Hairy black-hole solutions in generalized Proca theories,” Phys. Rev. D 96, 084049 (2017).
  * [16] R. Kase, M. Minamitsuji, S. Tsujikawa, and Y. L. Zhang, “Black hole perturbations in vector-tensor theories: The odd-mode analysis,” J. Cosmol. Astropart. Phys. 02 (2018) 048.
  * [17] N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer, and N. Weiner, “A theory of dark matter,” Phys. Rev. D 79, 015014 (2009).
  * [18] M. Pospelov and A. Ritz, “Astrophysical signatures of secluded dark matter,” Phys. Lett. B 671, 391 (2009).
  * [19] A. De Felice, L. Heisenberg, R. Kase, S. Mukohyama, S. Tsujikawa, and Y. l. Zhang, “Cosmology in generalized Proca theories,” J. Cosmol. Astropart. Phys. 06 (2016) 048.
  * [20] A. De Felice, L. Heisenberg, R. Kase, S. Mukohyama, S. Tsujikawa, and Y. l. Zhang, “Effective gravitational couplings for cosmological perturbations in generalized Proca theories,” Phys. Rev. D 94, 044024 (2016).
  * [21] L. C. Tu, J. Luo, and G. T. Gillies, “The mass of the photon,” Rep. Prog. Phys. 68, 77 (2005).
  * [22] A. De Felice, L. Heisenberg, R. Kase, S. Tsujikawa, Y. l. Zhang, and G. B. Zhao, “Screening fifth forces in generalized Proca theories,” Phys. Rev. D 93, 104016 (2016).
  * [23] V. Dzhunushaliev and V. Folomeev, “Proca tubes with the flux of the longitudinal chromoelectric field and the energy flux/momentum density,” Eur. Phys. J. C 80, 1043 (2020).
  * [24] V. Dzhunushaliev, V. Folomeev, T. Kozhamkulov, A. Makhmudov, and T. Ramazanov, “Non-Abelian Proca theories with extra fields: particlelike and flux tube solutions,” Phys. Scr. 95, 074013 (2020).
  * [25] V. Dzhunushaliev, V. Folomeev, and A. Tlemisov, “Linear energy density and the flux of an electric field in Proca tubes,” Symmetry 13, 640 (2021).
  * [26] V. Dzhunushaliev and V. Folomeev, “Axially symmetric particlelike solutions with the flux of a magnetic field in the non-Abelian Proca-Higgs theory,” Phys. Rev. D 104, 116027 (2021).
  * [27] V. Dzhunushaliev and V. Folomeev, “Proca balls with angular momentum or flux of electric field,” Phys. Rev. D 105, 016022 (2022).
  * [28] N.I.M. Gould, J.A. Scott, and Y. Hu, “A numerical evaluation of sparse direct solvers for the solution of large sparse symmetric linear systems of equations,” ACM Trans. Math. Softw. 33, 10 (2007).
  * [29] O. Schenk and K. Gartner, “Solving unsymmetric sparse systems of linear equations with PARDISO,” Future Gener. Comput. Syst. 20, 475 (2004).
  * [30] J. Kunz, I. Perapechka, and Y. Shnir, “Kerr black holes with synchronised scalar hair and boson stars in the Einstein-Friedberg-Lee-Sirlin model,” JHEP 07, 109 (2019).
  * [31] A. Hell, “The strong couplings of massive Yang-Mills theory,” JHEP 03, 167 (2022).
  * [32] M. Baker, P. Cea, V. Chelnokov, L. Cosmai, F. Cuteri, and A. Papa, “The confining color field in SU(3) gauge theory,” Eur. Phys. J. C 80, 514 (2020).