text:
# 1 Introduction

## 1 Introduction

Deformation quantization which is also known as Weyl–Wigner–Groenewold–Moyal quantization describes a quantum system in terms of c–number phase space variables. The correspondence with the operator quantization is due to symbol maps and star products. Operators are mapped into c–number functions however, their composition is given by star product which is noncommutative but associative. One of the most important ingredients of deformation quantization is Wigner function. It is essential to define phase space integrals corresponding to expectation values of operator formalism.

When classical coordinates are converted into noncommuting variables they mutate into operators. Thus after performing canonical quantization of noncommutative systems one should deal with two kind of operators. This can be avoided by formulating noncommutativity of coordinates as deformation of commuting ones. One can adopt all the machinery developed in deformation quantization.

An electron moving on plane in a uniform external magnetic field orthogonal to plane is known as the Landau problem. Different aspects of the Landau problem in noncommutative coordinates have been widely studied. Our aim is to derive Wigner functions for the Landau problem for two different noncommuting spaces: When only coordinates are noncommuting and then letting coordinates as well as momenta are noncommuting.

In Section 2 we briefly recall definitions of Wigner function. In Section 3 we study the Landau problem on the noncommutative plane by adopting a differential representation for quantum mechanical operators. This yields Wigner functions which are formal. However, guided by this formulation we define an effective Hamiltonian and obtain Wigner functions utilizing the usual deformation quantization. When both quantum mechanics and noncommutativity are formulated as deformations, it is plausible to consider all of phase space variables on the same footing. Indeed, in Section 4 deformation quantization is adopted by letting coordinates as well as momenta are noncommuting for the Landau problem and the related Wigner functions are calculated. Moreover, we show that noncommutativity relations can be chosen such that the resultant Wigner functions are not aware of noncommutativity of phase space variables. These Wigner functions can be used to calculate some physical quantities in noncommutative spaces.

## 2 Wigner Functions

Wigner functions can be defined equivalently in two different ways. One of them refers to the operator formulation of quantum mechanics: Let the Hamiltonian operator be $`\widehat{}`$ and define the eigenvalue problem

$$\widehat{}\psi (\stackrel{}{x})=\psi (\stackrel{}{x}).$$
(1)
If we denote its spectrum by $`_l`$ and the corresponding eigenstates by $`\psi _l(\stackrel{}{x}),`$ (diagonal) Wigner functions are defined as

$$W_l(\stackrel{}{x},\stackrel{}{\pi })=\frac{1}{2\pi }d^ds\psi _l^{}(\stackrel{}{x}\frac{\mathrm{}}{2}\stackrel{}{s})e^{i\stackrel{}{s}\stackrel{}{\pi }}\psi _l(\stackrel{}{x}+\frac{\mathrm{}}{2}\stackrel{}{s}),$$
(2)
where $`(\stackrel{}{x},\stackrel{}{\pi })`$ are $`2d`$ dimensional phase space variables.

Wigner functions can also be defined without referring to the operator formalism of quantum mechanics (1), by introducing the $`\mathrm{}`$–star product

$$_{\mathrm{}}\mathrm{exp}\frac{i\mathrm{}}{2}(\stackrel{}{}_\stackrel{}{x}\stackrel{}{}_\stackrel{}{\pi }\stackrel{}{}_\stackrel{}{\pi }\stackrel{}{}_\stackrel{}{x}).$$
(3)
Now, the Wigner functions are defined as $`_{\mathrm{}}`$–genvalue equation

$$_{\mathrm{}}W_l(\stackrel{}{x},\stackrel{}{\pi })=_lW_l(\stackrel{}{x},\stackrel{}{\pi }),$$
(4)
where $``$ is the classical Hamiltonian function corresponding to $`\widehat{}.`$

These two definitions are equivalent. Wigner function are used to define integrals in the deformation quantization corresponding to expectation values in operator formulation.

## 3 Wigner Functions for the Landau Problem in Noncommutative Coordinates

An electron moving in the plane which is perpendicular to the external uniform magnetic field $`B,`$ neglecting spin, is described by the Hamiltonian

$$H=\frac{1}{2m}(\stackrel{}{p}+\frac{e}{c}\stackrel{}{A})^2.$$
(5)
where we adopt the symmetric gauge

$$\stackrel{}{A}=(\frac{B}{2}y,\frac{B}{2}x).$$
(6)
Consider the coordinates satisfying

$$[\stackrel{~}{x},\stackrel{~}{y}]=i\theta ,$$
(7)
where $`\theta `$ is a constant parameter. Because of being operators even at the classical level, we denoted noncommuting coordinates as $`\stackrel{~}{x},\stackrel{~}{y}.`$ Noncommutativity can be imposed by treating the coordinates as commuting, but introducing the $`\theta `$–star product

$$_\theta \mathrm{exp}\frac{i\theta }{2}(\stackrel{}{}_x\stackrel{}{}_y\stackrel{}{}_y\stackrel{}{}_x).$$
(8)
Now, we deal with the commutative coordinates $`x`$ and $`y`$ but composition of functions is given by $`\theta `$–star product (8). For example, instead of the commutator (7) one introduces the Moyal bracket

$$x_\theta yy_\theta x=i\theta .$$
(9)
We quantize this system in terms of the standard canonical quantization by establishing the usual canonical commutation relations

$$[r_i,\widehat{p}_j]=i\mathrm{}\delta _{ij}$$
(10)
and adopting the usual differential representation

$$\widehat{p}_i=i\mathrm{}\frac{}{\delta r_i}.$$
(11)
Note that we use the same notation for the (quantum) operator and classical coordinates $`\stackrel{}{r}=(x,y)`$.

One of the receipts to define quantum mechanics in noncommutative coordinates is to interpret derivatives as momentum operators: $`\stackrel{}{}=\frac{i}{\mathrm{}}\widehat{\stackrel{}{p}}.`$ According to this receipt, in the symmetric gauge (6) the Hamiltonian operator $`\widehat{H}`$ corresponding to (5), acting on an arbitrary function $`\mathrm{\Psi }(\stackrel{}{r})`$ yields

$$\widehat{H}_\theta \mathrm{\Psi }(\stackrel{}{r})=\frac{1}{2m}\left[\left(\widehat{p}_x\frac{eB}{2c}y\right)^2+\left(\widehat{p}_y+\frac{eB}{2c}x\right)^2\right]_\theta \mathrm{\Psi }(\stackrel{}{r})\widehat{H}_{nc}\mathrm{\Psi }(\stackrel{}{r}),$$
(12)
where, in terms of $`\kappa =\frac{eB\theta }{4c\mathrm{}},`$ we defined

$$\widehat{H}_{nc}=\frac{1}{2m}\left[\left((1+\kappa )\widehat{p}_x\frac{eB}{2c}y\right)^2+\left((1+\kappa )\widehat{p}_y+\frac{eB}{2c}x\right)^2\right].$$
(13)
The eigenvalue equation

$$\widehat{H}_{nc}\varphi (\stackrel{}{r})=E\varphi (\stackrel{}{r}),$$

can be solved in polar coordinates and its eigenstates can be given in terms of the Laguerre polynomials in $`(x^2+y^2).`$ These solutions are not suitable to find the related Wigner functions although, formally we can write them by making use of the definition (2). Another method may be to use definition of Wigner functions in polar coordinates.

We would like to use the fact that effectively the Landau problem is one dimensional harmonic oscillator. The operator

$$a=\frac{1}{\sqrt{2m\mathrm{}\stackrel{~}{\omega }}}\left[\left(\widehat{p}_x\frac{m\stackrel{~}{\omega }}{2}y\right)i\left(\widehat{p}_y+\frac{m\stackrel{~}{\omega }}{2}x\right)\right],$$
(14)
and its hermitian conjugate $`a^{}`$ can be shown to satisfy

$$[a,a^{}]=1,$$
(15)
where $`\stackrel{~}{\omega }(1+\kappa )\omega `$ and $`\omega =\frac{eB}{mc}`$ is the cyclotron frequency.

The Hamiltonian becomes

$$\widehat{H}=\mathrm{}\stackrel{~}{\omega }(a^{}a+\frac{1}{2}),$$
(16)
which is equivalent to harmonic oscillator in one dimension whose frequency is $`\stackrel{~}{\omega }.`$ Thus, its spectrum is given by

$$E_n=\mathrm{}\stackrel{~}{\omega }(n+1/2),$$

where $`n=0,1,\mathrm{}.`$ In terms of the ground state

$$\varphi _0(\stackrel{}{r},\theta )=\frac{m\stackrel{~}{\omega }}{2\pi \mathrm{}}\mathrm{exp}[\frac{m\stackrel{~}{\omega }}{4\mathrm{}}(x^2+y^2)],$$

its eigenfunctions are given as

$$\varphi _n(\stackrel{}{r},\theta )=\frac{1}{\sqrt{n!}}(a^{})^n\varphi _0.$$

These states can be used to write the Wigner functions as

$$W_{L,n}(\stackrel{}{r},\stackrel{}{p},\theta )=\frac{1}{4\pi ^2}d^2v\varphi _n^{}(\stackrel{}{r}\frac{\mathrm{}}{2}\stackrel{}{v})e^{i\stackrel{}{v}\stackrel{}{p}}\varphi _n(\stackrel{}{r}+\frac{\mathrm{}}{2}\stackrel{}{v}).$$

However, as it is mentioned above these Wigner functions are only formal, few of the integrals can be calculated under some conditions.

Although we are using the specific representation (11), let us assume that we can choose the realization

$$a\psi (q)=\frac{1}{\sqrt{2}}(q+\frac{}{q})\psi (q),a^{}\psi (q)=\frac{1}{\sqrt{2}}(q\frac{}{q})\psi (q),$$
(17)
where the parameter $`q`$ is a c–number. Then, eigenfunctions of (16) coincide with the eigenfunctions of one dimensional harmonic oscillator, which lead to Wigner functions

$$W_{0,n}(P,q)=\frac{(1)^n}{\pi \mathrm{}}e^zL_n(2z),$$
(18)
where we used the variable

$$z=(P^2+q^2).$$

$`P`$ is the parameter appearing in the definition (2).

As we have already mentioned in Section 2, there exists another method of calculating Wigner functions: We can use the $`\mathrm{}`$–star product (3) and the definition (4). However, in this case we should define a c–number object corresponding to $`\widehat{H}_{nc}.`$ Let us define the effective c–number Hamiltonian

$$H_{eff}=\frac{(1+\kappa )^2}{2m}\stackrel{}{p}^2+\frac{m\omega ^2}{8}\stackrel{}{r}^2+\frac{\stackrel{~}{\omega }}{2}(xp_yyp_x).$$
(19)
We would like to emphasize the fact that it is not classical because $`\kappa `$ depends on $`\mathrm{}.`$ Now, to find the related Wigner functions we use the definition

$$H_{eff}_{\mathrm{}}W_1(\stackrel{}{p},\stackrel{}{r},\theta )=EW_1(\stackrel{}{p},\stackrel{}{r},\theta ).$$
(20)
This can be written as a differential equation only in one variable

$$\rho H_{eff}$$

as

$$\rho \frac{d^2W_1(\rho )}{d\rho ^2}+\frac{dW_1(\rho )}{d\rho }\frac{4}{\mathrm{}^2\stackrel{~}{\omega }^2}(\rho +E)W_1(\rho )=0.$$
(21)
To solve this equation let us change the variable as

$$R=\frac{4\rho }{\mathrm{}\stackrel{~}{\omega }}$$
(22)
and write

$$W_1=e^{R/2}L(R).$$
(23)
Thus (21) becomes

$$\left[R\frac{d^2}{dR^2}+(1R)\frac{d}{dR}(\frac{1}{2}\frac{E}{\mathrm{}\stackrel{~}{\omega }})\right]W_1=0,$$
(24)
whose solutions are Laguerre polynomials $`L_n(4R)`$ where $`n=0,1,\mathrm{}`$ and

$$E_n=\mathrm{}\stackrel{~}{\omega }(n+1/2).$$

Thus we have

$$W_{1,n}(\stackrel{}{p},\stackrel{}{r},\theta )=e^{2H_{eff}/\mathrm{}\stackrel{~}{\omega }}L_n(4H_{eff}/\mathrm{}\stackrel{~}{\omega }).$$
(25)
Indeed, this is not surprising: The system is equivalent to one dimensional harmonic oscillator, which also can be observed at the “classical level”. By performing the canonical transformation

$`p_1={\displaystyle \frac{(1+\kappa )}{\sqrt{m}}}p_x{\displaystyle \frac{\omega \sqrt{m}}{2}}y,q_1={\displaystyle \frac{(1+\kappa )}{\sqrt{m}}}p_y+{\displaystyle \frac{\omega \sqrt{m}}{2}}x,`$ (26)
$`p_2={\displaystyle \frac{(1+\kappa )}{\sqrt{m}}}p_y{\displaystyle \frac{\omega \sqrt{m}}{2}}x,q_2={\displaystyle \frac{(1+\kappa )}{\sqrt{m}}}p_x+{\displaystyle \frac{\omega \sqrt{m}}{2}}y,`$ (27)
the effective Hamiltonian becomes

$$H_{eff}=\frac{1}{2}(p_1^2+q_1^2).$$
(28)
The other set of variables $`(p_2,q_2)`$ leads to the degeneracy of the Landau levels after quantization.

As it may be expected, using the effective Hamiltonian (19) and assuming the representation (17) yield the same Wigner functions.

## 4 Wigner Functions for the Landau Problem in Noncommutative Phase Space

Quantum mechanics in noncommutative space can be defined in terms of Wigner functions found by making use of the $`_{\mathrm{}}`$ and $`_\theta `$ simultaneously. However, when $`_{\mathrm{}}`$–genvalue equations are considered it is plausible to treat $`\stackrel{}{r}`$ and $`\stackrel{}{p}`$ on the same footing. Thus, we would like to deform also momenta of the phase space by introducing

$$_k\mathrm{exp}\frac{ik\theta }{2}(\stackrel{}{}_{p_x}\stackrel{}{}_{p_y}\stackrel{}{}_{p_y}\stackrel{}{}_{p_x}),$$
(29)
where $`k`$ is any constant.

Now, consider

$$_{\mathrm{}}_\theta _{\theta ^\pm },$$

which we use to define the $``$–genvalue equation

$$H\mathrm{\Psi }=E\mathrm{\Psi },$$
(30)
where $`H`$ is the classical Hamiltonian (5) of the Landau problem in the symmetric gauge (6).

A similar approach for harmonic oscillator was given in .

$``$–genvalue equation can be written as

$`\{H{\displaystyle \frac{\mathrm{}^2}{8}}[{\displaystyle \frac{1}{m}}\stackrel{}{}_r\stackrel{}{}_r+{\displaystyle \frac{m\omega ^2}{4}}\stackrel{}{}_p\stackrel{}{}_p\omega (_{p_x}_y_{p_y}_x)]`$
$`{\displaystyle \frac{\theta ^2}{8}}[{\displaystyle \frac{m\omega ^2}{4}}\stackrel{}{}_r\stackrel{}{}_r+{\displaystyle \frac{k^2}{m}}\stackrel{}{}_p\stackrel{}{}_p+k\omega (_{p_x}_y_{p_y}_x]`$
$`+{\displaystyle \frac{\theta \mathrm{}}{4}}[{\displaystyle \frac{\omega }{2}}\stackrel{}{}_r\stackrel{}{}_r+k\stackrel{}{}_p\stackrel{}{}_p+({\displaystyle \frac{\omega ^2m}{4}}+{\displaystyle \frac{k}{m}})(_x_{p_y}_y_{p_x})`$ (31)
$`{\displaystyle \frac{i\mathrm{}}{2}}[{\displaystyle \frac{1}{m}}\stackrel{}{p}\stackrel{}{}_r{\displaystyle \frac{m\omega ^2}{4}}\stackrel{}{r}\stackrel{}{}_p{\displaystyle \frac{\omega }{2}}(y_x_xp_x_{p_y}x_y+p_y_{p_x})]`$
$`+{\displaystyle \frac{i\theta }{2}}[{\displaystyle \frac{m\omega ^2}{4}}(x_yy_x)+{\displaystyle \frac{k}{m}}(p_x_{p_y}p_y_{p_x})+{\displaystyle \frac{\omega }{2}}(p_y_y+y_{p_y}kx_{p_x}kp_x_x)]E\}W_2=0,`$
where $`\stackrel{}{}_r/\stackrel{}{r}`$ and $`\stackrel{}{}_p/\stackrel{}{p}.`$ Introducing the variable

$$\zeta H,$$
(32)
one can show that the imaginary part of (31) vanishes and it reads

$$4\alpha ^2(k,\theta ,\mathrm{})(\zeta \frac{d^2W_2(\zeta )}{d\zeta ^2}+\frac{dW_2(\zeta )}{d\zeta })+\zeta W_2(\zeta )=EW_2(\zeta ),$$
(33)
where we defined

$$\alpha ^2(k,\theta ,\mathrm{})=\frac{1}{16}\mathrm{}^2\omega ^2+\frac{\theta ^2}{16}(\frac{k^2}{m^2}+\frac{\omega ^4m^2}{16}+\frac{k\omega ^2}{2})\frac{1}{2}\mathrm{}\theta (\frac{\omega ^3m}{4}+\frac{k\omega }{m}).$$
(34)
In terms of the variable

$$\eta =\frac{\zeta }{\alpha (k,\theta ,\mathrm{})},$$

(33) can be written as

$$(\eta _{\eta }^{}{}_{}{}^{2}+_\eta \eta /4+\stackrel{~}{E})W_2(\eta )=0,$$

where $`\stackrel{~}{E}=\frac{E}{4\alpha (k,\theta ,\mathrm{})}.`$ By writing

$$W_2(\eta )=e^{\eta /2}L(\eta ),$$

we see that $`L(\eta )`$ satisfies the differential equation

$$[\eta _{\eta }^{}{}_{}{}^{2}+(1\eta )_\eta +\stackrel{~}{E}1/2]L(\eta )=0.$$

Thus $`L(\eta )`$ is given as Laguerre polynomials $`L(\eta )=L_n(4\eta )`$ and

$$E_{k,n}=4\alpha (k,\theta ,\mathrm{})(n+1/2),$$

where as usual $`n=0,1,\mathrm{}.`$ Therefore, Wigner functions of the Landau problem in noncommutative phase space are

$$W_{2,n}(\stackrel{}{p},\stackrel{}{r},k,\theta )=\frac{(1)^n}{\pi \mathrm{}}e^{H/2\alpha (k,\theta ,\mathrm{})}L_n(4H/\alpha (k,\theta ,\mathrm{})).$$
(35)
There exists a very interesting case: When

$$kk_s=\frac{1}{4}\omega ^2m^2$$

$`\alpha (k,\theta ,\mathrm{})`$ becomes independent of $`\theta :`$

$$\alpha ^2(k_s,\theta ,\mathrm{})=\frac{\mathrm{}^2\omega ^2}{16}$$

Thus, noncommutativity relations of phase space can be taken such that the resultant quantum mechanics is not aware of noncommutativity.

One can use the Wigner functions (35) to calculate physical quantities. One of the first applications is to calculate the quantities considered in other frameworks and compare the results. Nevertheless, this approach constitutes an alternative method of defining the Landau problem in noncommutative classical phase space.


format: text

text:
# “Nonbaryonic” Dark Matter as Baryonic Color Superconductor

## 1 Introduction

The presence of large amounts of non-luminous components in the Universe has been known for a long time. In spite of the recent advances in the field ( see e.g. recent summary ), the mystery of the dark matter/energy remains: we still do not know what is it. The main goal of this work is to argue that the dark matter could be nothing but well-known quarks which however are not in the “normal” hadronic phase, but rather in some “exotic”, the so-called color superconducting (CS) phase.

This is a novel phase in QCD when light quarks form the condensate in diquark channels, and it is analogous to Cooper pairs of electrons in ordinary superconductors described by BCS theory. There existence of CS phase in QCD represents our first crucial element for our scenario to work. The study of CS phase received a lot of attention last few years, see original papers, and recent reviews on the subject. It turns out that CS phase is realized when quarks are squeezed to the density which is few times nuclear density. It has been known that this regime may be realized in nature in neutron stars interiors and in the violent events associated with collapse of massive stars or collisions of neutron stars, so it is important for astrophysics. The goal of this work is to argue that such conditions may occur in early universe during the QCD phase transition. Therefore, it might be important for cosmology as well.

The force which squeezes quarks in neutron stars is gravity; the force which does a similar job in early universe during the QCD phase transition is a violent collapse of a bubble formed from the axion domain wall. If number of quarks trapped inside of the bubble (in the bulk) is sufficiently large, the collapse stops due to the internal Fermi pressure. In this case the system in the bulk may reach the critical density when it undergoes a phase transition to CS phase with the ground state being the diquark condensate. These configurations with large number of quarks in color superconducting phase, will be named the QCD balls. Therefore, an existence of the axion domain wall represents our second crucial element for our scenario to work. We should note at this point that the axion field was introduced into the theory to explain the lack of CP violation in the strong interactions. Later on the axion field became one of the favorite candidates for the cold dark matter, see original papers - and recent reviews  on the subject. In the present scenario the axion field plays the role of squeezer rather than dark matter itself. In principle, it can be replaced by some other, yet unknown fields with similar properties. However, to be more concrete in estimates below we shall use the specific properties of the axion field with known constraints on its coupling constant.

We do not address the problem of formation of QCD-ball in this letter. Instead we concentrate on the problem of stability of these objects. As we will show, once such a configuration is formed, it will be extremely stable soliton like particle. The source of the stability of the QCD-balls is related to the fact that its mass $`M_B`$ becomes smaller than the mass of a collection of free separated nucleons with the same baryon charge. The region of the absolute stability of the QCD-balls is determined by inequality $`m_N>M_BM_{B1}`$ which is satisfied in some region of $`B`$, i.e. $`B_{min}<B<B_{max}`$. The lower limit $`B_{min}`$ in this region determined by inequality $`m_N>M_BM_{B1}`$ when the system becomes unstable with respect to decay to the nucleons. The upper limit $`B_{max}`$ is determined by the region of applicability of our approach when the baryon density in the bulk becomes close to the nuclear density, and therefore, our calculation scheme (based on description in terms of quarks ) becomes unjustified at this point. Different approaches (based on consideration of hadronic rather than quark degrees of freedom) have to be used in this regime. It could happen that some metastable (or even stable) states may exist in this low-density regime. However, the corresponding analysis is beyond the scope of the present work and it shall not be considered here.

Therefore, if sufficiently large number of quarks (determined mainly by the axion properties) is trapped inside the axion bubble during its shrinking, it may result in formation of an absolutely stable QCD-ball with the ground state being a diquark condensate. Such QCD-balls, therefore, may serve as the cold dark matter candidate which amounts about 30% of the total matter/energy of the Universe, $`\mathrm{\Omega }_{DM}0.3`$.

Strictly speaking, the QCD-balls being the baryonic configurations, would behave like nonbaryonic dark matter. In particular, QCD-balls, in spite of their QCD origin, would not contribute to $`\mathrm{\Omega }_Bh^20.02`$ in nucleosynthesis calculations because the QCD-balls would complete the formation by the time when temperature reaches the relevant for nucleosynthesis region $`T1MeV`$. Once QCD-balls are formed, their baryon charge is accumulated in form of the diquark condensate, rather than in form of free baryons, and in such a form the baryon charge is not available for nucleosynthesis. Therefore, the observed relation $`\mathrm{\Omega }_B\mathrm{\Omega }_{DM}`$ within an order of magnitude finds its natural explanation in this scenario: both contributions to $`\mathrm{\Omega }`$ originated from the same physics at the same instant during the QCD phase transition. As is known, this fact is extremely difficult to explain in models that invoke a dark matter candidate not related to baryons.

Before we continue the description of our proposal we would like to make few comments on what have happened on the theoretical side during the last few years, which are crucial elements in our present discussions, and which were not available to researchers earlier.

First of all, there existence of the axion domain walls, related to the symmetry under discrete rotations of the so-called $`\theta `$ angle $`\theta \theta +2\pi n`$ has been known for a long time since . However, the structure of the domain wall considered in  had only one typical scale, $`m_a^11fermi`$. Therefore, the quarks, even if they were trapped inside the bubble at the very first moment, could easily penetrate through such domain wall configuration during the bubble evolution. In this case the axion domain wall ( without support of the fermi pressure from the bulk) would completely collapse. What was realized only quite recently, is the fact that the axion domain walls have actually sandwich substructure on the QCD scale $`\mathrm{\Lambda }_{QCD}^11fermi`$. Therefore, the fermions which are trapped inside the bubble at the very first instant, can not easily penetrate through the domain wall due to this QCD scale substructure, and will likely stay in the bulk, inside the bubble. In this case, the collapse of the axion domain wall stops due to the fermi pressure in the bulk. The arguments ( regarding there existence of the QCD scale substructure inside the axion domain walls) are based on analysis of QCD in the large $`N_c`$ limit with inclusion of the $`\eta ^{}`$ field111 Uniqueness of the $`\eta ^{}`$ field in this problem is related to the special structure of interaction of the axion field $`\theta (x)`$ and the singlet $`\eta ^{}(x)`$ field in low energy description of QCD when only a special combination $`[\theta (x)\eta ^{}(x)]`$ is allowed to enter the low energy QCD Lagrangian.  and independent analysis  of supersymmetric models where a similar $`\theta `$ vacuum structure occurs.

The second important element of our proposal not available earlier, is related to the recent advances, in understanding of CS phase. The fact that the color superconducting phase may exist at high baryon density was discussed a while ago, however it was not a widely accepted phenomenon until recent papers where a relatively large superconducting gap $`\mathrm{\Delta }100MeV`$ with a large critical temperature $`T_c0.6\mathrm{\Delta }`$ were advocated.

To conclude the Introduction we should remark here that the idea that some quark matter, such as strange quark “nuggets” may play a role of the dark matter, was suggested long ago, see also original papers and relatively recent review on the subject. The idea that soliton-like configurations may serve as a dark matter, is also not a new idea. Most noticeable example is being Q-balls. The idea that the dark matter may be just solitons containing large baryon (or even antibaryon) charge is, again, an old idea, see also . The new element of this proposal is the observation that one can accommodate all the nice properties (discussed previously \- ) without invoking any new fields and particles (apart from the axion). Rather, our QCD-balls formed from the ordinary quarks which however are not in the “normal” hadronic phase, but rather in color superconducting phase when squeezed quarks organize a single coherent state described by the diquark Bose–condensate, similar to the Cooper pair condensate in BCS theory in conventional superconductors.

In many respects ( in terms of phenomenology) the QCD balls are similar to strangelets- with few important differences, see below:
$`1.`$ In our proposal the first order QCD phase transition is not required for the formation of the QCD-balls. Axion domain walls of a large size (in comparison with a typical QCD scale) are able to form the large bubbles. These bubbles, filled by $`u,d,s`$ quarks, play the same role as the bubbles formed during the first order phase transition as discussed in.
$`2.`$ The Stability of strange quark matter at zero external pressure, as described in-, is highly model dependent result. In particular, the stability of strangelets is very sensitive to the magnitude of the bag constant within MIT bag model calculations. The idea which is advocated in the present work has a new element, the external pressure due to the axion domain walls. With this new element the stability of the system is very likely to occur in very wide region of the parametric space even in the models which would not support strangelets in the absence of the external pressure.
$`3.`$ The bulk of the QCD ball is in the superconducting phase. This property obviously influences the phenomenology of how the QCD balls interact with a normal matter. In particular, if the energy of a hadron which hits the QCD ball is smaller than the superconducting gap $`\mathrm{\Delta }`$, the hadron can not penetrate into the bulk and excite the internal degrees of freedom of the system, but rather it will be reflected (the so-called, Andreev reflection). Similar property is also true for the strangelets. The elastic cross -section of hadrons on the QCD balls is large, of the order of the geometrical size of QCD balls; the inelastic cross -section (when internal degrees of freedom are excited) is almost identically zero for small energies as mentioned above. Electromagnetic interactions of photons with QCD balls contain the standard fine structure constant $`\alpha `$ with an addition suppression due to the neutrality of the CFL dense quark matter. Such features for the interactions imply that if the QCD ball with small velocity $`v/c10^3`$ enters the Earth, it will not decay by exploding. Rather it will go through the Earth and exit on the opposite side of the Earth leaving behind the shock waves. It is tempting to interpret the recent seismic event with epilinear source as the process which involves the dark matter particle, similar to the QCD ball.
$`4.`$ There is a maximum size of the QCD-ball above which such an object can not be formed and can not be absolutely stable. This is due to the fact that for very large system the axion domain wall pressure becomes a negligible factor which can not stabilize the system.
$`5.`$The property on a maximum size mentioned above has a profound phenomenological consequence. Indeed, if one assumes that stable state as described in exists, the strangelets can collide with an ordinary neutron star which results in formation of a quark star. In such a case all neutron starts would be transformed into quark stars. In our proposal, when the maximal size of the QCD ball is determined by the external axion domain wall pressure, this transition (from neutron stars to quark stars) does not happen as a routine effect.
$`6.`$ If the size of the QCD ball slightly exceeds the maximum critical size, it becomes metastable, rather than stable configuration. Such QCD balls could also be interesting particles for the dark matter phenomenology, see footnote on page 14 and discussions in the Conclusion.

## 2 QCD-balls

Crucial for our scenario is the existence of a squeezer, axion domain wall which will be formed during the QCD phase transition. As is known, there are many types of the axion domain walls, depending on a model. We assume that the standard problem of the domain wall dominance is resolved in some way as discussed previously in the literature, see e.g.,, and we do not address this problem in the present paper222It is widely accepted that the domain walls in the so-called, N=1 axion model will be eaten up by the axion strings at a very high rate. That is true for the axion walls bounded by strings. However, if a domain wall is formed as a closed surface, the probability for such a wall to decay is extremely small. Therefore, such domain walls in $`N=1`$ model can play the same role in our scenario as stable domain walls in $`N1`$ models. Besides that, $`N=1`$ model has a nice property that the domain wall dominance problem is automatically resolved.. We also assume that the probability of formation of a closed bubble made from the axion domain wall is non-zero333We do not attempt to develop a quantitative theory of the formation of the QCD-balls in this work; It is sufficient for our following discussions that this probability is finite.. We also assume that quarks which are trapped in the bulk, can not easily escape the interior when the bubble is shrinking. In different words, the axion domain wall is not transparent due to the QCD sandwich structure of the wall as discussed in ,. The collapse is halted due to the Fermi pressure. Therefore, we assume that a large number of quarks remains in the bulk, inside the bubble when the system reaches the equilibrium.

### 2.1 Equilibrium

The equilibrium is reached when the Fermi pressure equals the surface tension and pressure due to the bag constant $`E_B`$. To put this condition on the quantitative level, we represent the total energy $`E`$ of a QCD-ball with the fixed baryon charge $`B`$, in the following way,

$`E=4\pi \sigma R^2+{\displaystyle \frac{g\mu ^4}{6\pi }}R^3+{\displaystyle \frac{4\pi }{3}}E_BR^3`$ (1)
$`B=gV{\displaystyle _0^\mu }{\displaystyle \frac{d^3p}{(2\pi )^3}}={\displaystyle \frac{2g}{9\pi }}\mu ^3R^3,\mu =\left({\displaystyle \frac{9\pi B}{2gR^3}}\right)^{\frac{1}{3}},`$
where we assume the quarks to be massless. We assume that the relativistic fermi gas is non-interacting in the first approximation, see corrections due to the interactions below. In this formula $`\mu `$ is the Fermi momentum of the system to be expressed in terms of the fixed baryon charge $`B`$ trapped in the bulk; $`R`$ is the size of the sysytem; $`g`$ is the degeneracy factor, $`g2N_cN_f=18`$ for massless degrees of freedom; $`E_B`$ is bag constant which describes the difference in vacuum energy between the interior and exterior. The bag constant is a phenomenological way to simulate the confinement. Finally, $`\sigma f_am_\pi f_\pi `$ is the axion domain wall tension with $`f_a(10^{10}10^{12})GeV`$ being constrained by the axion search experiment.

In what follows, it is convenient to introduce dimensionless scaling variable $`x`$, as follows, $`x\sqrt[3]{B}=R\sqrt[4]{E_B}`$ such that energy per quark $`ϵ_{tot}E/B`$ can be expressed in the following simple way in terms of dimensionless parameters $`x`$ and $`\sigma _0`$,

$`ϵ_{tot}(x){\displaystyle \frac{E}{B}}=E_B^{1/4}\left(4\pi \sigma _0x^2+{\displaystyle \frac{3}{4x}}\sqrt[3]{{\displaystyle \frac{9\pi }{2g}}}+{\displaystyle \frac{4\pi }{3}}x^3\right)`$ (2)
$`xR{\displaystyle \frac{E_B^{1/4}}{B^{1/3}}},\sigma _0{\displaystyle \frac{\sigma }{B^{1/3}E_B^{3/4}}}.`$
The minimization of this expression $`ϵ_{tot}(x)/x|_{x=x_0}=0`$ determines the stability radius $`x_0`$ which fixes the energy of the system at the equilibrium, $`ϵ_{tot}(x_0)`$. In particular, if one neglects $`\sigma _0`$ in eq. (2) originated from the axion domain wall tension, one reproduces the well known results, $`x_00.48,ϵ(x_0)1.9E_B^{1/4}`$. Such a relation means that if $`E_B`$ is relatively small such that the energy per quark is less than $`m_N/3`$, the configuration becomes an absolutely stable state of matter-.

In eqs.(1, 2) we have neglected many important contributions which can drastically change the results. We shall review the role of these contributions below. The main goal of this subsection is the incorporation of these contributions into eqs.(1, 2). First of all, in eq. (1) we neglected the quark-quark interaction on the Fermi surface, which brings the system into superconducting phase for relatively large baryon density . The corresponding contribution $`\mathrm{\Delta }E_{int}`$ to the total energy (1) is negative and at asymptoticaly large $`\mu `$ is equal to,

$`\mathrm{\Delta }E_{int}={\displaystyle \frac{3\mathrm{\Delta }^2\mu ^2}{\pi ^2}}({\displaystyle \frac{4\pi }{3}}R^3)`$ (3)
The negative sign of $`\mathrm{\Delta }E_{int}`$ is quite obvious: the formation of the diquark condensate due to the quark-quark interaction lowers the energy of the system. For appropriate treatment of this term one should express $`\mu `$ as a function of $`B,R`$ according to the relation (1) and substitute this into eq. (2). In principle, one should also take into account that the superconducting gap $`\mathrm{\Delta }(\mu )`$ also strongly varies with $`\mu `$ (and therefore, with $`R`$) in the relevant region of $`\mu `$. However, in what follows we shall ignore this dependence for numerical estimates and shall treat $`\mathrm{\Delta }100MeV`$ as constant. Our last remark regarding eq. (3). This formula was derived for very large $`\mu `$. Nevertheless for illustrative purposes we shall use the expression for $`\mathrm{\Delta }E_{int}`$ for small $`\mu `$ as well. We shall see that in the relevant region of densities the contribution $`\mathrm{\Delta }E_{int}`$ does not exceed $`15\%`$. This somewhat justifies the use of expression (3) for our numerical estimates which follow. With all these reservations in mind, we account the additional contribution to energy per quark, describing the quark-quark interaction on the Fermi surface by adding $`\mathrm{\Delta }ϵ_{tot}^{int}`$ into eq. (2) in the following way

$`\mathrm{\Delta }ϵ_{tot}^{int}=E_B^{1/4}\left(\sqrt[3]{{\displaystyle \frac{4}{\pi }}}{\displaystyle \frac{\mathrm{\Delta }^2}{\sqrt{E_B}}}x\right),`$ (4)
where we expressed everything in terms of dimensionless parameter $`\frac{\mathrm{\Delta }^2}{\sqrt{E_B}}`$ and dimensionless variable $`x`$.

Now we want to consider the modification of eq. (1) which is related to the actual variation of the bag “constant” $`E_B`$ with $`\mu `$. To explain the physical meaning of this effect, we remind the reader that the bag “constant” $`E_B`$ describes the differences of vacuum energies in the interior and exterior regions. It is a phenomenological way to simulate the confinement. The bag “constant” contribution goes with the positive sign to $`E`$, see eq.(1). The physical reason for this sign is obvious: the vacuum energy outside the bubble is lower than inside, thus the positive contribution to $`E`$, in contrast with the interaction term, $`\frac{3\mathrm{\Delta }^2\mu ^2}{\pi ^2}`$ discussed above.

Our main point is as follows: the contribution related to $`E_B`$ can be expressed formally in terms of the difference between the vacuum condensates calculated at zero (exterior) and non-zero (interior) baryon densities. The most important contribution to $`E_B`$ is due to the gluon condensate, such that $`E_B(\mu )\frac{b\alpha _s}{32\pi }G_{\mu \nu }^2_{\mu =0}\frac{b\alpha _s}{32\pi }G_{\mu \nu }^2_{\mu 0}`$ with $`b=\frac{11}{3}N_c\frac{2}{3}N_f`$ where we used the well-known expression for the conformal anomaly in QCD in the chiral limit. We do not know $`E_B(\mu )`$ as a function of $`\mu `$ for the relevant region of the baryon density. However we do know the behavior of this quantity for relatively small densities corresponding to the nuclear matter densities,

$`{\displaystyle \frac{\frac{\alpha _s}{\pi }G_{\mu \nu }^2_{\mu 0}}{\frac{\alpha _s}{\pi }G_{\mu \nu }^2_{\mu =0}}}1{\displaystyle \frac{(0.65GeV)\rho _N}{\frac{\alpha _s}{\pi }G_{\mu \nu }^2_{\mu =0}}}1{\displaystyle \frac{\rho _N}{(264MeV)^3}}`$ (5)
where $`\rho _N`$ is baryon density, and the magnitude for the gluon condensate is known to be, $`\frac{\alpha _s}{\pi }G_{\mu \nu }^2_{\mu =0}1.210^2GeV^4`$. As expected the gluon condensate (and therefore, the absolute value of the vacuum energy) decreases when the baryon density increase. Similar formulae are known for the chiral quark condensate where for the small densities one can derive the following relation $`\frac{\overline{q}q_{\mu 0}}{\overline{q}q_{\mu =0}}=1\frac{\sigma _N\rho _N}{m_\pi ^2f_\pi ^2}`$ with sigma term measured to be $`\sigma _N45MeV`$ see  for the details. One should emphasize here that the formula (5) describing the variation of the gluon vacuum condensate at small baryon densities $`\rho _N`$, is a direct consequence of the QCD low energy theorems. It is a firm result of QCD, not based on any model dependent considerations, and should be accepted as it is.

More specific information on the bag “constant” $`E_B`$ contribution as function of $`\mu `$ in the entire region of of $`\mu `$ can be calculated in some non-physical models such as QCD with two colors, $`N_c=2`$. Such a knowledge can not be literally used for our numerical estimates which follow, however it can be used for modeling the functional dependence of the vacuum energy.

Therefore, we want to model two properties discussed above in order to incorporate them into the corresponding eq. (2). First, the bag constant contribution must vanish when the baryon density in the bulk vanishes. This corresponds to the case when vacuum energy inside and outside of the bubble is the same, and therefore, it should be no additional vacuum energy contribution to the equation for the equilibrium. Secondly, the bag constant contribution should vary with density as we discussed above.

Our first parametrization is motivated by analysis of the vacuum condensates in QCD-like theories at finite baryon density as a function of $`\mu `$. If we assume a similar behavior in real QCD than we should replace the bag constant $`E_B`$ by the expression $`E_BE_B(1\frac{\mu _c^2}{\mu ^2})`$ for $`\mu \mu _c`$ and $`E_B0`$ for $`\mu \mu _c`$, where $`\mu _c`$ would correspond to a magnitude of the critical chemical potential at which the baryon density vanishes. In QCD, one expects that this is to happen at $`\mu _c330MeV`$.

As before, one should express the corresponding contribution to $`ϵ_{tot}`$ in terms of fixed baryon charge $`B`$ and radius $`R`$, such that the bag“ constant” contribution actually becomes a complicated function of $`B,R`$. In terms of dimensional parameter $`x`$ the corresponding contribution to (2) is accounted for by the following replacement,

$`E_B^{1/4}{\displaystyle \frac{4\pi }{3}}x^3E_B^{1/4}{\displaystyle \frac{4\pi }{3}}x^3\left(1({\displaystyle \frac{4}{\pi }})^{2/3}{\displaystyle \frac{\mu _c^2}{\sqrt{E_B}}}x^2\right)`$ (6)
Let us emphasize: we are not attempting to solve a difficult problem of evaluation of nonperturbative vacuum energy as a function of $`\mu `$ in QCD. Rather, we want to make some simple estimates to account for this effect in order to analyze the stability of QCD balls later in the text.

We want to be confident that the results on stability of QCD balls (to be discussed later) are not sensitive to the specific parameterization (6) motivated by the study of QCD with two colors. Therefore, we would like to have a different, independent parameterization of the same effect to be used in our stability analysis. We make use of eq.(5) which is valid for small densities $`\rho _N`$. This formula gives us an idea about typical variation of vacuum condensates when the baryon density changes. We assume that the vacuum energy difference in QCD (the bag “constant” contribution in eq. (2)) can be expressed in terms of different vacuum condensates with the typical scale for the variation given by eq.(5).

We want to implement the QCD property (5) into the MIT bag model. If the phenomenological numerical magnitude for the bag constant $`E_B`$ were closed to the numerical value for the vacuum energy $`\frac{b\alpha _s}{32\pi }G_{\mu \nu }^2(340MeV)^4`$ we could literally use eq (5), such that the bag constant contribution can be parameterized as follows, $`E_B(\rho _N)E_B\frac{\rho _N}{(264MeV)^3}`$. Unfortunately, these two are very different numerically, and we will introduce the corresponding correction factor $`r\sqrt[4]{\frac{b\alpha _s}{32\pi }G_{\mu \nu }^2/E_B}(340MeV)/(150MeV)2.25`$ in our implementation of QCD property (5) into the MIT bag model, see below.

Still, formula $`E_B(\rho _N)\rho _N`$ can not be used literally for our purposes because we need an expression for the bag “constant” contribution which goes to constant $`E_B`$ at large densities, $`E_B(\rho _N)E_B`$. A simple model which satisfies this requirement is to make the following replacement,

$`E_B(\rho _N)E_B{\displaystyle \frac{r^3\rho _N}{(264MeV)^3}}{\displaystyle \frac{E_B}{\left(1+\frac{(264MeV)^3}{r^3\rho _N}\right)}},`$ (7)
where we introduced the correction factor $`r`$ to match the scales. As before, one should express the bag “constant” contribution proportional to (7) in terms of a fixed baryon charge $`B`$ and radius $`R`$. We shall analyse the corresponding equation (2) with improvements (7) in the next subsection. To anticipate the events, one should mention that our two models (6, 7) describing the effect of the bag “constant” variation with baryon density lead to the similar results, see below.

The next approximation we have made in eqs.(1, 2) is related to the assumption of a thin-wall approximation for the domain wall. This may not be well justified assumption because the typical width of the domain wall and the size of QCD ball could be the same order of magnitude, such that thin-wall approximation is failed. However, we neglect these complications at this initial stage of study. Nevertheless, we do not expect that this effect can drastically change our qualitative results which follow.

We also neglected in eqs.(1, 2) all complications related to the finite magnitude of the quark masses, first of all $`m_s`$, which result in additional $`K`$ condensation along with diquark condensation in CFL phase. Finally, the expression for the energy $`E`$ with corrections (3,6), changes the simple relation (1) between baryon charge $`B`$ and chemical potential $`\mu `$ according to the standard thermodynamical relations, $`B=\frac{F}{\mu }`$, where $`F=E\mu B`$ is the free energy. However, we checked that these changes are relatively small( do not exceed $`5\%`$ in the relevant region of $`\mu `$). Therefore, in what follows, in order to avoid the technical complications in the qualitative analysis, we use a simple algebraic expression (1) which is formally valid only for noninteracting quarks, $`B\mu ^3`$, but numerically remains a good approximation in a large region of $`\mu `$. This allows us to use the dimensional variable $`x`$ which we introduced before for the non-interacting case. Let us repeat again: we do not attempt to solve the problem quantitatively with all uncertainties in parameters discussed above; rather, we want to give some qualitative arguments demonstrating that stability region might occur in the wide region of $`B`$ with realistic choice of parameters specified below.

With all these reservations regarding eqs.(1, 2) in mind we express the energy of a QCD-ball per baryon charge $`B`$ in units of $`\sqrt[4]{E_B}`$ , as follows

$`y(x)_{tot}E_B^{1/4}ϵ_{tot}(x)={\displaystyle \frac{4\pi }{3}}x^3\left(1({\displaystyle \frac{4}{\pi }})^{2/3}{\displaystyle \frac{\mu _c^2}{\sqrt{E_B}}}x^2\right)`$ (8)
$`+\left(4\pi \sigma _0x^2+{\displaystyle \frac{3}{4x}}\sqrt[3]{{\displaystyle \frac{\pi }{4}}}\sqrt[3]{{\displaystyle \frac{4}{\pi }}}{\displaystyle \frac{\mathrm{\Delta }^2}{\sqrt{E_B}}}x\right).`$
In this formula, in comparison with eq.(2), we took into account the effect describing the quark-quark interaction on the Fermi surface given by eq. (4) and the effect of the variation of the vacuum energy with baryon density, given by eq. (6).

The equilibrium condition $`ϵ_{tot}(x=x_0)/x=0`$ determines the radius $`x_0`$ of the QCD ball with baryon charge $`B`$. We shall analyze this condition in the next subsection; now we want to constraint $`x_0\overline{x}`$ to be considered. The constraint follows from the condition that the baryon density should be relatively large. In this case our treatment of the problem by using the quark degrees of freedom, eq.(8), rather than hadronic degrees of freedom, is justified. The baryon number density $`\rho _N`$ for the QCD ball configuration is given by444Our normalization for the baryon charge corresponds to $`B=1`$ for the quark, thus factor $`B/3`$ in eq. (9).,

$`\rho _N{\displaystyle \frac{B}{3V}}={\displaystyle \frac{E_B^{3/4}}{4\pi x^3}}n_0,n_0(108MeV)^3,`$ (9)
which gives upper limit $`\overline{x}`$ above which our approach is not justified. Numerically, with our choice of parameters, see below, $`\overline{x}0.6`$, and therefore, any solution $`x_0`$ of the equilibrium condition $`ϵ_{tot}(x=x_0)/x=0`$ must satisfy to the constraint $`x_0\overline{x}0.6`$.

### 2.2 Stability of QCD balls

As expected, the equation describing the equilibrium $`ϵ_{tot}(x=x_0)/x=0`$ has a nontrivial solution (minimum) in a large region of parametrical space deterimed by parameters $`E_B,\sigma ,\mathrm{\Delta },\mu _c,B`$. It is not our goal to have a complete analysis of this allowed region of solutions. Rather, we shall make a specific choice for all parameters except for the baryon number $`B`$ and analize the stability condition as a function of $`B`$. We shall also comment on results with $`\sigma =0`$ corresponding to pure QCD configuration without any involvement of the axion field (case considered previously in MIT bag model, \- ). The first step is to calculate the point $`x=x_0`$ which is determined by equation $`ϵ_{tot}(x=x_0)/x=0`$. The next step is to analyze the stability of the obtained configuration as a function of external parameters. Condition when the QCD-ball becomes an absolutely stable object can be derived from the following arguments. Total energy per quark $`ϵ_{tot}(x_0)`$ in eqs. (2, 8) is a combination of two factors: the first one, $`ϵ_{QCD}(x_0)`$, is due to the strong interactions; the second factor, $`ϵ_{axion}(x_0)`$ is mainly due to the axion domain wall tension555the QCD contribution to $`\sigma `$ due to the $`\eta ^{}`$ and pions is suppressed by a factor $`f_\pi ^2/f_a^21`$., i.e. $`ϵ_{tot}(x_0)=ϵ_{QCD}(x_0)+ϵ_{axion}(x_0)`$, with $`ϵ_{axion}(x_0)E_B^{1/4}\left(4\pi \sigma _0x_0^2\right)`$ and $`ϵ_{QCD}(x_0)`$ is determined by rest of terms in eq. (8). The absolute stability of the system implies that a nucleon can not leave a system because the energy of the configuration with baryon charge $`B`$ is smaller than the energy of configuration of charge $`B3`$ plus energy of a nucleon with baryon charge $`B=3`$ and energy of the axion emission. Such a situation is analogous to the three dimensional quantum mechanical problem with an effective potential being a step-function and the energy of the bound state is lower than the potential energy at the large distances. In this case a particle obviously can not leave the system.

We should emphasize here that the quarks can not leave the system due to the energetic conditions which take place after the QCD ball is formed. In different words, the stability occurs due to the differences in properties inside/outside of the QCD ball, and not due to the features of the original axion domain wall. The axion domain wall already had played its role during the formation period when a large number of quarks could not escape the system and were trapped in the bulk during the collapse of the wall. A similar situation when a configuration may become a stable one due to a difference in conditions (inside/outside the bulk) was discussed long ago as an example of a non-topological soliton in quantum field theory. We further comment on the similarities with non-topological solitons later in the text.

It is quite obvious that the axion domain wall with a typical correlation length $`m_a^1\mathrm{\Lambda }_{QCD}^1`$ can not produce nucleons by itself when it shrinks due to the nucleon emission. Instead, typically, the axion domain wall reduces its size by emitting the axions while the nucleon leaves the system. In this case the term $`ϵ_{axion}(x_0)E_B^{1/4}\left(4\pi \sigma _0x_0^2\right)`$ is responsible for the emission of axions rather than production of nucleons. As a result of this, this term should be ignored for the analysis of the stability. However, with exceedingly small probability the emitted axion, in principle, can be absorbed by the nucleon which leaves the system. In this case the energy, in principle, can be transformed from the axion domain wall to the produced nucleon and the term $`ϵ_{axion}(x_0)E_B^{1/4}\left(4\pi \sigma _0x_0^2\right)`$ should be accounted in the energy budget for analysis of the decay. We estimate in appendix that the probability for the corresponding absorption of the axion by the leaving nucleon is negligible. Therefore in what follows we neglect this process.

The relevant term which describes the emission of nucleons is the one related to the QCD physics i.e. $`ϵ_{QCD}(x_0)`$. Therefore, the condition when configuration becomes a sufficiently stable (with the life time exceeding the life time of the Universe, see Appendix for details) is determined from the following inequality

$`ϵ_{QCD}(x_0)<{\displaystyle \frac{m_N}{3}},{\displaystyle \frac{ϵ_{tot}(x)}{x}}|_{x=x_0}=0,x_0<\overline{x},`$ (10)
where the last condition follows from (9).

To analyse eq.(10) we shall accept the following magnitudes for the dimensional parameters:

$`\mathrm{\Delta }100MeV;\sigma 1.810^8GeV^3;`$ (11)
$`\mu _c330MeV;E_B(150MeV)^4.`$
Having these external parameters fixed, we left with the only one unknown number, the baryon charge $`B`$, which eneters $`\sigma _0`$ in our dimensionless parametrization (2,8). We shall treat $`\sigma _0`$ as a free parameter and our goal is to find the region of $`\sigma _0`$ when conditions (10) are satisfied. As we discussed above, we shall use two different models to account the effect of the variation of the bag constant contribution with density, see eqs. (6, 7).

Having defined our stability condition (10), external parameters (11) and two simple models accounting the effect of the variation of the bag constant, eqs. (6, 7), we reduce our problem to analysis of dimensionless functions, $`y_{QCD}^{(1)}(x)`$ and $`y_{QCD}^{(2)}(x)`$ defined as follows, see eqs. (6, 7, 8),

$`y_{tot}(x)y_{QCD}^{(1,2)}(x)+y_{axion}(x);y_{axion}(x)4\pi \sigma _0x^2`$ (12)
$`y_{QCD}^{(1)}{\displaystyle \frac{0.69}{x}}+4.2x^3{\displaystyle \frac{1}{1+6x^3}}0.48x,`$ (13)
$$y_{QCD}^{(2)}\frac{0.69}{x}+4.2x^3(15.68x^2)0.48x,$$
(14)
where three consequent terms describe: the fermi pressure, the bag constant contribution accounting the variation of the vacuum energy with the baryon density (7,6), and, finally, the quark-quark interaction on the fermi surface (4) correspondingly. Stability condition (10) in dimensionless variables becomes

$`y_{QCD}^{(1,2)}(x_0)<{\displaystyle \frac{m_N}{3\sqrt[4]{E_B}}}2.1,{\displaystyle \frac{y_{tot}(x)}{x}}|_{x=x_0}=0.`$ (15)
Before we discuss some specific numerical results which follow from analysis of eqs. (13 -15), we would like to list some general model-independent properties of the solutions. We believe that the properties listed below are quite common features of the QCD balls, which likely to remain untouched even in a more general treatment of the problem when many additional effects are included (some of these effects were mentioned above).

a). As we already mentioned, in the absence of the axion field, $`\sigma 0`$, the problem was extensively discussed earlier using MIT bag model,-. Our original remark here is: when a variation of the vacuum energy with density is taken into account, a stable solution disappears provided that a typical QCD scale for the vacuum variation (6,7) is used. The physical reason for that behavior is quite obvious: a density -dependent vacuum energy is not a sufficiently strong squeezer to equilibrate the fermi pressure. A typical scale for the variation should be reduced (in comparison with what we assumed in eqs.(6,7) ) by an order of magnitude, in order for the solution to reappear. Specifically, we checked that the equilibrium is possible for $`\sigma 0`$ if coefficient $`5.68`$ in (14) describing the vacuum energy variation is replaced by $`0.5`$. We demonstrate this effect in Fig.1 where we display the total energy of QCD ball per baryon charge (12) with $`\sigma 0`$ for three different values of parameter $`(\frac{4}{\pi })^{2/3}\frac{\mu _c^2}{\sqrt{E_B}}`$ describing the effect of variation of vacuum energy with the baryon density. For a typical choice of physical values (11) the relevant parameter is $`(\frac{4}{\pi })^{2/3}\frac{\mu _c^2}{\sqrt{E_B}}=5.68`$. In this case the minimum does not exist which implies that the stability can not be achieved as announced above. The minimum starts to reappear only when a typical scale for the variation of density is considerably reduced, see Fig.1 with the curve corresponding $`(\frac{4}{\pi })^{2/3}\frac{\mu _c^2}{\sqrt{E_B}}=0.3`$. Such a small value for the critical value $`\mu _c^2`$ does not look appealing from the physics point of view. Therefore, we incline to accept that there is no solution for such a configuration (strange quark nuggets,-) in QCD if no external pressure (such as gravity or axion domain wall) is applied. It is certainly not a very new result: special study on stranglets reveals a strong model dependence of the stability of strange quark matter. In particular, the Nambu Jona- Lasinio model does not support any kind of strangelets.

b). In general, one expects there existence of a minimal and maximal sizes for the QCD balls in the region of stability. The minimal charge $`B_{min}`$ corresponds to the maximum $`\sigma _0^{max}B_{min}^{1/3}`$. At this point the the stability requirement (10) is marginally satisfied. When $`B<B_{min}`$, $`\sigma _0`$ becomes too large such that nucleons can leave the system. On the other hand, the maximum possible charge, $`B_{max}`$, corresponds to the minimum value of $`\sigma _0^{min}B_{max}^{1/3}`$. For larger $`B`$, the baryon density (9) becomes too low to justify our approach (based on the quark degrees of freedom). At lower baryon densities some metastable states may form; they could decay to some heavy elements which might be of interests for astrophysics. However the corresponding study would require an analysis of the system in terms of nuclear degrees of freedom, which is beyond the scope of the present work. When $`\sigma _0`$ becomes even smaller, the problem is essentially equivalent to $`\sigma =0`$ studied earlier where stable solutions are not expected to occur.

Numerically, we analyzed two models (13, 14) which lead to the similar results. In particular, for model (13) the maximum possible tension, $`4\pi \sigma _0^{max}10`$ corresponds to the minimum baryon charge $`B_{min}`$. For such $`\sigma _0`$ the equilibrium is reached at $`x_00.32`$ wnen the energy per quark $`y_{QCD}^{(1)}(x_0)2.1`$ hits the upper energy bound of the stability region (15). When $`4\pi \sigma _0^{max}>10`$, the energy per quark becomes too high such that nucleon can escape and the system would decay. In physical units this solution corresponds to $`B_{min}10^{32}`$ and stabilization radius $`R_0=x_0\sqrt[3]{B}/\sqrt[4]{E_B}10^{11}GeV^1`$. Energy per quark for this configuration $`ϵ_{QCD}^{(1)}=y_{QCD}^{(1)}(x_0)\sqrt[4]{E_B}2.1\sqrt[4]{E_B}320MeV`$ is smaller than constituent quark mass, as it should be666We remind that we discuss the QCD part of energy only; the total energy of the configuration which includes the axion part is larger..

For the same model, the minimum possible tension when our approach is justified, $`4\pi \sigma _0^{min}2`$ corresponds to the maximum possible baryon charge $`B_{max}`$. According to the scaling $`B\sigma _0^3`$, the maximum baryon charge $`B_{max}=(\frac{\sigma _0^{max}}{\sigma _0^{min}})^3B_{min}10^{34}`$ is two orders of magnitude larger than $`B_{min}`$. In this case the equilibrium is reached at $`x_00.52`$ when the baryon density (9) is already relatively low, and close to the boundary when the quark based lore can not be trusted.

Our second model (14) gives quantitatively similar results, and it is not worthwhile to discuss numerical details here. The most important features of the solution for this model remain the same: there is a region between $`B_{max}`$ and $`B_{min}`$ when solutions are stable; at $`\sigma =0`$ solution does not exist at all provided that a typical QCD scale for the vacuum variation (6,7) is used.

However, one should take all these numerical estimates very cautiously because of a number approximations we have made in eqs. (1, 2). Nevertheless, in what follows, mainly for the illustrative purposes, we shall stick with these numerical estimates.

The quark number density $`n`$ in the region $`B_{min}<B<B_{max}`$ when our approach is justified is estimated as

$`n{\displaystyle \frac{B}{V}}={\displaystyle \frac{3E_B^{3/4}}{4\pi x^3}}(1.56.5)3n_0,`$ (16)
where we used the expression (1) for the baryon density. As we already mentioned the expression (1) is formally valid only for noninteracting quarks, but numerically remains a good approximation in a large region of $`\mu `$. In eq. (16) $`3n_03(108MeV)^3,`$ is the nuclear saturation density normalized with our convention ( $`B=1`$ for quarks), thus factor $`3`$ in front of the numerica value $`0.16(fm)^3(108MeV)^3`$. It is quite remarkable that the numerical value for $`n`$ is in the region where color superconductivity phase is likely to realize, and therefore, our treatment of the squeezed fermi system as the quark dense matter (rather than ordinary nuclear matter) is justified a posteriori.

Few remarks are in order regarding eq.(16). First of all, the estimates presented above demonstrate that we are in the region of the phase diagram where CFL phase is likeley to realize. Therefore, our original assumption is justified. Secondly, for large $`BB_{max}`$ our treatment of the system is not valid anymore, and a different type of QCD balls with an ordinay nuclear matter (instead of diquark condensate) in the bulk may be formed and could be even stable in some regions of parametrical space. Though this region of large $`BB_{max}`$ could be an interesting region from the phenomenological point of view, it shall not be discussed here777The corresponding proper treatment would require the knowledge of the dynamics of the interacting nuclear matter, which is not the subject of the present work. In principle such nuclear matter could be also stable.. However, even in this case when the QCD balls made of nuclear matter, rather than quark dense matter, we still expect that there should exist a maximum size above which the stability is not possible. This follows from our analysis that stability can not be achieved without the external pressure $`P_\sigma `$ due to the axion domain $`P_\sigma 2\sigma /R`$ which vanishes at very large $`R`$.

Another factor which also constraints the size of the balls is related to the suppression of large size closed axion domain walls during the formation stage. It is clear that the formation of the large size closed domain walls is suppressed according to the Kibble-Zurek mechanism,; however an explicit estimation for this effect is still missing.

As we mentioned in the Introduction, we do not address the problem of formation of QCD balls in this letter, it will be a subject of a different work. However we would like to mention some relevant elements of a possible scenario of how QCD-balls, in principle, can be formed after the QCD phase transition, at a temperature of order 150 MeV which is much higher than the critical temperature for quark pairing estimated to be $`0.6\mathrm{\Delta }`$. The main point is this: the axion domain wall with the QCD-scale substructure as discussed in is very selective with respect to the momentum of the particles; it is almost transparent for light $`\pi `$ mesons with large momentum $`km_\pi `$ such that the transmission coefficient is close to one. Therefore, the highly energetic pions can easily penetrate through the domain wall and leave the system.

At the same time, the transmission coefficient is close to zero for slow-moving particles such as baryons with $`km_\pi `$. Eventually, this “selective” feature of the domain wall may cool down the system considerably. Due to the domain wall pressure it may reach the critical density when it undergoes a phase transition to a color superconducting phase with the ground state being the quark condensate. At this point we assume that the baryon number trapped in the bulk is sufficiently large. If $`BB_{max}`$, the quarks will leave the system by forming nucleons until the upper limit $`B_{max}`$ is achieved. At this point the energy per unit baryon charge $`ϵ_{QCD}=\sqrt[4]{E_B}y_{QCD}(x_0)<\frac{m_N}{3}`$, becomes sufficiently small such that quarks can not leave the system. Some specific calculations are required before any statements regarding a possibility to form the QCD balls can be made. We do not see any fundamental obstacles which would prevent the formation of such objects. Terefore, at this moment we simply assume that this probability does not vanish.

### 2.3 QCD- balls versus Q-balls

In this subsection we would like to mention a striking resemblance of the QCD-balls (which is the subject of this letter) and non-topological solitons, as well as Q-balls which is a special case of a nontopological soliton configuration associated with some conserved global $`Q`$ charge. Both cases, QCD balls and Q-balls demonstrate a similar behavior for a soliton mass as function of $`Q`$. Namely, QCD balls as well as Q-balls may become very stable configurations for relatively large $`Q`$ charge. Therefore, an effective scalar field theory with some specific constraint on the effective potential (when $`Q`$ ball solution exists) is realized for QCD in high density regime by formation of the diquark scalar condensate which plays the role of the effective scalar field. The big difference, of course, that underlying theory for QCD-balls is well known, it is QCD with no free parameters. This is in contrast with the theory of Q-balls when the underlying theory is not known. Formal similarity becomes even more striking if one takes into account that the ground state of the CFL phase in QCD is determined by the diquark condensate with the following time dependence $`e^{i2\mu t}`$,

$`\mathrm{\Psi }_{La}^{i\alpha }\mathrm{\Psi }_{Lb}^{j\beta }^{}\mathrm{\Psi }_{Ra}^{i\alpha }\mathrm{\Psi }_{Rb}^{j\beta }^{}`$ $`(e^{i2\mu t})ϵ^{ij}ϵ^{\alpha \beta c}ϵ_{abc},`$ (17)
with $`\mathrm{\Psi }`$ being the original QCD quark fields, and $`\mu `$ being the chemical potential of the system, see formula (40) from ref. . As is known, such time-dependent phase is the starting point in construction of the Q balls. In the expression (17) we explicitly show the structure for the diquark condensate corresponding to CFL (color-flavor locking) phase with ($`\alpha `$, $`\beta `$, etc.) to be flavor, ($`a`$, $`b`$, etc.) color and ($`i`$, $`j`$, etc.) spinor indices correspondingly. Of course, there are many differences in phenomenology between Q balls and QCD-balls. For example, in CFL phase the baryon symmetry is spontaneously broken, and corresponding Goldstone massless boson carries the baryon charge. However, the evaporation of this massless particle into hadronic phase from the surface of the QCD-ball is not possible, because hadronic phase does not support such excitation. This is in contrast with phenomenology of Q-balls, where the theory is formulated in terms of one and the same scalar $`\varphi `$ field, such that evaporation of $`\varphi `$ particles from the surface of the Q-ball is possible if some conditions are met. In spite of many differences, the analogy with Q-balls is quite useful and can be used for analysis of different experimental bounds on QCD-balls, which is the subject of the next section.

## 3 Experimental bounds on masses and fluxes of QCD-balls

In this section we adopt the results of paper to constraint the free parameter (charge $`B`$) of the QCD-balls. In the paper  the authors re-analyzed the results of various experiments, originally not designed for the Q-ball searches, but nevertheless these experimental results were successfully used in  to bound different properties of the Q-balls. We actually repeat this analysis for a specific type of the QCD-balls when original quarks are in the CFL (color-flavor locking) phase.

As we mentioned earlier, at sufficiently large baryon density, the color superconductivity phenomenon takes place. However, there are many different phases (as a function of parameters like $`m_s`$, number of light flavors, etc.) associated with color superconductivity. In particular, for 3 degenerate flavors of light quarks, the CFL phase with nonzero value for the diquark condensate (17) is realized. Due to the fact that equal numbers of $`u,d,s`$ quarks condensed in the system, the electric charge of the ground state is zero, i.e. no electrons required to neutralize the system. This is quite important feature for the phenomenology of the QCD-balls we about to discuss. Nature is less symmetric, and other CS phases could be realized. In particular, for relatively large $`m_s`$, along with diquark condensate, the $`K`$ condensate may also be formed. In the limit of very large $`m_s`$, QCD becomes effectively a theory with two light quarks. In this case, the Cooper pairs are $`uddu`$ flavor singlets. This phase, the so-called 2SC (2 flavor super-conductor ) phase is a phase with non-zero electric charge. Electrons neutralize the system, however, all properties, such as interaction cross sections, the rate of energy loss of QCD balls in matter, are very different for QCD-balls with quarks in CFL or 2SC phase. In what follows, to avoid many complications, we limit ourself with analysis of QCD balls where quarks are in the most symmetric CFL phase, in which case the QCD-ball has zero electric charge.

We assume, in analogy with, that a typical cross section of a neutral QCD-ball with matter is determined by their geometrical size, $`\pi R_0^2`$. In this case, the only information we need to constraint the QCD-ball parameters, is its size and mass. We also assume that the QCD-balls is the main contributor toward the dark matter in the Galaxy. Their flux $`F`$ then should satisfy

$`F<F_{DM}{\displaystyle \frac{\rho _{DM}v}{4\pi M_B}}7.210^5{\displaystyle \frac{GeV}{M_B}}cm^2sec^1sr^1,`$ (18)
where $`\rho _{DM}`$ is the energy density of the dark matter in the Galaxy, $`\rho _{DM}\frac{0.3GeV}{cm^3}`$, and $`v310^3c`$ is the Virial velocity of the QCD-ball. We identify $`M_B`$ in the expression (18) with the total energy $`E`$ of the QCD ball at rest with given baryon charge $`B`$. The Gyrlyanda experiments at Lake Baikal reported that the flux of neutral soliton-like objects has the bound

$`F<3.910^{16}cm^2sec^1sr^1,`$ (19)
which translates to the following lower limit of the neutral QCD-ball mass $`M_B`$ and baryon charge $`B`$,

$`M_B^{exp}>210^{21}GeV,`$ (20)
$`B^{exp}({\displaystyle \frac{M_B}{\sigma ^{1/3}}})^{9/8}[{\displaystyle \frac{3}{2}}(8\pi c^2)^{\frac{1}{3}}]^{9/8}>1.610^{20}.`$
Similar constraints follow from the analysis of the Baksan experiment and analysis  of the Kamiokande Cherenkov detector, and we do not explicitly quote these results. These experimental bounds are well below the critical line of the stability of the QCD-balls.

## 4 Discussions and Future directions

Complete theory of formation of the QCD-balls is still lacking. Only such a theory would predict whether QCD-balls can be formed in sufficient number to become the dark matter. Such a theory of formation of the QCD balls would answer on questions like this: 1. What is the probability to form a closed axion domain wall with size $`\xi `$ during the QCD phase transition? 2. How many quarks are trapped inside the domain wall at the first instant? 3. How many quarks will leave the system and how many of them will stay inside the system while the bubble is shrinking? 4. What is the dependence of relevant parameters such as: size $`\xi (t)`$, baryon number density $`n(t)`$ and internal temperature $`T(t)`$ as function of time? 5. Do these parameters fall into appropriate region of the QCD phase diagram where the color superconductivity takes place? 6. What is the final density distribution of the QCD-balls as a function of their size $`R`$ after the formation period is complete? 7.Will the QCD balls survive the evaporation and boiling even if they formed? Clearly, we do not have answers on these, and many other important questions at the moment. All these interesting, but difficult questions are obviously beyond the scope of the present work, and shall not be discussed here. However, we want to make a short comment on issue 7 which was an important element in many previous studies.

The question on evaporation of quark nuggets was discussed earlier, see original papers - and recent review . The first study of this question is due to Alcock and Farhi who argued that only very large nuggets with $`B10^{52}`$ could survive the evaporation. This result would essentially eliminate the possibility of any quark nuggets surviving till the present epoch. However, Madsen et al. then point out that few important effects can considerably reduce the original estimation given in ref.. The first important effect is related to the deficiency of $`u`$ and $`d`$ quarks (in contrast with $`s`$ quark) in the surface area. This leads to the suppression of the evaporation rate such that $`B10^{46}`$ can be stable against evaporation. In this calculation the penetrability of the phase boundary was assumed to be near $`100\%`$. This assumption was questioned in  and  where it was demonstrated that nuggets with $`B10^{43}`$ ($`B10^{39}`$ according to ref. ) could survive the evaporation even if the first effect ( described above and which led to $`10^6`$ suppression, see  for details) is neglected. As discussed in - the limit on $`B`$ may be further reduced by reabsorption. All these effects taken together suggest that nuggets with $`B10^{30}`$ are not ruled out and can survive the evaporation.

Our original remark here is: along with the suppression effects discuseed above, we have two additional effects which may further reduce the evaporation rate.

Indeed, the core of the axion domain wall as discussed in, has a QCD sub-structure with a typical scale $`1GeV`$. It is quite obvious that this sub-structure certainly reduces the penetrability of particles from inside to outside, and therefore, it suppresses the evaporation rate. Also, the baryon charge in superconducting phase is in the form of the diquark condensate rather than in form of free quarks discussed in the previous analysis-. This fact may also considerably reduce the evaporation rate because it requires the breaking of the Cooper pair before the evaporation becomes possible. This effect certinly increases the effective binding energy and decreases the evaporation rate. It is difficult to make a precise estimate of these effects at the moment, due to the many compications discussed earlier as well as many additional difficulties mentioned above. However, we believe, it is fair to say that the QCD balls with $`B10^{32}`$ as discussed in the previous section, can safely survive the evaporation, and therefore, the possibility seems worth exploring.

Now, we wish to estimate the absolute value for the dark matter number density $`n_{DM}`$ assuming that the nonbaryonic dark matter is actually the QCD balls. In this case, at the QCD phase transition at $`TT_c`$ soon after the QCD balls are formed, $`n_{DM}`$ can be estimated as follows,

$`n_{DM}510^9{\displaystyle \frac{2\pi ^2}{45}}g_{}T_c^3{\displaystyle \frac{m_N}{M_B}},`$ (21)
where we used the known magnitudes for the baryon to photon ratio, $`n_B/n_\gamma 510^{10}`$, and the dark matter to baryon ratio, $`\mathrm{\Omega }_{DM}/\mathrm{\Omega }_B10`$. Numerically, for the baryon charge $`B10^{32}`$ and effective massless degrees of freedom, $`g_{}10`$ the estimation (21) leads to

$`rT_cn_{DM}^{1/3}T_c3.510^{13},r10cm,`$ (22)
where $`r`$ has the physical meaning of an average distance between QCD-balls after they formed. As expected, average distance $`r`$ is much smaller than the horizon radius $`R_H^{QCD}`$ at the QCD phase transition, $`r10^5R_H^{QCD}`$. It is quite remarkable that $`r`$ is much larger than the size of the QCD-ball, see eq.(2), such that QCD-balls become well separated soon after they formed. Besides that we expect that the QCD ball size should be related, through dynamics, to the correlation length $`\xi m_a^1`$ of the original axion field. We also expect that the spatial extend of a typical closed wall at the instant of formation has the same order of magnitude $`\xi `$ . Initial size of a closed wall $`\xi `$ eventually (after some shrinking as a result of tension, and after some expansion as a result of evolution of the Universe) determines the size of the QCD-balls. However, the dynamics of this transition is quite complicated, and we are not able to derive a relation between initial domain wall size distribution and QCD-ball size distribution at the later stage. Close numerical values for the QCD ball size and $`\xi m_a^1`$ also suggest that these parameters are related somehow. Therefore, it is at least possible, that the decay of the axion domain wall network may result in formation of the QCD-balls with their nice properties discussed in this work.

To conclude: we advocate the idea that the QCD-balls could be a viable cold dark matter candidate which is formed from the ordinary quarks during the QCD phase transition when the axion domain walls form. As we argued the system in the bulk may reach the critical density when it undergoes a phase transition to a color superconducting phase in which case the new state of matter representing the diquark condensate with a large baryon number $`B`$ becomes a stable soliton-like configuration. The scenario is no doubt lead to important consequences for cosmology and astrophysics, which are not explored yet. In particular, some unexplained events, such as Centauro events, or even the Tunguska-like events (when no fragments or chemical traces have ever been recovered), can be related to the very dense QCD balls. The recent detection of two seismic events with epilinear (in contrast with a typical epicentral ) sources may also be related to the very dense QCD balls. Also, the “missing” baryons in Galaxy Clusters may also be related to the QCD balls. Finally, the cuspy halo problem in dwarf galaxies might be related to the unstable cold dark matter , which, again, could be related to the QCD balls discussed in this work. Indeed, as we mentioned in the Introduction, if the QCD ball size exceeds the critical value, it becomes metastable (rather than stable) configuration. The life time of these metastable QCD balls could be very large. Therefore, they could serve as decaying dark matter particles suggested in.

Therefore, the “exotic”, dense color superconducting phase in QCD, might be much more common state of matter in the Universe than the “normal” hadronic phase we know. More than that: one can present some arguments to support the idea that the observed in nature asymmetry between baryons and antibaryons may also be originated from the same physics during the QCD phase transition. In this case the antimatter is hidden inside of the anti- QCD balls in the form of the diquark condensate similar to the QCD ball case. One could naively think that such a scenario is in contradiction with observations on absence of antimatter around us. However, such a conclusion would be very premature one due to the specific interaction features of the matter in hadronic phase with the matter in color superconducting phase. Namely, if the energy of the quark which hits the anti-QCD ball is smaller than the superconducting gap $`\mathrm{\Delta }`$, the quark can not penetrate into the bulk, break the Cooper pair and excite the internal degrees of freedom of the system. Rather it will be reflected . A similar property is commonly known as the “Andreev Reflection” in the literature on conventional superconductivity. Therefore, at low energies, the anti-QCD balls behave as QCD balls with respect to the interaction with environment, and there is no contradiction with known constraints on such kind of anti matter in our Universe.

In this case, without fine tuning of parameters, one can easily understand the relation between $`\mathrm{\Omega }_{DM}\mathrm{\Omega }_B`$ which both originated at the same instant. As is known this ratio is very difficult to understand if these quantities do not have the same origin.

In conclusion, qualitative as our arguments are, they suggest that the dark matter could be originated at the QCD scale.

## Acknowledgments

I am thankful to Robert Brandenberger, Michael Forbes, Sasha Kusenko, Jes Madsen, Shmuel Nussinov, Krishna Rajagopal, Edward Shuryak, Pierre Sikivie, Paul Steinhardt and Frank Wilczek for useful discussions, comments and remarks. This work was supported in art by the National Science and Engineering Research Council of Canada.

## Appendix

The main goal of this Appendix is to argue that the probability of the absorption of axion by emitted quark is extremely small. In this case our criteria of stability of the QCD balls, see eq.(10), which neglects the axion domain wall energy contribution, is justified. Such a treatment of the problem essentially implies that we impose a weaker condition of metastability (rather than a stronger condition of the absolute stability) on the QCD balls. In different words, we assume that the energy can not be transformed from the axion domain wall to the quark which is about to leave the system. In what follows we make some estimates which support this assumption. Indeed, as we shall see in a moment this probability is exceedingly small due to the very small nucleon-axion coupling constant, $`\frac{m_q}{f_a}(10^{13}10^{15})`$.

We start from the estimation of the probability to absorb the axion from the axion domain wall background field $`a(z)`$ by an elementary excitation $`|\psi _{in}>`$ with mass $`\mathrm{\Delta }`$ proportional to the gap. This elementary excitation carries the unit baryon chage in superconducting phase. We assume that the final state is represented by the wave function of constituent quark $`<\psi _{out}|`$ with mass $`m_N/3`$ ( hadronic phase). We take a simpe expression for the axion-quark interaction to be $`\frac{m_q}{f_a}\psi _{out}^{}a(z)\psi _{in}`$, where $`m_q`$ is the current mass quark of order few $`MeV`$, and $`\frac{a(z)}{f_a}1`$ is the axion domain wall background field. The precise expression for the domain wall profile function $`a(z)`$ is known, however, in our estimate we shall use a simple expression $`a(z)f_ae^{m_az}`$ in order to emphasize that the magnitude of the axion field vanishes at infinity and the typical scale where axion field varies is $`m_a^1`$. We also assume that the quark has a trajectory $`z=vt`$ with velocity $`v`$ close to the speed of light. In this case the time-dependent interaction takes the form $`m_qe^{m_at}`$, and the probability for the transition can be estimated from the dimensional arguments as follows,

$`W|m_q{\displaystyle _0^{\mathrm{}}}dte^{i\omega tm_at}+h.c.|^2|{\displaystyle \frac{m_qm_a}{m_a^2+\omega ^2}}|^2,`$ (23)
where $`\omega `$ is the energy difference between $`|\psi _{in}>`$ and $`|\psi _{out}>`$ states. It is impotant to note that, typically $`\omega (100200)MeV`$ is large, and therefore, the probability (23) is very small. We neglected many factors in estimate (23). In particular, we neglected the momentum dependence of $`|\psi _{in}>`$ and $`|\psi _{out}>`$ states; the mismatch between these momenta would bring an additional suppression to (23), and we ignore this effect at the moment. It is easy to understand the source of the suppression in eq.(23): the probability for a considerable excitation $`\omega `$ of the system by a smooth field with a typical correlation scale $`m_a^1`$ is very small.

In order to derive a total number of events of absorption $`W_{tot}`$ one should multiply the expression(23) by an additional factor describing a total number of elementary quark excitations close to the surface of the system such that they can leave the system without re-scattering. This requirement (to be close to the surface of the QCD ball) is important because the distance from the surface should not exceed the mean free path. Otherwise, the quark even if it absorbs the axion, would not be able to leave the system. Assuming the thermodynamical equilibrium at temperature $`T`$ soon after the formation of the QCD balls, we can estimate this factor as follows $`\frac{2\pi ^2}{45}g_{}T^3\mathrm{exp}(\frac{\mathrm{\Delta }}{T})4\pi R^2\xi ,`$ where $`\xi `$ is the mean free path which we estimate to be $`1/T`$. Our final expression for the total probability of absorption of the axion (while the temperature is of order $`T`$) is estimated to be

$`W_{tot}|{\displaystyle \frac{m_qm_a}{m_a^2+\omega ^2}}|^2{\displaystyle \frac{2\pi ^2}{45}}g_{}T^24\pi R^2\mathrm{exp}({\displaystyle \frac{\mathrm{\Delta }}{T}}),`$ (24)
where we neglected many additional suppression factors, such as factor $`1/6`$ describing the probability for the quark to move in the direction pointing off the center of the QCD ball. Numerically, even if we neglect the factor $`\mathrm{exp}(\frac{\mathrm{\Delta }}{T})`$ in eq.(24), the probability is already quite small, $`W_{tot}<10^3`$ for the typical values of $`g_{}10,T0.6\mathrm{\Delta }`$. When temperature becomes considerably smaller than $`T`$, the probability of absorption diminishes due to the small number of excitations, $`\mathrm{exp}(\frac{\mathrm{\Delta }}{T})`$. This late epoch of evolution can be ignored. Also, one should keep in mind that the quark excitations are not supported in the hadronic phase due to the confinement. Therefore, one should have three quarks (or quark and diquark pair from the condensate) to be organized in a color singlet state such that it can propagate in the hadronic phase. It definitely gives an additional suppression which we even did not try to estimate: the suppression factor (24) is already sufficiently strong for our purposes. Therefore, our treatment of the problem when we use a weaker condition of metastability, see eq. (10), rather than a stronger condition of the absolute stability, is justified.

## References