Cosmological Axion and Quark Nugget Dark Matter Model

We study a dark matter (DM) model offering a very natural explanation of two (naively unrelated) problems in cosmology: the observed relation $\Omega_{\rm DM}\sim\Omega_{\rm visible}$ and the observed asymmetry between matter and antimatter in the Universe, known as the"baryogenesis"problem. In this framework, both types of matter (dark and visible) have the same QCD origin, form at the same QCD epoch, and both proportional to one and the same dimensional parameter of the system, $\Lambda_{\rm QCD}$, which explains how these two, naively distinct, problems could be intimately related, and could be solved simultaneously within the same framework. More specifically, the DM in this model is composed by two different ingredients: the (well- studied) DM axions and (less-studied) the quark nuggets made of matter or antimatter. The focus of the present work is the quantitative analysis of the relation between these two distinct components contributing to the dark sector of the theory determined by $\Omega_{\rm DM}\equiv [\Omega_{\rm DM}(\rm nuggets)+ \Omega_{\rm DM}(\rm axion)]$. We argue that the nugget's DM component always traces the visible matter density, i.e. $\Omega_{\rm DM}(\rm nuggets)\sim\Omega_{\rm visible}$ and this feature is not sensitive to the parameters of the system such as the axion mass $m_a$ or the misalignment angle $\theta_0$. It should be contrasted with conventional axion production mechanism due to the misalignment when $\Omega_{\rm DM}(\rm axion)$ is highly sensitive to the axion mass $m_a$ and the initial misalignment angle $\theta_0$. We also discuss the constraints on this model related to the inflationary scale $H_I$, non-observation of the isocurvature perturbations, $r_T<0.12$, and also, varies axions search experiments.


I. INTRODUCTION AND MOTIVATION
The idea that the dark matter may take the form of composite objects of standard model quarks in a novel phase goes back to quark nuggets [1], strangelets [2], nuclearities [3], see also review [4] with large number of references on the original results. In the models [1][2][3][4] the presence of strange quark stabilizes the quark matter at sufficiently high densities allowing strangelets being formed in the early universe to remain stable over cosmological timescales. There were a number of problems with the original idea 1 and we refer to the review paper [4] for the details.
The quark nugget model advocated in [7] is conceptually similar, with the nuggets being composed of a high density colour superconducting (CS) phase. An additional stabilization factor in the quark nugget model is provided by the axion domain walls which are copiously produced during the QCD transition 2 . The only new additional crucial element in the proposal [7] is that the 1 In particular, the first order phase transition is a required feature of the system for the strangelet to be formed during the QCD phase transition. However it is known by now that the QCD transition is a crossover rather than the first order phase transition as the recent lattice results [5] unambiguously show. Furthermore, the strangelets will likely evaporate on the Hubble time-scale even if they had been formed [6]. 2 In this case the first order phase transition is not required for the nuggets to be formed as the axion domain wall plays the role of the squeezer. Furthermore, the argument related to the fast evaporation of the strangelets as mentioned in footnote 1 is not applicable for the quark nugget model [7] because the vacuum ground state energies inside (CS phase) and outside (hadronic phase) the nuggets are drastically different. Therefore these two systems can coexist only in the presence of the additional exter-nuggets could be made of matter as well as antimatter in this framework, see original papers [7,8]. This novel crucial element of the model [7,8] completely changes entire framework because the dark matter density Ω dark and the baryonic matter density Ω visible now become intimately related to each other and proportional to each other Ω dark ∼ Ω visible (1) as will be explained in next section II. In other words, the nature of dark matter and the problem of the asymmetry between matter and antimatter in the Universe, normally formulated as the so-called baryogenesis problem, become two sides of the same coin in this framework. To reiterate this claim: the relation (1) is very generic outcome of the framework, and it is not sensitive to any specific details of the model. The proposal [7,8] represents an alternative to baryogenesis scenario when the "baryogenesis" is actually a charge separation process in which the global baryon number of the Universe remains zero. In this model the unobserved antibaryons come to comprise the dark matter in the form of dense antinuggets in colour superconducting (CS) phase. The dense nuggets in CS phase also present in the system such that the total baryon charge remains zero at all times during the evolution of the Universe. The detail mechanism of the formation of the nuggets and antinuggets has been recently developed in refs. [9,10], and the present work can be considered as a continuation of these studies. We review the nal pressure provided by the axion domain wall, in contrast with strangelet models [1,4] which must be stable at zero external pressure.
basics elements of this proposal, its predictions and the observational consequences including presently available constraints, in next Section II. We highlight in few paragraphs below the main questions which will be addressed in the present work.
Our first paper [9] devoted to the formation of the nuggets was focused on the dynamics of a single nugget, its formation, and its time evolution. If the Universe were symmetric during the formation time, the equal number of nuggets and antinuggets would form. However, the coherent (on the enormous scale of the Universe) CP− odd axion field was present in the system. Precisely this asymmetry generates the preferential evolution of the system such that number of nuggets and antinuggets would not be identically the same, which was the main result of our paper [10]. As the total baryon charge is conserved in the system, the remaining baryon charge (not hidden in form of the nuggets and antinuggets) represents the visible baryonic matter contributing to Ω visible . This represents a qualitative explanation of generic consequence of this entire framework as represented by (1).
However, in the previous papers [9,10] our main objective was to understand a novel fundamental phenomenon on a qualitative level by simplifying the system itself and neglecting a large number of specific features of the system which may produce quantitative effects of order one, but cannot change the qualitative picture advocated in [9,10].
Indeed, in our previous studies we ignored a number of effects. In particular: a) We assumed that the energy per baryon charge is the same for hadronic and CS phases; b) We ignored the contribution of the axion domain wall surrounding the quark nugget in CS phase using a simple, order of magnitude estimate for the mass of the nugget, M nugget ∼ m p B; c)We neglected the contribution of the conventional propagating axions produced by misalignment mechanism (or by decay of the topological defects) by assuming that Ω dark is saturated by the nuggets. All these simplifications allowed us to relate Ω visible and Ω dark in terms of a single parameter c defined as a ratio of the hidden baryon charge in nuggets and antinuggets, see precise definition in eq. (4).
The main goal of the present work is to take into account all these complications and numerical factors and incorporate them into a new equation which would be generalization of eq. (4). We think that the timing for this task is quite appropriate due to a number of reasons. First, there is number of groups searching for the axion, see relevant references in next section. If the axion is found with definite mass m a it may not necessary saturate the entire dark matter density Ω DM . A finite portion of Ω DM could be related to the nuggets contribution which, according to refs. [9,10], always accompanies the conventional axion production. We refer to the original papers [11][12][13] and recent reviews [14][15][16][17][18][19][20][21][22] on the axion dynamics and recent axion searches.
Our goal is to relate these two complementary mechanisms of the axion dynamics: some portion of the ax-ion field will radiate the propagating axions, while another portion of the coherent axion field will lead to the nugget's formation. These two mechanisms accompany each other and we want to establish a quantitative relation between the two.
Secondly, it has been recently claimed that the nonobservation of the isocurvature perturbations and the tensor modes impose strong constraints on relation between the axion parameter f a and inflation scale H I . Our comment here is that such kind of constraints are normally imposed by assuming that the relic axions saturates the DM density. However, if the nuggets become the dominant contributor to Ω DM , the corresponding constraint drastically weakens because the axion's contribution to Ω DM becomes subdominant, and as a consequence, the isocurvature perturbations will be suppressed for the same axion mass m a and coupling f a . These few comments suggest that there is a number of relevant parameters in the theory which crucially depend on relative contribution of the relic axions versus nugget's portion to Ω DM . The goal of this work is to elaborate on those relations.
The paper is organized as follows. In section II we overview the basic idea and phenomenological consequences on the proposal [7,8] with emphasis how the "baryogenesis" is replaced by the charge separation effect [9,10]. In section III, we will introduce our notations and definitions of the parameters to be computed. In section IV, we will focus on the internal structure of the nuggets. In particular, we compute the energy per baryon charge of a single stabilized (anti)nugget at zero temperature. Then in section V, using the results of previous two sections, we produce a number of numerical plots and show that there is a large region of parametrical space in this model which is consistent with all presently available constraints. The corresponding constraints have been obtained from a variety of independent astrophysical, cosmological, satellite and ground based observations. Furthermore, we will show that this model is also consistent with known constraints from the axion search experiments. In section VI, we briefly summarize the main results of this paper and discuss the potential future directions. In particular, we speculate that unexpected correlations between the frequency of appearance of the solar flares and the intensity of the dark matter fluxes observed in the Sun, might be related to the antiquark nuggets studied in this work.

II. AXION QUARK NUGGET (AQN) MODEL
The AQN Model in the title of this section stands for the axion quark nugget model to emphasize on essential role of the axion field and avoid confusion with earlier models such as quark nuggets, strangelets, nuclearities mentioned in the Introduction.
The basic ideas of the original proposal [7,8] can be explained as follows: It is commonly assumed that the Universe began in a symmetric state with zero global baryonic charge and later (through some baryon number violating process, so-called "baryogenesis") evolved into a state with a net positive baryon number. As an alternative to this scenario we advocate a model in which the baryogenesis is actually a charge separation process in which the global baryon number of the universe remains zero. In this model the unobserved antibaryons come to comprise the dark matter in the form of dense nuggets of quarks and antiquarks in CS phase. The formation of the nuggets made of matter and antimatter occurs through the dynamics of shrinking axion domain walls, see original papers [9,10] with the technical details.
If the fundamental θ parameter of QCD were identically zero during the formation time, equal numbers of nuggets made of matter and antimatter would be formed. However, the fundamental CP violating processes associated with the θ term in QCD result in the preferential formation of antinuggets over nuggets. This source of strong CP violation is no longer available at the present epoch as a result of the axion dynamics when θ eventually relaxes to zero, see Fig. 1. As a result of these CP violating processes the number of nuggets and antinuggets being formed would be different. This difference is always of order of one effect irrespectively to the parameters of the theory, the axion mass m a or the initial misalignment angle θ 0 , as argued in [9,10]. As a result of this disparity between nuggets and antinuggets a similar disparity would also emerge between visible quarks and antiquarks. Precisely this disparity between visible baryons and antibaryons eventually lead (as a result of the annihilation processes) to the system when exclusively one species of visible baryons remain in the system, in agreement with observations.
We illustrate this interrelation between the axion dynamics and the nugget's formation on Fig. 1 when both processes occur during the same QCD epoch. The key element of the proposal is the presence of the coherent axion θ(t) field which continues to oscillate when the nugget's formation starts. The time evolution of the nuggets continues until low temperature T form when the axion field assumes its final destination θ = 0 which is observed today. However, the originally produced asymmetry between nuggets and antinuggets cannot be washed out as explained in [10].
This disparity between nuggets and antinuggets unambiguously implies that the total number of antibaryons will be less than the number of baryons in early universe plasma. This is precisely the reason why the resulting visible and dark matter densities must be the same order of magnitude (1) in this framework as they are both proportional to the same fundamental Λ QCD scale, and they both are originated at the same QCD epoch. If these processes are not fundamentally related, the two components Ω dark and Ω visible could easily exist at vastly different scales.
Another fundamental ratio is the baryon to entropy ratio at present time In our proposal (in contrast with conventional baryogenesis frameworks) this ratio is determined by the formation temperature T form 41 MeV at which the nuggets and antinuggets complete their formation. We note that T form ≈ Λ QCD . This temperature is determined by the observed ratio (2). The T form assumes a typical QCD value, as it should as there are no any small parameters in QCD, see Fig. 1.
Unlike conventional dark matter candidates, such as WIMPs (Weakly Interacting Massive Particles) the darkmatter/antimatter nuggets are strongly interacting and macroscopically large, nuclear density objects with a typical size (10 −5 − 10 −4 ) cm, and the baryon charge which ranges from B ∼ 10 23 to B ∼ 10 28 . However, they do not contradict any of the many known observational constraints on dark matter or antimatter for three main reasons [23]: • They carry very large baryon charge |B| > 10 23 which is determined by the size of the nugget ∼ m −1 a . As a result, the number density of the nuggets is very small ∼ B −1 . Their interaction with visible matter is thus highly suppressed. In particular, the quark nuggets essentially decouple from cosmic microwave background (CMB) photons, and therefore, they do not destroy conventional picture for the structure formation; • The nuggets has a huge mass M nugget ∼ m p B, therefore the effective interaction is very small σ/M nugget ∼ 10 −10 cm 2 /g, which is evidently well below the upper limit of the conventional DM constraint σ/M DM < 1cm 2 /g; • The quark nuggets have very large binding energy due to the large gap ∆ ∼ 100 MeV in the CS phase. Therefore, the strongly bound baryon charge is unavailable to participate in the big bang nucleosynthesis (BBN) at T ≈ 1MeV, long after the nuggets had been formed.
We emphasize that the weakness of the visible-dark matter interaction in this model is due to a small geometrical parameter σ/M ∼ B −1/3 which replaces the conventional requirement of sufficiently weak interactions for WIMPs. In this short overview of the model we also want to make few comments regarding the observational consequences of the proposal. It is known that the galactic spectrum contains several excesses of diffuse emission the origin of which is uncertain, the best known example being the strong galactic 511 keV line. Our comment here is that this model offers a potential explanation for several of these diffuse components (including 511 keV line and accompanied continuum of γ rays in 100 keV and few MeV ranges, as well as x-rays, and radio frequency This diagram illustrates the interrelation between the axion production due to the misalignment mechanism and the nugget's formation which starts before the axion field θ relaxes to zero. Possible cooling paths are denoted as path 1, 2 or 3. The phase diagram is in fact much more complicated as the dependence on the third essential parameter, the θ is not shown as it is largely unknown. It is assumed that the initial axion field starts at θ 0 , while the nuggets complete their formation in the CS region at T form ≈ 41 MeV, µ > µ c and θ ≈ 0. bands). For further details see the original works [24][25][26][27][28][29][30] with specific computations in different frequency bands in galactic radiation, and a short overview [31].
As we already mentioned the typical size of a nugget ∼ m −1 a is determined by the axion mass, yet to be discovered. Therefore, a typical baryon charge of the nuggets and their size distribution are also unknown characteristics of the model. However, it is important to emphasize that the allowed window for the axion mass 10 −6 eV ≤ m a ≤ 10 −2 eV corresponds to the range of the nugget's baryon charge B which largely overlaps with all presently available and independent constraints on such kind of dark matter masses and baryon charges see e.g. [31,32] for review 3 . To reiterate the same claim: the allowed window (3) for the baryon charge B is perfectly consistent with all presently available constraints. 3 The smallest nuggets with B ∼ (10 23 − 10 24 ) naively contradict to the constraints cited in [32]. However, the corresponding constraints are actually derived with the assumption that nuggets with a definite mass (smaller than 55g) saturate the dark matter density. In contrast, we assume that the peak of the nugget's distribution corresponds to a larger value of mass, B ≥ 10 25 , while the small nuggets represent a tiny portion of the total dark matter density. The same comment also applies to the larger masses excluded by Apollo data as reviewed in [31]. Large nuggets with B ∼ 10 28 may exist, but represent a small portion of the total dark matter density, and therefore, do not contradict the Apollo's constraints, see also some comments on the baryon charge distribution (34) in concluding section VI.
Furthermore, if one takes the recent conjecture [33,34] that the extreme ultra violet (EUV) and soft x-ray emission from the solar Corona are due to the antiquark nuggets hitting the Sun from outer space, then one can make one step further and identify the observed "nanoflares" with annihilation events of the antinuggets in the solar corona. Taken this identification literally one can translate the observed energy distribution of the "nanoflares" in solar corona in terms of the nugget's baryon charge distribution. This identification leads to an independent estimate on possible values for B which is perfectly consistent with the allowed range (3) for the baryon charge of the nuggets [33,34]. Therefore, in our analysis which follows we treat the window (3) as the solid constraint on the allowed magnitude of the nugget's baryon charge.

III. THE AQN MODEL. NOTATIONS AND DEFINITIONS
In this section, we will first briefly introduce our notations for the AQN model which is constituted by two parts: one portion is represented by the nuggets and antinuggets, while the other component is represented by free propagating axions which were produced by the axion field (through e.g. the misalignment mechanism). In order to quantitatively describe this model, we have to derive a general expression which relates these two complimentary contributions to the Ω DM as they both originated from one and the same axion field during the same QCD epoch as illustrated on Fig. 1.
To proceed with this task we first remind a simplified expression derived in our previous papers [9,10] where a number of effects have been neglected as mentioned in Section I. In the simplified treatment we characterized the system by a single crucial parameter "c" which describes the disparity between nuggets and antinuggets and defined as follows 4 : The most important consequence of this model is given by relation (1). It is expressed now in terms of numerical parameter c as eq. (4) states. As we already stated this equation represents a very generic consequence of the entire framework. The parameter 0 < c < 1 is obviously positive and cannot be greater than 1 as eq. (4) suggests. This is because the baryon charge hidden in the nuggets should be less than the baryon charge hidden in the antinuggets as the difference represents the baryon charge of the visible matter in hadronic phase we presently observe.
As we already mentioned in Introduction in our previous studies [9,10] the relation (4) is not accurate as we previously ignored a number of numerical factors. In particular, we assumed that the energy per baryon charge for hadronic and CS phases is the same; we ignored the contribution of the axion domain wall surrounding the quark nugget in CS phase; we neglected the contribution of the conventional propagating axions produced by misalignment mechanism (or by decay of the topological defects) by assuming that Ω dark is saturated by the nuggets. All these simplifications allowed us to derive a simple formula (4) in terms of a single parameter c, which is very important for qualitative, but not quantitative, arguments. Now we want to generalize relation (4) by accounting for all these (previously neglected) effects. The corresponding modifications obviously do not affect the basic qualitative claim (1), but may change some numerical factors, which is precisely the main objective of the present studies.
Traditionally, the axion is regarded as one of the leading candidates for DM, see e.g. recent reviews [20][21][22] and references therein. It is normally assumed that the initial coherent axion field radiates the propagating axions through the misalignment mechanism or as a result of decay of the topological defects.
However, these two distinct processes of axion's field relaxation will be always accompanied by another complimentary mechanism (related to the so-called N DW = 1 axion domain walls) as recently advocated in refs. [9,10]. The corresponding scenario is overviewed in Sections I and II, and shall not be elaborated here. The only important elements to be mentioned now are: a) this scenario assumes that the PQ symmetry breaking occurs before or during inflation, and b) this scenario leads to a very generic consequence (1) which holds irrespectively to any parameters of the system such as m a or initial misalignment angle θ 0 . This generic prediction of this new proposal should be contrasted with two other well known and well studied mechanisms of the axion field's relaxation when the corresponding contribution to the DM density is highly sensitive to the axion parameters, a . We now proceed with our definitions and notations of the AQN model which incorporates both contributions: the conventional axion production reviewed in [20][21][22] and the nugget's contribution studied in [9,10]. As usual, we use the ratio of different components' densities to critical density today to mark their proportions where label i stands for different species: b for baryon; '+' for nugget and '−' for antinugget; a for free axions. ρ cr is the critical density of the Universe today. Then we have where Ω DM is represented by two parts, the axion contribution Ω DM (axion) ≡ Ω a and the nugget's portion Ω DM (nugget) ≡ (Ω + + Ω − ). Next, to describe the difference between nugget and antinugget, we define the following parameters: and N ± is the number density of (anti)nuggets in the Universe today; E ± , B ± and ± = E ± /B ± are respectively the energy, baryon number and energy density per baryon charge for a single nugget or antinugget. In these relations both baryon charges B + and B − are defined to be positive numbers. The definition of c in (8) coincides with the definition given in (4).
As we already mentioned after eq. (4) we neglected a number of important numerical factors in our previous studies [9,10]. We can now formalize these effects in a very precise way using our definitions (7) and (8). In particular, ± is the energy per baryon charge in CS phase is not the same as in hadronic phase, i.e. ± = m p /3. Furthermore, E ± which was previously estimated as E + = E − = Bm p now includes the contribution from the surrounding axion domain wall, and obviously has much more complicated structure. In addition, we previously ignored the contribution from free propagating axions by assuming Ω a = 0. All these new elements will be incorporated into our analysis in section IV.
To simplify our analysis we want to make a technical assumption that N + N − in eq. (7). It does not affect any of our conclusions as we argue below. The basic justification for this assumption is as follows: the initial distributions of the nuggets and antinuggets is the same. It is just their evolution in the background of the coherent axion CP-odd field generates the asymmetry between them as argued in [10]. As time evolves the nuggets and antinuggets assume different size distribution, but their number densities, N + and N − should stay the same 5 . This information is coded in equations (7) and (8)  Substituting N + = N − in eq. (7) we arrive to desired expression This equation has obvious physical meaning and essentially states that the relative contribution of the nuggets 5 There is some loophole in this argument as smaller bubbles may completely collapse, while larger bubbles survive the evolution.
As the size-distribution for B − is different from B + the corresponding number densities for N + and N − could be also different as a result of different rates for collapsed nuggets. We expect that this effect is quite small. In any event, this effect can be easily accommodated by redefinition of the effective B eff + and B eff − accounting for the collapsed nuggets. and antinuggets to Ω DM is proportional to the corresponding baryon charges expressed in terms of parameter ∼ c, and it is also proportional to the difference of their energy densities parametrized by c . Now we want to derive an equation, similar to (4) which accounts for a number of the effects which were previously ignored, including the differences between nuggets and antinuggets. With this purpose we express the baryon number conservation in the following form where m p is the mass of a single baryon charge and can be approximated by the proton mass. The coefficient 3 in eq. (10) corresponds to our normalization of the baryon charge in the present work. This normalization is consistent with our definition of µ corresponding to the quark (rather than baryon) chemical potential. Furthermore, B ± in all our formulae are positive numbers as mentioned above. Therefore, the B ± parameters count for the number of quarks in the system, rather then their baryon charges. With these comments in mind and using eqs. (5), (7) and (10), we arrive to the following relations The coefficient c ∈ (0, 1) in these relations satisfies the same constraint as in our simplified treatment of the problem presented in (4). Therefore, Our next step is to rewrite the equations (6) and (7) in the following convenient form The last step to achieve our goal is to use eqs. (11) and (9) to arrive to the final expressions which will be used in our numerical studies in SectionV: The equation (14) is a generalization of our previous simplified expression (4). It accounts for a number of numerical effects mentioned previously. Formula (14) obviously reduces to our previous expression (4) in the limit when Ω a = 0 and nuggets and antinuggets have the same energies, i.e., + = − = m p /3 such that c = 1.
While the numerical estimates for parameters 3 ± /m p and c entering (14) will be discussed in details in next section IV, the rest of this section is devoted to a short overview of known estimates of the parameter Ω a which also enters eq. (14).
The corresponding computations of Ω a have been carried out in a number of papers. In what follows we limit ourselves by reviewing the estimates of Ω a resulted from the misalignment mechanism [35], while leaving out the contribution related to the decay of the topological defects 6 . We emphasize that we exclude the corresponding contributions related to the topological defects not because they are unimportant. Rather, we omit them because their role is largely unknown under present circumstances when the PQ symmetry is broken before/during inflation 7 . In addition, even in a different scenario when the PQ symmetry broken after inflation, the question whether it saturates the observed dark matter density remains controversial as mentioned in footnote 6. Thus, we leave out this contribution to simplify our notations and our analysis as the focus of the present work is the nugget's contribution to Ω DM rather than the direct axion production represented by Ω a . In other words, Ω a contribution is kept in our formulae for normalization purposes to illustrate the significance (or insignificance) of the nugget's contribution to Ω DM as a function of parameters. The contribution to the Ω a due to the decay of the topological objects can be always incorporated into our formulae once the uncertainties of this contribution are better understood.
For the vacuum misalignment production of free propagating axions, we adopt the general formula as presented in ref. [20] Ω a h 2 ≈ 2 × 10 4 f a 10 16 GeV with where H I is the inflationary Hubble scale and F anh (x) is the correction factor due to the anharmonic cosine part 6 There is a number of uncertainties and remaining discrepancies in the corresponding estimates. We shall not comment on these subtleties by referring to the original papers [36][37][38][39][40]. According to the most recent computations presented in ref. [40], the axion contribution to Ω DM as a result of decay of the topological objects can saturate the observed DM density today if the axion mass is in the range ma = (2.62 ± 0.34)10 −5 eV, while the earlier estimates suggest that the saturation occurs at a larger axion mass. One should also emphasize that the computations [36][37][38][39][40] have been performed with assumption that PQ symmetry was broken after inflation. 7 It was explained in [9] in great details that the so-called N DW = 1 domain walls must be present in the system even when PQ symmetry was broken before/during inflation, see item 5 in section III of [9]. It should be contrasted with conventional studies on the role of the topological defects when the PQ symmetry is assumed to be broken after the inflation, see also footnote 6 with related comment.
in axion potential [20,41]. The parameter θ a,i in (15) is the initial misalignment angle and H I /(2πf a ) is the backreaction contribution to this homogeneous field displacement due to the isocurvature perturbations. The parameter f a and the axion mass m a are not independent parameters as their product is fixed by the topological susceptibility of QCD, χ = f 2 a m 2 a . Using the recent value χ = 0.0216 fm −4 at zero temperature [42], we have This completes our short overview of Ω a contribution entering our basic formulae (14). The next section IV is mainly devoted to the estimates of 3 ± /m p and c which, along with Ω a also enter our basic equation (14). Then in section V, using the results of sections III and IV, we do numerical analysis to study the allowed window and constraints related to the phenomenological parameters of the AQN dark matter model.
We conclude this section with the following generic comment. The nugget's contribution given by Ω ± and the direct axion production represented by Ω a always accompany each other during relaxation of the dynamical axion field to zero at the QCD epoch. These contributions to Ω DM represent complementary mechanisms and cannot be formally separated (e.g. by variation of a free parameter of the system such as f a ) as the closed N DW = 1 domain wall bubbles, which are responsible for the nugget's formation, can be produced irrespectively whether the PQ scale is above or below the inflationary scale H I , as argued in [9].

IV. INTERNAL STRUCTURE OF THE NUGGETS
The main goal of this section is to estimate the energy per baryon charge of a stabilized quark nugget at zero temperature. We consider two drastically different models to accomplish this task. The main reason to consider two different models (constructed on fundamentally different principles) is to test the sensitivity (or non-sensitivity) to different phenomenological parameters effectively describing the strongly coupled QCD. For simplicity and without loss of generality, we assume that the CS phase assumes the simplest possible structure in the form of the colour flavour locked (CFL) phase without any additional complications such as possible meson condensation.
The first model largely follows the original work [7]. However, the difference with the previous analysis is that the first paper on the subject [7] was mostly dealing with fundamental and basic questions on principle possibility to stabilize the nuggets by the axion domain walls. The goal of the present studies is quite different as we want to produce some quantitative results on the parameters entering the basic equation (14).
The model [7] considers the equilibration between the Fermi pressure, the domain wall surface tension, the "bag constant" pressure ∼ E B , and, finally, the quark-quark interaction related to the CS gap. The energy of a stabilized nugget can be represented in the following form while the nugget's baryon number can be estimated as follows where our normalization corresponds to B = 1 per single quark degree of freedom in order to remain consistent with notations of Sec. III. Few comments are in order. First of all, the "squeezer" parameter µ 1 is estimated to be 330 MeV, the degeneracy factor for CFL phase is estimated as g 2N c N f 18. Parameters µ and R entering (18) and (19) are the chemical potential and radius of the nugget, respectively, while the bag constant E B (150 MeV) 4 and the CS gap ∆ 100 MeV assume their commonly accepted magnitudes, see original work [7] with some arguments supporting these numerical values.
The domain wall tension σ eff entering (18) requires some additional comments. First of all, the effective domain wall tension σ eff should not be confused with the conventional surface tension σ 8f 2 a m a which normally enters the computations [36][37][38][39][40] of the axion production due to the decay of the topological defects.
There are two main reasons for this important difference. First of all, the axion domain wall solution in our case interpolates between topologically distinct vacuum states in hadronic and CS phases, in contrast with conventional axion domain wall which interpolates between topologically distinct hadronic vacuum states. The chiral condensate may or may not be formed in CS phase. It strongly affects the topological susceptibility in CS phase which could be much smaller than in the conventional hadronic phase. The well known manifestation of this difference is the expected smallness of the η mass in CS phase in comparison with the hadronic phase.
One should emphasize the 2π periodicity of the axion (θ) and the Nambu-Goldstone fields (η ) still hold in presence of the chemical potential µ in dense matter CS phases [43]. Therefore, the topological reason for mere existence of the axion domain wall still persists, while the numerical value of the tension σ eff will deviate from its conventional expression σ computed in hadronic phase.
The second reason for strong deviation of the σ eff from conventional expression for σ is that formula σ 8f 2 a m a was derived assuming the thin-wall approximation when the domain wall is assumed to be almost flat, i.e. a typical curvature of the domain wall structure is much smaller than its width. This approximation is obviously badly violated because the axion domain wall width is of order m −1 a , while typical curvature is comparable with the width of the domain wall as these two parameters are related in our framework, R ∼ m −1 a . The physical consequence of this relation is that the axion field strongly overlaps within the nugget's volume. This effect is expected to drastically reduce the domain wall tension 8 .
To account for these complicated QCD effects we define σ eff ≡ κ · σ, with an unknown phenomenological parameter 0 < κ < 1 which accounts for the physics mentioned above. In particular, the violation of the thinwall approximation was modelled in [9] by introducing a suppression factor exp(−R 0 /R form ). The corresponding suppression could be quite strong and can be as small as 10 −5 assuming a typical formation radius R form ∼ 0.1R 0 as studied in [9]. In what follows we treat κ as a free phenomenological parameter.
Our goal now is to minimize the expression (18). To achieve this goal we introduce two dimensionless variable x and σ 0 as where we express σ in terms of topological susceptibility χ = f 2 a m 2 a = 0.02fm −4 . Then, the energy per quark ≡ E/B can be expressed as QCD (x), The equilibrium point can be found following the condition The solution can be well approximated from numerical computation as 8 The corresponding large modifications can be understood from a simple model when the domain wall is bended allowing a strong overlap between opposite sides of the wall. The effective domain wall tension obviously receives the modifications as a result of this bending geometry when the axion field configuration deviates from a simple well-known 1-dimensional solution.
This solution has accuracy up to 0.6% comparing to the exact numerical solution within the range (23b), see Appendix A with more technical details. This energy per baryon charge will be used in next section for our numerical analysis of this AQN dark matter model to test its consistency with all available constraints. Our second model to be considered has very different building principle, and essentially is based on the ideas of old constituent quark model being applied to the dense matter systems [44]. This model will serve as a complementary tool which allows us to test the sensitivity of our results to different phenomenological parameters effectively describing the strongly coupled dense QCD.
We assume the following form for the energy of the nuggets Here M q is the effective constituent quark's mass in CS phase in the bulk of the nugget. The M q should have a typical QCD scale and serves as a parameter of the model. The energy per quark degree of freedom is therefore where we substitute R B 1/3 = 1 µ 9π 2g 1 3 from Eq. (19). We should note that M q and µ are not completely free parameters according to various phenomenological models for dense phases. To be consistent with previous studies (applied to the neutron star physics where this model was used) we adopt the following numerical values for relevant parameters [44]: (M q , µ) (200, 400) MeV, and (160, 500) MeV. (26) In spite of the fact that the Model-2 has a fundamentally different building principle, we observed that these choices (26) for (M q , µ) parameters produce the results which numerically very close to the results (23) obtained for Model-1, see Appendix A with detail analysis. We conclude this "technical section" with the following comment. Both models discussed in this section and represented by eqs. (23) and (25) are based on fundamentally different principles. Nevertheless, both models lead to similar results for relevant parameters, and furthermore, demonstrate that the energy density in both models depends on a single dimensionless parameter κ −1 B 1/3 m a /m π , which is clearly a highly nontrivial feature. A more detailed discussions of these two models can be found in Appendix A.

V. THE AQN MODEL CONFRONTING THE OBSERVATIONS
The main purpose of this section is to analyze all possible constraints on the parameters of the AQN model. In previous sections we introduced a number of phenomenological parameters describing this model. Our goal here is to analyze the allowed region in parametrical space where these parameters may vary but remain consistent with all known observations/experiments. We use the result (23) or (25) from section IV to fit the parameters ± . We use (15) from section III to fit the density of the relic dark matter propagating axions into our AQN model. Finally, we implement constraints on physically observable parameters B, θ a,i , m a , H I from large number of independent experiments and observations to analyze the allowed window for most important parameter of the AQN model represented by the coefficient "c". This key parameter is defined by eq. (8) and describes the disparity between nuggets and antinuggets. As we argued in [10] this coefficient c must be order of one as a result of interaction with coherent CP odd axion field. This coefficient, in principle, is calculable from the first principles along with other parameters of the model as all fields, coupling constants and interactions are represented by the SM physics accompanied by the axion field θ(x) with a single additional fundamental parameter f a . However, such computations presently are not feasible as even the phase diagram shown on Fig.1 at θ = 0 is not yet understood.
We represent our numerical results in subsections V B and V C. However, first of all, in next subsection V A we overview the known constraints on relevant parameters of the AQN model.

A. Constraints on parameters of the AQN model
We start with eq. (12) which defines coefficient c. This should not be considered as a constraint on c as it essentially represents our convention that we define the visible matter as the baryons with positive baryon charge. Therefore, the absolute value of the baryon charge hidden in the antinuggets B − must be greater than the baryon charge B + hidden in the nuggets (as a result of global conservation of the baryon charge which is assumed to be zero at all times). This leads to the formal relation c < 1 which reflects our convention. The parameter c is obviously a positively defined parameter which is explicitly represented by eq. (8).
Another constraint 0.2 κ −1 B 1/3 m a /m π 0.95 follows from (23b). This constraint is related to our studies of the stability of the nuggets in CS phase, see section IV with more technical details.
The next item we want to discuss is the known constraints on the baryon charge B of the nuggets (or antinuggets) represented by (3). We already made few comments about the constraints related to the galactic obser-vations, Apollo data, and ancient Mica analysis as mentioned in section II, footnote 3 and reviewed in [31,32]. Now we want to make few comments regarding other constraints related to different observations and analysis.
First of all we want to mention the constraints related to limits from the total geothermal energy budget of the Earth [45] which is consistent with (3). It has also been suggested that the anita experiment may be sensitive to the radio band thermal emission generated by these objects as they pass through the antarctic ice [45]. These experiments may thus be capable of adding direct detection capability to the indirect evidence mentioned previously in section II.
It has been also suggested recently [46] that the interactions of the antinuggets with normal matter in the Earth and Sun will lead to annihilation and an associated neutrino flux. Furthermore, it has been claimed [46] that the antiquark nuggets cannot account for more than 20% of the dark matter flux based on constraints for the neutrino flux in 20-50 MeV range where sensitivity of the underground neutrino detectors such as Su-perK have their highest signal-to-noise ratio. However, the claim [46] was based on assumption that the annihilation of visible baryons with antiquark nuggets generate the neutrino spectrum similar to conventional baryonantibaryon annihilation spectrum when the large number of produced pions eventually decay to muons and consequently to highly energetic neutrinos in the 20-50 MeV energy range. This claim has been dismissed in ref. [47] by emphasizing that antinuggets cannot be treated as an usual antimatter in conventional hadronic phase as the quarks in nuggets (and antiquarks in antinuggets) belong to CS phase rather than to the hadronic phase we are familiar with.
One could also expect a set of strong constraints on the antinuggets by considering the radio observations or strong 511 keV line emission from nearby galaxies if one assumes that the antinuggets (almost) saturate the corresponding emissions in our Milky Way. Indeed, similar consideration in other models essentially rules out the dark matter explanation of 511 keV line observed in our galaxy. However, the corresponding analysis carried out in [48] for the radio radiation and in [49] for 511 keV line emission from nearby galaxies does not lead to any new constraints, in huge contrast with conventional WIMPs models which typically predict the radiation from nearby galaxies exceeding the observed values. We refer to the original papers [48,49] for the details and references.
Our last comment related to constraints on B as given by (3) is related to the recent arguments [33] advocating that that the extreme ultra violet (EUV) and soft x-ray emission from the solar corona might be due to the antiquark nuggets hitting the Sun from outer space. If one identifies the observed "nano-flares" with annihilation events of the antinuggets in the solar corona one can infer the nugget's baryon charge distribution from the (measured) energy distribution of the "nano-flares" in solar corona, see also some related comments in concluding section VI. This identification leads to an independent estimate [33] on possible values for B which is perfectly consistent with the previous studies, which we express as 10 23 B 10 28 .
In our numerical analysis which follows we treat the window (27) as the solid constraint of the allowed magnitude of the nugget's baryon charge. We next consider the classical window 9 for axion mass, see recent review paper [21]: By using the relation (17) this window for m a can be expressed in terms of the corresponding classical window for f a : 5.7 × 10 8 GeV f a 5.7 × 10 12 GeV. (29) One should emphasize that the constraint (28) or equivalently (29) is commonly accepted axion window, and it is not, by any means, originated from our analysis of the AQN model. Nevertheless, all our constraints depend on m a as it explicitly enters the equation (23b). From these discussions it is clear that the axion mass m a plays a dual role in our analysis because it enters the formulae related to the physics of the nuggets as equation (23b) states. It also enters the expression (15) for Ω a . It unambiguously implies that the remaining portion of the dark matter represented by the nugget's contribution (13) becomes also (implicitly) highly sensitive to m a through dependence of the axion portion of the dark matter represented by Ω a .
As we argue below, for a reasonable values of κ in range 10 −4 κ 10 −2 , the constraints (23b), (27) and (28) become mutually compatible which we consider as a highly nontrivial consistency check as all the parameters entering these relations have been constrained by very different physics related to independent observations, experiments and analysis.
The next constraint to consider is related to analysis of the inflationary scale H I and the related constraints on the tensor-to-scalar ratio and the isocurvature perturbations. The basic assumption of this work is that PQ symmetry breaking occurs before/during inflation, in which case 9 For the main purposes of this paper, we will only consider "the classical axion window", where the initial misalignment angle θ a,i is not fine-tuned. Note that while the upper bound is a very solid constraint as it is given by stellar physics (e.g. see review [17]), the lower bound on the axion mass in Eq. (28) should be treated as an order of magnitude estimate provided that θ a,i is not fine-tuned. If the fine tuning is allowed, the θ a,i may assume arbitrarily small value, in which case the corresponding lower bound on ma is shifted. The only exclusion interval in this case 6 × 10 −13 eV < ma < 2 × 10 −11 eV is obtained from black hole superradiance effects [50]. see e.g. [20] for review. This assumption plays a crucial role in our analysis [9,10] because the CP odd axion field must be coherent on enormous scale of the entire Universe to separate the baryon charges on these gigantic scales with the same sign of θ. Precisely this coherent axion field generates the disparity between the nuggets and antinuggets which eventually leads to the generic and fundamental prediction (1) of this entire framework. It is known that the Inflationary Hubble scale is tied to the value of the tensor-to-scalar ratio r T which measures H I . Assuming a simplest single field inflationary model, the non-observation of the tensor modes (r T < 0.12) imposes the upper limit for the inflation scale, see [20,51]: Important comment here that (31) is highly modeldependent result and varies from one inflationary model to another. It is presented here exclusively for illustrative purposes to provide some orientation with the relevant scales of the problem. The isocurvature perturbations related to the axion field provides another independent constraint on H I . We recall that the amplitude for the isocurvature power spectrum is determined by the following expression, see the original papers [52][53][54][55][56] and review article [20]: The corresponding isocurvature amplitude is strongly constrained by CMB analysis, A I /A s < 0.038, where A s is conventional amplitude for the scalar power. It is normally assumed that the non-observation of the isocurvature perturbation provides a strong constraint on axion properties in scenario where the PQ symmetry broken before/during inflation. Our original comment here is that the axion contribution represented by Ω a in eq. (32) to the dark matter density Ω DM could be numerically quite small in the AQN model as the nuggets in most cases play the dominant role by saturating the dark matter density. Such a scenario drastically alleviates some severe constraints on parameters in conventional analysis where one normally assumes that the axions saturate the dark matter density. We conclude this subsection with the following remark. The conventional analysis on the relation between dark matter axions, inflationary scale, isocurvature perturbations very often assumes that the axions saturate the dark matter density. It should be contrasted with our AQN model where the axions themselves with the same f a may contribute very little to Ω DM as the dominant contribution may come form the nuggets which always satisfy the relation Ω DM ∼ Ω visible according to (1) irrespectively to f a or initial misalignment angle θ a,i . It may alleviate some severe constraints on the parameters (such as H I , f a , r T , θ a,i ) which other models normally face.

B. Numerical plots
The goal of this subsection is to analyze the dependence of the internal (with respect to the AQN model) parameter "c" from external parameters of the system such as B, θ a,i , m a , H I which are well defined observables irrespectively to the specific features of the AQN model. As the parameter "c" cannot be negative or larger than one, the corresponding plots provide us with information on the typical values of the external parameters B, θ a,i , m a , H I when the AQN model is self-consistent with all presently available constraints.
We start our analysis by plotting on Fig. 2a the parameter "c" as a function of m a and B, where we fix specific values for parameters κ = 10 −4 and H I /2π = 5.7 × 10 8 GeV and θ a,i = 10 −3 to simplify the arguments and analysis. We also plotted the f a = H I 2π by red solid line to localize the physical parametric space and remove unphysical (within the AQN model) solutions. We also plotted (by green and blue dashed lines) the region in parametrical space where the condition (23b) is satisfied and our computations in CS phase are justified. For this specific choice of the parameters one can explicitly see that parameter c is constrained in a parallelogram with the range 0.4 c 0.6. This region of the parametrical space satisfies all internal and external constraints listed in previous section.
From the same plot one can also identify the allowed region of the nugget's baryon charge B for a given axion mass m a . One should emphasize that the dark matter density on Fig. 2a for θ a,i = 10 −3 is entire saturated by the nuggets as the direct axion production is strongly suppressed by small initial misalignment angle θ a,i = 10 −3 . To see the role of the direct axion production one can choose θ a,i = 10 0 as shown on Fig. 2b. In this case the direct axion production saturates the dark matter density Ω DM = Ω a at small axion masses m a 10 −5 eV, as shown by solid yellow line. When 10 −3 < θ a,i < 1 varies between these two values the corresponding allowed region for the nuggets's parametrical space (B, m a ) will be modified accordingly as the allowed region for the nuggets obviously shrinks when the direct axion production starts to play an essential role.
The key observation here is that there will be always a region (B, m a ) when the total dark matter density assumes its observational value through the parameter c which determines the nugget's contribution to Ω DM . The corresponding contribution varies to accommodate the related axion portion Ω a as the total dark matter density Ω DM on Fig. 2 is fixed and assumes its observational value.
On Fig. 3 we wish to demonstrate the sensitivity of the allowed region for region (B, m a ) with respect to H I /2π parameter, where the constraint (31) is applied.
We want to demonstrate the corresponding sensitivity to H I /2π parameter by showing that there is no any dependence on H I for sufficiently small H I /2π 10 10 GeV as shown on Fig. 3a. In all respects the plot is very much the same as the one shown in Fig. 2a. In both cases the dark matter density is dominated by the nuggets, and the the allowed region (B, m a ) is not sensitive to the H I /2π as long as the Hubble parameter is sufficiently small. However, when H I /2π becomes close to f a ∼ 10 11 GeV, the window for c is shifted as shown on Fig. 3b to accommodate the conventional contribution of the propagating axions. The main point is that there will be always a region (B, m a ) when the total dark matter density assumes its observational value, though the magnitude of c ∈ (0, 1) assumes somewhat different values, depending on external parameters of the system.
Our next task is to analyze the sensitivity of our results to the QCD parameters related to CS properties of the nuggets. To accomplish this goal we plot parameter c on Fig. 4 as a function of m a and B using (2) tot for Model-2 determined by eqs. (25) and (26). The corresponding plot for (M q , µ) = (200, 400) MeV is presented on Fig. 4a while for (160, 500) MeV is shown on Fig. 4b respectively.
The main conclusion is that the Model-2 (which is based on the fundamentally different building principles than Model-1) with varies parameters produces nevertheless quantitatively similar results as Model-1 analyzed previously and shown on Fig. 2a. A more detail comparison presented in Appendix A supports this claim. This conclusion essentially implies that our phenomenological results are not very sensitive to the specifics of the QCD parametrization of the system describing the dense CS phase of matter in strongly coupled regime. Therefore, we treat our results as the solid consequences of the AQN model.
As an additional note, the parameter c Ω as defined by eqs. (7) and (9) and which describes the mass-difference between the nuggets and antinuggets (in contrast with parameter c which describes the baryon charge-difference between the nuggets and antinuggets) is numerically very close to parameter c studied above. Specifically, one can show that for Model-1 (23) the parameter c Ω 1.17c within 15% accuracy for the region 0.4 c 0.6 which dominates the parametrically allowed region as discussed in the preceding paragraphs. Therefore, we do not show the plots for c Ω as a function of external parameters because they are very similar to the plots for c presented and discussed above.

C. No fine-tuning in the AQN scenario
As we mentioned in section III, the dark matter propagating axion itself may not saturate the total dark matter. In contrast, the nugget's formation always generates a large contribution Ω DM ∼ Ω visible and always accompanies the conventional axion production. This property of the AQN model is demonstrated on Fig. 5, where we plot Ω a /Ω DM as a function of m a and φ = θ 2 a,i + (H I /(2πf a )) 2 . The function φ(θ a,i , H I ) enters formula (15) for Ω a and counts together with the initial  homogeneous displacement contribution and the backreaction contribution to the free dark matter axions. The Fig. 5 explicitly shows that that m a and φ have to be highly fine tuned to make Ω a to saturate Ω DM exactly, shown as a bright green solid line. In other words, for a specific magnitude of m a there is a single value of φ when the dark matter density assumes its observable value. Once these two parameters, m a and φ slightly deviate from the appropriate values, the Ω a strongly deviates from Ω DM .
This conventional fine-tuning scenario should be contrasted with the results of the AQN model when Ω a may contribute very little to Ω DM . Nevertheless, the Ω DM assumes its observation value as a result of an additional nugget's contribution which always accompanies the axion production and always generates a contribution of order one, as we already emphasized. In other words, the AQNs play the role of the "remaining" DM density which, in fact, could be the dominating portion of the Ω DM . As we have seen in previous subsection V B, there . will always be contribution to the DM model constituted of free axions and the AQNs for the allowed parametric space. In other words, for a specific magnitude of m a there is a large window of φ corresponding to the different values of parameter c ∈ (0, 1) when the dark matter density assumes its observable value. Different colours on the Fig. 5 correspond to the different axion contribution represented by parameter Ω a /Ω DM . Therefore, the fine-tuning problem does not even occur in the AQN scenario as entire parametrical space on the left from the green solid curve Ω a = Ω DM corresponds to the observable dark matter density. The white region in Fig. 5 is excluded because the requirement Ω a ≤ Ω DM is violated.
As the final technical remark we also notice from Fig. 5 that for most part of the parameter space, we have Ω a 0.2Ω DM . To make the above statement more precise, we plot on Fig. 6 the ratio Ω a /Ω DM as a function of c and parameter κ −1 B 1/3 m a /m π determined by the QCD physics as given by eq. (23). The white region in Fig. 6 stands for the excluded region of parameters (c, κ −1 B 1/3 m a ). This plot shows that the parameter c cannot be very close to ∼1 for the allowed QCD window 0.2 κ −1 B 1/3 m a /m π 0.95. This property, in fact, can be understood analytically from eq.(4) or its generalized version eq. (14) where c → 1 implies Ω DM Ω visible which violates the observable relation Ω DM 5Ω visible . The main conclusion to be drawn from this plot is that the parametrical region where Ω DM assumes its observable value is very large and perfectly consistent with the QCD constraints related to parameters c and κ −1 B 1/3 m a /m π .
This conclusion is another manifestation of the basic consequence of the AQN framework when the relation (1) is not sensitive to any details of the system such as m a or θ a,i , but rather represents a direct outcome of this proposal. This fundamental result is essentially incorporated into the initial building principle of the entire framework, and cannot disappeared as a result of some additional technical details and modifications.

VI. CONCLUSION AND FUTURE DEVELOPMENT
This work is the first attempt to extract the key element of the proposal, the coefficient "c" describing the disparity between nuggets and antinuggets, as defined by eq.(4) or its generalized version eq. (14), from the observational constraints in the quantitative way. Precisely this asymmetry eventually determines the dark matter density (within this framework) we observe today as a result of the charge separation mechanism replacing the conventional "baryogenesis" scenario as overviewed in section I. This work, to a large extend, is motivated by the recent progress in the axion search experiments which, hopefully, in few years should cover almost entire interesting region of the allowed parametrical space for m a as recently reviewed in [19][20][21][22].
A distinct feature of the AQN model is that the fundamental relation (1) is always satisfied in this framework irrespectively to the parameters of the system, such as the axion mass m a or misalignment angle θ a,i . Therefore, the discovery of the axion with mass m a which gives only a fraction of the observed dark matter density Ω a < Ω DM is perfectly consistent with our scenario as the remaining portion of the dark matter is generated by the nugget's contribution (6) which always accompanies the axion production within this proposal.
The corresponding numerical results are presented in section V B which explicitly demonstrate that this model is perfectly consistent with varies constraints listed in section V A. One should emphasize that the corresponding constraints have been derived from a variety of astrophysical, cosmological, satellite and ground based observations. Furthermore, there is a number of frequency bands where some excess of emission was observed, but not explained by conventional astrophysical sources. This model may explain some portion, or even entire excess of the observed radiation in these frequency bands, see short review [31]. This model is also consistent with known constraints from the axion search experiments. Finally, in section V C we argued that the corresponding results are not very sensitive to varies QCD parameters describing the dense CS phase of the AQNs.
Therefore, there is no much sense to repeat here, in conclusion, the main results of the sections V A, V B, V C. Rather, we would like to discuss the implications of these results on the future related studies.
In particular, it has been recently argued that the dark matter axions can generate peculiar effects in the socalled Josephson junctions [57][58][59]. If a small measured peak of unknown origin is identified with the dark matter axions as suggested in [57][58][59], one can infer that the axion mass should be m a = (104.0 ± 1.3) · 10 −6 eV, ρ a 0.1 GeV cm 3 , (33) where we also quote the estimated axion dark matter density ρ a which was required to interpret the measured signal in the Josephson junctions as a result of the DM axion [59]. This value of ρ a represents approximately a quarter of the expected dark matter density in the halo. One should emphasize that the estimate (33) for the axion mass was extracted from analysis of four different experiments observed by four different groups pointing to the same mass value (33). Therefore, it is very unlikely that (33) is a statistical fluke.
If one literally accepts the measured signal in the Josephson junctions as a result of the DM axion (33) one could immediately infer, for example, from Fig.5 that for m a 10 −4 eV and Ω a /Ω DM 0.3 as eq. (33) would suggest, the misalignment angle θ a,i should be order of one provided that (H I /2πf a ) 1. Furthermore, from Fig.  2b for the axion mass m a 10 −4 eV and θ a,i 1 one could estimate the average baryon charge of the AQNs which is estimated to be order of B ∼ 10 25 . This estimate is very encouraging as it is perfectly consistent with our previous phenomenological analysis reviewed in section II.
One should emphasize that such kind of estimates and self-consistency checks are very important for subsequent development and future studies. For example, in ref. [33] it has been argued that the Sun could serve as an ideal lab to study the AQN model. The main point is that the antinuggets deposit a huge amount of energy in the corona as a result of annihilation events with the solar material. It may explain the so-called "the Solar Corona Mystery" when the temperature of the corona is about 10 6 K, i.e., being a few 100 times hotter than the solar surface temperature. It may also explain the extreme UV and soft x-ray emissions from corona, which is very hard to explain using conventional astrophysical processes. It has been also suggested in ref. [33] that the observed brightening-like events called "nanoflares" in the Sun can be identified with the annihilation events of the antinuggets, in which case the observed energy distribution of the nanoflares must coincide with the baryon charge distribution studied in the present work. In other words, the energy distribution of the nanoflares and the baryon charge distribution of the nuggets is one and the same function within the AQN model [33].
This statement can be formally expressed as follows where dN is the number of the nanoflare events per unit time with energy between W and W + dW which occur as a result of complete annihilation of the antinuggets carrying the baryon charges between B and B + dB. These two distributions are tightly linked as these two entities are related to the same AQN objects when the annihilation of the antinugget's baryon charge generates the energy which is interpreted as the observed nanoflare event within AQN model. The modelling and observations of the nanoflares in corona suggest α ∼ 2, see [33] with the details and references. One should emphasize that this interpretation is consistent with constraints on B, see (3), refs. [31,32] for review and footnote 3 with some additional remarks. Therefore, further studies of the nanoflares in solar corona may shed some light on the nature of the dark matter.
In fact this interpretation may receive further support by future analysis of some specific correlations studied in [60]. If subsequent development along the lines advocated in [60] suggesting that the frequency of the flares observed in the sun is directly related to the dark matter particles (called "invisible matter" in ref. [60]), it would be a major breakthrough not only in our understanding of the solar corona, but also in our understanding of the nature of the dark matter.
The first theoretical estimates providing a specific mechanism on how the correlations observed in [60] could be, in principle, explained have been recently suggested in [34]. The basic idea of [34] is that the nuggets entering the solar corona will inevitably generate the shock waves as the typical velocities of the nuggets are much higher that the speed of sound in solar atmosphere. The shock waves due to the AQNs may serve as the triggers igniting the large flares which are shown to be correlated with "invisible matter" of ref. [60]. This mechanism would also relate naively unrelated entities such as the axion, its mass, the baryon charge distribution B of the nuggets (34), the flare's intensity and their frequency of appearance in the Sun.
Our next comment about possible future studies is related to recent activities of the axion search experiments. We want to specifically mention some present and future axion search experiments such as ADMX and ADMX-HF [61], IAXO [62], CAST [63], ORPHEUS [64], MADMAX [65], ORGAN [66]. These, and possibly many other experiments should eventually discover the QCD axion irrespectively to the assumption on validity of the axionic Josephson effect leading to (33). The only original comment we would like to make is that the discovery of the axion with mass m a which may generate only a small fraction of the observed dark matter density Ω a Ω DM nevertheless would be a major discovery because the remaining portion of the dark matter could be generated by the nugget's contribution (6) which always accompanies the conventional axion production within the AQN scenario, and always satisfies the generic relation (1). This would conclude a long and fascinating journey of searches for this unique and amazing particle conjectured almost 40 years ago.
We want to make one more comment on the axion search experiments when the observable is sensitive to the axion amplitude θ itself, in contrast with conventional proposals which are normally sensitive to the derivatives ∼ ∂ µ θ. The basic idea of the proposal [67] is that the θ becomes a physically observable parameter even in the abelian Maxwell QED if the system is defined in a topologically nontrivial sector. This can be easily achieved by placing the system into the background of an external magnetic field. The phenomenon in all respects is very similar to well known Witten's effect when the θ parameter becomes a physical observable in the presence of the magnetic monopole. This novel phenomena was coined as the Topological Casimir Effect (TCE), and there are some specific ideas how to design and fabricate the corresponding apparatus (the so-called aKWISP project), see talk by Cantatore [68] with relevant information.
One more possible direction for future studies from our "wish list" is a development of the QCD-based technique related to nuggets evolution, cooling rates, evaporation rates, viscosity, transmission/reflection coefficients, etc., in an unfriendly environment with non-vanishing T, µ, θ. The problem here is that the axion mass scale m a and the Hubble scale H(T ) at the QCD epoch are drastically different from the QCD scale itself ∼ Λ QCD describing the CS phase and the nugget's structure. In similar circumstances some researches normally change the scales (in which case the relevant parameters obviously assume some unphysical values) to attack the problem numerically. After the numerical computations are done, one can return to the physically relevant values for the parameters by using plausible arguments with some specific assumptions (which may or may not be correct) on the scaling features of the parameters when they assume their physical values.
Our goal in the nearest future is to develop some numerical methods and approaches to study the real time evolution of nuggets with real physical parameters which assume drastically different scales. If the project turns out to be successful it would be a major technical step forward relating the analysis of refs [9,10] which was devoted to the formation period of the nuggets at high temperature and present work which is dealing with the present epoch of the cold Universe, billion years after the nuggets had been formed.  This constraint does not imply that very large nuggets cannot exist. In fact the quark nuggets can perfectly coexist with conventional nuclear matter making an absolutely stable object, when the hadronic nuclear matter is surrounding the dense the quark nugget. Such objects could be captured by stars or planets and stay in the cores of the astronomical object indefinitely 11 . They can also accrete the visible hadronic matter during a long Hubble evolution since the formation time at T 41 MeV as shown on Fig.1. However, an antinugget made of antimatter will annihilate its baryon charges when it is in contact with visible matter. As the main phenomenological manifestation of the AQN model is precisely such dition (22). Large value of κ −1 B 1/3 ma/mπ > 0.95 corresponds to the small chemical potential µ 330 MeV. At this relatively small value of µ, our treatment of the dense matter as the colour superconductor is not justified. The reason is that µ 330 MeV corresponds to the conventional stable nuclear matter which should be treated differently as the relevant degrees of freedom are hadrons, rather than quarks. 11 In fact, there are many arguments suggesting that the cores of the neutron stars could be in CS phase.
annihilation events of the antinuggets with visible matter in the present work we present our plots for the antinuggets and their baryon charges B.
The absolute stability of the AQNs is determined by the condition 3 QCD mp ≤ 1, while the metastability (with very large life time exceeding the life time of the Universe) this condition may not be strictly satisfied. The corresponding plots are shown on Fig. 8. The main lesson from these studies is that the Model-2 is absolutely stable in the entire parametrical region 0.2 κ −1 B 1/3 m a /m π 0.95, while Model-1 is absolutely stable for 0.6 κ −1 B 1/3 m a /m π 0.95, and it becomes metastable for 0.2 κ −1 B 1/3 m a /m π 0.6. This metastability (in contrast with absolute stability) should not be a point of concern as discussed in the original paper [7] because it corresponds to a very long life time of the nuggets. This metastability region can be ignored for any practical purposes in our discussions in the present work.
For illustrative purposes we also present on Fig. 9 the plot for the domain wall contribution DW / tot as a function of the same dimensionless parameter κ −1 B 1/3 m a /m π . One can explicitly see from this plot that the energy related to the axion field represents a considerable portion of the nugget's total energy representing approximately 1/3 of the total energy for the Model-1 (red curve). However, this contribution is very distinct from the conventional propagating axions (represented by Ω a in this work) which are produced as a result of misalignment mechanism. The energy of the axion field represented by DW / tot cannot be easily released to space as the axions describing the axion domain wall are not on-shell axions. The corresponding energy in form of the propagating axions can be only released to the space when the nuggets get completely annihilated and the axions released, for example in the solar corona as discussed in [33,34]. It would be a major discovery if these axions can be observed by CAST [63] or IAXO [62] type instrument.