DAMA/LIBRA annual modulation and Axion Quark Nugget Dark Matter Model

The DAMA/LIBRA experiment shows $9.5 \sigma$ evidence for an annual modulation in the $(1-6)~ {\rm keV}$ energy range, strongly suggesting that the observed modulation has the dark matter origin. However, the conventional interpretation in terms of WIMP-nucleon interaction is excluded by other experiments. We propose an alternative source of modulation in the form of neutrons, which have been liberated from surrounding material. Our computations are based on the so-called Axion Quark Nugget (AQN) dark matter model, which was originally invented long ago to explain the similarity between the dark and visible cosmological matter densities, i.e. $ \Omega_{\rm dark} \sim \Omega_{\rm visible}$. In our proposal the annual modulation is shown to be generated in keV energy range which is consistent with DL observation in $(1-6)~ {\rm keV}$ range. This keV energy scale in our proposal is mostly determined by spectral properties of the neutrinos emitted by the AQN dark matter particles, while the absence of the modulation with energies above 6 keV is explained by a sharp cutoff in the neutrino's energy spectrum at $\sim 15$ MeV. This proposal can be directly tested by COSINE-100, ANAIS-112, CYGNO and other similar experiments. It can be also tested by studying the correlations between the signals from these experiments and the signatures from drastically different detectors designed for studies of the infrasonic or seismic events using such instruments as Distributed Acoustic Sensing (DAS).

However the annual modulation observed by DL is excluded by other direct detection experiments if interpreted in terms of the WIMP-nuclei interactions. This motivated a number of alternative explanations for the DL signal. In the present work we argue that the modulation observed by DL is due to the neutrons surrounding the detector. In this respect our proposal is similar to the previous suggestions [7][8][9][10] where the authors agued that the induced neutrons (which have been liberated from material surrounding the detector) may be responsible for the observed annual modulation.
Our proposal is drastically different from previous suggestions [7][8][9][10] in one but crucial aspect: the neutrons in our case are induced by neutrinos emitted by the Axion Quark Nugget (AQN) dark matter particles. Therefore, the annual modulation observed by DL has truly genuine DM origin, though it is manifested indirectly in our framework through the following chain: AQN → (neutrinos) → (surrounding neutrons) → DL. (1) In this framework the modulation should obviously show a proper period of 1 year with proper phase as the neutrinos from (1) are originated from dark matter nuggets, and the corresponding time variation will be obviously transferred to the modification of the neutron flux, which eventually generates the modulation observed by DL.
One should emphasize that the emission features of the neutrinos emitted by AQNs such as the intensity and spectrum (which ultimately determines the (1 − 6) keV energy recoil for the observed DL annual modulation) have been computed in the AQN model long ago for completely different purposes, and we will use exactly the same original parameters of the model without any intension to modify them to fit the DL observations.
We overview the basic ideas of the AQN model in the next section II. One should emphasize that this model is consistent with all available cosmological, astrophysical, satellite and ground based constraints, where AQNs could leave a detectable electromagnetic signature. While the model was initially invented to explain the observed relation Ω dark ∼ Ω visible , it may also explain a number of other (naively unrelated) phenomena, such as the excess of galactic emission in different frequency bands. The AQN model may also resolve other, naively unrelated astrophysical mysteries. It includes, but not limited: the so-called "Primordial Lithium Puzzle" [11], the so-called "The Solar Corona Mystery" [12,13], the recent EDGES observations [14], and unexpected annual modulation in x-rays in (2−6) keV energy band observed by XMM -Newton observatory [15], among many others. These cosmological puzzles could be resolved within AQN framework with the same set of physical parameters to be used in the present work to explain the annual modulation observed by DL, without fitting or modifications of any of them.
Our main results are formulated in section III where we estimate the intensity of modulation due to the neutrino flux emitted by dark matter nuggets when they traverse through the earth. This section is separated to 4 different subsections, III A, III B, III C and III D according to 4 different elements of the proposal (1). We use precisely the same set of parameters obtained previously for a different purpose, as reviewed in the Appendix A. We further estimate the neutrino-induced neutron's flux in 2 surrounding rocks which according to the proposal (1) is the source of the observed DL modulation signal.
In section IV we comment on DL arguments [1][2][3][4] suggesting the irrelevance of the induced neutron flux (due to muons and neutrinos). We explicitly show why DL arguments are not applicable for our proposal (1) with its truly genuine DM nature, though manifested indirectly.
In subsection V A we comment on a number of experiments which exclude the DL results, while in subsection V B we make few comments on recent results by COSINE-100 [16][17][18] and ANAIS-112 [19], and recent proposal CYGNO [20] which have been largely motivated by the DL observations. Finally, in subsections V C and V D we suggest few novel tests which could unambiguously support or rule out the proposal (1).

II. AXION QUARK NUGGET DM MODEL
In this section we overview the AQN model. Specifically, in subsection II A we list most important and very generic features of the framework, which do not depend on specific parameters of the model. In subsection II B we overview the observational constraints on a single fundamental parameter of this framework, the average baryon charge of the nugget B which enters all the observables. Finally, in subsection II C we highlight a number of other model-dependent properties, such as ionization features of the AQNs, their survival pattern, their annihilation rate and some other questions which are relevant for the present studies.

A. Generic features of the AQN model
The idea that the dark matter may take the form of composite objects of standard model quarks in a novel phase goes back to quark nuggets [21], strangelets [22], nuclearities [23], see also review [24] with large number of references on the original results. In the models [21][22][23][24] the presence of strange quark stabilizes the quark matter at sufficiently high densities allowing strangelets being formed in the early universe to remain stable over cosmological timescales.
The axion quark nuggets (AQN) model advocated in [25] is conceptually similar, with the nuggets being composed of a high density colour superconducting (CS) phase. As with other high mass dark matter candidates [21][22][23][24] these objects are "cosmologically dark" not through the weakness of their interactions but due to their small cross-section to mass ratio. As a result, the corresponding constraints on this type of dark matter place a lower bound on their mass, rather than coupling constant.
There are several additional elements in AQN model in comparison with the older well-known and well-studied constructions [21][22][23][24]: • First, there is an additional stabilization factor provided by the axion domain walls (with a QCD substructure) which are copiously produced during the QCD transition and which help to alleviate a number of the problems inherent in the older models 1 .
• The core of the AQN is in CS phase, which implies that the two systems (CS and hadronic) can coexist only in the presence of the external pressure which is provided by the axion domain wall. It should be contrasted with original models [21][22][23][24] which must be stable at zero external pressure.
• Another crucial additional element in the proposal is that the nuggets could be made of matter as well as antimatter in this framework.
The direct consequence of this feature is that the dark matter density Ω dark and the baryonic matter density Ω visible will automatically assume the same order of magnitude Ω dark ∼ Ω visible without any fine tunings, and irrespectively to any specific details of the model, such as the axion mass m a or size of the nuggets R. Precisely this fundamental consequence of the model was the main motivation for its construction. The presence of a large amount of antimatter in the form of high density AQNs leads to many observable consequences as a result of possible (but, in general, very rare) annihilation events between antiquarks from AQNs and baryons from visible Universe. We highlight below the basic results, accomplishments and constraints of this model.
Let us recapitulate the original motivation for this model: it is commonly assumed that the Universe began in a symmetric state with zero global baryonic charge and later (through some baryon number violating process, non-equilibrium dynamics, and CP violation effects, realizing three famous Sakharov's criteria) evolved into a state with a net positive baryon number.
As an alternative to this scenario we advocate a model in which "baryogenesis" is actually a charge separation (rather than charge generation) process in which the global baryon number of the universe remains zero at all times. In this model the unobserved antibaryons come to comprise the dark matter in the form of dense nuggets of antiquarks and gluons in CS phase. The result of this "charge separation" process is two populations of AQN 1 In particular, in the original proposal the first order phase transition was the required feature of the construction. However it is known that the QCD transition is a crossover rather than the first order phase transition. It should be contrasted with AQN framework when the first order phase transition is not required as the axion domain wall plays the role of the squeezer. Furthermore, it had been argued that the nuggets [21][22][23][24] are likely to evaporate on the Hubble time-scale even if they were formed. In the AQN framework the fast evaporation arguments do not apply because the vacuum ground state energies in the CS and hadronic phases are drastically different.
carrying positive and negative baryon number. In other words, the AQN 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, the number of nuggets and antinuggets will be different 2 . This difference is always an order of one effect irrespectively to the parameters of the theory, the axion mass m a or the initial misalignment angle θ 0 . We refer to the original papers [43][44][45][46] devoted to the specific questions related to the nugget's formation, generation of the baryon asymmetry, and survival pattern of the nuggets during the evolution in early Universe with its unfriendly environment. The only comment we would like to make here is that the disparity between nuggets Ω N and antinuggets ΩN generated due to CP violation unambiguously implies that the baryon contribution Ω B must be the same order of magnitude as ΩN and Ω N because all these contributions are proportional to one and the same fundamental dimensional parameter Λ QCD . If these processes are not fundamentally related the two components Ω dark and Ω visible could easily assume drastically different scales. This represents a precise mechanism of how the "baryogenesis" can be replaced by the "charge separation" processes in the AQN framework. The remaining antibaryons in the early universe plasma then annihilate away leaving only the baryons whose antimatter counterparts are bound in the excess of antiquark nuggets and are thus unavailable for fast annihilation. All asymmetry effects are order of one which eventually results in similarities for all components, visible and dark, i.e.
as they are both proportional to the same fundamental Λ QCD scale, and they both are originated at the same QCD epoch. In particular, the observed matter to dark matter ratio Ω dark 5 · Ω visible corresponds to a scenario in which the number of antinuggets is larger than the number of nuggets by a factor of (ΩN /Ω N ) ∼ 3/2 at the end of nugget formation. It is important to emphasize that the AQN mechanism is not sensitive to the axion mass m a and it is capable to saturate (2) itself without any other additional contributors. It should be contrasted with conventional axion production mechanisms when the corresponding contribution scales as Ω axion ∼ m −7/6 a . This scaling unambiguously implies that the axion mass must be fine-tuned m a 10 −5 eV to saturate the DM density today while larger axion mass will contribute very little to Ω dark . The relative role between the direct axion contribution Ω axion and the AQN contribution to Ω dark as a function of mass m a has been evaluated in [45], see Fig. 5 in that paper.
Unlike conventional dark matter candidates, such as WIMPs (Weakly interacting Massive Particles) the darkmatter/antimatter nuggets are strongly interacting and macroscopically large. However, they do not contradict any of the many known observational constraints on dark matter or antimatter due to the following main reasons [47]: • They are absolutely stable configurations on cosmological scale as the pressure due to the axion domain wall (with the QCD substructure) is equilibrated by the Fermi pressure. Furthermore, it has been shown that the AQNs survive an unfriendly environment of early Universe, before and after BBN epoch [46]. The majority of the AQNs also survive such violent events as the galaxy formation and star formation; • They carry a huge (anti)baryon charge |B| 10 25 , and so have an extremely tiny number density; • The nuggets have nuclear densities, so their effective cross section σ ∼ R 2 is small relatively to its mass, σ/M ∼ 10 −10 cm 2 /g. This key ratio is well below the typical astrophysical and cosmological limits which are on the order of σ/M < 1 cm 2 /g ; • They have a large binding energy such that the baryon charge locked in the nuggets is not available to participate in big bang nucleosynthesis (BBN) at T ∼ 0.1 MeV, and the basic BBN picture holds, with possible small corrections of order ∼ 10 −10 which may in fact resolve the primordial lithium puzzle [11]; • The nuggets completely decouple from photons as a result of small σ/M ∼ B −1/3 1 ratio, such that conventional picture for structure formation holds.
• The nuggets do not modify conventional CMB analysis, with possible small corrections which, in fact, may resolve a tension [14] between standard prediction and EDGES observation (stronger than anticipated 21 cm absorption features).
To reiterate: the weakness of the visible-dark matter interaction is achieved in this model due to the small geometrical factor σ/M ∼ B −1/3 rather than due to a weak coupling of a new fundamental field to standard model particles. In other words, this small effective interaction ∼ σ/M ∼ B −1/3 replaces a conventional requirement of sufficiently weak interactions of the visible matter with WIMPs.

B. Mass distribution constraints
One should emphasize that the AQN construction does not specify the size of the nuggets. In general, for a given 4 axion mass m a there is always a range of B when a stable solution exists [25]. The average size of the nugget within this stability range scales as R ∝ m −1 a while the baryon charge of the AQN itself scales as B ∝ R 3 ∝ m −3 a . In the AQN framework we treat B as a fundamental parameter to be constrained observationally. It is clear that larger B values produce weaker observational signals as σ/M ∼ B −1/3 . Furthermore, any such consequences assume the largest values where the densities of both visible and dark matter are sufficiently high such as in the core of the galaxy, the early universe, or the stars and planets. In other words, the nuggets behave like conventional cold dark matter in environments where the density of the visible matter is small, while they become interacting and radiating objects (i.e. effectively become visible matter) if they enter an environment of sufficiently large density.
As mentioned above the flux of AQN on the earth's surface is scaled by a factor of B −1 and is thus suppressed for large nuggets, see eq. (8) below. For this reason the experiments most relevant to AQN detection are not the conventional high sensitivity dark matter searches but detectors with the largest possible search area. For example it has been proposed [48] that large scale cosmic ray detectors such as the Auger observatory of Telescope Array may be sensitive to the flux of AQN in an interesting mass range. An obvious challenging problem with such studies is that the conventional cosmic ray detectors are designed to analyze the time delays measured in µs as cosmic rays are assumed to be moving with the speed of light, while AQNs move with non-relativistic velocity v AQN 10 −3 c. It obviously requires the latency time to be measured in ms range in order to study the correlated signals from AQN. The modern cosmic ray detectors are not designed to analyze such long time correlations, and normally such signals would be rejected as a background.
The strongest direct detection limit is set by the Ice-Cube Observatory's non-detection of a non-relativistic magnetic monopole [49]. While the magnetic monopoles and the AQNs interact with material of the detector in very different way, in both cases the interaction leads to the electromagnetic and hadronic cascades along the trajectory of AQN (or magnetic monopole) which must be observed by the detector if such event occurs. A non-observation of any such cascades puts the following limit on the flux of heavy non-relativistic particles passing through the detector, where we assume 100 % efficiency of the observation of the AQNs passing through IceCube Observatory, see Appendix A in ref. [50].
Similar limits are also derivable from the Antarctic Impulsive Transient Antenna (ANITA) [51] though this result depends on the details of radio band emissivity of the AQN. In the same work the author also derives the constraint arising from potential contribution of the AQN annihilation events to earth's energy budget which re-quires |B| > 2.6 × 10 24 [51], which is also consistent with (3). There is also a constraint on the flux of heavy dark matter with mass M < 55g based on the non-detection of etching tracks in ancient mica [52]. It slightly touches the lower bound of the allowed window (3), but does not strongly constraint entire window (5) because the dominant portion of the AQNs lies well above its lower limit (3) assuming the mass distribution (6) as discussed below.
The authors of ref. [53] use the Apollo data to constraint the abundance of quark nuggets in the region of 10 kg to one ton. It has been argued that the contribution of such very heavy nuggets must be at least an order of magnitude less than would saturate the dark matter in the solar neighbourhood [53]. Assuming that the AQNs saturate the dark matter the constraint [53] can be reinterpreted that at least 90% of the AQNs must have the masses below 10 kg. This constraint can be approximately expressed in terms of the baryon charge as follows: Therefore, indirect observational constraints (3) and (4) suggest that if the AQNs exist and saturate the dark matter density today, the dominant portion of them must reside in the window: 3 · 10 24 B 10 28 [constraints from observations]. (5) Completely different and independent observations also suggest that the galactic spectrum contains several excesses of diffuse emission the origin of which is not well established, and remains to be debated. The best known example is the strong galactic 511 keV line. If the nuggets have a baryon number in the B ∼ 10 25 range they could offer a potential explanation for several of these diffuse components. It is very nontrivial consistency check that the required B to explain these excesses of the galactic diffuse emission belongs to the same mass range as reviewed above. For further details see the original works [54][55][56][57][58][59] with explicit computations of the galactic radiation excesses for varies frequencies, including excesses of the diffuse x-and γ-rays. In all these cases photon emission originates from the outer layer of the nuggets known as the electrosphere, and all intensities in different frequency bands are expressed in terms of a single parameter B such that all relatives intensities are unambiguously fixed because they are determined by the Standard Model (SM) physics. 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 [12,13,60].
In this scenario the AQN composed of antiquarks fully annihilate within the so-called transition region (TR) providing a total annihilation energy which is very close to the observed EUV luminosity of 10 27 erg/s. The EUV emission is assumed to be powered by impulsive heating events known as nanoflares conjectured by Parker long ago. The physical origin of the nanoflares remains to be unknown. If our identification of these nanoflares with the AQN annihilating events is correct, we may extract the baryon charge distribution dN/dB for the AQNs because the energy distribution dN/dE of the nanoflares has been previously modelled by the solar physics people to fit the observations, i.e.
where slop parameter α slightly varies for different solar models 3 . It is a highly nontrivial consistency check that the typical nanoflare energy range E nano (10 22 − 6 · 10 25 ) erg corresponds (within AQN framework) to the baryon charge window 3 · 10 24 B 2 · 10 28 which strongly overlaps with all presently available constraints (5) on the AQN sizes as reviewed above. Precisely this highly nontrivial consistency check on size distribution along with (essentially model-independent) computation of the total luminosity of the EUV radiation ∼ 10 27 erg/s from the solar corona, which is consistent with observations, gives us a strong confidence for the plausibility of the identification (6) between the nanoflares conjectured by Parker long ago and the AQN annihilation events.
Encouraged by these consistency checks we adopted the AQN size distribution (6) with parameter α being extracted from the heating corona studies for all our subsequent Monte Carlo simulations including the computations of the AQN flux on Earth as given by (8) to be discussed in great details in next section.

C. Few additional comments
Here we make few additional comments on basic features of the AQNs which are important for understanding of this framework in general and for specific estimates (1) relevant for DAMA/LIBRA signal in particular.
We start this overview with mentioning of the electric charge which AQNs may carry while propagating in a media. It is normally assumed that all types of nuggets including the old version models [21][22][23][24] are neutral at 3 One should comment here that the algebraic behaviour for the distribution (6) is a very generic feature of the percolation theory within AQN framework as recently argued in [46]. However, a numerical estimate of parameter α from the theoretical side requires a deep understanding of the QCD phase transition dynamics at θ = 0, when the axion domain walls are formed. This knowledge is not likely to be available any time soon as the QCD lattice simulations cannot study a system with θ = 0, which represents a well-known sign problem in the lattice community.
zero temperature T = 0 because the electrosphere made of leptons will be always formed even when the quark nuggets themselves are electrically charged (for example, due to the differences in quark's masses). However, the neutrality will be lost due to the ionization at T = 0, in which case the nuggets will esquire the positive charge due to the ionized electrons while the anti-matter nuggets will esquire the negative electric charge due to the ionized positrons. The corresponding charge Q for AQNs can be estimated as follows [11,12]: where n(z) is the density of the positrons in electrosphere which has been computed in the mean field approximation. In this estimate it is assumed that the positrons with p 2 /(2m e ) < T will be stripped off the electrosphere as a result of high temperature T . These loosely bound positrons are localized mostly at outer regions of electrosphere at distances z > z 0 = (2m e T ) −1/2 which motivates our cutoff in estimate (7). Numerically Q ∼ 10 8 represents very tiny portion Q/B in comparison with the baryon charge B ∼ 10 25 even for relatively high temperature T 100 eV in the solar corona. These objects behave in all respects as neutral objects (for example the cosmic magnetic field does not affect the AQN's trajectory) because Q/M e/m p . Nevertheless, a nonvanishing charge Q may play a very important role in some circumstances, such as propagation of AQNs in highly ionized plasma in solar corona [12,13,60]. The non-vanishing charge Q may also suppress the primordial lithium abundance at T 20 keV due to the strong attraction between the negatively charged AQNs and positively charged lithium nuclei [11].
Our next comment is related to the survival pattern of the AQNs. The comprehensive studies on this matter can be found in the original work [46]. The only comment we would like to make here is that the dominant portion of the AQNs will survive the evolution of the Universe. However, very small portion of the AQNs which is gravitationally captured by stars and planets may drastically decrease their baryon charges and may even experience of complete annihilation 4 . The typical size L of the region of the medium with density ρ where the complete annihilation occurs is estimated as σρL Bm p , where σ 6 is an effective cross section which could be much larger than naive πR 2 due to the long-range Coulomb interaction when Q = 0 according to (7). This effect plays a very important role in the solar corona [12,13,60]. All nuggets which are gravitationally captured by the Sun will be completely annihilated.
One may wonder what happens to the axions from AQNs which are now liberated and become propagating axions with average energy E 1.3m a . Small portion of these axions will be converted to photons in the background of the solar magnetic field, and in principle can be observed on Earth. The effect is very small though even if one takes into account the resonance condition due to the plasma effects in the solar corona [61].
In case of the AQN traversing the Earth's interior only some small portion of the baryon charge will be annihilated. The full scale Monte Carlo simulations suggest that on average approximately (10-20)% of the total baryon charge will be lost as a result of traversing of the AQNs through the Earth's interior, see column for ∆B/B in table III in ref. [50]. The average amount of the lost baryon charge depends, of course, on the size distribution and parameter α defined by eq. (6). To avoid confusion, let us emphasize one more time that all these AQNs which get completely annihilated or they lost a finite portion of their baryon charges represent a very tiny portion of all AQNs during entire evolution of the Universe as estimated in [46].
Final comment we would like to make in this subsection is related to the the spallation which represents a very common process when heavy nuclei lose their baryon charge as a result of interaction with a medium. In contrast with conventional nuclei the spallation cannot play any essential role in AQN survival pattern due to several reasons. First of all, the gap ∆ in CS phase is typically in 100 MeV range in contrast with conventional nuclear physics where the binding energy normally is in few MeV range. The most important distinct feature, however, is as follows. When a large amount of energy is injected into a heavy nucleus, the spallation takes place and large number of neutrons may be liberated leading to the decreasing of the baryon charge of the heavy original nucleus. Such process cannot occur with AQN because all particles which get excited due to the energy injection in CS phase are the coloured objects. Therefore, these elementary excitations cannot enter the hadronic vacuum where the normal baryons live and must stay inside of the AQN. Therefore, the AQNs do not suffer from spallation processes as heavy nuclei do. This is direct manifestation of the same feature of the AQN construction already mentioned in footnote 1 which states that the confined hadronic and CS phases have different vacua. The same feature precludes transformation of the entire star/planet into a new phase if the AQNs made of matter stop in the deep interior of the star/planet, see also footnote 4.

III. DL MODULATION BY AQNS
This section is separated to four different subsections, III A, III B, III C and III D according to four different elements of the proposal (1).

A. AQN flux on Earth
We start our task with the 1-st element from proposal (1) by estimating the AQN hit rate per unit area on earth surface assuming that ρ DM is entirely saturated by the nuggets. The relevant rate has been studied previously in [50] for a different problem of computing the axion flux produced by the AQNs. Now we estimate the AQN hitting rate assuming conventional dark matter density ρ DM 0.3GeVcm −3 surrounding the Earth. Assuming the conventional halo model one arrives to the following result [50]: . (8) The averaging over all types of AQN-trajectories with different masses M N m p |B|, with different incident angles and different initial velocities and different size distribution does not modify much this estimate. The result (8) suggests that the AQNs hit the Earth's surface with a frequency approximately once a day per 100 2 km 2 area. The hitting rate for large size objects is suppressed by the distribution function f (B) ∝ B −α as given by (6).
The estimate (8) explicitly shows that conventional DM detectors are too small in size to detect AQNs directly as the corresponding flux is many orders of magnitude smaller than the one due to the conventional WIMPs. However, some modern cosmic ray detectors, such as Pierre Auger observatory, in principle are capable to study small flux of order (8) as suggested in [48]. One should also mention that a smaller size IceCube detector imposes a direct constraint on the average baryon charge of the nugget B ≥ 3 · 10 24 , see Appendix A in ref. [50].
From (8) one can derive the total hit rate for entire earth's surface which is given by [50]: After the nugget hits the surface it continues to propagate by annihilating the material along its path. The trajectory of the AQN is a straight line as only small portion of the momentum (and the baryon charge) will be lost in this journey. The energy produced due to the annihilation events will be isotropically dissipating (in the rest frame of the nugget) along the propagation. The rate (9) includes all types of the AQN's trajectories inside the earth interior: the trajectories when AQNs hit the surface with incident angle close to 0 o (in which case the AQN crosses the earth core and exits from opposite site of earth) as well as the trajectories when AQNs just touch the surface with incident angle close to 90 o (when AQNs leave the system without much annihilation events in deep underground). The result of summation over all these trajectories can be expressed in terms of the average mass (energy) loss ∆m AQN per AQN. The same information can be also expressed in terms of the average baryon charge loss per nugget ∆B as these two are directly related: ∆m AQN ≈ m p ∆B . The corresponding MC simulations with estimates for ∆m AQN have been carried out in [50], see table III in that paper. This information will be very important for our analysis in next section III B as it provides a normalization for the total neutrino flux being emitted by the AQNs when they traverse the earth's interior.

B. Neutrino production from AQNs
The 2-nd element of the proposal (1) requires the estimation of the neutrino intensity due to the AQN annihilation processes. Before we proceed with corresponding estimates we want to make a short detour related to the axion production due to the same AQN annihilation events. We also need to know the basic features of the neutrino spectrum emitted by AQNs. The lessons from that studies can be used for estimations of the neutrino flux from AQNs, which is the main subject of this subsection.

Detour on the axion production
It has been noticed in [62] that the large number of axions will be produced because the axion domain wall 5 will start to shrink during the AQN annihilation events and emit the propagating axions which can be observed. The corresponding spectrum has been computed in [61] where it has been shown that the emitted axions will have a typical velocities v a 0.6c in contrast with conventional galactic axions characterized by small velocities ∼ 10 −3 c such that these two different production mechanisms can be easily discriminated.
The average number of the emitted axions N a as a result of the AQN annihilation events in the earth's interior can be estimated as follows [63]: 5 As mentioned in Sect. II the axion domain wall plays an important role of a squeezer stabilizing the nugget. This energy has been accumulated and stored during the QCD epoch at the moment of formation. The corresponding energy accounts for a considerable portion of the nugget's total energy which is parametrized by κa ≈ 1/3, see below eq. (10) where coefficient κ a determines the relative amount of annihilating energy (per unit baryon charge) transferred to the axion production. The computation of the coefficient κ a 1/3 as well as E a 1.3m a is straightforward exercise [61] as it represents a conventional quantum field theory problem for weakly interacting axion field.
The energy flux of the axions (being averaged over all emission angles and summed over all trajectories) measured on earth surface is estimated as [63]: which has a proper dimensionality [GeV · cm −2 · s −1 ] for the energy flux. The axion flux for this mechanism is estimated as which has a proper dimensionality [cm −2 · s −1 ] for the axion flux.
We are now ready to estimate the neutrino emission from the AQNs, which is the main topic of this subsection. One could follow the same logic of computations as highlighted above, with the only difference is that instead of emission of the axions one should study the neutrinos, which are similar to axions as they can easily propagate through entire earth interior. This is because the relevant cross section with surrounding material is very small in both cases, and formula (11) can be applied for estimation of the neutrino flux on the earth surface with corresponding replacements κ a → κ ν and E a → E ν , while factor Ṅ in eq. (12) of course remains the same.

On neutrino spectrum in CS phase
If we had a conventional hadronic phase inside the nuggets the corresponding computations of neutrino emission would be a very simple exercise. Indeed, we know a typical yield of the pseudo Nambu Goldstone (NG) bosons per annihilation event of a single baryon charge (such aspp). We also know a typical decay pattern for all NG bosons such as K → µν and π → µν with consequent µ decays. We also know with very high precision the branching ratios for the non-leptonic decays of the NG bosons such as K → 3π and η → 3π with consequent decays to neutrinos. It would allow us to compute the total number of neutrinos per single annihilation event. This would also allow us to compute the energy spectrum of neutrinos 6 emitted by the AQNs. 8 Unfortunately, we do not have this luxury of knowing all these key features in the CS phase. Therefore, we cannot predict the spectrum of neutrinos emitted by AQNs. The only solid and robust information which is available today is the typical scale of the lightest NG mass in CS phase, which is normally estimated in (20 − 30)MeV range. This scale for NG masses can be translated to an estimation for a typical neutrino's energy scale in AQN framework. Assuming that one half of the lightest NG mass goes to the neutrino's energy we expect that E ν 15 MeV. While a typical energy scale for the NG mass and the corresponding neutrino's energies in the AQN framework is established with a reasonable accuracy because it is based on well established theory of CS phases [66,67], the computation of the corresponding neutrino spectrum is a much harder problem. The basic problem is of course, the lack of knowledge of CS phase (in contrast with the confined phase where all NG masses and relevant branching ratios are measured with very high accuracy). An additional uncertainty also comes from the lack of understanding of the fermion excitations of the CS phases which may be, in fact, the dominant contributors to the neutrino fluxes. This is because the NG bosons which are produced as a result of annihilation events in CS phase cannot leave the system and consequently decay to emit neutrino, as they must stay inside the nuggets 7 . It should be contrasted with hadronic phase when pions and Kaons (produced as a result of pp annihilation) decay to muons and neutrinos in vacuum. Therefore, the NG bosons in CS phases are likely to be absorbed by fermion excitations (if kinematically allowed), which consequently decay to neutrinos.
Indeed, in unpaired quark matter neutrino emissivity is dominated by the direct Urca processes 8 such as d → u+e − +ν e . In case of antinuggets it should be replaced by anti-quarks with emission of neutrino ν e and positron e + with the energy determined by the energy of the fermion excitation, which itself assumes the energy of order of the low conventional (for confined phase) pattern, in which case a large number of neutrinos will be produced in the (20-50) MeV range, where the sensitivity of underground neutrino detectors such as Super-Kamiokande have their highest signal-to-noise ratio. The basic claim of ref. [65] is that annihilation processes involving an antiquark nugget in CS phase proceed in drastically different way than assumed in [64] when the lightest pseudo Nambu-Goldstone meson has mass in the (20-30) MeV range. As a result of this crucial difference the neutrino's energies will be in 15 MeV range, well below the present day constraints. 7 This is because all excitations in CS phase are the colour objects, and cannot propagate in hadronic vacuum. 8 Another direct Urca process e − + u → d + νe which for antinuggets would correspond to emission of the antineutrino is likely to be strongly suppressed as it requires the presence of the positrons with sufficiently high energy above the Fermi surface. For low temperature the corresponding density for positron excitations is exponentially small ∼ exp(−E e + /T ). This argument suggests that κν κν as eq. (13) states.
NG mass, see Appendix A 3 with more comments on this matter. If this process indeed becomes the dominant mechanism of the neutrino emission from AQNs than one should expect that the number of neutrinos per single event of annihilation should greatly exceed the number of anti-neutrino per single event of annihilation, which we assume to be the case. Formally, this case can be expressed as follows: where the coefficients κ ν and κν describe the number of neutrinos and antineutrinos produced per single annihilation event, similar to parameter κ a 1/3 entering expression (10) and describing the axion emission due to the AQN annihilation events.
We conclude this subsection with the following generic comment. Our system is the strongly coupled gauge theory, the QCD. It should be contrasted with conventional weakly coupled gauge theories when all computations are under complete theoretical control. In our system it is very hard to predict realistic spectra and intensities for neutrino and anti-neutrino fluxes in the 15 MeV energy range due to variety of possible phases, high sensitivity to the parameters and large number of possible decay channels producing neutrinos and antineutrons, as we discussed above. Such an analysis could be coined as "the nuclear physics of CS phases". The complexity and uncertainties of such studies (though it is entirely based on the Standard Model physics) is the main reason to introduce phenomenological parameters κ ν and κν in eq. (13) which will be treated in what follows as unknown parameters and can only be constrained by experiment. However, the basic scale of the problem is fixed by the lowest NG mass in CS phase with a reasonable accuracy, and it is given by eq. (13), i.e. E ν 15 MeV. Precisely this basic neutrino energy scale determines the maximum recoil energy around 6 keV in DL signal.

Neutrino flux from AQNs
In what follows we need the expression for the neutrino flux similar to our formula for the axion flux (12), Similar expression is also valid for antineutrinos: one should replace κ ν → κν and dN ν → dNν in (14). Now we are ready for the numerical estimates. We use ∆m AQN and Ṅ from subsection III A to arrive to the following estimate for the neutrino flux in terms of unknown parameter κ ν dN ν dtdA 0.6 · 10 6 · κ ν · ∆B B and similar expression for antineutrinos obtained from (15) by replacing κ ν → κν and dN ν → dNν. In formula (15) the dimensionless ratio ∆B / B counts the relative portion of baryon charge being annihilated in the interior while AQNs traversing the earth. This parameter depends on the nugget's size distribution as reviewed in Sect.II B. Numerically, it is close to 10% for models with large average charge B 10 26 and around 30% for models with smaller average charge B 10 25 [50]. In what follows to simply things we want to ignore all these numerical factors, and represent the total neutrino and antineutrino fluxes produced by AQN mechanism over entire energy range 0 E ν,ν 15 MeV as follows It is instructive to quote few known constraints on neutrino and antineutrino fluxes in this energy band E ν,ν ∼ 15 MeV in order to compare them with the fluxes produced by the AQN mechanism as expressed by eq.(16). The largest flux relevant for this energy band comes from solar 8 B which is around Φ νe 5 · 10 6 (cm −2 s −1 ), while solar "hep" component is close to Φ νe 8 · 10 3 (cm −2 s −1 ) [68]. The atmospheric and diffuse supernova neutrino backgrounds are at least two orders of magnitude smaller than "hep" component, and can be safely ignored in our discussions.
To conclude this subsection we would like to emphasize that the neutrino flux (16) generated by AQNs could be the same order of magnitude as the dominant 8 B solar contributor Φ νe 5 · 10 6 (cm −2 s −1 ) in this energy band. It is not presently ruled out by any experiment. The key point here is that the AQN-induced neutrino flux (16) is the subject of the annual modulation as it has inherent DM origin. Therefore, it could be the source of the observed DL modulation signal.

C. Neutrino-induced Neutrons
The 3-rd element of proposal (1) is the liberation of the neutrons from surrounding rocks due to the coupling with neutrinos. The intensity of these liberated neutrons is the subject of annual modulation as the corresponding neutron intensity is proportional to the neutrino flux (16), which itself is directly proportional to the DM flux in form of the AQNs according to (9).
The idea that the surrounding neutrons might be the origin of the DL modulations has been suggested previ-ously in a number of papers, see [7][8][9][10]. This proposal, of course, remains a subject of debates as the DL collaboration rejects the idea that surrounding neutrons could play any essential role in their observations [2]. We will make some comments within AQN scenario in next section IV to address this question.
In the present subsection we elaborate on the 3-rd element of our proposal (1) by estimating the neutron flux which is induced by neutrinos. The crucial difference with previous proposals [7][8][9][10] is that the neutrinoinduced neutron flux (17) manifests itself with proper annual DM modulation [5,6] characterized by the expected phase t 0 0.4yr corresponding approximately to June 1 for the Standard Halo Model (SHM).
We use conventional expressions from [2,10] for cross section σ and the number density n 10 29 m −3 of the target to make straightforward comparison with the previous estimates: where we used the AQN induced neutrino flux (16) and cross-section σ 10 −41 cm 2 for the neutrino-induced neutron spallation for 208 Pb target. The effective volume V entering eq. (17) will be discussed later in the text. The AQN-induced neutron rate production per unit volume can be represented as follows This rate is approximately one order of magnitude lower if κ ν 1 in comparison with corresponding estimates adopted in [2,10] for the neutron's rate induced by the solar neutrinos. It could be the same order of magnitude as used in [2,10] if κ ν 10, see estimate (13). We return to significance of the estimate (18) later in the text. The only comment we would like to make here is that the typical kinetic energy of the neutrons liberated by this mechanism will be in the 10 2 keV range, see estimate (19) below.
Indeed, the first thing to notice is that the typical neutrino energy (13) is well above the neutron emission threshold E ν > 7.37 MeV for 208 Pb target such that neutron spallation is kinematically allowed. Furthermore, if one assumes that momentum resulting from the neutron spallation is mostly transferred to the liberated neutron, such that p n (p ν − p ν ), the kinetic energy of the neutron can be estimated as follows One should also add that there is a sharp cutoff of order ∼ 10 2 keV for the neutron's energy (19) produced by this mechanism. It is determined by the neutrino energy (13) in 15 MeV range, which itself is kinematically bound from above by corresponding NG masses in CS phase as reviewed in Appendix A. The presence of such sharp cutoff will be important element in our arguments supporting the proposal (1). What do the neutrons, characterized by the rate (18) and energy (19), do if they enter the DL detector? This is the subject of the next subsection.

D. DL modulation (1) due to the neutrons
The 4-th element of proposal (1) represents the most controversial portion of our analysis. We shall try to argue that neutrons characterized by the flux (18) and energy (19) may serve as the source of the observed DL annual modulation. The corresponding computations are very hard to carry out as they are inevitably based on nuclear physics of large number of very complicated systems. Fortunately, there are many specific experiments and tests which can be in principle performed, to be discussed in sections IV and V. These tests can support or rule out the proposal (1).
The starting point is the standard formula for energy transfer ∆E as a result of elastic 2 → 2 scattering when a target of mass m 2 at rest, struck by a particle with mass m 1 with kinetic energy E n : The entire section of ref. [7] was devoted for explanation of numerous uses and misuses of this formula in different circumstances. We agree with most of the comments, careful explanation of misconceptions, detail analysis of neutron-nuclear interaction given by Ralston in [7]. We defer our specific comments within present context until next section IV. Now we want to make numerical estimate for the recoil energy using expression (20) and identifying m 1 with the neutron characterized by the kinetic energy (19), while m 2 23m 1 with the lightest 11 Na nucleon from DL detector consisting 25 radio-pure NaI crystal scintillators: The significance of this estimate is hard to overstate as it unambiguously shows that the recoil energy cannot exceed the value (21) which is amazingly close to 6 keV cutoff observed by DL. Furthermore, the scale (19) was not invented for the specific proposal (1), in contrast with many different WIMP-based suggestions to fit the observed DL modulations. Rather, this scale is entirely determined by the cutoff in neutrino energy (13) which itself unambiguously fixed by typical NG masses in CS phase. All these scales have been known for quite sometime, as reviewed in Appendix A.
The key element of our proposal (1) is that the flux of the induced neutrons (18) with kinetic energy (19), which eventually generates the signal in DL detector with recoil energy (21), is the subject of annual modulation because all the intensities are proportional to DM velocity v AQN as eq. (9) states. Therefore, to describe the corresponding modulations one can use conventional formula [5,6]: where V = 220 km s −1 is the Sun's oribital speed around the galactic center, V ⊕ = 29.8 km s −1 is the Earth's orbital speed around the Sun, ω = 2π yr −1 is the angular frequency of the annual modulation, and |b| ≤ 1 is the geometrical factor associated with the direction of v AQN relative to the orbital plane of Earth, t 0 0.4yr corresponding to June 1. Hence, it is natural to expect the modulation must be of order O(V ⊕ /V ) ∼ 10%, as incoming flux of the AQN particle explicitly proportional to v AQN (t) according to equation (9). The corresponding value Ṅ (t) enters the expression for the AQN-induced neutrino flux (14), eventually generating the neutrons (18).
The very hard and challenging question remains to be answered if this neutron intensity (17) is sufficient for explanation of the observed DL modulation. To answer this question one has to understand the numerical value for the effective volume V entering formula (17). As we already mentioned this is very complicated problem of nuclear physics as emphasized and nicely explained in [7].
Therefore, instead of theoretical speculations about the value for effective volume entering formula (17) we reverse the problem and estimate the volume V which would match the DL modulation. We suggest several experiments how this proposal (1) and large rate (17) with effective volume V can be tested in next two sections IV and V.
The DL modulation amplitude in terms of (counts per day) cpd [4] reads: DL modulation rate = (0.0103 ± 0.0008) cpd kg keV . (23) This rate must be multiplied to 4 keV for (2-6) keV energy range and 250 kg to get total modulation rate DL total modulation 10 counts day .
On other hand, assuming that 10% of neutrons (18) is the subject of annual modulation with proper phase t 0 as explained in (22)  Few comments are in order. First of all, parameter κ ν was introduced as a number of neutrinos being produced as a result of annihilation of a single baryon charge, see (13). It could be as small as κ ν ∼ 1, but it could be as large as κ ν ∼ 10 being consistent with presently available constraints as mentioned at the end of section III B. The basic reason for this uncertainty is that our system is the strongly coupled gauge theory, to be contrasted with conventional weakly coupled gauge theories when all computations are under complete theoretical control. In our system it is very hard to predict realistic spectra and intensities for neutrino and anti-neutrino fluxes in the 15 MeV energy range due to variety of possible phases, high sensitivity to the parameters and large number of possible decay channels producing neutrinos and antineutrons, as we discussed above. This deficiency in our computational power should not be treated as the weakness of the proposal. Instead, it should be considered as a consequence of the complexity of the system.
The observed rate (24) matches the AQN induced modulation (25) if κ ν 10 and L ∼ 10 m. This required length L ∼ 10 is definitely much greater than the neutron's mean free path λ 2.6 m which was extracted from the studies on the muon-induced background [72] and adopted by refs. [2,10] in the context of the present work on DL modulations. In next section IV we comment on the consistency (25) with DL observations, while in section V we make few comments on relation to other experiments. In the same section V we also suggest possible tests (such as the measuring of the spatial directions of moving neutrons along with time modulation) which could support or rule out the proposal (1).

IV. COMMENTS ON DL ARGUMENTS.
The DL collaboration, of course, discussed the possibility that their signal is associated with neutron flux (induced by muons or neutrinos or both). In fact, entire paper [2] was devoted to analysis of a possible role of neutrino-induced and muon-induced neutrons. These possibilities were discarded in [1][2][3][4] based on the following arguments: 1. Modulation phase arguments; 2. Energy range arguments; 3. Intensity arguments. We want to make few comments on each of the items from this list. We start with the modulation arguments, while the energy and intensity arguments will follow.
1. Quote from [2]: " ... It is worth noting that neutrons, muons and solar neutrinos are not competing background when DM annual modulation signature is investigated since in no case they can mimic this signature...". In the proposed scenario (1) this argument obviously does not apply because both the neutrino flux (16) and the neutron flux (17) with typical energy (19) are automatically the subject of the annual modulation (22) with proper phase t 0 0.4yr corresponding to June 1. This is because the source of the modulation in this framework has truly genuine DM origin represented by AQNs.
2. The DL has carried out the comprehensive studies on dependence of the annual modulation as a function of the energy interval. The claim is that the modulation is not observed above the energy ∼ 6 keV. In particular the modulation amplitude for the energy above 6 keV for the whole data sets (DAMA/NaI, DAMA/LIBRA-phase-1, DAMA/LIBRA-phase-2) is shown to be consistent with zero, see Fig 11 in ref. [4]. This property is indeed very hard to understand in terms of the conventional physics advocated in [7][8][9][10].
At the same time this unique feature of the system (characterized by a sharp cutoff at ∼ 6 keV) automatically emerges in our framework. Indeed, the neutrino flux (16) with typical neutrino energy E ν ∼ 15 MeV is determined in our framework by the NG masses in the CS phase, see (A3). The corresponding neutrino-induced neutron flux (17) is characterized by the typical energy (19) formulated in terms of E ν . The sharp cutoff for the recoil energy (21) in this framework (which falls into the proper ∼ 6 keV energy range) is determined by (19) which essentially determined by the NG masses in CS phase as reviewed in Appendix A. One should emphasize that all these energy scales have not been specifically invented in this work to explain the observed DL modulations with (1-6) keV energy; rather the relevant energy scales have been established long ago in unrelated studies for different purposes in a different context.
3. The neutrino flux (16) originated from AQNs in this framework is lower than the background solar neutrino flux Φ νe 5·10 6 cm −2 s −1 for this energy band (originating from solar 8 B) at least by factor of 5 for κ ν 10, as mentioned at the end of section III B. The corresponding neutrino-induced neutron flux (17) is also must be lower in comparison with intensity of neutrons induced by the solar neutrinos. However, the key point is that this subdominant neutron component is originated from AQNs, and therefore, is the subject of conventional DM annual modulation (22) with proper phase t 0 .
Is the corresponding neutrino-induced neutron's intensity is sufficient 9 to explain DL modulation? The DL argued that the answer is "No" [1][2][3][4]. However, the DL arguments on the neutron's intensity were challenged in [7]. We have nothing new to add to these extensive discussions on possible role of neutrino-induced neutrons. Instead of speculations about this very complex nuclear physics system with complicated resonance structure, we suggest to test the proposal (1) by measuring the modulation, intensity and the directionality in coordinate space of the neutrons which are responsible for the recoil (21).
The subdominant flux of the AQN-induced neutrons can be, in principle, discriminated from the dominant components, including the solar neutrino-induced neutrons if the direction of the neutron's momentum and the modulation are measured. This is because the solar neutrinos are propagating from a single direction in the sky, while AQN-induced neutrinos (which have truly genuine DM origin) are randomly distributed in space. This 12 topic on possible tests of the proposal (1) represents the subject of the next section V.

V. RELATION TO OTHER EXPERIMENTS. POSSIBLE FUTURE TESTS.
This section is separated to three different subsections. First, in subsection V A we make few comments on the previous experiments which exclude the DL signal if interpreted in terms of WIMP-nuclei interactions. We continue with more recent analysis in subsection V B where we make few comments on some recent experiments designed to reproduce (or rule out) the DL signals. Finally, in subsections V C and V D we offer few novel possible tests which can support or refute the proposal (1) explaining the DL signal in terms of the AQN-induced neutrons. In particular, in subsection V D

A. Previous experiments
We start with few comments on the experiments which exclude the DL results. The corresponding collaborations include but not limited to: CDEX [73], CDMS-II [74,75], EDELWEISS-II [76], LUX [77], SuperCDMS [78], XENON10 [79] and XENON100 [80], CoGeNT [81]. The main claim of these collaborations can be formulated as follows: if DL modulation is interpreted in terms of WIMP-nuclei interactions with given σ and given m W IM P then DL signal is excluded with very high level of confidence [82].
From the perspective of the proposal (1) it could be a number of reasons of why DL observes the signal while other collaborations do not. First of all, most of the collaborations (with few exceptions such as CoGeNT [81] and CDMS-II [75]) did not carry out some dedicated studies on the time modulation, which was the crucial ingredient in DL arguments. From the AQN perspective the time modulation is the key element when the subdominant neutron flux can manifest itself if proper time modulation studies are performed.
Another reason (for negative results) could be related to different neutron shields used by different collaborations. We refer to paper by Ralston [7] where the subject on complex behaviour of neutrons in complicated environment is nicely presented. This analysis obviously shows that even minor differences in neutron shields may have dramatic effects and drastically change the impact of neutron's background.
Yet, one more reason, probably the most important one in context of the present work, is as follows. The AQN-induced neutrons with energy (19) are scattering off Na in DL experiment generating recoil energies which fall into the (2-6) keV bin according to (21). The recoil energies for the heavier targets such as xenon or germanium in different experiments could be below threshold because the energies of the time modulated neutrons (19) are bound from above with a cutoff being determined by energies of the AQN-induced neutrinos. The energies of these neutrinos are also bound from above and cannot exceed (13) as they are determined by NG masses in CS phase. As we already mentioned previously, all these energy scales in proposal (1) have not been invented to fit the DL signals. Rather the relevant energy scales have been established long ago in unrelated studies for different purposes in a different context.

B. Recent activities
Now we want to make few comments on recent dedicated experiments which were specifically designed to test DL annual modulation signal. It includes COSINE-100 [16][17][18], ANAIS-112 [19] and CYGNO [20] experiments. We also want to mention other type of experiments [83][84][85] which were not originally designed to test DL annual modulation signal. However, their capabilities to measure the directionality could play a decisive role in detecting of the DM signals. To be more specific, we choose to mention these experiments due to the following reasons: The COSINE-100 experiment was mostly motivated by DL annual modulation. The aim of the collaboration is to reproduce (or refute) the signal and to search for possible origin for the modulation, if observed. The COSINE-100 collaboration uses the same target medium (sodium iodide) which is the crucial element in the context of the present proposal (1) as recoil energies fall into the (2-6) keV bin in our scenario according to (21). Presently the COSINE-100 data is consistent with both a null hypothesis and DL (2-6) keV best fit value with 68% confidence level. More data are obviously needed. It is important that COSINE-100 is planning to measure the neutron's intensity and neutron's modulation [17], in which case the COSINE-100 would know if the possible modulation is due to the neutrons 10 . It is obviously the key element of the proposal (1) based on the subdominant AQN-induced component of neutrons which, however, manifests itself by annual modulation.
The ANAIS-112 collaboration also uses the same target medium as DL and the COSINE-100. The ANAIS-112 has recently published the first results on annual modulation [19]. Their best fits are incompatible at 2.5σ with DL signal. The goal is to reach the sensitivity at 3σ level in five years. As the ANAIS-112 collaboration uses the same target material our comments from the previous paragraph in context of the present proposal (1) also apply to ANAIS-112 experiment especially as ANAIS-112 and COSINE-100 agreed to combine their data.

13
The CYGNO proposal [20] is different from the COSINE-100 and the ANAIS-112 experiments due to the capability to measure the directionality which is the key element of the CYGNO [20] proposal. It is important that it will be located at the same site (LNGS) where DL is located. Therefore, the neutron flux must be the same, including the subdominant AQN-induced component (25) which is the subject of annual modulation. The collaboration is planning to measure (initially) the neutron flux and its modulation without neutron shielding 11 . Such measurements may play a key role in supporting (or ruling out) the proposal (1) because the AQN-induced neutrons are responsible for the recoil (21). The collaboration is also planning in future to reach the neutrino floor by measuring the neutrinos and their directions. In particular the CYGNO could discriminate the neutrinos from the sun by identifying their directions. Furthermore, the CYGNO instrument will be capable to determine nuclear recoil directions, which would allow the collaboration to discriminate WIMP-like DM from AQN-induced events (1). Such measurements, if successful, would obviously play an important role in supporting (or ruling out) the proposal (1). In addition to that, the dominant solar neutrinos in the energy range E ν 12 MeV could be discriminated from subdominant AQNinduced neutrinos (A3). Furthermore, neutrinos in the energy band E ν ≥ 12 MeV cannot be originated from the Sun at all as a result of 8 B threshold. The discovery of such neutrinos and measuring of their annual modulation with intensity in the range (16) would be enormous support for the proposal (1) as the flux of the atmospheric neutrino is at least three orders of magnitude lower than (16).
We also want to mention several other experiments [83][84][85] which are capable to measure the directionality. The idea is to use the carbon nanotube arrays or the graphene layers to measure the directionality of the DM signals. As we mentioned above, the measurements of the directionality could play a decisive role in detecting of the DM signal.

C. Possible future tests
We already mentioned in previous subsection V B few possible tests which can support or rule out the proposal (1) with existing or planning experiments: COSINE-100 [16], ANAIS-112 [19], CYGNO [20]. In this subsection we want to mention a specific for the AQN framework phenomenon when the intensity of the AQN-induced neutrinos may be amplified by very large factor (up to 10 4 ) which greatly increases the chance for discovery such AQN-induced neutrons which always accompany the neutrinos according to section III C. Therefore, this neutrino amplification factor will be obviously accompanied by the corresponding amplification of the neutrino-induced neutron flux (17) and (18).
The idea was originally formulated for the axions in [63], and the effect was coined as the "local flash". The computations can be easily generalized for the neutrinos in straightforward way, see Appendix B with some technical details. It can be explained as follows.
If the AQN hits the surface at distance d R ⊕ from the detector the short-lasting flash occurs with amplification factor A ν (d) measuring the relative short lasting spark of the neutrino flux with respect to the neutrino flux (16) computed by averaging over entire earth surface over long period of time. The amplification A ν (d) is highly sensitive to distance d and can be approximated as follows, see (B4) for derivation: One should note that the correction to the neutrino flux (B3) due to the traversing of a nearby AQN depends on unknown parameter κ ν . However, the relative amplification A ν (d) with respect to the averaged neutrino flux (16) does not depend on κ ν as eq. (26) states. In formula (26) Ṅ is determined by eq. (9), while ∆t 2R ⊕ / v AQN is the time for the AQN to cross the earth averaged over entire ensemble of AQN's trajectories traversing the earth.
As one can see from (26) a huge amplification may indeed occur for d R ⊕ . However, the probability for such event to happen is very tiny, and can be estimated as Event rate 1 see Appendix B with details. The "local flash" lasts for a short period of time which can be estimated as follows We summarize in Table I few choices of time duration τ and the event rate as a function of amplification factor A ν . In particular, it would be a daily short lasting "flash" when the intensity of the subdominant AQNinduced neutrino component (16) is amplified by factor A ν (d) 10 2 such that it becomes the dominant one for a short period of time lasting for about 1 second. Important lesson to be learnt from these estimates is as follows. A subdominant neutrino flux induced by AQNs may become a dominant portion of the neutrinos, overpassing the solar neutrino flux in this energy band for very short period of time. Needless to say, that this AQNinduced neutrino flux (16) and the corresponding neutron flux (17) are also the subject of annual (22) and daily

D. Search for correlations with other AQN-induced phenomena
The flux (8) suggests that the AQNs hit the Earth's surface with a frequency approximately once a day per 100 2 km 2 area, which is precisely the source for a shortlasting amplification in the neutrino production discussed in previous subsection V C. It is important to emphasize that such events occur along with other processes which always accompany the AQNs propagation through the atmosphere and the Earth's underground. Therefore, there will be always some correlations between the amplifications in the neutrino flux and other associated phenomena in the vicinity of the area where AQN event occurs.
For example, if the DM detector (sensitive to neutrinoinduced neutrons such as DL) and an axion search detector are localized close to each other on a distance d νa ∼ d there will be the axion signal due to the amplification A a and the neutrino-induced neutron signal due to the amplification A ν . These signals must be correlated in form of two almost simultaneous short lasting sparks between these two signals in two different detectors. The observation of these cross-correlated signals (collected during a long period of time and by measuring the directionality in DM detector to discriminate the background) would unambiguously support the proposal (1) on the nature of the observed DL modulation signal. Similar crosscorrelation between different synchronized axion stations from a Global Network as suggested in [86] and a nearby neutrino detector would also strongly support the proposal (1).
In particular, the position of the Center for Axion and Precision Physics Research located at Daejeon and COSINE -100 detector located at the Yangyang Underground Laboratory in South Korea obviously satisfy the criterial d νa 0.1R ⊕ when strongly correlated amplifications for A ν and A a may occur in both detectors almost simultaneously with time delay of order of few seconds.
Another correlation which is worthwhile to study can be explained as follows. It has been recently argued that the AQN propagating in the Earth's atmosphere and underground emits the infrasound and weak seismic waves [87]. In fact one such event, according to [87], occurred on July 31-st 2008. It was properly recorded by the dedicated Elginfield Infrasound Array (ELFO) near London, Ontario, Canada and corresponds to relatively large nugget 13 with B 10 27 if interpreted as AQN event [87]. The infrasound detection was accompanied by non-observation of any meteors by an all-sky camera network. The impulses were also observed seismically as ground coupled acoustic waves around South Western Ontario and Northern Michigan. The estimates [87] for the infrasonic frequency ν 5 Hz and overpressure δp ∼ 0.3 Pa are consistent with the ELFO record. It has been also proposed in [87] a detection strategy for a systematic study to search for such events originating from much smaller and much more common AQNs with typical B 10 25 by using Distributed Acoustic Sensing (DAS).
Our original remark here is that the amplification in the neutrino flux (and corresponding enhancement in the neutrino-induced neutrons) as discussed in previous subsection V C must be accompanied by infrasonic and weak seismic waves which can be studied by the DAS instruments implemented in networks of optical-fiber telecommunication cables. The observation for such correlations, if successful, is absolutely unique to AQN framework and would obviously play an important role in supporting of the proposal.
In particular, the DL (and future CYGNO detector) is located at the LNGS site with large number of seismic detectors located in the same area. Proposal is to search for the correlations between enhanced flux of neutrons and weak seismic and infrasound events with delay measured in fraction of a second depending on a precise localization of the seismic detectors. The measurements of the directionality by CYGNO would be the key element in establishing such a correlation. 13 The frequency of appearance for such large nuggets is very tiny according to (6). This explains why such events occur approximtaely once every 10 years instead of observations them once a day which would be the case for much more common nuggets with B ∼ 10 25 .

VI. CONCLUDING COMMENTS
The main results of our work can be summarized as follows: 1.We argued that the annual modulation observed by DL might be explained as a result of the AQN-induced neutrons through the chain (1). In this framework the annual modulation has truly genuine DM origin, though it is manifested indirectly.
2.We also argued that the recoil energy must have a sharp cutoff at ∼ 6 keV consistent with the observed DL signal.
3. We proposed specific tests which can support or rule our the proposal (1), see subsection V C.
4. We proposed to study a specific correlation between the amplification in the AQN-induced neutron flux and the impulses of the infrasound and weak seismic waves, which if found, would strongly support our proposal, see subsection V D.
Why should we consider this AQN model seriously? There is a number of reasons. Originally, this model was invented to explain the observed relation Ω DM ∼ Ω visible and the baryon asymmetry of the Universe as two sides of the same coin, when the baryogenesis framework is replaced by a "charge separation" framework, as reviewed in Section II. After many years since its original formulation this model remains to be consistent with all available cosmological, astrophysical, satellite and ground based constraints, where AQNs could leave a detectable electromagnetic signature. Furthermore, it is shown that the AQNs can be formed and can survive the unfriendly environment during the evolution of the early Universe, such that they entitled to serve as the DM candidates. Finally, the same AQN framework may also explain a number of other (naively unrelated) observed phenomena such as excess of the galactic diffuse emission in different frequency bands, the so-called "Primordial Lithium Puzzle" and "The Solar Corona Mystery", the seasonal variations observed by XMM-Newton observatory, to name just a few, see Section I for the references.
We want to emphasize that all these cosmological puzzles mentioned in Section I could be resolved within the AQN framework with the same set of physical parameters being used in the present work on explanation of DL modulation signal.
The observation of the subdominant AQN-induced neutrons by measuring the directionality and modulation as discussed in Section V B would be a direct manifestation of the AQN dark matter model. The observation of variety of amplifications as discussed in Section V C would be also a strong support for the proposal (1). Finally, the recording of the correlations between the AQN induced neutrino amplification and the impulses of the infrasound and weak seismic waves would strongly support our proposal as argued in subsection V D. We finish this work on this positive and optimistic note.

Nambu-Goldstone modes in CS phase
There are many possible CS phases due to the generation of a gap ∆ through different channels with slightly different properties. While the relevant physics is a part of the standard model, QCD with no free parameters, the corresponding phase diagram is still a matter of debate as it strongly depends on the precise numerical value of the gap ∆, see review articles [66,67]. For our purposes though the key characteristics are very much the same for all phases. Therefore, we limit ourself below to reviewing the most developed, the so called CFL (colour flavour locking) phase. The spontaneous breaking of chiral symmetry in colour-superconductors gives rise to lowenergy pseudo-Nambu-Goldstone (NG) modes with similar quantum numbers to the mesons (pions, kaons, etc.). These objects, however, are collective excitations of the CS state rather than vacuum excitations as is the case for the conventional confined hadronic phase. The finite quark masses explicitly break chiral symmetry, giving rise to these "pseudo"-Nambu-Goldstone modes on the order of 20 MeV, in huge contrast with the hadronic confined phase where the lightest mass meson has m π 140 MeV.
To be more precise, we consider large µ limit for which the masses and other relevant parameters in the CFL phase can be explicitly computed [66,67]: As one can see from (A1) the NG bosons are much lighter than in vacuum. This is because their masses are proportional to m 2 q rather than to m q , as at zero chemical potential. As a result, the lightest NG meson, the kaon, has a mass in the range of 10 to 20 MeV depending on precise value of ∆ and µ MeV, see for example [67].
Another important difference between the NG modes in dense matter and in vacuum is in the dispersion relations for the NG modes which assume the following form, see for example [66]: such that the rest energy of the lightest NG mesons does not exceed the 10-20 MeV range. In fact E K 0 may even vanish, in which case the K 0 field forms a condensate (the so-called CF L K 0 -phase). In these formula v N G deviates from speed of light c due to the explicit violation of the Lorentz invariance in the system such that the dispersion relations for all quasiparticles are drastically different from their counterparts in vacuum.
One should comment here that the dispersion relations for the NG modes within the anti-nuggets (which is most relevant for our purposes) can be obtained from (A2) by replacing µ → −µ such that the lightest NG states become theK 0 and K − for nuggets made of antimatter. This comment is important for identification of the neutrino and anti-neutrino spectra to be discussed in next subsection A 2.

Neutrino emission from NG bosons
The neutrino emission from CFL phase quark matter has been studied previously in a number of papers mostly in context of the physics of neutron stars, see the original papers [88][89][90] and the review article [66].
In this subsection we will overview the basic results of [65] on neutrino and antineutrino fluxes from the sun. The main goal of ref. [65] was to argue that the Super-Kamiokande stringent constraint on anti-neutrino flux Φν e < (1.4 − 1.9)cm −2 s −1 at large energies E > 19.3 MeV [69] (which played the key role in analysis of [64]) does not apply to our case on neutrino and antineutrino production by AQNs because the typical energies of neutrinos and antineutrinos will be much lower.
Indeed, the muons cannot be produced at all in the CFL phase for purely kinematical reasons. Therefore, the energetic antineutrinos which are normally produced in the µ ± → e ± ν e ν µ decay channels are not produced in the CS matter. This is the crucial point of the arguments presented in [65].
In the simplest possible scenario (which was adopted in ref. [65]) the majority of neutrinos will be emitted by the lightest NG bosons, in which case the energy of the emitted neutrinos does not normally exceed 15 MeV for the CFL phase because the lightest NG bosons do not normally exceed mass in 30 MeV range as mentioned in previous subsection A 1. This is very basic and very generic feature of the CS phase. We postpone the important discussions on basic features of the neutrinos emitted by the quarks in CS phases to subsection A 3. Below we list few features on the neutrino spectrum if it is saturated exclusively by NG bosons: 1. A specific choice of µ and ∆ determines the basic mass scales for the light NG modes, which eventually determine the ν spectrum; 2. It is normally assumed that a required value for µ for CFL to be realized is µ 400 MeV. The corresponding value for ∆ is estimated in this case as ∆ 100 MeV in the CFL phase, see review [66]. Such large value of µ can be indeed reached in the AQN's core as recent numerical studies suggest, see Fig. 2 in ref [46]; 3. The neutrino spectrum is qualitatively different from that of the antineutrinos because the annihilation occurs not in vacuum but in dense CS state with µ 400 MeV.
4. Larger chemical potentials µ generally lead to even lighter NG modes while the masses increase with the size of the gap ∆ as eqs. (A1) and (A2) state; We conclude this subsection with the following generic comment. The main goal of [65] was to argue (in simplified setting) that the neutrino's energies in the AQN framework are typically in the 15 MeV range (in contrast with 20-50 MeV range as assumed in [64]) such that the stringent constraint from SuperK [69] does not apply.
3. Fermion excitations and neutrino's production in CS phases.
As we already mentioned the neutrino emission from CS phases quark matter have been studied previously in a number of papers mostly in context of the physics of neutron stars, see the review article [66]. We cannot literally apply the machinery developed in previous studies because all excitations (such as NG bosons) in our case are produced not in a thermally equilibrium system when the density of of the excitations is unambiguously determined by the temperature. Instead, all excitations in our scenario are produced as a result of rare annihilation events. These annihilation events excite the NG bosons as well as fermion excitations, which may decay by emitting neutrinos. The coefficients κ ν and κν entering (16) precisely correspond to this mechanism of the neutrino production.
In this subsection we want to overview some features of the fermion excitations of the CS phases which may be the dominant contributors to the neutrino fluxes. This is because the NG bosons which are produced as a result of annihilation events in CS phase cannot leave the system and consequently decay to emit neutrino, as they must stay inside of the nuggets. It should be contrasted with hadronic phase when pions and Kaons (produced as a result e.g. pp annihilation) decay to muons and neutrinos in vacuum. Therefore, the NG bosons in CS phases are likely to be absorbed by fermion excitations (if kinematically allowed), which consequently decay to neutrinos. We start our overview by mentioning some CS phases which support low energy fermion excitations. Detail discussions can be found in review article [66]. First of all, the so-called 2SC phase (when two out of three colours and flavours are paired and condensed) supports unpaired modes which could be light and couple to NG bosons. Another phase which also supports the light fermion excitations is the so-called CF L−K 0 phase when K 0 energy vanishes according (A2). For antinuggets this phase corresponds toK 0 condensation. In both cases (CFL and CFL-K 0 ) the gap of the excitations decreases with increases of p F such that the gaps for p and n become lower. Indeed, ∆ p,n = ∆ − m 2 s 2p F such that p, n fermion excitations could become completely ungapped. If this happens these modes may become the dominant producers of the neutrinos.
Indeed, in unpaired quark matter neutrino emissivity is dominated by the direct Urca processes such as d → u+e − +ν e . In case of antinuggets it should be replaced by anti-quarks with emission of neutrino ν e and positron e + with the energy determined by the energy of the fermion excitation, which itself assumes the energy of order of the lightest NG mode according to analysis of the previous section A 2.
If this process indeed becomes the dominant mechanism of the neutrino emission from AQNs than one should expect that κ ν κν, E ν 15 MeV, κ ν 1, which is assumed to be the case as Eq. (13) states.

Appendix B: Local flashes
In this Appendix we generalize our axion studies [63] to include the neutrinos into consideration, similar to what we have done in Section III B. The main topic for the present studies is the enhancement effect and great amplification of the axion density which was coined as a "local flash" in [63]. It occurs on rare occasions when an AQN hits (or exits) the Earth surface in vicinity of an axion search detector. In this Appendix we generalize the arguments of Section III B to estimate a similar local flash for neutrinos.
We follow the same logic of ref [63] and consider a case when an AQN is moving in a distance d close enough to the detector, as shown in Fig. 1