The Axion Quark Nugget Dark Matter Model: Size Distribution and Survival Pattern

We consider the formation and evolution of Axion Quark Nugget dark matter particles in the early universe. The goal of this work is to estimate the mass distribution of these objects and assess their ability to form and survive to the present day. We argue that this model allows a broad range of parameter space in which the AQN may account for the observed dark matter mass density, naturally explains a similarity between the"dark"and"visible"components, i.e. $\Omega_{\rm dark}\sim \Omega_{\rm visible}$, and also offer an explanation for a number of other long standing puzzles such as"Primordial Lithium Puzzle"and"the Solar Corona Mystery"among many other cosmological puzzles.


I. INTRODUCTION
In this paper we describe a scenario in which dark matter (DM) consists of macroscopically large, nuclear density, composite objects known as axion quark nuggets (AQNs) [1]. In this model the "nuggets" are composed of large numbers of standard model quarks bound in a nonhadronic high-density color superconducting (CS) phase. As with other high-mass dark matter candidates (such as Witten's quark nuggets [2]; see Ref. [3] for a review) these objects are "cosmologically dark" not through the weakness of their interactions, but due to their small cross section-tomass ratio which scales all observable consequences. As such, constraints on this type of dark matter place a lower bound on their mass distribution, rather than the coupling constant.
There are two additional elements in the AQN model in comparison with the older well-known construction [2,3]. First, there is an additional stabilization factor provided by the axion domain walls which are copiously produced during the QCD transition and which help to alleviate a number of the problems inherent in the older models. 1 Another crucial additional element in the proposal is that the nuggets could be made of matter as well as antimatter in this framework. This novel key element of the model [1] completely changes the 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 , irrespective of any specific details of the model, such as the axion mass or size of the nuggets. It was precisely this fundamental consequence of the model that was the main motivation for its construction.
The presence of a large amount of antimatter in the form of high-density AQNs leads to a large number of observable consequences within of this model as a result of annihilation events between antiquarks from AQNs and visible baryons. We refer to Sec. II for a short overview of the basic results, accomplishments, and constraints of this model. 1 In particular, the first-order phase transition was a required feature of the system for the original nuggets 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. Furthermore, the nuggets [2,3] will likely evaporate on the Hubble time scale even if they had been formed. In case of the AQNs, the first-order phase transition is not required as the axion domain wall plays the role of the squeezer. Furthermore, the argument related to the fast evaporation of the nuggets is not applicable for the AQNs 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 external pressure provided by the axion domain wall, in contrast with original models [2,3] which must be stable at zero external pressure.
The only comment we would like to make here is that some long-standing problems may find their natural resolutions within the AQN framework. The first of these is the "primordial lithium puzzle" which has persisted for at least two decades. It has been recently shown [4] that this longstanding mystery might be naturally resolved within the AQN scenario. Another example is the 70-year-old mystery (since 1939) known in the community as "the solar corona mystery." It has been recently suggested that this mystery may also find its natural resolution within the AQN scenario [5,6] as a result of the annihilation of AQNs in the solar corona. These two examples show the very broad application potential of this model. Furthermore, the corresponding quantitative results are highly sensitive to the size distribution of the AQNs, and their ability to survive in unfriendly environment such as the solar corona or high-temperature plasma during big bang nuclear synthesis (BBN).
The main goal of the present work is to focus on these two specific questions about the model which have been previously ignored (mostly due to oversimplified settings), with the main goal of a qualitative (order-of-magnitude) estimate rather than a quantitative description. We are now in a position to fill this gap and address the hard questions on the size distribution and survival pattern during the long evolution of the Universe.
The central result is a demonstration that AQNs of a sufficient size will survive the high-density plasma of the early Universe from their formation at the QCD transition until the present day. However, before turning to the details of formation and evolution we will briefly review the relevant properties of the AQN dark matter model, and its basic predictions, results, and accomplishments in Sec. II. In Secs. III and IV we discuss the size distribution during the formation period, while Secs. V-VIII are devoted to an analysis of the survival features of the AQNs at high (before the BBN epoch) and low (after the BBN epoch) temperatures. In Sec. IX we analyze the present-day observational constraints on the mass distribution.

II. THE AQN DARK MATTER MODEL
A. The basic predictions, results, and accomplishments The AQN dark matter model was originally introduced to resolve two important outstanding problems in cosmology: the nature of dark matter and the mechanism of baryogenesis. The connection between these seemingly unrelated questions is motivated by the apparently coincidental similarity of the visible and dark matter energy densities, If the dark matter is in fact a new fundamental particle then there is no a priori reason for this similarity, and visible and dark matter could have formed with widely different energy densities. However, if these two forms of matter share a common origin the ratio (1) may have a physical explanation rather than being a tuned parameter of the theory. In the case of visible matter the energy density is fixed at the QCD transition when baryons form and acquire their observed mass, 2 so we may ask if the dark matter density may also form at this time.
The AQN proposal represents an alternative to the baryogenesis scenario when the "baryogenesis" is replaced by 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 a CS phase. Dense nuggets in a CS phase are also present in the system such that the total baryon charge remains zero at all times during the evolution of the Universe. The detailed mechanism of the formation of the nuggets and antinuggets has been recently developed in Refs. [7][8][9]. The only comment we would like to make here is that the energy per baryon charge is approximately the same for nuggets in a CS phase and visible matter in a hadronic phase, as both types of matter are formed during the same QCD transition and both are proportional to the same fundamental dimensional parameter ∼Λ QCD . Therefore, the relation (1) is a natural outcome of the framework rather than a consequence of a fine-tuning.
In the context of the AQN model the dark matter is formed by the action of a collapsing network of axion domain walls formed at the QCD transition. These processes will contribute to direct axion production through the misalignment mechanism, domain-wall decay, as well as the nuggets' formation; see Fig. 1. This process will be described in detail in Sec. III, and we shall not elaborate on this topic here.
The result of this "charge separation" process is two populations of AQNs carrying positive and negative baryon number. That is, the AQNs may be formed of either matter or antimatter. However, due to the global CP-violating processes associated with θ 0 ≠ 0 during the early formation stage (see Fig. 1), the number of nuggets and antinuggets formed will be different. This difference is always an order one effect irrespective of the parameters of the theory, the axion mass m a , or the initial misalignment angle θ 0 , as argued in Refs. [7,8]. The disparity between nuggets Ω N and antinuggets ΩN unambiguously implies that the baryon contribution Ω B must be of the same order of magnitude as all contributions are proportional to one and the same dimensional parameter Λ QCD .
If we assume that all other dark matter components (including conventional axion production, as shown on Fig. 1) are subdominant players, one can relate the baryon charge hidden inside the nuggets with the visible baryon charge. Indeed, the observations suggest that Ω dark is 5 times greater than Ω B , which in our framework implies that 2 Prior to the QCD transition, the quarks and leptons carry masses generated via the Higgs mechanism, but these are 3 orders of magnitude smaller than the baryon masses and thus represent a negligible fraction of the present-day energy density.
This approximate relation represents a direct consequence of baryon charge conservation, If direct axion production is not negligible, the ratio (2) will obviously be modified. However, all numerical coefficients entering Eq. (2) will always be order of one, unless the axion mass m a is fine-tuned to saturate the present value for Ω dark . We refer to the original paper [9] and specifically Fig. 5 in that paper for a more precise relation between these two distinct contributions. One should note that the direct axion production contribution Ω axion to dark matter Ω dark is highly sensitive to the axion mass m a in contrast with the nuggets' contribution (1) which holds irrespective of the value of the axion mass m a . In particular, Ω axion becomes negligible for sufficiently large axion mass because it scales as Ω axion ∼ m −7=6 a ; see footnote 4 for recent numerical estimates. These estimates suggest that if m a ≥ 10 −4 eV, direct axion production is numerically small, Ω axion ≪ Ω dark , and therefore, the ratio (2) is approximately valid.
We may also reformulate Eq. (2) in terms of the number densities and average mass of the various components, where the AQNs masses are related to their baryon number by M N ≈ MN ≈ m p jBj. The resulting AQNs will be macroscopically large (typically with radii above 10 −5 cm) and of roughly nuclear density, resulting in masses above roughly a gram. The density of the color superconducting phase is not precisely known and depends on the exact phase of matter realized in the AQN. Furthermore, the axion domain wall surrounding the nugget also contributes to its mass; see Ref. [9] for quantitative relations. For the present work we will simply adopt a typical nuclear baryon number density of order 10 40 cm −3 for all our estimates, such that a nugget with jBj ∼ 10 25 has a typical radius of R ∼ 10 −5 cm. As a result, the effective interaction is very small, σ N =M N ∼ 10 −10 cm 2 =g, where σ N ∼ R 2 assumes a typical geometrical cross section. This estimate is well below the upper limit of the conventional DM constraint σ=M DM < 1 cm 2 =g. This is the main reason why, despite being made from strongly interacting particles, the AQNs nevertheless behave as cold DM from the cosmological perspective.
Another fundamental ratio is the baryon-to-entropy ratio at present time, In the AQN 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; see Fig. 1.
We refer to Ref. [7] for relevant estimates of the parameter η within the AQN framework. We note that T form ∼ Λ QCD assumes a typical QCD value, as it should. This is because there are no small parameters in QCD and all observables must be expressed in terms of a single fundamental parameter, which is Λ QCD . One should add here that the numerical smallness of the factor (5) is a result of an exponential sensitivity of η to the formation temperature as η ∼ expð−m p =T form Þ, with the proton's mass being a numerically large parameter when it is written in terms of the QCD critical temperature m p ≃ 5.5T c . It is the purpose of this work to investigate the mass distribution of AQNs generated by the collapse of the axion domain-wall network and assess their survival pattern within the high-density plasma of the early Universe. The efficiency of the formation process is reflected in the observed baryon-to-photon ratio (5), which implies that only a small fraction of the primordial baryonic matter is successfully bound when the formation is completed at T ≈ 41 MeV into AQNs and thus protected from further annihilation.
To conclude this short overview of the basic features of the AQN framework, we would like to mention that the AQN dark matter model has recently been applied to a variety of situations in the early Universe. In particular, it has been suggested that the anomalously strong 21-cm absorption feature reported by the EDGES Collaboration [10] may be driven by an additional component to the largewavelength end of the radiation spectrum at early times [11]. Such a component may be produced by thermal FIG. 1. This diagram illustrates the interrelation between axion production due to the misalignment mechanism and the nuggets' formation which starts before the axion field θ relaxes to zero. Adopted from Ref. [9]. emission from a population of AQNs which primarily emit well below the cosmic microwave background (CMB) peak and which are not expected to be in thermal equilibrium with the photons [12].
It has also been suggested that the presence of partially ionized AQNs at temperatures T ≃ 20 keV soon after the BBN epoch may result in the preferential capture (and eventual annihilation) of the highly charged heavy nuclei with Z ≥ 3 produced during BBN. This proposal offers a resolution [4] to the long-standing "primordial lithium puzzle." While both of these phenomena-which occurred at very early times in the evolution of the Universe (at redshifts z ≃ 17 and T ≃ 20 keV, respectively)-are potentially very interesting and important, the underlying physics and astrophysical backgrounds are not sufficiently well understood to impose strong constraints on the AQN size distribution during these earlier times. Instead, the strongest constraints come from current, more readily observable, and better understood environments and will be the subject of Sec. II B of this work.

B. Mass distribution constraints
As stated above, the nuggets are not fundamentally weakly interacting but are effectively "dark" due to their large mass and consequent low number density. For example, the flux of nuggets that should be observed on or near Earth is and thus direct-detection experiments can impose lower limits on the value of hBi for the distribution of AQNs. Limits may also be obtained from astrophysical and cosmological observations. In this case, any observable consequences will be scaled by the matter-AQN interaction rate along a given line of sight, where R ∼ B 1=3 is the typical size of a nugget which determines the effective cross section of interaction between DM and visible matter. Thus, as with direct detection, astrophysical constraints impose a lower bound on the value of hBi.
The relevant constraints come from a variety of both direct-detection and astrophysical observations, which we summarize here. Further details are available in the original papers and references therein.

Direct detection
As mentioned above, the flux of AQNs at the Earth's surface is scaled by a factor of B −1 and is thus suppressed for large nuggets. For this reason, the experiments most relevant to AQN detection are not the conventional highsensitivity dark matter searches, but rather detectors with the largest possible search area. For example, it has been proposed that large-scale cosmic-ray detectors such as the Auger Observatory or Telescope Array may be sensitive to the flux of AQNs in an interesting mass range; however, this sensitivity is strongly limited by the relatively low velocity of the AQNs, v N ∼ 10 −3 c [13].
The strongest direct-detection limit is likely set by the IceCube Observatory's nondetection of a nonrelativistic magnetic monopole [14]. While magnetic monopoles and AQNs interact with the material of the detector in very different ways, in both cases the interaction leads to electromagnetic and hadronic cascades along the trajectory of AQNs (or magnetic monopole) which must be observed by the detector if such events occur. The nonobservation of any such cascades would put a limit on the flux of heavy nonrelativistic particles passing through the detector, which limits the AQN flux to Φ N ≲ 1 km −2 yr −1 which is mainly determined by the size of the IceCube Observatory. Similar limits are also derivable from the Antarctic Impulsive Transient Antenna [15], though this result depends on the details of the radio-band emissivity of AQNs. There is also a constraint on the flux of heavy dark matter with mass M < 55 g based on the nondetection of etching tracks in ancient mica [16].
If we take the local dark matter mass density to be ρ ≈ 0.3 GeV=cm 3 and assume that this value is saturated by AQNs, we may translate the flux constraint Φ N ≲ 1 km −2 yr −1 into a lower limit with 3.5σ confidence on the mean baryon number of the nugget distribution of hBi > 3 × 10 24 ðdirect detection constraintÞ; ð8Þ where we assume 100% efficiency of the observation of the AQNs passing through the IceCube Observatory. If this efficiency is lower, the limit (8) will also be weakened.

Indirect detection
We also consider the constraints arising from possible dark matter interactions within the Solar System. These include a limit from the potential contributions to Earth's energy budget which require jBj > 2.6 × 10 24 [15], which is consistent with Eq. (8). It has also been suggested that there is a strong limit on any annihilating AQN population due to the low flux of high-energy neutrinos from the Dun [17]. However, the composite nature of AQNs when the bulk of quark matter is in a CS phase is likely to result in the majority of neutrino emission occurring at relatively low energies where they will be lost in the much larger conventional solar neutrino background [18].

Galactic observations
It is known that the spectrum from the Galactic center (where the dark and visible matter densities assume high values) contains several excesses of diffuse emission the origin of which is unknown, the best-known example being the strong galactic 511 keV line [19].
The emission spectrum of AQNs has been studied in a variety of environments and over a range of energy scales. As discussed above, the expected diffuse emission due to AQNs scales as the product of the dark and visible matter densities (ρ DM · ρ B ) so that the strongest consequences will be from relatively high-density regions such as the Galactic center. The dependence on the visible matter density also means that the AQN model predicts a dark matter contribution to the galactic spectrum that is less spherically symmetric than is expected for either decaying (∼ρ DM ) or self-annihilating (∼ρ 2 DM ) dark matter models, consistent with observations [19]. The emission spectrum associated with an AQN population will be distinct from that of more conventional dark matter candidates in that the most energetic annihilations will be at the nuclear (∼100 MeV) scale, implying that any signal will be limited to sub-GeV energies.
The potential AQN contribution to the Galactic spectrum has been analyzed across a broad range of frequencies, from a radio-band thermal contribution to x-ray and γ-ray photons produced by more energetic annihilation events. In each case the predicted emission was consistent with observations and had the possibility of improving the global fit of galactic emission models. All of these emissions from different frequency bands are expressed in terms of the same integral (7), and therefore the relative intensities are unambiguously and completely determined by the internal structure of the nuggets which is described by conventional nuclear physics and basic QED. For further details, see the original studies in Refs. [20][21][22][23][24][25][26] which contain specific computations in different frequency bands for galactic radiation, and the short overview in Ref. [27].
To summarize this subsection: the most significant potential dark matter signal in this context is the galactic 511 keV line (plus related continuum emission in the 100 keV range) resulting from low-energy electron-positron annihilations through 1 S 0 and 3 S 1 positronium formation with consequent decays. These emission features of the galactic spectrum have proven difficult to explain with conventional astrophysical sources. At the same time, the AQN model with the constraints from direct detection discussed in Sec. II B 1 could offer a potential explanation for the entire observed 511 keV emission feature (including the width, morphology, 3γ continuum spectrum, etc.).
If further contributions from conventional astrophysical sources are discovered the constraint (8) from Sec. II B 1 may be tightened to higher values of hBi. As the line of sight through the galactic center samples the emission from a large number of individual AQNs, this measurement is sensitive only to the average baryon number hBi of the antimatter AQNs and does not provide any information on their size distribution, in contrast with the solar observations discussed in the next subsection, which are highly sensitive to the size distribution.

Solar corona observations
Yet another AQN-related effect might be intimately linked to the so-called "solar corona heating mystery." The renowned (since 1939) puzzle is that the corona has a temperature T ≃ 10 6 K which is 100 times hotter than the surface temperature of the Sun, and conventional astrophysical sources fail to explain the extreme UV (EUV) and soft x-ray radiation from the corona 2000 km above the photosphere. Our comment here is that this puzzle might find its natural resolution within the AQN framework, as recently argued in Refs. [5,6,28].
In this scenario the AQNs composed of antiquarks fully annihilate within the so-called transition region (TR), providing a total annihilation energy of order where b ⊙ is the Sun's gravitational capture parameter, The estimate (9) is determined by the local dark matter density, independent of the mass distribution. This value is suggestively close to the observed EUV luminosity of 10 27 erg=s. The EUVemission is believed to be powered by impulsive heating events known as nanoflares, the origin of which is unknown. If these nanoflares are in fact AQN annihilation events [which was precisely the main conjecture formulated in Ref. [5] and formally expressed by Eq. (11); see below] we may extract an upper limit on the mass distribution for the AQNs because the energy distribution of the nanoflares has been previously studied in the context of solar physics. The main reason for our ability to study the AQN mass distribution is due to the fact that the nanoflare distribution has been modeled using magnetohydrodynamics (MHD) simulations in plasma and solar physics studies to match the solar observations with simulations. We can use the corresponding results to constrain the AQN mass distribution.
A few comments on nanoflares and the AQN distribution are in order. First of all, the majority of nanoflares (and, therefore, individual AQN annihilation events) must be below the resolution of solar telescopes and they must interact sufficiently with the corona to deposit the majority of the available energy in the TR rather than at lower radii. According to Ref. [29], the resolution limit for flares is E res ≈ 3 × 10 24 erg ≈ 2 × 10 27 m p c 2 which implies that the majority of AQNs must have a baryon number below B ∼ 10 27 or, alternatively, that only a small fraction of their mass is able to annihilate in the corona, though the later possibility is disfavored by the analysis of Ref. [6].
Second, various analyses of coronal heating based on MHD have considered the nanoflares to be a lowenergy continuation of the higher-energy class of solar flares, despite the fact that they have significantly different spacial and temporal distributions. 3 Under this assumption, a number of energy distributions have been considered generally with power-law fits consistent with the better observed population of flares at larger energies. The analysis of Ref. [30] favors a power law with slope α ≈ 2.5, where α is defined as follows: where dN is the number of nanoflare events per unit time with energy between E and E þ dE. According to the conjecture formulated in Ref. [5], this distribution coincides with the baryon charge distribution dN=dB which is the topic of the present work. These two distributions are tightly linked as these two entities are related to the same AQN objects according to our interpretation of the observed nanoflare events. Third, any population with a slope shallower than α ¼ 2 will experience too few nanoflares to dominate the total heating budget. However, an alternate analysis [31] considered a broken power law in which a shallow slope (α ≈ 1.2) transitions to a steeper slope (α ≃ 2.5) at large energies. In this case, the heating contribution will be peaked at the energy where the break in the spectrum occurs, i.e., the position of the knee. In the model off Ref. [31] it occurs at E ≃ 10 24 erg, which is slightly below the instrumental resolution energy E res ≈ 3 × 10 24 erg. In many respects we consider this model to be preferable from the AQN perspective because it explicitly shows that nanoflares and flares have different natures, in agreement with indirect evidence pointing to their distinct origins; see footnote 3.
While the mechanism of energy release in the AQN model is substantially different from conventional flare studies, the nanoflare models of Refs. [30,31] provide a useful parametrization of the distribution of AQN masses which may be consistent with the observed degree of coronal heating and EUV emission.
With this set of observational constraints in mind, we now turn to the main purpose of this work: a study of the formation and subsequent evolution of the AQN population.

III. FORMATION OF THE AQNS
This section should be viewed as an introduction to the domain-wall formation mechanism, its basic ideas (such as percolation and the formation of closed surfaces), basic generic results, and main assumptions. We also give an overview of some results from our previous studies [7,8], which represent the starting point of the quantitative approach which is the subject of the present work. We also report some new numerical results (supporting the entire framework) at the very end of this section.
It is known that axion domain walls can form in the early Universe [32,33]. When the Universe cools down to T osc ∼ 1 GeV, the axion mass effectively turns on, the axion potential gets tilted, and the axion field starts to oscillate; see During this time when T c < T < T osc , the axions are emitted and may contribute to the dark matter density. The conventional mechanism of emission is the misalignment mechanism [34][35][36]. The axions may also be radiated due to the decay of topological defects [37][38][39][40][41][42]. In both cases the corresponding contribution is highly sensitive to the axion mass m a as the corresponding contribution to the dark matter scales as Ω axion ∼ m −7=6 a . There are a number of uncertainties and remaining discrepancies in the corresponding estimates. We do not comment on these subtleties here and refer to the original papers. 4 The formation of AQNs always accompanies these two distinct contributions to Ω DM . However, in comparison with the misalignment mechanism [34][35][36] and the decay of topological defects [37][38][39][40][41], the contribution of the nuggets to dark matter is always an order-one effect and is not sensitive to the axion mass m a or the misalignment angle θ 0 , as overviewed in the Introduction and expressed by Eq. (1). The axion field plays a dual role in this framework: it is responsible for the direct production of the propagating axions with a contribution which scales as 3 Nanoflares appear to occur uniformly across the solar surface, while larger flares are strongly correlated with active regions. Furthermore, the EUV intensity (which represents the nanoflare activity) shows very modest variation within the solar cycle. It is in drastic contrast with large flares which demonstrate a huge variation (a factor of 10 2 or more) in frequency of appearance within the solar cycle; see the detailed discussions in Ref. [6]. 4 According to the most recent computations presented in Ref. [41], the axion contribution to Ω DM as a result of the decay of topological objects can saturate the observed DM density today if the axion mass is in the range m a ¼ ð2.62AE0.34Þ10 −5 eV, while earlier estimates suggest that the saturation occurs at a larger axion mass. This could be due to some other complications with conventional computations, as argued in Ref. [42]. One should also emphasize that the computations in Refs. [37][38][39][40][41][42] were performed with the assumption that Pecci-Quinn symmetry was broken after inflation. Ω axion ∼ m −7=6 a , and it plays a key role in the formation of AQNs, as discussed in Refs. [7,8].
In the AQN model, we assume the preinflation scenario in which the Pecci-Quinn phase transition occurs before inflation [7]. Normally, in this case no topological defects can form as there is a single vacuum state which occupies the entire observable Universe; see footnote 4 for clarification. This argument is absolutely correct for N DW ≠ 1 axion domain walls which require the presence of different physical vacua with the same energy. However, N DW ¼ 1 axion domain walls are special in the sense that the axion field θ interpolates between one and the same physical vacuum but corresponding to different topological k branches with θ → θ þ 2πk. As explained in Ref. [7], different k branches of the same vacuum must be present at each point in space to provide the 2π periodicity of the vacuum energy [43,44]. Inflation cannot separate these k branches. As a consequence, N DW ¼ 1 axion domain walls can form even in this preinflation scenario with θ interpolating between the k ¼ 0 branch (θ ¼ 0) and k ¼ 1 branch (θ ¼ 2π); see Ref. [7] for details.
The key point is that a finite portion (a few percent) of N DW ¼ 1 walls are formed as closed surfaces. Such behavior has been observed in numerous numerical simulations [37,45]; see also Sec. IV for comments and more details. In previous studies this contribution to Ω axion (due to the closed axion domain walls) has been ignored because these closed surfaces (representing only a few percent of the total area) collapse as a result of the wall tension and do not play any significant role in the dynamics of the system. However, in the AQN framework this small but finite portion of the closed surfaces plays a key role. This is because the collapse of the closed N DW ¼ 1 bubbles will be halted due to the Fermi pressure created by the accumulated fermions [7]. As a result, the closed N DW ¼ 1 bubbles will eventually become stable nuggets and serve as dark matter candidates.
The dynamics of nuggets with size RðtÞ is governed by the following equation [7]: where σ eff ¼ κ · 8f 2 a m a ðtÞ describes the effective domainwall tension (which does not coincide with the well-known expression σ ¼ 8f 2 a m a computed in the thin-wall approximation), and ΔP is the pressure difference inside and outside the nugget.
An important element that we want to discuss here (in addition to our previous studies) is related to the viscosity term ∼η which enters Eq. (12) and effectively describes the friction for the domain-wall bubble oscillating in hightemperature plasma. It is precisely this term that describes the slow change of the nuggets' size before the formation is completed and the nuggets assume their final form at Fig. 1. For small oscillations the solution of Eq. (12) can be approximated as which shows the physical meaning of the frequency ω ∼ R −1 form and damping time τ. The parameter τ describes the time at which the formation is completed. This is a highly nontrivial parameter as it represents a combination of very different scales. Indeed, the viscosity η along any path shown in Fig. 1 is always assumes Λ QCD scale (of course it is, not known exactly in different phases); the axion scale appears as it enters through σ eff and finally, the cosmological scale enters as the formation effectively starts at T ≃ T c ≃ 170 MeV and must end at T ≃ T form which represents a very long cosmological journey with typical time scale t ∼ T −2 ∼ 10 −4 seconds.
It is a highly nontrivial observation that all of these drastically different scales nevertheless lead to a consistent picture. Indeed, a typical time for a single oscillation is ω −1 ∼ 10 −14 s for an axion mass m a ∼ 10 −4 eV, while the number of oscillations is very large and of order ωτ ∼ 10 10 according to Eq. (13); see also Appendix A. Therefore, the complete formation of the nuggets occurs on a time scale of 10 −4 s, which is precisely the cosmological scale when the temperature drops to 41 MeV. This scale is known from completely different arguments related to the estimate of the baryon-to-photon ratio (5).
Unfortunately, we could not numerically test this amazing "conspiracy of scales" in our original studies [7,8]. This is because the factor ωτ is very large in comparison with the other scales of the problem. It is very hard to deal with very large (or very small) factors in numerical computations. 5 This is precisely the reason why in the numerical analyses in Refs. [7,8] the viscosity term ∼η was artificially enlarged ∼10 8 times to make Eq. (12) numerically solvable, which is a conventional technical trick; see footnote 5.
It is one of the goals of the present work to overcome this technical difficulty by adopting a new numerical methodcoined as the envelope-following method-which can solve our system successfully while allowing the viscosity term to keep its real physical magnitude η ∼ Λ 3 QCD when the parameter ωτ ∼ 10 10 assumes its very large physical values. We describe the method and present the numerical analysis in Appendix A. Here we only summarize the basic results 5 Of course, our case is by no means special in this respect: it is a common problem in any numerical study when some parameters assume parametrically large/small values. This is obviously the case in any numerical studies related to axion physics because of the drastic separation of scales; see, e.g., Refs. [40][41][42].
of these studies which confirm the main features of the AQN model; see Fig. 2: (1) The nugget completes its evolution by oscillating a large number of times ωτ ∼ 10 10 before it assumes its final configuration with size R form at T form ≈ 40 MeV. Therefore, the "conspiracy of scales" phenomenon mentioned above has been explicitly tested. (2) The chemical potential inside the nugget indeed assumes a sufficiently large value μ form ≳ 450 MeV during this long evolution. This magnitude is consistent with the formation of a CS phase. Therefore, the original assumption about the CS phase used in the construction of the nugget is justified a posteriori.

IV. BARYON CHARGE DISTRIBUTION
The main goal of this section is to calculate the baryon charge distribution of nuggets in the AQN scenario and compare it to the observational constraints listed in Sec. II B. We start in Sec. IVA by describing the basic idea of the computations. In Secs. IV B and IV C we study initial size and temperature distributions, respectively. Finally, in Sec. IV D we present our main results on the nuggets' distribution dN=dB.

A. Basic idea of the computations
In the present section we need a relation between the initial size of the nugget R 0 formed at the temperature T 0 and its total baryon charge when the stage of formation is completed. The desired relation reads [see Appendix A for a detailed analysis regarding Eq. (14)]. This relation implies that the total baryon charge B of a stable nugget is completely determined by the initial size R 0 and the initial temperature T 0 of the closed domain wall such that B ∝ ðR 0 T 0 Þ 3 . Equation (14) tells us that closed domain walls with different initial radii and temperatures will eventually carry different baryonic charges B. Since the closed domain walls can form with different initial radii and at different temperatures (T c ≲ T 0 ≲ T osc ), we may map these initial conditions onto a baryon charge distribution of the nuggets dN=dB.
According to Eq. (14), the baryon charge distribution (dN=dB) of stable nuggets can be obtained from the initial size (R 0 ) and initial temperature (T 0 ) distributions of the closed axion domain walls which form between T osc and T c (the initial stage of the nuggets). We start with the following equation: where dN is the number of closed domain walls with initial radii in the range ðR 0 ; R 0 þ dR 0 Þ and initial temperatures in the range ðT 0 ; T 0 þ dT 0 Þ; fðR 0 ; T 0 Þ is a two-parameter distribution function which represents the probability density of a closed domain wall with R 0 and T 0 in the above ranges. The factor N 0 is the total number of closed bubbles that form in the early Universe when T a ≲ T 0 ≲ T c , while P is a normalization factor that normalizes the probability density fðR 0 ; T 0 Þ to one, i.e., The main goal of this section is to develop a technique which allows to compute fðR 0 ; T 0 Þ.
To simplify our analysis we assume that all initial closed domain walls will eventually become stable nuggets. We clarify this assumption later in the text when we compare the prediction of our construction with observational constraints. As the next step, we use Eqs. (14) and (15) to represent the number of stable nuggets with baryon charge less than B as follows: where K · R 3 0 T 3 0 ≤ B constrains the parameter space of the integration. From Eq. (17), we can further calculate the baryon charge distribution dNðBÞ=dB which is the main topic of this section. Obviously, the distribution fðR 0 ; T 0 Þ which depends on T 0 and R 0 in a very nontrivial way plays a crucial rule in our calculations of the dNðBÞ=dB distribution. The study of the function fðR 0 ; T 0 Þ can be approximately separated into two distinct pieces: one part describes the R 0 dependence, while the T 0 distribution can be incorporated separately. The next two subsections are devoted to an analysis of these two different elements of the main problem.

B. Initial size distribution
As discussed in Sec. III, in the AQN model the N DW ¼ 1 domain walls are topological defects with the axion field θ interpolating between the k ¼ 0 (θ ¼ 0) and k ¼ 1 (θ ¼ 2π) branches. Although the k ¼ 0 and k ¼ 1 branches correspond to the same unique physical vacuum, they effectively act as two different vacua with the same energy. The domain walls can interpolate between these (physically identical but topologically distinct) vacua, similar to a model with a VðθÞ ∼ cos θ potential, when θ ¼ 0 and θ ¼ 2π correspond to one and the same physical vacuum. Therefore, the N DW ¼ 1 axion domain walls in this scenario can be treated as Z 2 domain walls which greatly simplifies the computations.
The closed Z 2 domain walls have been observed in simulations of Z 2 -wall systems [45]. In our case, this means that closed N DW ¼ 1 axion domain walls can form, which are the sources of stable nuggets as we discussed in Sec. III. Furthermore, this analogy will provide us with more useful information about the initial size distribution of these closed bubbles. The authors of Ref. [45] pointed out that the probability of forming a closed Z 2 domain wall with initial radius R 0 ≫ ξ (where ξ is the correlation length of the topological defects) is exponentially suppressed, ∼ exp ð−R 2 0 =ξ 2 Þ. The procedure in Ref. [45] to derive this relation is briefly reiterated below.
To simulate the Z 2 system in three dimensions, we first divide a big cubic volume into many small cubic cells, each of which has length ξ. Then, to each cell a number a number þ1 or −1 is assigned at random with equal probability p ¼ 0.5. This is the simulation of the phenomenon that different patches (with volume ∼ξ 3 ) of the space during the phase transition will settle randomly with equal probability in one of the two vacua (θ ¼ 0 and θ ¼ 2π in the case of N DW ¼ 1 axion domain walls). The domain walls lie on the boundaries between cells of opposite sign. Two neighboring cells are connected if they have the same sign. Many connected cells can form a cluster with the same sign. The size s of a cluster is defined as the number of cells in the cluster. We then can look for the size distribution of þ1 clusters. (Of course, the size of −1 clusters will follow the same distribution.) It turns out that this is a typical problem of percolation theory, which deals with the statistics of clusters at different values of p. See Refs. [46,47] for a review of percolation theory. 6 In our case, where p ¼ 0.5 in three dimensions, the size distribution of the finite clusters is known from percolation theory [46], where n s is the number density of finite clusters as a function of the cluster size s (the number of cells inside a cluster). Although the distribution (18) is derived for large clusters s ≫ 1 [46], it turns out that this relation can be extrapolated down to s ¼ 1 as a very good approximation [48]. As a consequence, we adopt Eq. (18) for the whole spectrum s ≥ 1 for further calculations. The two coefficients τ and λ are p dependent. According to the Ref. [48], λ has a typical value ∼10 and τ ranges from 1.5 to 2.2 based on the three-dimensional lattice simulations. Discussing the exact values of τ and λ at p ¼ 0.5 is beyond the scope of this work. Instead, we simply adopt λ ¼ 10 and τ ¼ 2 for further calculations. 7 However, as we will see, the shape of the baryon charge distribution dNðBÞ=dB of nuggets is not sensitive to the precise numerical values of τ and λ.
The result (18) can be translated into the language of domain walls straightforwardly. The probability of forming a closed bubble with radius R 0 decreases exponentially when R 0 increases, which can be formally expressed as To derive this distribution as a function of R 0 from Eq. (18), we use the relations s ≃ R 3 0 =ξ 3 and n s ¼ 1 V dN ds where we get rid of the simulation volume V (a constant) in Eq. (19). The parameter ξ is the correlation length of topological defects as mentioned above, which is also set as the length of a 6 In percolation theory, there is a percolation threshold p c at which an infinite cluster first appears on an infinite lattice. p c ¼ 0.31 in three dimensions for a cubic lattice. In our case where the probability of a cell picking þ1 is p ¼ 0.5, we have one infinite þ1 cluster (p > p c ) and one infinite −1 cluster (1 − p > p c ). In the language of domain walls, this can be interpreted as the system being dominated by one infinite wall of very complicated topology [45]. In addition to this infinite domain wall, there are some closed domain walls (finite clusters) and they satisfy the size distribution (18). The structure and dynamics of the infinite domain wall are less important for our present work which focuses on the closed domain walls. 7 λ can also be calculated using the relation λ −1 ≃ jp − p c j −1=σ p , where λ −1 is the crossover size (see, e.g., Refs. [46,49,50]). This relation is valid for jp − p c j ≪ 1. The parameter σ p ¼ 0.45 in three dimensions [47]. We then get λ ≈ 0.025 when jp − p c j ≪ 1 is satisfied. In addition, τ ¼ −1=9 for p > p c is obtained in a field-theoretical formulation of the percolation problem [51]. However, the exact values of λ and τ are not important for us, since they do not affect the slope of the distribution dNðBÞ=dB, as we will see in Sec. IV D. single cell. The smallest cluster is a cell (s ≥ 1), implying that the lowest bound of the radius of closed bubbles is R 0 ≳ ξ. Since the relation (18) is applicable for all finite clusters s ≥ 1 as mentioned above, we adopt Eq. (19) as the size distribution of all closed bubbles R 0 ≳ ξ.
It is very instructive to consider an oversimplified case where there is no initial temperature distribution. This can be realized if all of the closed bubbles form at the same moment (at the same temperature). In this case the distribution fðR 0 ; T 0 Þ does not depend on T 0 and, according to Eq. (19), can be written as fðR 0 Þ ¼ dN=dR 0 ∝ ξ −1 ðR 0 =ξÞ 2−3τ exp ½−λðR 0 =ξÞ 2 . Using Eq. (14), this R 0 dependence can be translated into a dN dB distribution: where B min ≡ K · ξ 3 T 3 0 . In this oversimplified model where there is no T 0 distribution, we find that dN=dB is greatly suppressed by the exponential factor ∼ exp ½−λðB=B min Þ 2=3 . This essentially would imply that the distribution is strongly peaked at B ≈ B min , while larger bubbles are strongly suppressed.
As we discuss in the next subsection, the T 0 dependence drastically and qualitatively changes this simplified picture.
The key element is that the closed bubbles initially form at different temperatures between T osc and T c , as discussed above. The correlation length ξ ∼ m −1 a , which is inversely proportional to the axion mass m a , drastically changes during this evolution because of the dramatic changes of the axion mass in this interval.
These profound changes completely modify the basic features of the distribution function fðR 0 ; T 0 Þ, which is the subject of the following subsection. As we shall see below, the baryon charge distribution satisfies a power law dN=dB ∝ B −α when T 0 dependence is properly incorporated, rather than following the exponential behavior (20). This power law is consistent with the parametrization (11) which has been postulated to fit the observations. Furthermore, the power-law behavior dN=dB ∝ B −α (as we discuss below) is not very sensitive to the parameters of the coefficients τ and λ, and therefore represents a very robust consequence of the framework.
C. Initial temperature distribution and the correlation length ξðTÞ As we discussed in Sec. III, the closed axion domain walls could form anywhere between T osc and T c ; see Fig. 1 for the phase diagram corresponding to this evolution. It is hard to calculate the exact T 0 distribution. It is known, though, that normally the temperature dependence enters implicitly through the correlation length ξðTÞ which is highly sensitive to the temperature.
To account for the corresponding modifications we adopt the conventional assumption that the correlation length is a few times the domain-wall width, ξðTÞ ∼ m −1 a ðTÞ. The axion mass is known to be a temperature-dependent function before it reaches its asymptotic value near T c because it is proportional to the topological susceptibility. At sufficiently high temperature T ≫ T c one can use the instanton liquid model [52,53] to estimate the power law m a ðTÞ ∝ T −β . When the temperature is close to T ≃ T c one should use the lattice results to account for the proper temperature scaling of the axion mass.
The recent lattice QCD result shows 8 that m a ðTÞ ∝ T −β with β ¼ 3.925 just above T c [54]. We then can approximate the correlation length in the entire interval as where ξ min ≡ ξðT 0 ¼ T c Þ is the minimal correlation length.
The same ξ min also serves as the minimal radius that closed bubbles could have because R ≳ ξ.
In what follows we also assume the following simple model to account for the temperature variation of the dN=dT 0 distribution 9 : where δ is a free parameter that can be adjusted to shape different T 0 distributions. This parametrization has the advantage of producing a simple final expression for the baryon number distribution while still capturing the essentials of the temperature dependence. The constant 1=T c has no special physical meaning but is introduced to balance the units of the right-hand side and the left-hand side of the relation. Perhaps the simplest case is δ ¼ 0, in which case T 0 is uniformly distributed, i.e., the probability of forming nuggets is uniform between T osc and T c . One should emphasize that δ ¼ 0 case is still not reduced to the oversimplified example mentioned at the end of the previous subsection. This is because the temperature dependence explicitly enters through Eq. (22), but it also enters implicitly through the temperature dependence of the correlation length ξðTÞ in Eq. (19). 8 Reference [54] did not show the value of β explicitly, but it provided the related data in its Supplement Information. We get β ¼ 3.925 by fitting the data provided. The lattice results are, in fact, consistent with analytical models [52,53]; see also Appendix A for additional details. 9 One subtlety is that the effect of the expansion of the Universe between T osc and T c is also included in the model (22), since N is defined as the number of closed domain walls rather than the number density.
For positive δ > 0, the nuggets tend to form close to the point T osc , while for negative δ < 0, nuggets tend to form when the tilt becomes much more pronounced close to the QCD transition temperature T c . A sufficiently large numerical value of jδj > 1 with any sign corresponds to a very sharp (almost explosive for jδj ≫ 1) increase of the probability for axion bubble formation at T ≃ T osc or at T ≃ T c depending on the sign of δ. At the same time, jδj ∼ 0 corresponds to a very smooth behavior in the entire temperature interval (22). We, of course, do not know any properties of the distribution (22) in strongly coupled QCD when θ ≠ 0. Therefore, we proceed with our computations with arbitrary δ and make comments on the obtained properties of the baryon distribution dN=dB as a function of the unknown parameter δ in next Sec. IV D.
Combining the T 0 distribution (22) with the R 0 distribution (19), and substituting Eq. (21) into Eq. (19), we arrive at the following two-parameter distribution function: Notice that here we use "¼" rather than "∝". This is because we have an extra factor P in Eq. (15) which serves as the normalization factor, and the constant multipliers in fðR 0 ; T 0 Þ can be collected and included in P.
With this expression for fðR 0 ; T 0 Þ and the basic Eq. (17), we can now proceed with our calculation of the baryon charge distribution dN=dB. The corresponding results will be discussed in the next subsection.

D. The dN=dB distribution: results
Substituting Eq. (23) into Eq. (17), one can explicitly compute the function NðBÞ and the distribution dN=dB. In what follows it is convenient to introduce the following dimensionless variables: the baryon charge b ¼ B=B min of the nugget measured from its minimum value B min ¼ Kξ 3 min T 3 c ; the relative size r ¼ R 0 =ξ min of the nugget measured from its minimum size ξ min ; and the relative temperature u ¼ T 0 =T c during formation evaluation in units of T c .
In terms of these dimensionless variables the desired distribution dN=dB can be represented as follows: (see Appendix B for technical details).
One can easily estimate the integral (24) by observing that it is saturated for very large b ≫ 1 by u sat of order when the exponential factor in Eq. (24) assumes a value of order one. Substituting the expression back into Eq. (24), one arrives at the following asymptotical behavior for the distribution: where the final result is expressed in terms of the physical baryon charge B rather than in terms of the dimensionless parameter b. The parameter α here is defined in precisely the same way as it was defined in the observational fitting formula (11). The exponent α entering Eq. (26) can be approximated in the limit B ≫ B min as follows: where in the last step we ignored the factors of order one in comparison with the known (and very large) value of β ≃ 4 to simplify the qualitative discussions below. The approximate analytical formula (26) at very large B ≫ B min is in perfect agreement with the numerical analysis presented in Appendix B.
The behavior (26) is an amazingly simple and profoundly important result. Indeed, it shows that the exponential suppression is replaced by the algebraic decay (26) which is consistent with the observational fitting formula (11). The "technical" explanation of why this happens is that the integral (24) is saturated by u sat when the exponential factor in Eq. (24) assumes a value of order one. In terms of the physical parameters, it is related to the fact that the exponential suppression (23) due to the large size R 0 is effectively removed by a strong temperature dependence with a very large beta function β. Integration over the entire temperature interval eventually leads to the algebraic decay (26).
Another important property of Eq. (26) is that the final result for the slope (27) is not very sensitive to the parameters λ and τ. The total normalization factor of course is very sensitive to these parameters, as discussed in Appendix B. It is also not very sensitive to the well-known parameter β ≈ 4 as long as it is relatively large. The slope α is mostly determined by δ which may have any sign and effectively describes the temperature interval where the bubbles are produced with the highest efficiency. The fitting models (11) (based on observations that were discussed in Sec. II) can be reproduced with a negative δ < 0. As we previously mentioned, a negative δ corresponds to a preference for bubble formation close to T c where the axion potential tilt becomes much more pronounced. Furthermore, a model with α ≃ 2 corresponds to δ ≃ −3 (strongly peaked at T ≃ T c ), while another model with α ≃ 1.2 corresponds to a more smooth distribution of dN=dB over the entire temperature interval, with δ ≃ −1 corresponding to a mild preference for bubble formation at T ≃ T c .
The last comment we want to make is about the largest possible size of nuggets. According to percolation theory, there is no upper limit on the size of finite clusters (closed domain walls). However, the shape of large clusters may not be perfectly spherical (in three dimensions) while our computations are based on the assumption of exact spherical symmetry of the formed bubbles. Furthermore, the radius for nonsymmetric bubbles is defined in an average sense for large closed clusters; see, e.g., Ref. [47] for more details. The deviation from the ideal spherical shape makes the large collapsing closed domain walls fragment with high probability into smaller pieces, and thus could significantly suppress the possibility of forming large nuggets. 10 The detailed calculations of the suppression effect from the irregular shape for large clusters is hard to carry out and also well beyond the scope of the present work. However, we may introduce a cutoff B cut to roughly account for this extra suppression. Above B cut , no nuggets can form from the collapse of closed axion domain walls. This parameter turns out to be useful when we later calculate the total number of nuggets.
We conclude this section with the following remark. The main result of our analysis is expressed as Eq. (26) with the slope (27). This formula represents the baryon charge distribution immediately after the formation period is complete when the baryon-to-photon ratio η assumes its present value (5). This "primordial" distribution of the nuggets is the subject of a long evolution in hot plasma which may modify the properties of dN=dB. This problem of the "survival" of the primordial nuggets is the subject of the next section.

V. SURVIVAL OF THE PRIMORDIAL DISTRIBUTION
After the AQNs have formed at T ≈ 40 MeV, the process of "charge separation" is essentially complete and the plasma surrounding the nuggets contains exclusively protons, neutrons, electrons, and positrons. A nugget composed of matter will gradually collect electrons into its electrosphere as the plasma cools, but apart from this it will essentially remain in its initial form. The surface layer of electrons contributes negligibly to the total mass so that the distribution of nugget masses remains essentially identical to the primordial distribution discussed above. However, the AQNs composed of antimatter, which are present in larger numbers, will be subject to a much more complicated evolution. The details of this process will be laid out below, but we first give a general overview of the evolution of the antimatter AQN mass distribution.
Initially, the plasma surrounding the AQNs is dominated by electrons and positrons which are roughly as abundant as the photons, i.e., n e ≃ n e þ ≃ n γ ∼ T 3 . During this phase the electrosphere captures positrons in those states for which the binding energy is above the plasma temperature and expands similarly to the case of a matter nugget. However, once the temperature drops below the electron mass T ≤ m e the electrosphere can no longer capture free positrons at a rate sufficient to compensate for annihilations. Below this temperature the electrosphere will begin to capture free protons which, if they stay bound to the nugget for a sufficient period of time, will eventually annihilate with the central quark matter.
The process of capturing protons becomes much more pronounced after the temperature drops to T ≈ 20 keV when the dominant portion of the positrons in the plasma get annihilated, while the number densities of electrons and protons become equal, i.e., n e ≈ n p ∼ ηT 3 . However, even at this temperature (as we discuss below) only a very tiny portion of the AQNs' baryon charge will be annihilated, such that the mass distribution still remains essentially unaffected by the unfriendly environment of the hot plasma.
Finally, after recombination at T ≤ 1 eV the surrounding matter is largely neutral and at much lower densities. During this time, matter (primarily in the form of neutral hydrogen) will continue to collide with the antimatter AQNs with some probability of annihilation but at a relatively low rate. The rare events of annihilation during the present time lead to a number of observable effects, as reviewed in Sec. II B.
At each phase of evolution the scattering rate of baryons on the nugget (and thus the probability of an annihilation) scales with the cross section of the nugget. This is at least approximately true even in the case where long-range electrical effects must be considered as the nuggets' electrical charge is itself a surface effect. As such, any change in the mass distribution should be expected to show a ΔM=M ∼ σ=M ∼ B −1=3 behavior.
The following sections will trace the evolution of the mass distribution from formation to the present day. Specifically, in Sec. VI we study the evolution of the nuggets in very hot plasma before the BBN epoch. In Sec. VII we analyze the AQN evolution before recombination, while in Sec. VIII we study the evolution of the nuggets after recombination including the period of galaxy formation. Finally, in Sec. IX we study the evolution of the AQNs in the present-day Universe. We will demonstrate that, for a range of physically interesting parameters, the initial population of AQNs will survive until the present day as a population consistent with all observational constraints and with the parameter space allowed for the axion mass and the AQNs' baryon charge B as discussed in Sec. II B. 10 Ref. [55] presented a similar argument when the author discussed the possibility of domain-wall membranes (e.g., closed domain walls) collapsing into black holes.

VI. PRE-BBN EVOLUTION
The AQN form and settle into a stable color superconducting phase at a temperature of approximately T form ≈ 40 MeV; see Fig. 1. Once this transition is complete, the AQNs will cease accreting mass and annihilation with the free baryons in the plasma will become the dominant process. 11 However, annihilation between an energetic free baryon and the quark content of the AQNs is a highly nontrivial process, as we discuss below.
We start with the estimates of the collision rate in the pre-BBN epoch. The corresponding rate between an AQN and the baryons of the surrounding plasma is where the baryon number density in plasma n B can be approximated as n B ∼ ηT 3 . The total number of collisions during this time period is saturated by the highest temperature T form ≃ 40 MeV and can be estimated as follows: where we have used the relation t ∼ T −2 to change to a temperature integration.
While the number of collisions (29) is comparable to the total baryon charge B of a nugget, the probability of annihilation is quiet small. Instead, the most likely interaction of any incident matter with the nugget is total reflection due to a number of reasons: the sharp boundary between the hadronic and CS phases such that only a very small fraction κðTÞ ≪ 1 of collisions represented by Eq. (29) will result in an annihilation. We refer to Appendix C for order-of-magnitude estimates supporting the main claim that κðTÞ ≪ 1. A more precise value is not essential for the arguments that follow.
So long as the electrons and positrons remain relativistic (and thus are present in numbers comparable to the photons) all long-range interactions are effectively screened and the cross section appearing in Eq. (28) is purely the physical size of the AQNs. As such, the estimate (29) holds until much lower temperatures when the positrons have fully annihilated (which approximately occurs at T ≈ 20 keV) and longer-range interactions become possible. The estimate (29) implies that the number of annihilation events does not modify the primordial spectrum of the AQNs discussed in Sec. IV D because the relative number of annihilation events is very small, i.e., ðκN col Þ=B ∼ κ ≪ 1.
While the baryon charge annihilation events are strongly suppressed by the factor κ ≪ 1 the e þ e − annihilation events involving particles from the AQNs' electrosphere are much more numerous and unsuppressed. One may therefore wonder if the energy injected by these annihilation events may impact the conventional thermal history of the Universe. The answer is "no," as simple estimates for the extra injected energy (due to the annihilation events with AQNs) show. Indeed, the relative injection energy due to AQNs at temperature T in comparison with the average thermal energy ðTn e Þ of the plasma can be estimated as follows: (see the Appendix in Ref. [12]). The basic reason for this tiny rate is the same as discussed before: the cross section is proportional to R 2 ∼ B 2=3 , while the number density of the nuggets is proportional to η=hBi which results in a strong suppression rate [Eq. (30)]. It is clear that such a small amount of energy injected into the system will be quickly equilibrated within the system such that the standard pre-BBN cosmology remains intact. In other words, the conventional equation of state and conventional evolution of the system is unaffected by the presence of AQNs.

VII. POST-BBN EVOLUTION
At temperatures below T γ ≃ m e the electrons and positrons begin to annihilate, causing their density to fall as ∼e −m e =T until a major portion of the positrons get completely annihilated, while the number densities for electrons (n e ) and protons (n B ) become approximately equal at T Ã ≃ 20 keV. This is the consequence of the same "charge separation" effect (replacing "baryogenesis" in the AQN framework) when more antimatter than matter is hidden in the form of dense nuggets, as reviewed in Sec. II.
This regime when the AQNs are present in the plasma at T Ã ≃ 20 keV has been recently discussed in Ref. [4] in quite a different context and for very different purposes. To be more specific, it has been shown that the primordial abundance of Li and Be nuclei will be depleted in comparison with conventional BBN computations. 12 This effect represents the resolution of the "primordial Li puzzle" within the AQN framework. 11 The plasma already possesses the required baryon asymmetry at this time so only the antimatter AQNs will be subject to annihilation, while the AQNs made of matter experience only elastic scattering. 12 The effect is due to the exponentially strong enhancement of the capture probability (and subsequent annihilation of Li and Be ions) by antinuggets. Technically, the effect for heavy ions with large Z ≥ 3 occurs due to the very strong enhancement factor ∼ exp Z; see Ref. [4] for details.
The main goal of the present work is very different, though the plasma regime surrounding the AQNs is the same with T ≲ T Ã . In the present paper we study the survival pattern of the nuggets themselves, in contrast with the studies in Ref. [4] where the main question was the analysis of the relative densities δn Z =n Z of primordial nuclei with charge Z as a result of the presence of AQNs in plasma.
We start our analysis by highlighting the basic features of the AQN electrosphere in the regime T ≲ T Ã using simple qualitative arguments. Later in the text we will support these arguments by providing some analytical formulas. At T ≈ T Ã when the external positron density essentially vanishes the boundary conditions for the AQNs' electrosphere fundamentally change, resulting in a new charge distribution. 13 Below T Ã some fraction of the electrosphere positrons will be replaced by protons to fit with the new long-distance, proton-dominated boundary condition. The exact proton-to-positron ratio of the electrosphere will be determined by the rates at which captured protons are annihilated by the nuggets and positrons are annihilated by external electrons, and the rate at which beta processes can replace near-surface positrons.
Further from the quark surface the positrons are more weakly bound and the thermal behavior becomes important for the distribution. In this regime the density as a function of height scales as [24,25], where the approximate value ofz is taken from matching this solution to the numerical solution from higher densities.
The main observation here is that a T ≠ 0 environment leads to the ionization of the loosely bound positrons such that the antinuggets will be in a negatively charged configuration with charge −Q estimated as follows: where we assume that some loosely bound positrons will be stripped off the electrosphere as a result of the nonzero temperature. 14 This negative charge of the antinugget implies that the protons from the plasma might be captured by the nugget by screening the charge (32). This obviously implies that the effective cross section for capturing the protons ∼4πR 2 cap ðTÞ will be drastically larger than 4πR 2 from our previous estimates (28)-(29) when the electrosphere is entirely made of positrons, not protons.
In principle, the distribution of protons surrounding the nugget should be determined through a Thomas-Fermi computation similar to that performed in Ref. [25] but allowing for the presence of protons as well as positrons and using the early Universe plasma density as the r → ∞ boundary condition. However, for present order-of-magnitude estimates we will assume a simple power-law scaling with exponent p, with the normalization n 0 set to match the total charge given in Eq. (32). This assumption is consistent with our numerical studies [25] of the electrosphere with p ≃ 6 for positrons. It is also consistent with the conventional Thomas-Fermi model at T ¼ 0; see Ref. [4] for references and details. We keep the parameter p arbitrary to demonstrate that our main claim is not very sensitive to our assumption on the numerical value of p. With these assumptions the baryon number density n 0 in close vicinity of the nugget can be estimated as follows: One should note that the behavior of the proton cloud may deviate significantly from Eq. (33) at very small and very large radii; however, we simply want to determine the approximate scale over which electromagnetic effects will act. In this context the simple form of Eq. (33) should be sufficient. This behavior will continue until the proton density matches that of the surrounding plasma, which gives us a radius for the overdensity of protons surrounding the nugget, This increase in the effective scattering length of the nuggets will boost the number of interactions and may result in an increased annihilation rate. Most importantly, these captured protons will spend an extended amount of time near the surface of the AQNs, giving them an increased opportunity to annihilate. Again, we stress that this increased rate of proton capture is effective only after the positrons are fully annihilated and the protons are the only positive charge carriers in the plasma, which happens at T ≃ T Ã ≃ 20 keV. 13 One should emphasize that the presence of the electrosphere itself is a very generic phenomenon, and its main features are determined by the boundary conditions deep inside the nugget where the lepton's chemical potential is fixed as a result of beta equilibrium, similar to the analysis in the context of strange stars; see Ref. [3] for a review. 14 We estimated the position of the cutoff z > z 0 ¼ ð2m e TÞ −1=2 in Ref. [4]. We emphasize that all of our estimates that follow are not very sensitive to the cutoff scale z 0 .
In particular, for p ≃ 6 the effective capture distance R cap is of order which of course drastically changes the collision rate, as will be estimated below. The scaling (36) holds as long as the thermal equilibrium between the nuggets and surrounding plasma is maintained. Equation (36) breaks down at sufficiently large R cap when the power-law scaling (33) is replaced by an exponential behavior due to Debye screening. Numerically, Debye screening becomes operational at R ≃ 10 R cap at T ≃ 20 keV; see Appendix A in Ref. [4].
We may now perform the same estimates as in Eq. (28) but using the larger capture cross section of Eq. (36), i.e., The total number of collisions during this time is saturated by the highest temperature T ≃ T Ã such that the integral can be estimated as follows: where we use p ¼ 6 for numerical estimates. Equation (38) represents a full analog of the estimate (29) obtained for the pre-BBN epoch.
While the number of scatterings occurring in this lowtemperature regime is slightly below the number occurring just after nugget formation as estimated in Eq. (28), we expect the evolution of the AQN baryon number to be dominated by these low-energy collisions in which the AQNs and baryonic matter temporarily form a bound state. This is because, as argued above, collisions at tens of MeV are highly likely to result in elastic scattering, while electromagnetically bound protons have a much larger opportunity to overlap with the quark modes of the color superconductor and eventually annihilate.
The similarity in the total number of collisions in the estimates (29) and (38) at T form ≃ 40 MeV and T Ã ≃ 20 keV, respectively, can be easily understood from the following simple observations. The baryon number density in plasma scales as T 3 , the proton's velocity in plasma scales as T 1=2 , and the cosmic time scales as T −2 . All of these factors result in drastic changes in the rate between T form ≃ 40 MeV and T Ã ≃ 20 keV with an approximate suppression factor ðT Ã =T form Þ 3=2 ∼ 10 −5 . However, this substantial suppression in the temperature is mostly compensated by an enhancement in effective cross section ðR cap =RÞ 2 ∼ 10 4 ; see the estimate (36). These two effects work in opposite directions which explains why our estimates for the collision rates (38) and (29) are numerically close.
We summarize this section with the following comment. The total number of collisions (38) is still much smaller than the typical baryon charge hBi ∼ 10 25 of the nuggets, such that the majority of the nuggets will survive the post-BBN epoch as only a small portion of the collisions will eventually lead to annihilation events. Therefore, from these estimates we conclude that the post-BBN epoch does not modify the primordial spectrum of the AQNs.
At this point a thoughtful and careful reader may wonder how it could happen that the number of annihilation events estimated above is sufficiently small that all nuggets with B > 10 24 can easily survive the unfriendly hot and dense environment of the early Universe according to the estimates (38) and (29). At the same time, it has been argued recently in Refs. [5,6,28] that nuggets of all sizes will experience complete annihilation in the solar corona when the AQNs enter the solar atmosphere. How can these two claims be consistent? We refer the readers to Appendix D which specifically addresses this question, where we emphasize a number of crucial differences between the two cases.
The only comment we would like to make here is that the drastic enhancement of the rate of annihilation in the solar corona is due to the propagation of AQNs with supersonic speed (above the escape velocity v > 600 km=s at the solar surface) in the ionized plasma with a very large Mach number M ¼ v=c s ≃ 10, where c s is the speed of sound in the solar atmosphere. It is well known that a moving body with such a large Mach number will inevitably generate a shock wave and an accompanying temperature discontinuity with turbulence in the vicinity of a moving body. As a result of this complicated nonequilibrium dynamics, the effective cross section may be drastically increased in the course of shock-wave propagation due to the capture (with subsequent annihilation) of a large number of ions from the plasma. These features of AQNs in the solar corona should be contrasted with the AQNs' adiabatic evolution in the plasma of the early Universe with its relatively slow evolution.

VIII. POST-RECOMBINATION EVOLUTION
The integral (38) is saturated by the highest possible value T ≃ T Ã ≃ 20 keV because the largest collision rate occurs precisely at that time. As the temperature slowly decreases due to the Universe's expansion one should not expect any dramatic changes until the recombination epoch at T ≃ 0.3 eV. At this point the Universe becomes neutral and the scattering cross section of matter with the AQNs no longer receives a boost from electromagnetic effects analogous to Eq. (36). Therefore, these generic arguments suggest that the collision rate diminishes even faster after recombination, implying that the size distribution of the AQNs essentially does not change during that epoch. These generic arguments obviously do not apply to violent environments during galaxy or star formation (which occur during this epoch), which will be analyzed later in this section.
In spite of this relatively dilute environment, the rare events of annihilation of the AQNs with surrounding baryons still occur even at such low density. The corresponding radiation due to the annihilation processes, while negligible in comparison with the dominant CMB radiation, may nevertheless leave some imprints which could be observed today, as argued in Ref. [12]. This is due to some specific features of the spectrum characterizing the AQN annihilation events: the low-energy tail of the radiation due to the annihilation processes with the nuggets has a spectrum ∼ ln ν which should be contrasted with conventional CMB blackbody radiation characterized by ν 2 behavior at low ν ≪ T; see Ref. [12] for the details.
As we mentioned above, the violent environment during galaxy or star formation requires a special treatment because it may potentially change the generic argument that this epoch is essentially irrelevant to the AQNs' survival pattern. In what follows we compare the environment during galaxy and star formation with corresponding features in the solar corona (where it has been argued that complete annihilation of AQNs occurs) in the context of the AQN annihilation rate. The outcome of this comparison and corresponding conclusion will be formulated at the very end of this section.
The complete annihilation of AQNs in the solar corona-as discussed in Refs. [5,6] and reviewed in Appendix D-is the direct consequence of a few factors: (1) The relatively large density in the transition region, n ∼ 10 11 cm −3 . (2) The very large velocity of the nuggets on the solar km=s. This is of course a result of strong gravitational forces ∼M ⊙ localized over a relatively small distance R ⊙ .
(3) The large velocity v greatly exceeds the speed of sound c s ≃ ffiffiffiffiffiffiffiffiffiffiffiffi T=m p p . This implies that the Mach number is very large, M ≡ v=c s ≫ 1, such that shock waves inevitably form. (4) The high ionization of the plasma due to the high temperature T ≃ 10 6 K in the transition region. The combination of these factors leads to complete annihilation of AQNs, as reviewed in Appendix D. While individual (violent) conditions from the list above may emerge during galaxy or star formation, the combination of all four elements does not occur, in general, during this epoch. Therefore, we do not expect any considerable modification of the size distribution of the nuggets after the recombination.
Indeed, the baryon density during the structure formation epoch does not exceed n B ∼ 1 cm −3 . Furthermore, the typical velocities of particles in the gas are of the same order of magnitude as the speed of sound, i.e., v ∼ c s ∼ 10 2 km=s such that one should use the conventional formula to estimate the collision rate without any additional enhancement factors related to the Mach number M, i.e., The total number of collisions N col during the Hubble time H −1 at redshift z ∼ 10 is of order which represents a tiny portion of the average baryon charge hBi ∼ 10 25 of a nugget. Furthermore, even if in some small regions the relative velocities of the AQNs and baryons exceed the speed of sound c s this does not lead to shock-wave formation, similar to our discussions about the solar corona (reviewed in Appendix D). This is because the shock-wave phenomenon is based on an effective description of the system when the hydrodynamical description is justified, which implies that the typical distance between the particles ∼n −1=3 must be much smaller than the size of a moving body ∼R. This approximation is obviously badly violated for the AQNs during the structure formation epoch. Therefore, according to our estimate N col ≪ hBi ∼ 10 25 , and we conclude that the size distribution of AQNs is not modified during the era of structure formation. A similar conclusion also holds for another violent epoch-star formation-which could also potentially modify the AQN size distribution. The corresponding analysis can be separated into two different stages: the final stage of formation when typical parameters are similar to our analysis of the Sun, and the initial stage of star formation characterized by the density ranging from n ∼ 10 15 to n ∼ 10 0 cm −3 depending on the size of the infall cloud, which ranges from r ∼ 10 −1 to r ∼ 10 6 AU; see Ref. [56].
We start our estimates from the final stage of formation when the stars assume their final form. In this case, all of the nuggets that will be captured by a star will be completely annihilated, similar to our studies of the Sun. However, the portion of the AQNs that will be captured by the stars is very tiny in comparison with the total number of nuggets. The corresponding portion of the affected nuggets can be estimated in the terms of the capture impact parameter b cap which is typically only a few times the star's size R ⋆ , where v ∼ 10 −3 c is the typical velocity of the nuggets far away from the star. The rate of the total mass annihilation dM ann =dt of all nuggets captured (and consequently annihilated) by the star can be estimated as follows: where we used the solar parameters for the numerical estimates and assumed that ρ DM is saturated by AQNs. 15 The upper limit of the total mass annihilated by a single star can be estimated by multiplying Eq. (41) by the total lifetime of star, which can be approximated as H −1 , i.e., The estimate (42) should be compared with the total mass of the star M ⋆ ∼ M ⊙ ∼ 10 30 kg. As the stars represent only a fraction of the total baryonic matter of the Universe, and there is 5 times as much DM as baryonic matter, one can infer from Eq. (42) that the annihilated portion of DM (due being capturing by stars) represents only a small portion (≲10 −10 ) of the total dark matter content of the Universe. We now turn to the estimates of the AQN annihilation pattern during the initial stage of star formation. In this case, the DM nuggets passing through the infall cloud experience annihilation events. The corresponding total annihilated baryon charge for a single nugget (as a result of this passage) can be estimated as follows: which represents a tiny portion of the average baryon charge hBi ∼ 10 25 . In the estimate (43) we used the most generic configurations for when AQNs enter a large region with L ∼ 10 6 AU characterized by n B ∼ 1 cm −3 , while the passage of nuggets through a small, well-localized region L ∼ 10 −1 AU with high density n B ∼ 10 15 cm −3 is highly unlikely as it represents a very small portion of the total AQN flux. But even in this case the total number of annihilation events N col ∼ 10 17 remains small in comparison with the average baryon charge hBi ∼ 10 25 .
We conclude this section with the following comment. Violent events such as galaxy formation or star formation after recombination do not drastically modify the size distribution of the nuggets, similar to our previous analysis devoted to the different epochs of the Universe's evolution. The basic reason for this conclusion is that the nuggets-which can undergo complete annihilation-represent only small portion of the entire population according to Eq. (42), while the majority of the nuggets will lose a very tiny portion of their baryon charge during the Hubble time according to Eq. (39).

IX. PRESENT-DAY MASS DISTRIBUTION
After recombination, the size distribution of AQNs is essentially fixed until the present day. 16 We may formulate the change in baryon number as a result of post-formation annihilation as where f 1 and f 2 are the fractions of collisions in the positron-dominated and proton-dominated phases which result in the annihilation of a unit of baryon charge from the AQN with the collision rates taken from Eqs. (29) and (38), respectively. If we assume the original nugget distribution defined by B min , B cut , and α as discussed in Sec. IV D, we may use Eq. (44) to translate the distribution formed above T ≈ 40 MeV to a present-day mass distribution. Note that since R ∝ B 1=3 the first term in Eq. (44) scales as B 2=3 according to Eq. (29), while the second scales as B 5=9 according to Eq. (38).
We may then write the late-time baryon number distribution based on the result (26) obtained from percolation theory as follows 17 : where B is the present-day baryon number of the AQN and we have ignored the small portion ΔB related to the annihilation processes as discussed above. For a given environment the constant N 0 may be fixed if AQNs are assumed to provide all of the dark matter mass. 18 We also assume that one and the same α describes the mass distribution for all values of B. It may or may not be a correct assumption as some nanoflare models fit the solar 15 We comment here that in the case of the Sun the energy released as a result of annihilation events with the rate (41) represents approximately 10 −7 of the total solar luminosity radiated from the solar corona in the form of the EUV and x rays, which represents the resolution of the solar corona heating problem within the AQN scenario, as suggested in Refs. [5,6]. 16 Obviously, the distribution may develop localized anisotropy in regions of particularly high matter density such as stars or planets, but these local effects are not the subject of the present work. 17 Equation (26) formally holds only for asymptotically large B ≫ B min , but in fact it remains valid for almost the entire region of B with the exception of a small region in the vicinity of B min ; see Fig. 4 in Appendix B. 18 As we reviewed in Sec. II, conventional axion production due to the misalignment mechanism and domain-wall decay always accompanies the nuggets' formation and contributes to the total DM density; see Fig. 1. However, the relative contribution of these well-studied mechanisms to ρ DM strongly depends on the axion mass m a and it is likely to be negligible for sufficiently large axion mass; see Ref. [9] for the corresponding plots.
corona heating with different exponents α for small and large values of B; see below for additional comments.
With the assumptions just formulated, one can represent ρ DM as follows: where dN has the physical meaning of the AQN number density per unit volume 19 per baryon charge interval dB, while m p BðdN=dBÞ has the physical meaning of the mass density per unit volume hidden in the form of nuggets with baryon charge from B to B þ dB. The relation (46) allows us to solve for the normalization factor in Eq. (45) for different values of the exponent α: where we have again assumed that ΔB ≪ B represents a small fraction of the initial baryon number as argued to be the case in previous Secs. VI and VII. Note that α ¼ 2 marks the slope at which the distribution transitions from being mass dominated by the bottom end to the higher end. Distributions with α > 2 are largely defined by the mass scale B min , while the typical mass scale for shallower sloped distributions with α < 2 is dominated by B cut .
Given the profile of Eq. (45), we may now ask what range of parameters are consistent with the observational constraints discussed in Sec. II B. As argued above, we expect that f 1 ≪ f 2 so that the majority of annihilations occur below T ∼ 20 keV; if this is the case, then the parameter space of the AQN model may be largely defined by four variables: B min , B cut , α, and ΔB ∼ f 2 . As discussed in Sec. IV D, these parameters are not theoretically well constrained. Indeed, while the theoretical analysis predicts the generic power-law behavior (26), the numerical value for the exponent α is expressed in terms of the parameter δ according to Eq. (27) which itself describes some features of the bubble's formation during the QCD transition at T ∼ T c , which are basically unknown in strongly coupled QCD.
The observational constraints are less trivial and much more interesting. First of all, there are constraints on α which come primarily from solar data. As we discussed in Sec. II B 4, the solar corona measurements are sensitive to the full distribution rather than simply the average value hBi. If the nuggets are to explain the solar heating problem as argued in Refs. [5,6,28], then we require that the majority of heat input comes from lower-energy unobservable events and thus must have α > 2. This option is consistent with the analysis of the nanoflare distribution performed in Ref. [30] where the authors claimed that the best fit to the data is achieved with α ≃ 2.5, while numerous attempts to reproduce the data with α < 2 were unsuccessful.
Another option advocated in Ref. [31] is that the nanoflare distribution consists of two different exponents: the events below E ≤ 10 24 erg are described by α ≃ 1.2, while the higher-energy events with E ≥ 10 24 erg are described by α ≃ 2.5. In terms of the baryon charge distribution (45), the model in Ref. [31] corresponds to the nugget mass distribution with α ≃ 1.2 and a cutoff scale near B cut ∼ 10 27 . The higher-mass nuggets with B cut ≥ 10 27 are described by α ≃ 2.5.
Assuming that the distribution (45) can be approximately used in close vicinity of B min one can relate the average baryon number hBi and B min as so that for any α > 2 the average hBi and B min are at the same scale. One should not literally use the relation (48) as the distribution (45)  For the AQN population to survive as high-mass dark matter candidates we require B min > ΔB which, from Eq. (38), suggests a lower bound below ∼10 24 , although (as argued above) this scale considers the number of collisions between the background plasma with an AQN rather than the number of annihilation events ΔB. The only robust constraint on B min comes from the observational side [Eq. (8)] which can be expressed in terms of hBi as discussed above.
These considerations lead us to two general classes of AQN distributions which will be consistent with all known direct and indirect constraints. These models are essentially equivalent to a variety of nanoflare distributions [30,31] expressed in terms of the baryon charge number with the only additional condition. The nanoflare models must satisfy the condition B min > 10 24 in order to be consistent with the independent observational constraint (8). There are plenty of models in Refs. [30,31] that satisfy this condition and saturate the required energy budget for the corona heating, and we shall not elaborate on this topic here.
We consider this phenomenon to be a highly nontrivial self-consistency check of the AQN framework when the 19 Note that the definitions of N and N 0 here are slightly different from those in Sec. IV where N and N 0 are the total numbers rather than number densities. allowed window for the baryon charge B and the fitted energy spectrum for nanoflares overlap, with the identification E ≃ 2m p c 2 B. On the one hand, this window represents a cumulative constraint from a number of astrophysical, cosmological, satellite, and ground-based observations and experiments (as reviewed in Sec. II), which are consistent with analytical results based on percolation theory discussed in Secs. III and IV. On the other hand, this window largely overlaps with constraints originated from completely independent physics: solar corona heating when the nanoflare distribution (11) with energy E is identified with the AQN distribution (45) with baryon charge B.

X. CONCLUSION
The main results of this work can be formulated as follows.
(1) We used an approach coined as the envelopefollowing method to overcome a common numerical problem with the drastic separation of scales in the system. In our case the scales are the QCD scale ∼Λ QCD , the axion scale ∼m a , and the cosmological time scale t 0 ∼ 10 −4 s. The results support our original assumptions that the chemical potential inside the nugget indeed assumes a sufficiently large value μ form ≳ 450 MeV during this long cosmological evolution. This magnitude is consistent with the formation of a CS phase; see Fig. 2. (2) The nuggets complete their formation precisely in the region of T form ≈ 40 MeV where they should (see Fig. 2), as this corresponds to the temperature where the baryon-to-photon ratio η assumes its present value (5). (3) Items 1 and 2 represent a highly nontrivial consistency check of the AQN framework when three drastically different scales (Λ QCD , axion mass m a , and cosmological time scale t 0 ) "conspire" to produce a self-consistent picture. (4) We argued that the nuggets' distribution must have the algebraic behavior (26) as a direct consequence of a generic feature of the percolation theory. The exponent α cannot be predicted theoretically, but according to Eq. (27) it can be expressed in terms of another parameter which is sensitive to axion domain-wall formation during the QCD epoch. (5) We argued that the nuggets survive both the pre-BBN and post-BBN evolutions (as estimations in Secs. VI and VII show) as long as they are sufficiently large to satisfy the observational constraint (8). The essential reason for this is that the fraction of annihilated AQNs scales with the cross section-to-mass ratio ∼B −1=3 and thus nuggets of sufficiently large B will remain largely unaltered by annihilation events after their formation. (6) We argued that the present-day baryon charge distribution (45) is consistent with the nanoflare distribution (11) which was fitted to describe the solar corona observations. This represents a highly nontrivial consistency check of the proposal [5,6,28] that the AQNs made of antimatter are the nanoflares postulated long ago. It is the central claim of this work that there exists a larger amount of the allowed parameter space across which the AQN model is consistent with all available cosmological, astrophysical, satellite, and ground-based constraints. While the model was invented long ago to explain the observed relation Ω dark ∼ Ω visible , it may also explain a number of other observed phenomena, such as the excesses of galactic emission in different frequency bands as reviewed in Sec. II. This model also offers a resolution of the so-called "primordial lithium puzzle" and the 70-year old "solar corona mystery" (as mentioned in the Introduction) when one uses precisely the same parameters and the same distribution (45) advocated in the present work. However, all of these manifestations of the AQN model are indirect in nature. Therefore, one can always find tons of alternative explanations for the same phenomena.
In this respect the recent proposal advocated in Refs. [57][58][59] to search for axions (which will inevitably be produced as a result of the annihilation processes of the antimatter nuggets with surrounding matter) is a direct manifestation of the AQN model. These axions will be emitted when AQNs disintegrate in the solar corona. Axions will also be emitted when the nuggets hit the Earth and continue to propagate deep underground and loose baryon charge, accompanied by the emission of axions. In fact, the observation of these axions with very distinct spectral properties in comparison with conventional galactic axions would be the smoking gun supporting the entire AQN framework. We finish this work on this positive and optimistic note.

ACKNOWLEDGMENTS
This work was supported in part by the Natural Sciences and Engineering Research Council of Canada.
APPENDIX A: Nugget evolution from T 0 to T form In this Appendix we discuss the time evolution of a nugget from T 0 to T form both analytically and numerically. We start with the Lagrangian that dominates the nugget evolution, There is a slight difference between this expression and the Lagrangian adopted in Refs. [7,8]. Here we replace the domain-wall tension σ ¼ 8f 2 a m a with the effective domain-wall tension σ eff ¼ κ · σ. The phenomenological parameter κ accounts for the difference between the domainwall tension of a nugget σ eff and that of a planar domain wall σ [9]. In general, the effective domain-wall tension σ eff is smaller than σ with 0 < κ < 1. 20 ΔP is the pressure difference inside and outside the nugget, which is [7] ΔP ¼ P The first term on the right-hand side of Eq. (A2) is the Fermi pressure inside the nugget. The second term is the contribution from the MIT bag model with the famous "bag constant" E B ∼ ð150 MeVÞ 4 . Θ is the unit step function which implies that this term turns on at large chemical potential μ > μ 1 when the nugget is in the CS phase, while it vanishes at small chemical potential μ < μ 1 when the nugget is in the hadronic phase. The parameter μ 1 is estimated to be ∼330 MeV [1] when the baryon density is close to the nuclear matter density. The third term is the pressure from quark-gluon plasma (QGP) at high temperature outside the nugget. The parameter g out ≃ ð 7 8 4N c N f þ 2ðN 2 c − 1ÞÞ is the degeneracy factor of the QGP phase.
From the Lagrangian (A1) we obtain the equation of motion [Eq. (12)], Following the same procedure as in Ref. [7], we insert the QCD viscosity term ∼4η _ R R to effectively describe the friction for the domain-wall bubble oscillating in an unfriendly environment. The difference from the equation of motion in Ref. [7] is that here we have an extra term ∼ _ σ eff _ R. This term occurs because the tension σ eff itself is a function of time since σ eff ¼ κ · 8f 2 a m a ðtÞ. We treat the axion mass m a more precisely for it is a time (temperature)dependent function rather than a constant. The parameter f a and the function m a ðtÞ can be obtained from the axion model. Then, with an appropriate choice of κ, we can determine how σ eff evolves with the cosmological time.
According to Refs. [7,8], the baryon charge accumulated on the wall is where R is the radius of the nugget, and g in ¼ 2N c N f ≃ 12 and μ are the degeneracy factor and the chemical potential of the baryon charge in the vicinity of the wall, respectively. The accumulated baryon charge B wall is assumed to be a constant [7], which can be expressed as In principle, we can get the time evolution of the nugget by numerically solving the two ordinary differential equations (A3) and (A5). However, before we do the numerical calculations, we can get some profound analytical results from the above equations. The nugget starts its evolution at T 0 with the initial chemical potential of the baryon charge on the wall being approximately zero, μ 0 ≃ 0. From Eq. (A4), we can get the initial baryon charge accumulated on the wall, Then the nugget completes its formation at T form when the nugget stops oscillating with _ RðtÞ ≃ 0,̈RðtÞ ≃ 0, _ μðtÞ ≃ 0. All features of the nugget (radius, chemical potential, etc.) should remain almost constant after the formation point (T ¼ T form ) until the very end (t → ∞) (T → 0). Thus, with R form ≃ RðT ¼ 0Þ and μ form ≃ μðT ¼ 0Þ we get According to Eq. (A5), B wall is conserved during evolution. Therefore, by equating Eq. (A6) with Eq. (A7) we arrive at Also, with all of the derivative terms vanishing after T form , from Eq. (A3) we get 20 There are two main reasons for the difference between σ eff and σ, which are discussed in detail in Ref. [9]. We briefly summarize the two reasons here. The first reason is that nuggets with baryon charge accumulated inside will finally become stable in the CS phase. Thus, in our case the axion domain-wall solution interpolates between topologically distinct vacuum states in the hadronic (outside the nugget) and CS (inside) phases, in contrast to a conventional axion domain wall which interpolates between distinct hadronic vacuum states. The chiral condensate may or may not be formed in the CS phase, which could make the topological susceptibility in the CS phase much smaller than in the conventional hadronic phase. The second reason is that σ ¼ 8f 2 a m a is derived using the thin-wall approximation, which could be badly violated in the case of the closed domain wall when the radius and the width of the wall are of the same order of magnitude. This effect is expected to drastically reduce the domain-wall tension.
where the pressure difference ΔPðT ¼ 0Þ can be obtained from Eq. (A2), We notice that R form and μ form are completely solvable from Eqs. (A8) and (A9), and they are determined by the initial values R 0 and T 0 . An important feature of nugget evolution is that baryon charges not only occur on the wall of the nugget, but they also accumulate in the bulk of the nugget. The chemical potential μ on the wall gradually increases due to nugget contraction. As a consequence, the chemical potential in the bulk of the nugget increases, maintaining equilibrium with the chemical potential on the wall, which causes the accumulation of baryon charges in the bulk of the nugget. As explained in Ref. [7], the net flux of baryons entering and leaving the nugget ΔΦ≡Φ in −Φ out is negligibly small, but the sum of these two fluxes hΦi ≡ 1 2 ðΦ in þ Φ out Þ is very large and the nugget can entirely refill its interior with fresh particles within a few oscillations. The high exchange rate hΦi implies that the entire nugget could quickly reach chemical equilibrium. The result is that the initially baryonically neutral nugget evolves into one completely filled with quarks (or antiquarks for antinuggets). Therefore, we expect that the nugget will become stable at T form with the entire nugget in equilibrium, with the same chemical potential μ form . Thus, the total baryon charge carried by the stable nugget is Then, using Eq. (A8) we can express B as where K ≡ π 2 27 ffiffi 6 p g in is a constant introduced for convenience. This relation is applied in Sec. IV to calculate the baryon charge distribution of nuggets.
Next, we are going to numerically solve the differential equations (A3) and (A5) to get the nugget evolution. To numerically solve the two equations, we first need to know how the effective domain-wall tension σ eff ðtÞ ¼ κ · 8f 2 a m a ðtÞ evolves as a function of time. One of the most updated results for the axion mass m a ðTÞ is based on high-temperature lattice QCD [54]. The topological susceptibility of QCD, χðTÞ, is plotted in Fig. 2 in Ref. [54] as a function of the cosmological temperature T. The data points of the figure are also provided in Table 9 in the Supplementary Information of the same paper, fitting which we get the expression for χðTÞ as where Θ is the unit step function. Then we can get the axion mass using the relation Equations (A13) and (A14) explicitly show that before the QCD transition the axion mass increases rapidly with the exponent β ¼ 7.85=2 ¼ 3.925 as the cosmological temperature decreases (see Ref. [53] for a similar result). Then the axion acquires its asymptotic mass near the QCD transition and remains constant after that. The cosmological time-temperature relationship is also useful in our numerical calculations, which in the radiationdominated era is well known as where g ⋆ ðTÞ is the number of effective degrees of freedom of all relativistic particles at temperature T. Since the major part of the nugget's evolution occurs after the QCD transition, we treat g ⋆ ðTÞ as a constant for simplicity with g ⋆ ¼ 17.25 (see, e.g., Ref. [60]) as in the hadronic phase.
In the present work, we are not going to solve the nugget evolution with the parameters (such as κ, f a , T 0 , and R 0 ) taking many different values. Instead, we solve it with these parameters by taking a group of reasonable values as an example. They are taken as κ ¼ 0.04, f a ¼ 10 10 GeV, T 0 ¼ 200 MeV, and R 0 ¼ 6 × 10 −4 cm. Of course the group of parameters can vary a lot, but it is not the subject of the present work to numerically study how varying these parameters will affect the nugget evolution.
As we discussed in Sec. III, when we try to solve the two differential equations we immediately meet the multiplescale problem ωτ ∼ 10 10 which implies that the nugget will not settle down until after billions of oscillations. The enormous number of oscillations make the code extremely time consuming, which makes it almost impossible to completely solve the system numerically. To make the numerical solution feasible, the QCD viscosity term η was artificially enlarged by 8 or 9 orders in our previous calculations in Refs. [7,8].
In this paper, we adopt a different numerical method coined as the "envelope-following method" to solve the system, with the viscosity term keeping its physical magnitude η ∼ Λ 3 QCD when the parameter ωτ ∼ 10 10 assumes its very large physical value. We believe that this will make our numerical result for the nugget evolution more trustworthy. The motivation for using the envelope-following method and how it works are briefly explained below.
We notice that although a nugget oscillates very fast during evolution, the amplitude of oscillation decreases very slowly for each given cycle. The peaks of oscillations in fact form a "smooth" line which we call an envelope. We realize that if we can find a way to numerically solve the envelope, then it is unnecessary to know the full details of all oscillations. The envelope-following method turns out to be very beneficial for our study of the nugget oscillations. The method is efficient in solving highly oscillatory ordinary differential equations, which was illustrated in Ref. [61]. We briefly summarize the basic idea here.
We start with the initial conditions R ¼ R 0 , _ R ¼ 0 (and μ ¼ μ 0 ), which corresponds to the first peak of R oscillations and we label this peak as point a. Then we solve the differential equations until we get the next peak b of R oscillations, which should be slightly smaller than the first peak. This step is very fast since we solve the equations for just one oscillation. Joining points a and b, we get a secant line. This secant line is then used to project the solution to point a 0 which is many oscillations away. Starting with a 0 as the new peak, we solve the differential equations until we get the next peak b 0 , etc. We repeat the above procedure of drawing the secant line, projecting the solution, and finding the next peak. After several projections, we get the upper envelope of R oscillations. Using the same method, we can find the lower envelope of R oscillations and also the envelopes of μ oscillations. We should point out that, although the details of the oscillations are not important to us, we can recover them locally if we substitute the corresponding envelope information into the differential equations as the initial conditions.
We plot the numerical result for the nugget evolution solved using the envelope-following method in Fig. 2 with the group of parameters chosen above. In Fig. 2, we see that in this case the nugget completes its evolution at ∼40 MeV. The formation chemical potential is ∼450 MeV, well above the threshold of the CS phase μ 1 ¼ 330 MeV. Also, we see that the theoretical analysis of R form and μ form (denoted by the dashed blue and dashed orange lines, respectively) from Eqs. (A8) and (A9) matches the numerical result pretty well, which verifies the validity of the two equations as well as the relation (A12) (B ∝ R 3 0 T 3 0 ).

APPENDIX B: Calculations of NðBÞ and dN=dB
In this Appendix we show the details of calculating Eqs. (17) and (23), to support the results in Sec. IV D. We first rewrite Eq. (17) as with the limits of integration written explicitly, which can be explained as follows. The integral (B1) is performed over the region KR 3 0 T 3 0 ≤ B with the constraints T c ≲ T 0 ≲ T osc and R 0 ≳ ξðT 0 Þ from the model of T 0 and R 0 distributions. We show the region of integration in Fig. 3, where the parameter space T c ≲ T 0 ≲ T osc and R 0 ≳ ξðT 0 Þ is represented by the colored region. The green lines are the contour lines of B with KR 3 0 T 3 0 ¼ B for different values of B. Then the region of integration is the area enclosed by the solid black lines and one of the green lines (to the left of the green line), from which we can obtain the limits of integration. The lower limit of R 0 is R lower ¼ ξðT 0 Þ; the upper limit of R 0 is on the green line R upper ¼ ½B=ðKT 3 0 Þ 1=3 ; the lower limit of T 0 is T c . The upper limit of T 0 is a little complicated: it could either be the intersection of the line ξðT 0 Þ and the green line which is T upper ¼ T c · ½B=ðKξ 3 min T 3 c Þ 1=3ðβþ1Þ , or simply T upper ¼ T osc , depending on the value of B. However, we are not likely to have the chance to use the latter case with T upper ¼ T osc as the upper limit of T 0 , which is explained as follows.
If we want the upper limit of T 0 in the integration to be T osc , then B has to be larger than B cross ¼ Kξ 3 ðT osc ÞT 3 osc , the value of B at the crossing point where the line ξðT 0 Þ intersects the horizontal line T 0 ¼ T osc . We should compare B cross with the minimal baryon charge B min ¼ Kξ 3 min T 3 c which corresponds to the closed domain wall initially forming at T 0 ¼ T c with R 0 ¼ ξ min . We get where we approximate it using T osc =T c ≃ 10 and β ≃ 3.925. We see that the range is 15 orders of magnitude wide, which is large enough for us to match the baryon charge distribution of nuggets with the energy distribution of solar nanoflares. Therefore, we choose T 0 as T upper ¼ T c · ½B=ðKξ 3 min T 3 c Þ 1=3ðβþ1Þ for the upper limit in Eq. (B1). Next, we are going to calculate Eq. (B1). Using the definitions r ¼ R 0 =ξ min and u ¼ T 0 =T c , we rewrite Eq. (23) in a more concise way, Substituting fðr; uÞ into Eq. (B1) and using the definition b ¼ B=B min , we arrive at from which we get This can be further simplified using where m ≡ 3ðβ þ 1Þðτ − 1Þ þ βδ and n ≡ 2ðβ þ 1Þ; Γðs; xÞ ¼ R ∞ x t s−1 e −t dt is the incomplete gamma function. To obtain the last approximate equality, we neglect the term Γð− 1þm n ; λb 2=3 Þ since it is far smaller than the term Γð− 1þm n ; λÞ for b ≫ 1. The condition b ≫ 1 is satisfied for a wide range of B values, which are generally several orders larger than B min . Substituting Eq. (B6) into Eq. (B5), we arrive at We see that dN=dB follows a power-law distribution, which verifies the relation (26). The finite-cluster parameters τ (contained in m) and λ that we discuss in Sec. IV B only affect the relative magnitude of dN=dB, but not the slope of the power-law distribution −α.
The parameter β describing the relation between the axion mass and cosmological temperature is well calculated [54] (see also Appendix A for more details). The other parameter δ from the model of the T 0 distribution [Eq. (22)] is relatively adjustable, which can result in different slopes for the power-law distribution dN=dB. This parameter (which may have any sign) describes the distribution of bubble formation. As we explained in the main text, a positive δ corresponds to the preference for bubble formation close to T osc , while a negative δ corresponds to the preference for bubble formation close to T c with a much stronger tilt of the axion potential.
We plot the baryon charge distribution of nuggets in Fig. 4. We choose τ ¼ 2, λ ¼ 10, and β ¼ 3.925 for both panels. The difference between them is the value of δ which is highly underdetermined. In Fig. 4(a), we choose δ ≈ −1 to make α ¼ 1.2. The solid black line is the plot of Eq. (B5), which represents the exact result for the distribution. As a comparison, we also plot the approximate relation (B8) represented by the dashed green line, which is straight in the log-log scale. We see that the approximate relation (B8) matches the exact result (B5) pretty well after the turning point where the condition b ≫ 1 becomes valid. In Fig. 4(b), we consider the case δ ≈ −4 which corresponds to α ¼ 2. All other ingredients are the same as in the first panel.
APPENDIX C: On the estimate of the parameter κðTÞ during pre-BBN evolution The main goal of this Appendix is to make an order-ofmagnitude estimate of the parameter κðTÞ which is an important element of our analysis in Sec. VI. In simple quantum-mechanical terms, the parameter κðTÞ is defined as the transmission coefficient for the baryon to enter and annihilate inside the nugget, while the reflection coefficient ½1 − κðTÞ describes the probability for the baryon to reflect off of the sharp surface of the nugget. One should emphasize that the system cannot be formulated in the simple terms normally used for a one-particle quantummechanical description. Instead, a proper study of this phenomenon would require a nonperturbative quantumfield-theoretic description as many-body effects play a crucial role in the computations. This is because the key element of the analysis must be a description of physics at the interface between the hadronic phase (as the proton's wave function is formulated in terms of the quarks) and a CS phase characterized by the diquark vacuum condensate. One should emphasize that the diquark condensate is not a local object, but rather a complicated coherent superposition of quarks similar to a Cooper pair in conventional superconductors.
In spite of the complicated annihilation pattern of a baryon in the hadronic phase with diquarks in the CS phase as highlighted above, one can easily carry out a simple, order-of-magnitude estimate demonstrating that κðTÞ ≪ 1, which is precisely the main goal of this Appendix. To simplify things we separate the suppression factor κ into two pieces: κ ∼ κ 1 · κ 2 , where κ 1 is defined as the dynamical suppression factor, while κ 2 is defined as the kinematical suppression factor (see below).
The dynamical suppression factor κ 1 can be easily understood from the internal structure of the axion domain wall represented by the heavy η 0 field which accompanies the axion field in the AQN model. The corresponding very sharp QCD structure has a width ∼m −1 η 0 ≪ Λ −1 QCD ; see Ref. [62]. Therefore, one should expect some suppression due to the sharp potential, where we assume that the typical energy E of the three incoming quarks making up the proton is of order T, while the typical strength of the "potential" is order of Λ QCD . We emphasize that we should formulate the problem in terms of quarks (rather than in terms of a single proton's wave function) because the annihilation pattern should include three antiquarks from a different CS phase. The suppression parameter κ 1 ðTÞ represents the probability for simultaneous transmission of three quarks. We note that the factor E=Λ QCD ≪ 1 in Eq. (C1) for a sharp surface is a very generic quantum-mechanical feature, which is not even sensitive to the sign of the interaction. It can be checked with the simple quantum-mechanical problem of the scattering of a low-energy particle with E → 0 on a δðxÞ potential. There is another suppression factor κ 2 ðTÞ related to the strong mismatch between the wave functions of the quarks in the hadronic phase and antiquarks in the CS phase. The corresponding suppression always has a factor 1=N! for N constituents, which is N ¼ 3 for a proton. 21 κ 2 ðTÞ also depends on the overlapping features of the wave functions from drastically different phases (hadronic vs CS phase). Therefore, we can represent κ 2 ðTÞ as follows: Collecting the estimates (C1) and (C2) together, one should expect that 21 We use the generic factor N instead of 3 in Eq. (C2) on purpose to emphasize that the annihilation process of N different constituents must find their counterparts with precisely matching wave functions for successful annihilation. ðC3Þ which represents our final estimate for the suppression factor κðTÞ used in Sec. VI. The order-of-magnitude estimate (C3) is sufficient for qualitative arguments, suggesting that the nuggets easily survive the dense and hot environment after formation.
It is interesting to note that the solitons and antisolitons [which can be represented as coherent superpositions of a large number (N → ∞) of constituents] in condensed matter physics do not normally annihilate, but rather experience an elastic scattering, which can be thought as the manifestation of the factor 1=N! in Eq. (C2). It is also known that the antiproton does not easily annihilate with a large nuclei (considered to be in a nuclear matter phase), but could have a lifetime as long as ∼20 fm=c instead of the conventional ∼fm=c scale [63]. This suppression of annihilation can be attributed to our "overlapping" factor in Eq. (C2).

APPENDIX D: AQNs in the corona: Turbulence and effective cross section with plasma
The main goal of this Appendix is to explain the crucial differences between our analyses presented in Secs. V-VII where we argued that very few annihilation events may occur in the early Universe. At the same time, it has been argued in Refs. [5,6,28] that nuggets of all sizes will undergo complete annihilation in the solar corona when AQNs enter the solar atmosphere. There is no contradiction between these two claims.
Indeed, our estimates (29) and (38) indicate that ΔB ≪ B during the entire evolution of the Universe from soon after the nuggets' formation at T ≈ 40 MeV until the present time. We consider two different regimes: during the hottest period of evolution with T Ã < T < 40 MeV (when the nuggets' electrosphere is made of positrons and the annihilation of baryon charge is highly unlikely), and when the temperature drops to below T Ã ≈ 20 keV (when the electrosphere contains a significant number of protons and baryon charge annihilation may occur as estimated in Sec. VII).
This number of collisions in plasma at T Ã ≈ 20 keV should be compared with the number of collisions when an AQN enters the solar corona. The corresponding estimate goes as follows; see Refs. [5,6,28] for the details. First, we estimate the ionized charge of the nugget in terms of the internal temperature T I [similarly to our estimate Eq. (32)], where T P is the temperature of the surrounding plasma. The difference between T P and T I is enormous (as emphasized in Ref. [6]) and is related to the fact that the nugget propagates in the solar corona with velocity v ≳ 600 km=s (the escape velocity of the Sun) which greatly exceeds the speed of sound c s in the corona, i.e., the Mach number M ≡ v=c s ∼ 10. It is well known that a moving body with such a large Mach number will inevitably generate shock waves. It is also known that a shock wave generates a discontinuity in temperature, which for large Mach numbers M ≫ 1 can be approximated as follows: where the temperature T 1 ≃ T P is identified with the temperature of the surrounding unperturbed plasma, while the high temperature T 2 can be thought of as the internal temperature of the nuggets T I . This very high internal temperature T I =T P ∼ M 2 that develops due to the shock wave makes a huge difference in comparison with the analysis in Sec. VII for plasma in thermal equilibrium at T ≃ T Ã . The effective cross section of AQNs with surrounding protons can be estimated as πR 2 eff , where R eff corresponds to the distance where protons from the plasma with energy ∼T P can be captured by the nugget, Combining Eqs. (D1) with (D3), one can estimate the enhancement of the effective cross section in comparison with the naive estimate πR 2 as follows: This enhancement factor, of course, is different for different nugget velocities as the Mach number M varies with v. This enhancement factor obviously also changes with time and solar altitude as a result of the motion with the friction and radiation. One should also emphasize that there will be very efficient energy and momentum exchange and heat transfer to the surrounding plasma as a result of nonequilibrium dynamics in the form of the turbulence that develops in the vicinity of the shock-wave front. It should be contrasted with the analysis in Sec. VII when the system is in perfect thermal equilibrium. The numerical simulations performed in Ref. [6] suggest that most of the nuggets lose their entire baryon charge due to the annihilation ΔB ≃ B in the vicinity of the transition region at an altitude ∼2000 km where it is known that drastic changes in temperature and pressure occur. From the AQN perspective such unusual features of the transition regions as well as the "solar corona heating puzzle" are understood in terms of dark matter nuggets which continuously hit the solar atmosphere with a flux that has the correct magnitude to saturate the EUV and soft x-ray radiation from the corona.