Gravitationally bound axions and how one can discover them

As recently advocated in \cite{Fischer:2018niu}, there is a fundamentally new mechanism for the axion production in the Sun and Earth. However, the role of very slow axions in previous studies were neglected because of its negligible contribution to the total axion production by this new mechanism. In the present work we specifically focus on analysis of the non-relativistic axions which will be trapped by the Sun and Earth due to the gravitational forces. The corresponding emission rate of these low energy axions (below the escape velocity) is very tiny. However, these axions will be accumulated by the Sun and Earth during their life-times, i.e. 4.5 billion of years, which greatly enhances the discovery potential. The computations are based on the so-called Axion Quark Nugget (AQN) Dark Matter Model. This model was originally invented as a natural explanation of the observed ratio $\Omega_{\rm dark} \sim \Omega_{\rm visible}$ when the DM and visible matter densities assume the same order of magnitude values, irrespectively to the axion mass $m_a$ or initial misalignment angle $\theta_0$.This model, without adjustment of any parameters, gives a very reasonable intensity of the extreme UV (EUV) radiation from the solar corona as a result of the AQN annihilation events with the solar material. This extra energy released in corona represents a resolution, within AQN framework, a long standing puzzle known in the literature as the"solar corona heating mystery". The same annihilation events also produce the axions. The flux of these axions is unambiguously fixed in this model and expressed in terms of the EUV luminosity from solar corona. We make few comments on the potential discovery of these gravitationally bound axions.


Introduction
The Peccei-Quinn (PQ) mechanism accompanied by the axions remains the most compelling resolution of the strong CP problem, see original papers [2,3,4] and recent reviews [5,6,7,8,9,10,12,11,13] on the subject. We refer to the review articles for the discussions and analysis on the recent activities in the field of the axion searches by a numerous number of different groups using very different instruments.
For the purposes of the present work it is sufficient to mention that the conventional dark matter galactic axions are produced due to the misalignment mechanism [14] when the cosmological field θ(t) oscillates and emits cold axions before it settles down at its final destination θ final = 0. Another mechanism is due to the decay of the topological objects [15,16,17,18,19]. There is a number of uncertainties and remaining discrepancies in the corresponding estimates. We shall not comment on these subtleties 1 by referring to the original papers [15,16,17,18,19]. It is important that in both cases the produced axions are non-relativistic particles with typical v axion /c ∼ 10 −3 , and their contribution to the dark matter density scales as Ω axion ∼ m −7 /6 It has been suggested in recent work [1] that there is a fundamentally novel mechanism of the axion production in the Sun. This mechanism is deeply rooted to the so-called axion quark nugget (AQN) dark matter model when the stability of the nuggets is supported by the axion domain wall. The most important consequence of the new production mechanism is that the emitted axions (from the axion domain wall when the nugget gets annihilated) will be released with relativistic velocities with typical values v AQN axion 0.5c. These features should be contrasted with conventional galactic non-relativistic axions v axion ∼ 10 −3 c and solar ultra-relativistic axions with typical energies E = 4.2 keV.
The computations in ref. [1] of the spectral properties of the axions produced by this novel mechanism were based on the approximation which is known to be badly violated for low-energy part of the spectrum with v c. This part of the spectrum represents very tiny portion of the produced axions. Therefore, it had been ignored in the original studies [1]. However, the upgraded CAST instrument will be highly sensitive to the low energy part of the spectrum. Therefore, it is highly desirable to develop a new computational technique which would allow to carry out the computations in the region of small velocities v c.
Furthermore, the low-energy axions produced in the Sun might be trapped by strong gravitational force such that v ≤ v trapped will be trapped by the Sun indefinitely. The v trapped is numerically the same as the free fall velocity, While the portion of these low energy axions is tiny as we shall estimate below, these trapped axions may play an important role in physics as they will be accumulated around the Sun during entire life time of the solar system, i.e. around 4.5 billion years. The effects related to the trapped axions are not new, and discussed previously in the literature [22]. The goal here is to present some numerical estimates for our specific AQN model when the axions which are produced as a result of the annihilation events in the solar atmosphere and will be indefinitely bounded to the Sun. Therefore, the main goal of the present studies is to develop a new technique to generalize 2 the results of ref. [1] to perform the self-consistent computations of the axion spectrum in the regime when the axion velocities are small v c.
The paper is organized as follows. In next section 2 we overview the AQN model by paying special attention to the astrophysical and cosmological consequences of this specific dark matter model. In section 3 we highlight the basic arguments of ref. [23] advocating the idea that the annihilation events of the antinuggets with the solar material can be interpreted as the nanoflares conjectured by Parker long ago. The recent numerical simulations [24] strongly support the original estimates [23] and explicitly show that the disintegration of the AQNs occurs precisely at the altitude about 2000 km above the photosphere, in the vicinity of the so-called transition region (TR). Injection of this huge amount of energy as a result of the AQN annihilation events represents the resolution of the so-called " the solar corona heating puzzle" within AQN framework. The disintegration of the nuggets inevitably produces the axions in the vicinity of the TR. The computations of the intensity and spectral properties of these axions with relativistic velocities v ∼ c have been carried out in ref. [1]. In Section 4 we develop a new technique which allows to generalize these computations for low energy portion of the spectrum when v c. We use the corresponding results in Section 5 to discuss the physics of the trapped axions and we highlight the basic ideas how to discover them. We conclude in Section 6 with few thoughts on the future developments.

Axion Quark Nugget (AQN) dark matter model
The axion field plays a key role in the construction. Therefore, we would like to make a short overview of this model with emphasize on the role of the axion field and related astrophysical consequences of this proposal.
The AQN model was invented long ago [25] (though a specific formation mechanism of the nuggets was developed in much more recent papers [26,27,28]) as a natural explanation of the observed ratio Ω dark ∼ Ω visible . Indeed, the similarity between dark matter Ω dark and the visible matter Ω visible densities strongly suggests that both types of matter have been formed during the same cosmological epoch, which must be the QCD transition as the baryon mass m p ∼ Λ QCD represents the visible portion of the matter Ω visible .
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 [29], strangelets [30], nuclearities [31], see also review [32] with large number of references on the original results. The AQN model in the title of this section stands for the axion quark nugget model [25] to emphasize on essential role of the axion field in the construction and to avoid confusion with earlier models [29,30,31,32] mentioned above. The AQN model is drastically different from previous similar proposals in two key aspects: 1. There is an additional stabilization factor in the AQN model provided by the axion domain walls which are copiously produced during the QCD transition; 2. The AQN could be made of matter as well as antimatter in this framework as a result of separation of charges, see recent papers [26,27,28] with large number of technical details.
The basic idea of the AQN proposal can be explained as follows: It is commonly assumed that the Universe began in a symmetric state with zero global baryonic charge and later (through some baryon number violating process, the so-called baryogenesis) evolved into a state with a net positive baryon number. As an alternative to this scenario we advocate a model in which "baryogenesis" is actually a charge separation process when the global baryon number of the Universe remains zero. In this model the unobserved antibaryons come to comprise the dark matter in the form of dense nuggets of quarks and antiquarks in colour superconducting (CS) phase. The formation of the nuggets made of matter and antimatter occurs through the dynamics of shrinking axion domain walls, see original papers [26,27,28] with many technical details.
The nuggets, after they formed, can be viewed as the strongly interacting and macroscopically large objects with a typical nuclear density and with a typical size R ∼ (10 −5 − 10 −4 )cm determined by the axion mass m a as these two parameters are linked, R ∼ m −1 a . This relation between the size of nugget R and the axion mass m a is a result of the equilibration between the axion domain wall pressure and the Fermi pressure of the dense quark matter in CS phase. One can easily estimate a typical baryon charge B of such macroscopically large objects as the typical density of matter in CS phase is only few times the nuclear density. However, it is important to emphasize that there are strong constraints on the allowed window for the axion mass, which can be represented as follows 10 −6 eV ≤ m a ≤ 10 −2 eV. This axion window corresponds to the range of the nugget's baryon charge B which largely overlaps with all presently available and independent constraints on such kind of dark matter masses and baryon charges 10 23 ≤ |B| ≤ 10 28 , see e.g. [34,33] for review. The corresponding mass M of the nuggets can be estimated as M ∼ m p B, where m p is the proton mass. This model is perfectly consistent with all known astrophysical, cosmological, satellite and ground based constraints within the parametrical range for the mass M and the baryon charge B mentioned above (2). It is also consistent with known constraints from the axion search experiments. Furthermore, there is a number of frequency bands of radiation from the galactic centre where some excess of emission was observed, but not explained by conventional astrophysical sources. Our comment here is that this model may explain some portion, or even entire excess of the observed radiation in these frequency bands, see short review [33] and additional references at the end of this section.
The key element of the construction is the coherent CP-odd axion field θ which is not vanishing during the QCD transition in early Universe. As a result of the CP violating processes the number of nuggets and anti-nuggets being formed would be different. This difference is always of order of one effect [26,27,28] irrespectively to the parameters of the theory, the axion mass m a or the initial misalignment angle θ 0 . As a result of this disparity between nuggets and anti nuggets a similar disparity would also emerge between visible quarks and antiquarks. This is precisely the reason why the resulting visible and dark matter densities must be the same order of magnitude [26,27,28] as they are both proportional to the same fundamental Λ QCD scale, and they both are originated at the same QCD epoch. If these processes are not fundamentally related the two components Ω dark and Ω visible could easily exist at vastly different scales.
Unlike conventional dark matter candidates, such as WIMPs (Weakly interacting Massive Particles) the darkmatter/antimatter nuggets are strongly interacting and macroscopically large objects, as we already mentioned. However, they do not contradict any of the many known observational constraints on dark matter or antimatter in the Universe due to the following main reasons [35]: They carry very large baryon charge |B| 10 23 , and so their number density is very small ∼ B −1 . As a result of this unique feature, their interaction with visible matter is highly inefficient, and therefore, the nuggets are perfectly qualify as DM candidates. Furthermore, the quark nuggets have very large binding energy due to the large gap ∆ ∼ 100 MeV in CS phases. Therefore, the baryon charge is so strongly bounded in the core of the nugget that it is not available to participate in big bang nucleosynthesis (bbn) at T ≈ 1 MeV, long after the nuggets had been formed.
It should be noted that the galactic spectrum contains several excesses of diffuse emission the origin of which is unknown, the best known example being the strong galactic 511 keV line. If the nuggets have the average baryon number in the B ∼ 10 25 range they could offer a potential explanation for several of these diffuse components. The parameter B ∼ 10 25 was fixed in this proposal by assuming that this mechanism saturates the observed 511 keV line [36,37], which resulted from annihilation of the electrons from visible matter and positrons from anti-nuggets. The most relevant for the present purposes application is the proposal of ref. [23] that the AQN dark matter particles might be responsible for the heating of the solar corona. As the annihilation processes of the AQNs in the solar corona play crucial role in our present studies on the axion production (which is a direct consequence of these annihilation processes), we choose to overview the most important specific elements of this proposal in next section 3 to separate them from the basic generic ideas of the AQN framework presented above.

AQNs as the corona's heaters
It has been known for quite some time that the total intensity of the observed EUV and soft x-ray radiation (averaged over time) assumes the following value, which represents (since 1939) the renowned "the solar corona heating puzzle", see e.g. a general review [38] on the subject and also Ref. [22] with analysis of some specific features related to present work. The observation (4) implies that the corona has the temperature T 10 6 K which is ∼ 10 2 times hotter than the surface temperature of the Sun, and conventional astrophysical sources fail to explain the EUV and soft x ray radiation from corona. One should comment here that "the solar corona heating puzzle" includes a number of elements which are hard to explain using a conventional framework. In particular, the hot corona with drastically smaller mass density cannot be in equilibrium with the ∼ 300 times cooler solar surface [39]. In order to maintain the quiet Sun high temperature corona, some non-thermally supplied energy must be dissipated in the upper atmosphere [40]. It must occur over entire solar surface where typical magnetic field is B ∼ 1G, which is much smaller than in the sunspot regions (with B ∼ 10 2 G) occupying in less than 1% of the surface. The supply of energy must also take place somehow during quiet periods of the solar cycles when sunspots and/or flares may not be present in the system for months. There are many other problems which are nicely stated in review article [41] as follows "everything above the photosphere...would not be there at all".
A drastically different scenario is suggested in ref. [23] when the energy deposition is originated from outside the system, in contrast with previously considered proposals when the energy is originated from the deep dense regions of the sun. We want to argue that the observed peculiar behaviour might be intimately related to this fundamentally distinct scenario when the extra source of the energy is associated with the dark matter nuggets continuously entering the sun from outer space. A large amount of energy is available in the proposal as result of huge energy deposition of such dark matter constituents before being disintegrated.
Our goal here is to overview the proposal [23]. We start with simple estimates. The impact parameter for capture and crash of the nuggets by the Sun can be estimated as where v 10 −3 c is a typical velocity of the nuggets. Assuming that ρ DM ∼ 0.3 GeVcm −3 and using the capture impact parameter (5), one can estimate the total energy flux due to the complete annihilation of the nuggets, where we substitute constant v 10 −3 c to simplify numerical analysis. This order of magnitude estimate is very suggestive as it roughly coincides with the observed total EUV energy output from corona (4) representing ∼ (10 −7 − 10 −6 ) portion of the total solar luminosity. Precisely this "accidental numerical coincidence" was the main motivation to put forward the idea [23] that (6) represents a new source of energy feeding the EUV and soft x-ray radiation.
The main assumption made in [23] is that a finite portion of annihilation events have occurred before the antinuggets entered the dense regions of the Sun. This assumption has been recently supported by numerical Monte Carlo simulations [24] which explicitly show that indeed, the dominant energy injection occurs in vicinity of the transition region at the altitude ∼ 2000 km. These annihilation events supply the energy source of the observed EUV and x-ray radiation from the corona and the choromosphere. The crucial observation made in [23] and confirmed in [24] is that while the total energy due to the annihilation of the anti-nuggets is indeed very small as it represents ∼ 10 −6 fraction of the solar luminosity according to (4), nevertheless the anti-nuggets produce the EUV and x-ray spectrum because the most of the annihilation events occur in vicinity of the transition region at the altitude ∼ 2000 km characterized by the temperature T ∼ 10 6 K. Such spectrum observed in corona and the chromosphere is hard to explain by any conventional astrophysical processes as mentioned at the beginning of this section.
The basic claim of [23,24] is that the annihilation events of the anti-nuggets, which generate huge amount of energy (6) can be identified with the "nanoflares" conjectured by Parker long ago [42]. In most studies the term "nanoflare" describes a generic burst-like event for any impulsive energy release on a small scale, without specifying its cause. In other words, in most studies the hydrodynamic consequences of impulsive heating (due to the nanoflares) have been used without discussing their nature, see review papers [38,43]. The novel element of refs. [23,24] is that the the nature of the nanoflares was specified as annihilation events of the dark matter particles within AQN framework, i.e. nanoflares ≡ AQN annihilation events.
The main arguments supporting the identification (7) are: 1. In order to reproduce the measured radiation loss, the observed range of nanoflares needs to be extrapolated from sub-resolution events with energy 3.7 · 10 20 erg to the observed events interpolating between (3.1 · 10 24 − 1.3 · 10 26 ) erg, table 1 in ref. [44]. This energy window corresponds to the (anti)baryon charge of the nugget 10 23 ≤ |B| ≤ 4·10 28 which largely overlaps with allowed window (2) for AQNs reviewed in section 2. This is a highly nontrivial consistency check for the proposal (7) as the window (2) comes from a number of different and independent constraints extracted from astrophysical, cosmological, satellite and ground based observations. The window (2) is also consistent with known constraints from the axion search experiments within the AQN framework. Therefore, the overlap between these two fundamentally different entities represents a highly nontrivial consistency check of the proposal (7).
2. Our next argument goes as follows. The nanoflares are distributed very uniformly in quiet regions, in contrast with micro-flares and flares which are much more energetic and occur exclusively in active areas. In fact, it was reported 1.1 × 10 6 events per hour over the whole Sun for SoHO/EIT observations [45,46]. It is perfectly consistent with our identification (7) as the anti-nugget annihilation events should be present in all areas irrespectively to the activity of the Sun. At the same time the flares are originated in the active zones, and therefore cannot be uniformly distributed.
3. The observed Doppler shifts (corresponding to velocities 250 − 310 km/s) and the observed line width in OV of ±140 km/s far exceed the thermal ion velocity which is around 11 km/s, see Fig.5 in ref. [47]. These observed features can be easily understood within the AQN scenario. Indeed, the typical velocities of the nuggets entering the solar corona is about ∼ 300 km/s. Therefore, it is perfectly consistent with observations of the very large Doppler shifts and related broadenings of the line widths. Typical time-scales of the nanoflare events, of order of (10 1 − 10 2 ) seconds are also consistent with results of refs. [23,24] 4. It has been observed [48] that "the pre-flare enhancement propagates from the higher levels of the corona into the lower corona and chromosphere." It is perfectly consistent with our proposal as the dark matter AQNs enter the solar atmosphere from outer space. Therefore, they first enter the higher levels of the corona where they generate the shock wave, before they reach chromosphere in τ ∼ (10 − 10 2 ) seconds. 5. One should emphasize that the isotropic EUV radiation is very different in all respects from highly anisotropic distribution of sunspots and flares. It was emphasized in [23,24] on this qualitative difference to argue that the flares are originated at the sunspots area with locally large magnetic field B ∼ (10 2 −10 3 ) G, while the EUV emission (which is observed even in very quiet regions where magnetic field is relatively small in range B ∼ 1G) is isotropic and covers entire solar surface. A typical variation of the EUV radiation during the solar cycles (between their minimum and maximum values) is very modest in comparison with drastic variation of the flare occurrences. The differences in variation is of order of 10 2 , see e.g. [49]. It shows once again that the nature of nanoflares and large flares must be very different. This is consistent with our proposal (7) because the nanoflares are identified with AQN annihilation events while flares occur as a result of magnetic reconnection in active regions where magnetic field is large and plays the dominant role in dynamics. 6. Last but not least: the AQN resolution of the solar corona puzzle also resolves another mystery [49] where it was claimed that a number of highly unusual phenomena observed in solar atmosphere might be related to the gravitational lensing of "invisible" streaming matter towards the Sun which is correlated with positions of the planets. This is really a weird correlation because one should not expect any connections between the flare occurrences, the intensity of the EUV radiation, and the position of the planets. Nevertheless, the analysis [49] obviously demonstrates that this naive expectation is not quite correct. At the same time, such "weird" correlations naturally occur within AQN framework. This is because the dark matter AQNs, being the "invisible streaming matter" (in terminology of ref. [49]) can play the role of the triggers sparking the large flares [50]. Therefore, the observation of the correlation between the EUV intensity, the frequency of the flares and positions of the planets can be considered as an additional supporting argument of the dark matter explanation of the observed EUV irradiation (4), because both effects are originated from the same dark matter AQNs. As a direct consequence of this relation we expect that the intensity of the the axion emission from the Sun (which always accompanies the EUV emission) will be also correlated with the position of the planets.

Axions from AQNs: Intensity and the Spectrum
In Section 2 we explained that the axion field is the key element in the AQN construction. In Section 3 we argued that the AQNs may serve as the heaters of the corona. In recent paper [1] we argued that these annihilation events in corona will inevitably produce the axions. The computations of the spectral features of these axions in ref. [1] were based on an approximation which does not admit analysis of the spectrum in non-relativistic domain at v c. At the same time, as we discuss in Introduction this non-relativistic domain plays a key role in gravitational trapping. The corresponding computations is the main subject of this section. First, in subsection 4.1 we highlight the basic results from ref. [1] by providing a self-contained text for the convenience of the readers. In next subsections 4.2 and 4.3 we explain the computational framework and present the results of the computations, referring to Appendix A for the technical details.

Intensity
The axions play a key role in construction of the AQNs as they provide an additional pressure to stabilize the nuggets, see section 2 for review. The corresponding axion contribution into the total nugget's energy density has been computed in [28]. Depending on parameters the axion's contribution to the nugget's mass represents about 1/3 of the total mass. It can be translated in terms of the axion luminosity from the Sun as follows where L (AQN) is given by (6). The corresponding axion flux measured on Earth can be computed as follows [1] Φ(solar axions) ∼ L (axion) 4π E a D 2 ∼ 0.3 · 10 17 1 cm 2 s 10 −5 eV m a , D 150 · 10 6 km, where we assume that the axion's energy when the antinuggets get annihilated is slightly relativistic E a 1.2m a , but never becomes very relativistic. The corresponding energy flux is [1] m a Φ(solar axions) ∼ 3 · 10 11 eV cm 2 s .
These estimates should be compared with conventional cold dark matter galactic axion contribution assuming the axions saturate the observed DM density: Similar estimates can be also carried out for Earth. In this case as explained in [1] the observations of the E&M showers due to the nuggets entering the Earth's atmosphere (before hitting the Earth's surface) require very large area detectors. The nuggets will continue to radiate E&M energy in the deep underground. However, this radiation by obvious reasons cannot be recovered and observed. At the same time the observation of the axions (which have been produced as a result of the annihilation events in the very deep underground) is possible, and in fact very promising. Indeed, the corresponding axion flux can be estimated as follows [1] where ∆B/B is the portion of the AQNs which get annihilated in the Earth's core. Interestingly, the axion flux (12) which is generated due to the AQN annihilation events is much larger than the flux (10) generated due to the AQN annihilation events in the solar corona and measured on Earth. At the same time, the axion flux (12) is the same order of magnitude as the conventional cold dark matter galactic axion contribution (11). This is because the parameter ∆B/B ∼ 1 is expected to be order of one, as a finite portion of the AQNs will get annihilated in the Earth's core. However, the wave lengths of the axions produced due to AQN annihilations, are much shorter due to their relativistic velocities v ∼ 0.5c, in contrast with conventional galactic isotropic axions with v ∼ 10 −3 c. Therefore, these two distinct contributions can be easily discriminated.

Spectral properties. General Comments
The basic idea of the computation of the spectrum is as follows. Consider an AQN loosing its mass due to the annihilation with surrounding material, such that the axion contribution quickly becomes the dominant portion of total mass. One should comment here that the axion domain wall in the equilibrium does not emit any axions as a result of pure kinematical constraint: the domain wall axions are off-shell axions in the equilibrium. The time dependent perturbation obviously changes this equilibrium configuration. In other words, the configuration becomes unstable because the total energy of the system is no longer at its minimum. To retrieve its ground state, an AQN will therefore intend to lower its domain wall mass by radiating the axions. To summarize: the emission of axions is an inevitable consequence during the annihilation of antinuggets in simply for the reason to maintain the AQN stability. Now, we want to identify a precise mechanism which produces the on-shell freely propagating axions emitted by the axion domain wall. In this section we overview the basic idea of the computational technique to be used. To address this question, we consider the general form of a domain wall solution: where R 0 is the radius of the AQN, φ w is the classical solution of the domain wall, while χ describes the excitations due to the time dependent perturbation. We should note that, φ w is clearly off-shell classical solution, while χ describes the on-shell propagating axions. Thus, whenever the domain wall is excited, namely χ 0, freely propagating axions will be produced and emitted by the excitation modes. Few comments are in order before we proceed in subsection 4.3 with description of the technical details and corresponding results. First, if the domain wall can be considered to be infinitely large in xy direction such that the profile function depends on a single variable z the computations can be carried out easily as the classical profile function φ w (z) is known exactly. This is precisely the procedure which has been adopted in previous paper [1]. The corresponding technique is justified when a typical size L x ∼ L y m −1 a along x, y is much larger than the width of the domain wall of order m −1 a . If the wave length of the emitted axion is small, i.e. λ a ∼ m −1 a the axions cannot carry any information about the finite size of the system and the approximation is marginally justified (λ a stands for the de Brogile wavelength of the emitted axion). This is precisely the approximation, the so-called "thin wall approximation" which has been adopted in computations [1]. This approximate treatment is marginally justified for relativistic axions with v ∼ c, and we expect that accounting for the finite size of the system cannot drastically change the results in the relativistic domain v ∼ c. This will be explicitly confirmed below by present computations accounting for finite size of the system. Secondly, it is quiet obvious that the "thin wall approximation" is badly broken for non-relativistic axions with v c when λ a m −1 a and a new technique must be developed for proper analysis. The basic idea of computation accounting for finite size of the system R goes as follows. Suppose an AQN is traveling in vacuum where no annihilation event takes place, we expect the solution stays in its ground state φ(R 0 ) = φ w (R 0 ) which corresponds to the minimum energy state. Since there is no excitation (i.e. χ = 0), no free axion can be produced. However, the scenario drastically changes when some baryon charge from the AQN get annihilated. Due to these annihilation processes, the AQN starts to loose a small amount of its mass, and consequentially its size shrinks from R 0 to a slightly smaller radius R new = R 0 − ∆R. Note that its quantum state φ(R 0 ) = φ w (R 0 ) is no longer the ground state, because a lower energy state φ w (R new ) becomes available. Then, we may write the current state of the domain wall as φ(R 0 ) = φ w (R new )+φ w (R new )∆R, so the domain wall now has a nonzero exciting mode χ = φ w (R new )∆R and free axions can be produced during oscillations of the domain wall. To reiterate: the annihilation of antinuggets with surrounding matter forces the domain wall to oscillate. These domain wall oscillations generate excitation modes which ultimately lead to radiation of the propagating axions.
Our last comment deals with terminology and notations. The results for the spectrum obtained using the "thin-wall approximation" is coined as 1D spectrum. As we mentioned above it is marginally justified when λ a ∼ m −1 a , and it admits mathematically exact treatment which was previously presented in [1]. In the present work we mostly deals with 3D computations when a finite size of the system plays a key role, which is always the case for λ a m −1 a . The potential pitfall is that some technical simplifications are required to treat the 3D case. Consequentially, the obtained results might be sensitive to these technical simplifications. In order to we characterize the sensitivity to our technical simplifications we introduce a tunable parameter δ which varies from 0 to 1, so δ will serve as a probe to test the robustness of the obtained result. As we shall see below, the obtained results are not very sensitive to the choice of δ. Therefore we conclude that our 3D results are robust and reliable.
In what follows will express the normalized spectrum as a function of the speed of emitted axion v a /c defined in the nugget's frame, defined as follows where Φ tot axions is the axion flux inserted to eq. (14) for normalization purposes. It assumes the magnitude Φ(solar axions) given by eq. (10) for the solar axions, and the value Φ(Earth axions) given by (12) for the axions emitted from the Earth's core.

Spectral properties. Results
We follow the procedure described above in subsection 4.2 and present the axion field in time dependent background as follows where φ w,δ (r − R 0 ) satisfies the classical equation of motion while χ(t, r) describes the time-dependent excitations. As exact solution accounting for the finite size of the nugget is not known we parameterize different simplified solutions by parameter δ. We consider parameter δ as a probe-we treat a mild sensitivity to δ as a signal that our approach is reliable and robust. The next step is to expand the action S[φ] by keeping the quadratic terms only, where L 2 [δ] is the second order linear differential operator which depends on classical solution φ w,δ (r − R 0 ) parameterized by δ, see Appendix A for the technical details. The next step, as usual, is to expand the fluctuations χ in terms of complete basis and compute the coefficients a plm in this expansion. The result for the total radiated energy E rad is given by eq(A.24) from Appendix A. It can be presented it in the following form where the coefficients a plm can be explicitly computed and are given by (A.22). The expression for the radiated energy (17) allows us to compute the desired spectrum ρ(v a ) defined by (14). The results of the computations are presented on Fig. 1a for a specific choice of parameter δ = 0. The low energy portion of the spectrum with 0 ≤ v/c ≤ 0.01 is shown on Fig.1b. Few comments are in order. Parameter δ in our procedure was introduced as a probe to test our computational scheme. With the purpose of the test we performed similar computations for δ 0. The corresponding results are plotted in Appendix on Fig. A.2a and Fig. A.2b for δ = 0.5 and Fig. A.2c and Fig. A.2d for δ = 1. One can explicitly see that the results for the spectrum are not very sensitive to parameter δ. To reiterate: the basic qualitative results are not very sensitive to choice of parameter 0 ≤ δ ≤ 1. Therefore, our results can be considered to be very reliable and robust.
One next comment goes as follows. It is very instructive to compare our 3D computations with 1D computations presented in [1]. We had anticipated before the 3D computations have been carried out that the results in the relativistic domain v a /c 0.5 should not be drastically modified in comparison with simplified treatment in [1]. We can now confirm that this is indeed the case. At the same time we had expected the drastic modification of the spectrum in the non-relativistic regime v a /c ≤ 0.01. Indeed, the 3D spectrum in this domain behaves as ρ(v a ) ∼ v 3 a in contrast with linear dependence in simplified treatment in ref. [1]. This difference in behaviour at small v a /c 1 can be attributed to the phase volume suppression ∼ d 3 k in 3D case for λ a m −1 a .

Gravitationally trapped axions
In the previous section we computed the portion of the axions which have sufficiently low velocities (below escape velocity) such that they will be orbiting the Sun as long as it exists. This portion of the non-relativistic axions is extremely tiny as shown on Fig. 1b. Nevertheless, the effect could be drastically enhanced as we discuss below due to accumulation of these axions during entire life of the solar system, i.e. for ∼4.5 billion of years.
The condition for the axions to be bounded after being radiated is determined by the trapped velocity v trapped , defined as such that all axions with v ≤ v trapped will be trapped by the Sun and the axions with v ≤ v trapped ⊕ will be trapped by the Earth. The effect of the trapped axions is not new, and discussed previously in the literature [22]. The goal here is to present some numerical estimates for our specific AQN model when the axions which are produced as a result of the annihilation events can be trapped in the solar atmosphere. These estimates will play a key role in our discussions on the discovery potential of these axions.

Solar corona background. Non-resonance case.
According to Fig. 1b these highly non-relativistic axions represent a very tiny portion of the produced axions. The energy which is accumulated in the solar atmosphere per unit time as a result of trapping these axions can be estimated as follows dE dt (trapped axions) 1.6 · 10 27 · ξ · erg s 10 17 · ξ 10 −10 · erg s (19) where we used the expression (8) for the rate of the energy transfer to the axions. We also introduced the suppression factor ξ to account for the small fraction of the trapped axions with v ≤ v trapped . For numerical estimates in formula (19) we use suppression factor ξ ∼ 10 −10 computed in previous section and presented on Fig. 1b. The axions (19) could not leave the system during entire life time of the Sun, i.e. 4.5 billion years 10 17 s. Therefore, the total energy accumulated by the Sun and related to AQN annihilation events radiating the slow velocity axions can be estimated as follows which of course represents a very tiny fraction of the solar mass M 2 · 10 30 kg. The energy (20) corresponds to the following total number of the axions accumulated by the Sun during its lifetime: If we assume that the majority of these axions are localized within 2 solar radius R , we arrive to the following estimate for the average axion energy density inside this volume which is 3 orders of magnitude larger than the present average dark matter density today ρ DM 0.3 GeV cm 3 . One should comment here that this enhancement of the DM density in the vicinity of the Sun obviously not in contradiction with most precise observational upper limits on solar system (SS) -bound DM, which is normally expressed as ρ SS < 2 · 10 5 GeV cm 3 , see e.g. [51]. It is also interesting to note that some authors [52,53] previously argued that the DM in the SS might be greatly enhanced (on the level of 10 3 ) as a result of capturing of DM particles from the Galactic halo due to the 3 body interaction (the Sun, a planet and DM particle). Other authors [54,55] estimated that the effect of capturing is small. We refer to these original papers for the discussions and details. The only comment we would like to make here is that the effect estimated in eq. (23) is fundamentally distinct in nature in comparison with previously discussed effect [52,53,54,55]. The novel effect which is the subject of this work is entirely rooted to the AQN model when the nuggets get disintegrated when enter the solar atmosphere. The corresponding annihilation events produce the low velocities axions with v ≤ v trapped . These axions which behave as DM particles surrounding the Sun have no relation to the effect discussed in [52,53,54,55]. Now we want to estimate the number density n axions of these axions assuming, as before, that the majority of the axions are localized within 2 solar radius R .
Can these axions be observed? These axions cannot decay as the axion life time τ(a → 2γ) is very long. However, these axions can be converted to photons in the background of external magnetic field. The corresponding probability of this conversion is determined by the formula [56,57]: where L is a typical distance where the magnetic field B is present. For non-relativistic axions one can approximate q ± ±ω. Furthermore, for our present analysis we assume that typical B ∼ 300 G in the solar atmosphere, while L is very large such that sin 2 qL 2 can be approximated as 1 2 . Therefore, probability of the conversion can be approximated as follows where z = m u /m d 0.56 and parameter E/N = 0 for KSVZ model, and E/N = 8/3 for DFSZ model. For simple numerical analysis we take E/N = 0 to arrive to the following estimate The number of the produced photons (as a result of the conversion from the axions) per unit volume with the frequency ω = m a is estimated as follows dN(a → γ) dV n axions · P a→γ 10 −18 B 300 G 2 ξ 10 −10 10 −5 eV m a 1 cm 3 (28) where n axions is estimated in (24). These converted photons obviously can leave the system. The total number of photons leaving the system through area ∼ 4π(2R ) 2 per unit time is given by These photons are very monochromatic with ω = m a with accuracy of order 10 −3 . Potentially, it gives us some chance to observe them on Earth. The corresponding count of photons dF(a → γ) arriving from the Sun with monochromatic frequency ω = m a (due to the axion-photon conversion) is estimated as This count, of course, is extremely low. However, these estimates were based on rate (27) corresponding a → γ conversion in vacuum. As it is known since [56] the rate could be drastically enhanced if the system is placed in a media with non-vanishing plasma frequency ω p exactly matching the axion mass, i.e. ω p = m a , which represents the topic for the next subsection.

Solar Corona background. Resonance conversion in solar plasma
We start with numerical estimation for the plasma frequency ω p in the solar corona where the most axions are released as a result of the AQN's annihilation events, ω p ≡ 4παn m e 3.5 · 10 −6 · n 10 10 cm −3 1 2 eV. (31) The numerical similarity between ω p and the expected value for the axion mass m a from allowed window m a ∈ (10 −6 − 10 −3 ) eV represents the basic motivation for analysis in this subsection. In other words, our goal here is to study possible observational consequences of the resonance case when the condition ω p = m a could occur in the corona, which is explicit manifestation of the so-called "level-crossing effect" as formulated in ref. [56]. If the condition ω p = m a is fulfilled the corresponding resonance a → γ conversion in media is determined by formula [56]: where we adopted the notations for ∆ M from [56] and expressed the fundamental PQ mass scale M from [56] in terms of the original definition for g aγ . Of course we do not expect that this condition can be exactly satisfied in reality in nature. Furthermore, the oscillation length l deg = π/∆ M is very long, much longer than the size of the system such that P a→γ never becomes of order one effect. However, our goal here is different, and we present formula (32) exclusively for illustrative purposes to illuminate the role of the distance scale where the conversion occurs. With this purpose we expand the resonance expression (32) assuming that ∆ M L 1 and compare with non-resonance case (27) to arrive where we consider special case θ = π/2 to simplify the arguments. Formula (33) illustrates very important point: small conversion rate in non-resonance case (27) corresponds to very short distance ∼ m −1 a where this conversion takes place. Indeed, the first brackets in (33) identically coincides with formula (27) describing the conversion in non-resonance case. Precisely this first term describes a huge suppression factor.
For our present studies it is important to emphasize that the same formula (33) also explicitly shows that this suppressed conversion (27) can be greatly enhanced with the second factor ∼ (m a L) 2 if one can increase the coherence length L by maintaining ω p = m a . If the coherence can be maintained on much larger scale than m −1 a such that (m a L) 2 1 the effect of conversion P a→γ will be strongly enhanced in comparison with (27) by this large factor (m a L) 2 1 entering formula (33). It is clear that one should not expect that the effect could be of order one as one cannot maintain the coherence on the huge scale l deg = π/∆ M . However, some enhancement in comparison with (27) still can be achieved.
The same conclusion also follows from the following expression which was derived using the perturbation theory by treating the inhomogeneities of the magnetic field and plasma density as small perturbations [56] where we neglected ∆ vac || ∼ B 2 which is numerically very small for relatively weak typical solar magnetic field. The conversion rate given by eq. (34) generates the enhancement proportional to the large length L if the phases maintain the coherence and the cancellations between different phases do not occur due to the fast fluctuations. The requirement that the coherence is maintained up to the scale L is determined by the following condition If this condition is fulfilled then the conversion rate given by eq. (34) reduces to our previous expression (33) with enhancement factor ∼ L 2 , i.e.
This supports our previous conclusion that the enhancement factor (m a L) 2 is a result of constructive interference. The corresponding length scale L is determined by condition (35). Now we want to address the following question: What is the typical length scale L where the condition (35) can be satisfied in solar atmosphere? We limit our analysis with the trapped axions which have non-relativistic velocities with ω m a as discussed in previous section 5.1. These axions are distributed in the entire solar atmosphere. Therefore, there is always an extended region in corona or chromosphere where the electron density n is such that the plasma frequency (31) equals the axion mass, i.e. m a = ω p . The only question remains to be answered is what are the typical length scales where the average value n for the electron density varies 3 .
To estimate the corresponding scale L we notice that a typical variation of the density (and the plasma frequency) n by factor ∼ 10 occurs when the altitude changes by ∼ 10 3 km. Assuming a linear extrapolation (excluding very fast changes in the transition region) one should expect that the variation of the density δn/n ∼ 1 occurs on the scale of order l 0 ∼ 10 2 km. This estimate implies that the relative variation (mismatch) of the density on the coherence scale L of the axion/photon oscillation must not exceed λ/L to be consistent with (35). In other words, the scale L where coherence (35) can be maintained must satisfy the following condition Precisely at this coherence scale L the mismatch in plasma frequency is sufficiently small as δn/n ∼ λ/L. At the same time δn/n ∼ (λ/L) · (l 0 /L) ∼ 1 becomes order of one at much larger scales ∼ l 0 where coherence, of course, cannot be maintained. One should emphasize that this very rough estimate assumes a linear extrapolation. This assumption may or may not be justified in reality in solar atmosphere. One should emphasize here that any variation of the magnetic field entering (34) do not modify our estimate for the coherence length (37). This is because the estimate (37) is sensitive to the phase variation (rather than to the amplitude changes ∼ ∆ M (z)) determined by a steady and systematic variation of the plasma frequency ω p ∼ √ n in corona. If one literally accepts the estimate (37) the corresponding enhancement factor can be approximated as follows, (Lm a ) 2 ∼ 4 · 10 6 m a 10 −5 eV 1.
Our previous (non-resonance) estimate (30) should be multiplied the enhancement factor (38) for case if the resonance conditions can be satisfied and the linear extrapolation is justified. The rate of conversion with this factor becomes While the count (39) is still low, some hope is that this is an unique monochromatic line. Furthermore, the intensity of this line must be correlated with the EUV emission. In addition, during the flares the magnetic field B might be very large in the solar system which provides some enhancement factor and possible correlations with the flares. Finally, this monochromatic line can be easily discriminated from the noise as it should appear only along the line-of-sight in the direction of the Sun. Needless to say that a strong magnetic field in a detector is not required for the observation of these photons on Earth because the axion-photon conversion occurs in the solar atmosphere rather than on Earth. In a sense we use entire Sun as a one big helioscope where the trapped axions have been accumulated during 4.5 billion years and where the axions can be converted to photons in the entire solar atmosphere.

Axions from the Earth's Underground
According to ref. [1] the axion flux due to the AQN annihilation events in the very deep underground is given by (12). We integrate this rate over entire surface to arrive where we used the expression (8) for the rate of the energy transfer to the axions. We also introduced the suppression factor ξ ⊕ to account for the small fraction of the trapped axions with v ≤ v trapped . For numerical estimates in formula (19) we use suppression factor ξ ⊕ ∼ ξ · (v ⊕ /v ) 4 ∼ 10 −17 computed in previous section and given by (18). This is of course very tiny rate even when ∆B/B ∼ 1 as we expect. The axions (40) could not leave the system during entire life time of the Earth, i.e. 4.5 billion years 10 17 s. Therefore, the total energy accumulated by the Earth and related to AQN annihilation events radiating the slow velocity axions can be estimated as follows This energy can be expressed in terms of extra Earth's mass ∆M ⊕ accumulated by the Earth and represented by the trapped axions which of course represents a very tiny fraction of the Earth mass M ⊕ 5.9 · 10 24 kg. The energy (41) corresponds to the following total number of the axions accumulated by the Earth during its life-time: If we assume that the majority of these axions are localized within radius R ⊕ , we arrive to the following estimate for the average axion energy density inside this volume which is amazingly close to the average dark matter density today ρ DM 0.3 GeV cm 3 . The eq. (44) should be viewed as the order of magnitude estimate at the very best. The main uncertainty here is that the trapped axions are not distributed uniformly, as assumed in (44). Instead, they are obviously distributed in a highly nontrivial way determined by the position of the nugget when emission occurs (in deep underground) and the direction of the velocity at the moment of emission. Though the estimate (44) is rough, it is also very promising as it suggests that the galactic axion density and the axion density produced by the AQN mechanism could be the same order of magnitude. In addition, the contribution of the galactic axions to the dark matter density scales as m −7/6 a as mentioned in Introduction. It should be contrasted with estimate (44) which is not sensitive to the value of the axion mass m a as a result of very generic feature expressed as (3) of the AQN framework.
Furthermore, the distinct feature of the AQN trapped axions is very large wave length λ a = (m a v a ) −1 as the typical trapped velocity v a ≤ v trapped ⊕ is much smaller than a typical galactic DM velocity ∼ 10 −3 c according to (18). This unique feature of the trapped axions might be the "smoking gun" leading to their discovery. We conclude this section on this optimistic note.

Conclusion and future directions
This work represents a natural generalization of the previous studies [1] to properly account for the production of the low energy axions when the AQNs get annihilated in the Sun or Earth and emit axions with v ≤ v trapped in the solar in the deep Earth's underground. This portion of the non-relativistic axions is extremely tiny as shown on Fig. 1b. However, the effect is drastically enhanced as argued in Section 5 due to accumulation of these axions during entire life of the solar system, i.e. for ∼4.5 billion of years. The corresponding estimates represent the main results of the present studies 4 .
We shall not repeat and discuss here a large number of estimates presented in Section 5. Instead, we focus on a single formula (44) describing the energy density of the trapped axions ρ axions ⊕ . We think this estimate has a huge discovery potential because ρ axions ⊕ is relatively large and comparable with the average galactic dark matter density today ρ DM 0.3 GeV cm 3 . What is more important is that the spectral features of the trapped axions are very distinct from conventional galactic axions because the typical trapped velocity v a ≤ v trapped ⊕ is much smaller than a typical galactic DM velocity ∼ 10 −3 c according to (18). Therefore, the typical wave length λ a = (m a v a ) −1 of these axions is much longer in comparison with galactic axions. This unique feature makes the trapped axions are very distinct from conventional galactic axions. These axions obviously can be easily discriminated from anything else. The discovery of such axions would be a "smoking gun" for the entire AQN proposal unifying the DM and baryogenesis (separation of charges) problems.  , δ) vs v a /c in 3D case. The key element here is weak sensitivity to parameter δ supporting our claim that the computational scheme developed in this work is reliable and robust.
where R eff and δ R are functions of a tunable parameter R trans One can explicitly check this approximate solution (A.3) is continuous and first order differentiable. Also, it is precisely the exact solution in the near-field limit r ∼ R 0 and the far-field limit r m −1 a . Hence, the only unknown part the solution is the "transition" regime between these two limits, where we introduce a tunable parameter R trans to account for this type of error source. We will keep this parameter in the following calculations, so it serves as a probe to test whether the final result is sensitive to our crude approach in the transition regime. As we will see, the final result is robust as it is not sensitive to the tuning of R trans . Lastly, instead of using R trans directly, it is more convenient to define a simple parameter δ ≡ m a δ R which roughly varies from 0 to 1. As we will see, δ is the only parameter entering the final result.
We are now ready to compute the excitations χ(t, z) in the time dependent background. These excitations will be eventually identified as the axions emitted by the axions DW. To achieve this task we expand φ(t, z) = φ w (z − R 0 ) + χ(t, z), which gives where L 2 is a linear differential operator of the second order, The corresponding equation of motion is therefore To look for the initial conditions, we now want to describe the emission of axions in one cycle of oscillation. As mentioned in Sec. 4.2, annihilation of baryon charge results in oscillations of domain wall. Assuming the oscillation is approximately adiabatic, it is sufficient to only analyze the first half of an oscillation -say, the "contraction period"where the domain wall shrinks from R 0 to a slightly smaller size R 0 − ∆R. We assumed the rest half of the cycle, the "expansion period", is just the time-reversed and produces an equivalent contribution. We may write down such initial conditions as where t osc denotes the period of one full oscillation. The excitation modes in condition (A.8b) is unknown and depends on the conversion rate from excitation modes to freely propagating axions. In terms of χ, the initial conditions (A.8) imply χ(0, r) = 0 (A.9a) where we introduce a free parameter η(θ, ϕ) which may be interpreted as the "amplitude of efficiency" of the conversion rate from excitations to free axions, so η must vary between 0 to 1. However, η here may be also interpreted as a correction term like δ in the approximate solution (A.3) within the transition regime R 0 r m −1 a , so η can be greater than 1 in general. Nonetheless, we will expect η ∼ 1 and will treat it as a normalization factor regarding to the luminosity. And in general, η can be expanded by partial waves (A.10) If we assume a good spherical symmetry preserves during the most period of the annihilation process of AQN, then η 00 will be the dominant contribution and η 10 be the next order correction.
To solve for the excitation mode, it is convenient to write χ in terms of some normalized basis. The expansion for free wave is conventionally where j l (x) is the spherical Bessel function, and we have implicitly used two orthogonalities ∞ 0 dr r 2 j l (pr) j l (qr) = π 2p 2 δ(p − q), Note that L 2 is diagonal in basis of χ plm where K (l) p,q is a coefficient defined as for simplicity of calculation. In Appendix B we can show K (l) p,q δ(p − q) = δ(p − q). Then Eq. (A.7) is simplified to which clearly has solution a plm (t) = b plm sin E a t, (A. 16) following the initial condition (A.9a), where b plm is an time-independent coefficient to be determined. To find b plm , we should impose the second initial condition (A.9b) which implies Note that only the last term in the curly bracket is dominant because R trans m −1 a largely suppresses the first two terms. 7  where we have also drop R 0 in the hyperbolic secant function because it is of order R trans . This integral can be evaluated precisely if we expand the hyperbolic secant as where 2 F 1 (a, b, c; z) is the Gauss hypergeometric function, and f (a, b, c; z) is defined to be the regularized version of 2 F 1 (a, b, c, z) in a conventional way, see Refs. [59,60] and recent article [61]. As discussed in Sec. 5, we are especially interested in the non-relativistic domain, in this limit we have Then, the total radiation energy E rad of the domain wall is obviously (A.24) More generally, assuming now the flux is produced within a "cavity of radiation" V rad , the density of radiation energy (per unit volume) is therefore E rad /V rad . Then the net flux Φ rad going through the boundary of the cavity is clearly A few comments should be made regarding to the magnitude of R rad . First, V rad is defined as the cavity where radiation happens, so R rad ∆R in 1D case where thin-wall approximation is assumed. However, in 3D R rad may extend to order of R 0 or even m −1 a . More generally, it is reasonable to conjecture R rad can depend on the angular momentum l. It is clear that to compute or even estimate the order of R rad is very difficult. Thus, we should not bother the details of R rad , but rather treat it as a tunable normalization parameter and maybe absorb it into η lm if applicable.
Lastly, we also express the spectra as a function of flux velocity Up to this point, we note that K (0) p,q is always finite for any positive p, q > 0. Now, if we multiply both sides by δ(p − q) where we know the integral is quickly dominant by the the term L 0 dr · 1, so that the fraction in the limit L → ∞ gives trivial result. Now, we want to show that K (l) p,q δ(p − q) = K (l+1) p,q δ(p − q) for all l = 0, 1, 2.... First, let us see that Note that in the last line if we set p = q, the second term must vanish. It is because the numerator is obviously finite, while the denominator tends to infinity in the large L limit as indicated by Eq. A.12a. Thus, we conclude for all l = 0, 1, 2, 3... as expected.