Gravitationally trapped axions on Earth

TL;DR

This paper proposes a new mechanism for axion production on Earth via gravitational trapping, based on the AQN dark matter model, leading to a potential axion halo around Earth that could be detected by existing experiments.

Contribution

It introduces a novel Earth-based axion production mechanism using the AQN model, independent of axion mass or initial misalignment, and estimates the resulting axion flux and density profile.

Findings

Axion flux on Earth's surface estimated through simulations.

Trapped axions form a density halo around Earth.

Potential observability with current axion detection experiments.

Abstract

We advocate for the idea that there is a fundamentally new mechanism for axion production on Earth, as recently suggested in Fischer et al. (2018) and Liang & Zhitnitsky (2018). We specifically focus on production of axions within Earth, with low velocities such that they will be trapped in the gravitational field. Our computations are based on the so-called Axion Quark Nugget (AQN) dark matter model, which was originally invented to explain the similarity of the dark and visible cosmological matter densities. This occurs in the model irrespective of the axion mass $m_\mathrm{a}$ or initial misalignment angle $\theta_0$. Annihilation of antimatter AQNs with visible matter inevitably produce axions when AQNs hit Earth. The emission rate of axions with velocities below escape velocity is very tiny compared to the overall emission, however these axions will be accumulated over the 4.5 billion year life time of the Earth, which greatly enhances the discovery potential. We perform numerical simulations with a realistically modeled incoming AQN velocity and mass distribution, and explore how AQNs interact as they travel through the interior of the Earth. We use this to estimate the axion flux on the surface of the Earth, the velocity-spectral features of trapped axions, the typical annihilation pattern of AQN, and the density profile of the axion halo around the Earth. Knowledge of these properties is necessary to make predictions for the observability of trapped axions using CAST, ADMX, MADMAX, CULTASK, ORPHEUS, ARIADNE, CASPEr, ABRACADABRA, QUAX, DM Radio.

Figures

Tables

Table 1. Table 1 : Values of the mean baryon charge ⟨ B ⟩ delimited-⟨⟩ 𝐵 \langle B\rangle for different parameters of the AQN mass-distribution function.

$(B_{\min },\alpha )$	2.5	2.0	(1.2, 2.5)
${10}^{23}$	$2.99\times {10}^{23}$	$1.15\times {10}^{24}$	$4.25\times {10}^{25}$
$3\times {10}^{24}$	$8.84\times {10}^{24}$	$2.43\times {10}^{25}$	$1.05\times {10}^{26}$

Table 2. Table 2 : Our model for the density structure of the Earth. We consider the Earth to be made from 5 distinct layers, each of which we model as a constant-density shell, with density ρ i subscript 𝜌 𝑖 \rho_{i} . We show the outer radius for each shell, r i subscript 𝑟 𝑖 r_{i} , as well as the thickness of the shell.

Label	Layer	Thickness [km]	$r_i$ [km]	${\rho }_i$ [$g{\mathrm{cm}}^{-3}$]
1	Inner core	1221	1221	12.95
2	Outer core	2259	3480	11.05
3	Lower mantle	2171	5651	5.00
4	Upper mantle	720	6371	3.90
5	Crust	30	6401	2.55

Table 3. Table 3 : Summary of some results. Our baseline model is B min = 3 × 10 24 subscript 𝐵 min 3 superscript 10 24 B_{\rm min}=3\times 10^{24} , ϵ = 1 italic-ϵ 1 \epsilon=1 unless otherwise specified. For each model we show: ⟨ Δ ​ m AQN ⟩ delimited-⟨⟩ Δ subscript 𝑚 AQN \langle\Delta m_{\rm AQN}\rangle , the change in mass of a typical antimatter AQN as a result of interaction with the Earth; ⟨ Δ ​ B ⟩ / ⟨ B ⟩ delimited-⟨⟩ Δ 𝐵 delimited-⟨⟩ 𝐵 {\langle\Delta B\rangle}/{\langle B\rangle} , the typical loss of baryon charge of a typical AQN (also the fractional mass loss); ρ a ​ ( R ⊕ ) subscript 𝜌 a subscript 𝑅 direct-sum \rho_{\mathrm{a}}(R_{\oplus}) , the magnitude of the axion-halo density profile at the surface of the Earth after 4.5 Gyrs of accumulation.

$⟨B⟩$	$\alpha$	Other parameters	$⟨\Delta m_{\mathrm{AQN}}⟩$ [kg]	$⟨\Delta B⟩/⟨B⟩$	${\rho }_a\left(R_{\oplus }\right)\left[\mathrm{GeV}{\mathrm{cm}}^{-3}\right]$
$8.84\times {10}^{24}$	2.5	–	$5.01\times {10}^{-3}$	$33.9\%$	$2.36\times {10}^{-4}$
$8.84\times {10}^{24}$	2.5	$𝝁=-\mu \widehat{𝐳}$	$4.96\times {10}^{-3}$	$33.5\%$	$2.29\times {10}^{-4}$
$8.84\times {10}^{24}$	2.5	Solar gravitation666We consider an additional $42.1\mathrm{km}{s}^{-1}$ AQN impact velocity from AQN falling to Earth from infinity in the gravitational well of the Sun, this gives an additional velocity to the AQN that is always in the direction of travel. See Appendix F for more details	$5.01\times {10}^{-3}$	$33.9\%$	$2.74\times {10}^{-4}$
$8.84\times {10}^{24}$	2.5	$ϵ=3$	$1.24\times {10}^{-2}$	$84.0\%$	$3.73\times {10}^{-4}$
$2.43\times {10}^{25}$	2.0	–	$8.34\times {10}^{-3}$	$20.5\%$	$1.52\times {10}^{-4}$
$2.43\times {10}^{25}$	2.0	$𝝁=-\mu \widehat{𝐳}$	$8.15\times {10}^{-3}$	$20.1\%$	$1.43\times {10}^{-4}$
$4.25\times {10}^{25}$	(1.2, 2.5)	$B_{\min }={10}^{23}$	$1.08\times {10}^{-2}$	$15.1\%$	$1.17\times {10}^{-4}$
$1.05\times {10}^{26}$	(1.2, 2.5)	–	$2.57\times {10}^{-2}$	$14.6\%$	$1.14\times {10}^{-4}$
$1.05\times {10}^{26}$	(1.2, 2.5)	$𝝁=-\mu \widehat{𝐳}$	$2.48\times {10}^{-2}$	$14.1\%$	$1.07\times {10}^{-4}$

Table 4. Table 4 : Probability of seeing zero event in IceCube over 10 years

$\eta$	$⟨B⟩$	Hit number	Prob($k$=0) [%]	Significance
	$1.6\times {10}^{24}$	8.148	0.02893	$3.44\sigma$
1	$5.0\times {10}^{24}$	2.607	7.37260	$1.45\sigma$
	$1.0\times {10}^{25}$	1.304	27.15253	$0.61\sigma$
	$4.7\times {10}^{24}$	8.321	0.02432	$3.49\sigma$
3	$5.0\times {10}^{24}$	7.822	0.04007	$3.35\sigma$
	$1.0\times {10}^{25}$	3.911	2.00185	$2.05\sigma$

Equations

$$N_{\oplus ,4.5Gyr}^{(\mathrm{trap})}$$

$$\sim 10^{40}\left(\frac{{\xi }_{\oplus }}{10^{-17}}\right)\left(\frac{\Delta B}{B}\right)\left(\frac{10^{-5}\mathrm{eV}}{m_{\mathrm{a}}}\right),$$

$$E_{\oplus ,4.5Gyr}^{(\mathrm{trap})}$$

$$\sim 10^{35}\left(\frac{{\xi }_{\oplus }}{10^{-17}}\right)\left(\frac{\Delta B}{B}\right)\mathrm{eV}.$$

$$\Delta M_{\oplus }^{(\mathrm{trap})}\sim 0.1\left(\frac{{\xi }_{\oplus }}{10^{-17}}\right)\left(\frac{\Delta B}{B}\right)\mathrm{kg},$$

$${\rho }_{a,\oplus }^{(\mathrm{avg})}\sim \frac{E_{\oplus ,4.5Gyr}^{(\mathrm{trap})}}{4\pi R_{\oplus }^3/3}\sim 0.1\left(\frac{{\xi }_{\oplus }}{10^{-17}}\right)\left(\frac{\Delta B}{B}\right)\frac{\mathrm{GeV}}{\mathrm{c}{\mathrm{m}}^3}.$$

$$B_{\text{universe}}=0$$

$$|B{|}_{\text{DM}}$$

$${\Omega }_{\mathrm{DM}}\sim {\Omega }_{\mathrm{visible}}.$$

$$|B_{\mathrm{visible}}|:|B_{\mathrm{nuggets}}|:|B_{\mathrm{antinuggets}}|\simeq 1:2:3$$

$$\Phi =n_{\mathrm{N}}{v}_{\mathrm{N}}\approx \frac{{\rho }_{\mathrm{DM}}{v}_{\mathrm{N}}}{M_{\mathrm{N}}}\approx 1{\mathrm{km}}^{-2}{\mathrm{yr}}^{-1}{\left(\frac{⟨B⟩}{10^{24}}\right)}^{-1},$$

$$⟨B⟩>3\times 10^{24}[\text{non-detection constraint}],$$

$$\Phi \sim R^2\int \mathrm{d}\Omega \mathrm{d}l[n_{\mathrm{visible}}(l)\cdot n_{\mathrm{DM}}(l)]\sim \frac{1}{⟨B{⟩}^{1/3}},$$

$$\mathrm{d}N\sim B^{-\alpha }\mathrm{d}B\sim W^{-\alpha }\mathrm{d}W$$

$$⟨B⟩={\int }_{B_{\min }}^{B_{\max }}\mathrm{d}BBf(B),f(B)\propto B^{-\alpha }$$

$$\alpha =2.5,2.0,\mathrm{or}\{\begin{array}{rlr}& 1.2& B\lesssim 3\times 10^{26}\\ & 2.5& B\gtrsim 3\times 10^{26}.\end{array}$$

$$\varphi (R_0)={\varphi }_w(R_0)+\chi$$

$$\frac{\mathrm{d}E}{\mathrm{d}s}=\frac{1}{v}\frac{\mathrm{d}E}{\mathrm{d}t}=-\sigma \rho v^2$$

$$\frac{\mathrm{d}E}{\mathrm{d}t}=mv\frac{\mathrm{d}v}{\mathrm{d}t}+\frac{1}{2}{v}^2\frac{\mathrm{d}m}{\mathrm{d}t}=mv\frac{\mathrm{d}v}{\mathrm{d}t}-\frac{1}{2}\sigma \rho v^3,$$

$$\dot{m}=-\sigma \rho v.$$

$$m\frac{\mathrm{d}v}{\mathrm{d}t}=-\frac{1}{2}\sigma \rho v^2\equiv -\frac{1}{2}\pi (ϵR{)}^2\rho v^2.$$

$$\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}=\mathbf{v};r=|\mathbf{r}|;v=|\mathbf{v}|;$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-\frac{1}{2}\pi (ϵR{)}^2\frac{\rho (r)}{m}{v}^2\hat{\mathbf{v}}-\frac{G{M}_{\mathrm{eff}}(r)}{r^2}\hat{\mathbf{r}},$$

$$R={\left(\frac{3m}{4\pi {\rho }_{\mathrm{n}}}\right)}^{1/3}\simeq 1.045\times 10^{-13}{B}^{1/3}\mathrm{cm},$$

$$M_{\mathrm{eff}}(r)=\sum _{i=1}^j\frac{4\pi }{3}{\rho }_i(r_i^3-r_{i-1}^3)+\frac{4\pi }{3}{\rho }_{j+1}(r^3-r_j^3),$$

$$(r_j\le r<r_{j+1},\mathrm{with}j=1,\mathrm{2...5}\mathrm{and}{r}_0\equiv 0).$$

$${\rho }_{\mathrm{rest}}({\tilde{v}}_{\mathrm{a}})$$

$$\simeq \frac{{\tilde{v}}_{\mathrm{a}}^3}{N(\delta )}{\left(\frac{{\tilde{E}}_{\mathrm{a}}}{m_{\mathrm{a}}}\right)}^6|H_0(0,\delta ){|}^2[1+\mathcal{O}({\tilde{v}}_{\mathrm{a}}^2)],$$

$${\tilde{\mathbf{v}}}_{\mathrm{a}}={\mathbf{v}}_{\mathrm{a}}-{\mathbf{v}}_{\mathrm{AQN}};\tilde{\mathbf{p}}=\mathbf{p}-m_{\mathrm{a}}{\mathbf{v}}_{\mathrm{AQN}};$$

$${\tilde{v}}_{\mathrm{a}}=\sqrt{v_{\mathrm{a}}^2+v_{\mathrm{AQN}}^2-2{\mathbf{v}}_{\mathrm{a}}\cdot {\mathbf{v}}_{\mathrm{AQN}}};$$

$$\tilde{p}=\sqrt{p^2+m_{\mathrm{a}}^2{v}_{\mathrm{AQN}}^2-2m_{\mathrm{a}}\mathbf{p}\cdot {\mathbf{v}}_{\mathrm{AQN}}};$$

$${\tilde{E}}_{\mathrm{a}}=\sqrt{E_{\mathrm{a}}^2+m_{\mathrm{a}}^2{v}_{\mathrm{AQN}}^2-2m_{\mathrm{a}}\mathbf{p}\cdot {\mathbf{v}}_{\mathrm{AQN}}}.$$

$$1={\int }_{B^3}{\mathrm{d}}^3{\tilde{\mathbf{v}}}_{\mathrm{a}}\frac{{\rho }_{\mathrm{rest}}({\tilde{v}}_{\mathrm{a}})}{4\pi {\tilde{v}}_{\mathrm{a}}^2}={\int }_{B^3}{\mathrm{d}}^3{\mathbf{v}}_{\mathrm{a}}\frac{{\rho }_{\mathrm{rest}}[{\tilde{v}}_{\mathrm{a}}(v_{\mathrm{a}})]}{4\pi ({\tilde{v}}_{\mathrm{a}}(v_{\mathrm{a}}){)}^2},$$

$${\rho }_{\mathrm{obs}}(v_{\mathrm{a}})=\int \mathrm{d}\Omega v_{\mathrm{a}}^2\frac{{\rho }_{\mathrm{rest}}({\tilde{v}}_{\mathrm{a}})}{4\pi {\tilde{v}}_{\mathrm{a}}^2},$$

$${\rho }_{\mathrm{obs}}(v_{\mathrm{a}})=\frac{1}{2}{\int }_{-1}^1\mathrm{d}u\frac{v_{\mathrm{a}}^2{\tilde{v}}_{\mathrm{a}}}{N(\delta )}{\left(\frac{{\tilde{E}}_{\mathrm{a}}}{m_{\mathrm{a}}}\right)}^6|H_0(\tilde{p},\delta ){|}^2,$$

Full text

I Introduction

The Peccei-Quinn mechanism, accompanied by axions, remains the most compelling resolution of the strong charge-parity (CP) problem [orignal papers: 3, 4, 5, 6, 7, 8, 9][recent reviews: 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. In this model, conventional dark-matter (DM) axions are produced either by the misalignment mechanism [20, 21, 22], when the cosmological field $\theta(t)$ oscillates and emits cold axions before it settles at a minimum, or via the decay of topological objects [23, 24, 25, 26, 27, 28, 29]. There are uncertainties in estimates of the axion abundance from these two channels,111 According to recent computations [29] the axion contribution to $\Omega_{\rm DM}$ as a result of decay of topological objects can saturate the observed DM density today if the axion mass is in the range $m_{\mathrm{a}}=(2.62\pm 0.34)10^{-5}{\rm eV}$, while earlier estimates suggest that saturation occurs at a larger axion mass. There are additional uncertainties in this result, and we refer to the original studies [28] on this matter. One should also emphasize that the computations [23, 24, 25, 26, 27, 28, 29] have been performed with assumption that Peccei-Quinn symmetry was broken after inflation., and we refer the reader to [23, 24, 25, 26, 27, 28, 29] for a full discussion. In both of the traditional production mechanisms axions are produced as non-relativistic particles with typical $v_{\rm axion}/c\sim 10^{-3}$, and their contribution to the cosmological DM density scales as $\Omega_{\rm axion}\sim m_{\mathrm{a}}^{-7/6}$. This scaling implies that the axion mass must be tuned to $m_{\mathrm{a}}\simeq 10^{-5}$ eV in order to produce the observed cosmological DM density today. Higher axion masses will contribute very little to $\Omega_{\rm DM}$ and lower axion masses will over-saturate the density, resulting in a closed Universe which would be in strong conflict with cosmological data [30]. Cavity-type experiments have the potential to discover these non-relativistic axions.

Axions may also be produced via the Primakoff effect in stellar plasmas at high temperature [31]. These axions are ultra-relativistic; with a typical average energy of axions emitted by the Sun of $\langle E\rangle=4.2$ keV [32]. Searches for Solar axions are based on helioscope instruments like CAST (CERN Axion Search Telescope) [32].

Recent works [1, 2] have suggested a fundamentally novel mechanism for axion production in planets and stars, with a mechanism rooted in the so-called axion quark nugget (AQN) DM model [33], similar to the original quark nugget model by Witten [34] (see [35] for a review). AQN are “cosmologically dark”, not because of the weakness of their interactions, but due to their small cross-section-to-mass ratio, which scales down many observable consequences of an otherwise strongly-interacting DM candidate. There are two additional elements in our AQN model compared to [34, 35]. First, there is an additional stabilization factor for the nuggets provided by axion domain walls that are produced in great quantity during the quantum chromodynamic (QCD) transition and that help to alleviate a number of problems with the original nugget model.222In particular, a first-order phase transition is not required as the axion domain wall plays the role of the squeezer. Another problem with [34, 35] is that nuggets will likely evaporate on a Hubble time-scale even. For the AQN model this argument is not applicable because the vacuum-ground-state energies inside (color-superconducting phase) and outside (hadronic phase) the nugget are drastically different. Therefore, these two systems can coexist only in the presence of an external pressure, provided by the axion domain wall. This contrasts with the original model [34, 35], which must be stable at zero external pressure. Another feature of AQN is that nuggets can be made of matter as well as antimatter during the QCD transition. This element of the model completely changes the AQN framework [33] because the DM density, $\Omega_{\rm DM}$, and the baryonic matter density, $\Omega_{\rm visible}$, automatically become similar to each other, i.e. $\Omega_{\rm DM}\sim\Omega_{\rm visible}$, because they have the same QCD origin in that the same generating mechanism ensures they are of the same order of magnitude without any fine tuning as they both proportional to one and the same fundamental $\Lambda_{\rm QCD}$ scale. The existence of both matter and antimatter AQNs explains the asymmetry between visible components made of matter and antimatter. In other words, in the AQN framework the conventional idea of baryogenesis is replaced by the charge separation mechanism when the net baryon charge remains zero in the entire Universe. In short, the AQN model explains two fundamental questions at the same time: the nature of dark matter and the asymmetry between matter and antimatter. This is irrespective of specific details of the model, such as the axion mass $m_{a}$ or misalignment angle $\theta_{0}$.

The presence of bound axions in the construction of AQNs plays a crucial role for the present work: Annihilation events will always happen when an AQN propagates in a medium, and this inevitably produces propagating axions along with the conventional products of baryon-antibaryon annihilation. Emitted axions will be released with a velocity spectrum, with a typical value $v_{\rm axion}^{\rm AQN}\simeq 0.5c$. This can be contrasted with conventional galactic non-relativistic axions $v_{\rm axion}\sim 10^{-3}c$ and Solar ultra-relativistic axions with typical energies $\langle E\rangle=4.2$ keV.

Some of the axions emitted when AQNs annihilate with Earth will have very small velocities, below the escape velocity of Earth ($v\leq v_{\rm esc}\simeq 11\,\mathrm{km}\,\mathrm{s}^{-1}$). These axions will be gravitationally trapped and represent the main subject of this work. The idea that axions can be trapped in a gravitational field dates back to 2003 paper by DiLella and Zioutas [36]. The present work is a specific analysis of the gravitationally bound axions within the AQN framework, when the axions are produced in the interior of the Earth, as suggested in [1, 2]. The emission spectrum of axions from the AQN production mechanism is studied in Ref. [2]. An important conclusion of this work is that a tiny fraction $\sim 10^{-17}$ of the number of axions emitted within the Earth will be low-velocity ($v\leq v_{\rm esc}$), and thus will be trapped by the Earth. Over 4.5 billion years of accumulation these bound axions may become dense enough for detection by instruments such as ADMX and ADMX-HF [37], upgraded CAST [38], ORPHEUS [39], MADMAX [40], ORGAN [41, 42], CULTASK [43], ARIADNE [44], CASPEr [45], ABRACADABRA [46], QUAX [47], X3 [48], DM Radio [49], TOORAD [50], axion plasma haloscope[51].

In [2] a number of order-of-magnitude estimates were made regarding properties of the trapped-axions. The goal of the present work is to upgrade these estimates by performing detailed numerical simulations of axion production within the Earth. Questions that we are interested in answering are:

What fraction of an AQN is annihilated when it passes through the interior of the Earth? 2. 2.

What is the distribution of where annihilation events occur within the Earth? 3. 3.

How is the trapped-axion density distributed around the Earth? 4. 4.

What is the trapped-axion flux through the surface of the Earth, where they might be detected? 5. 5.

What are the axion spectral features on the surface of the Earth? 6. 6.

How do answers to the above questions depend on the AQN mass distribution and other model-dependent features?

Assuming that AQNs saturate the cosmological DM density, and defining $\Delta B/B$ as the fraction of a typical AQN that is annihilated when passing through the interior of the Earth, one can estimate the total number of gravitationally trapped axions accumulated by the Earth over 4.5 billion years as follows [2]:

[TABLE]

where $\xi_{\oplus}\sim 10^{-17}$ is the fraction of axions trapped in the gravitational field of the Earth compared to the total axion flux, estimated in [1, 2]. The corresponding accumulated energy is

[TABLE]

The trapped energy (2) can be expressed in terms of the trapped-axion mass $\Delta M_{\oplus}$ accumulated by the Earth

[TABLE]

which represents a tiny fraction of the mass of the Earth: $M_{\oplus}\simeq 5.9\times 10^{24}$ kg.

The key unknown element in estimates (1) and (2) is the parameter $\Delta B/B$: the fraction of the mass of a typical AQN that is annihilated within Earth. This parameter is difficult to estimate without numerical simulations because it depends on the trajectory of the AQN through Earth, which itself is affected by the the incoming AQN velocity distribution, the gravity of the Earth, the details of the density distribution within the Earth, as well as the axion-emission velocity spectrum.

If we make the assumption that axions are uniformally distributed within the Earth, we can estimate the average axion energy density within the Earth [2]:

[TABLE]

However, it is clear that the distribution of axions will not be uniform and instead, may be non-trivial, in a way determined by the position of the AQN when emission occurs and the direction of AQN velocity at the moment of emission. Evidently, the crude estimate of (4) is inadequate for planning experiments. To address this (and items 3, 4 and 5 from the list above) we use Monte-Carlo simulations to allow us to evaluate (4) precisely, as a function of a position from the Earth. This information is crucial for planning axion search experiments, which require a value for $\rho_{a,\oplus}(r=R_{\oplus})$ along with the spectral properties of the axions on the surface $r=R_{\oplus}$. Specifically, we simulate the entire process responsible for the generation of the axion halo around Earth. We start by simulating the propagation of AQNs through Earth, including the AQN annihilation processes. We then simulate the orbits of the low-velocity axions that are emitted by the AQN and calculate the accumulation over 4.5 Gyr. This axion-simulation process depends on the position of the AQN within the Earth where annihilation occurs and the axions are emitted, the velocity of the AQN at the moment of emission, and the velocities and directions of the emitted axions. This procedure allows us to compute $\rho_{a,\oplus}(r=R_{\oplus})$.

The present work is a natural continuation of [1, 2]. To avoid repetition we refer the readers to [2] for the derivation details of axion emission spectrum in the rest frame of the AQN and to [1] for an overview of the AQN model, including motivation, consequences and predictions. We note that the AQN model is consistent with all available cosmological, astrophysical, satellite and ground-based constraints. While the model was invented to explain the observed relation $\Omega_{\rm DM}\sim\Omega_{\rm visible}$, it may also explain a number of other observed phenomena, such as the excess of diffuse galactic emission in different frequency bands (including the 511 keV line) as reviewed in [52]. AQNs may also offer a resolution to the so-called “Primordial Lithium Puzzle” [53] and to the “The Solar Corona Mystery” [54, 55] and may also explain the recent EDGES observation [56], which is in some tension with the standard cosmological model. It is the same set of physical AQN parameters that address the aforementioned phenomena that we will adopt for the present work. What differentiates the present work is that we are advocating a direct search for the gravitationally trapped axions: a direct manifestation of the AQN model. An observation of axions with very distinct spectral properties compared to those predicted from conventional galactic axions would be a smoking gun for the AQN framework and this may then answer a fundamental question on the nature of DM.

This paper is organized as follows: In Sec. II, we present a brief overview of AQN model: its accomplishments, generic consequences and observational constraints. We also overview the axion emission spectrum, which plays an important role for the present study. In Sec. III, we develop analytical equations to describe the annihilation of AQNs impacting the Earth and the emission spectrum of axions in laboratory frame. The corresponding numerical simulations are described in Sec. IV, with results presented in Sec. V. Finally, we conclude with some thoughts on possible future developments in Sec. VI.

II Overview of the AQN model

II.1 Basic and generic features of the AQN framework

The original AQN model was developed by [33] as a natural explanation to two of the most fundamental problems in cosmology: namely, the nature of DM and the asymmetry between matter and antimatter. To be more specific, “baryogenesis” in this model is replaced by “baryon charge separation” such that the DM density $\Omega_{\rm DM}$ and the baryonic matter density $\Omega_{\rm visible}$ become similar to each other. Indeed, the conservation of the baryon charge implies

[TABLE]

where $B_{\text{universe}}=0$ is the total number of baryons in the universe, $|B|_{\text{DM}}$ counts total number of baryons and antibaryons hidden in both the nuggets and antinuggets that make up the DM. $B_{\text{visible}}$ is the total number of residual “visible” baryons (regular matter). The energy per baryon charge is approximately the same for AQNs and visible matter, because both types of matter are formed during the same QCD transition, and both are proportional to the same dimensional parameter of proton mass, $m_{\mathrm{p}}$, which implies that

[TABLE]

In other words, the nature of DM and the problem of the Universal asymmetry between matter and antimatter become two sides of the same coin in this framework. As has been argued in refs. [57, 58, 59] the relation (6) is a very generic outcome of the AQN framework, and it is not sensitive to specific details of the model.

To be more specific, from the observed ratio $\Omega_{\rm DM}\simeq 5\cdot\Omega_{\rm visible}$, one can infer that number of antinuggets is larger than the number of nuggets by a factor of $\sim 3/2$ at the conclusion of AQN formation. This results in a Universal matter budget containing baryons, quark nuggets, and antiquark nuggets in an approximate ratio,

[TABLE]

with a net baryonic charge of zero.

The axion field, which is a crucial element in the present work, plays an essential role during AQN formation in two ways: First, the existence of a coherent ${\cal CP}$-odd axion field $\theta$ is the fundamental mechanism to break the global symmetry of baryon charge between nuggets and antinuggets. This difference is always an order unity effect, irrespective of the parameters of the theory, as argued in [57, 58]. This is precisely the reason why the resulting visible and DM densities must be the same order of magnitude (6) in this framework, as they are both proportional to the same fundamental $\Lambda_{\rm QCD}$ scale, and they both originated at the same QCD epoch. If these processes are not fundamentally related, the two components $\Omega_{\rm DM}$ and $\Omega_{\rm visible}$ could easily exist at vastly different scales. Second, the so-called ${N}_{\rm DW}=1$ axion domain walls provide sufficient pressure to accumulate baryon charge inside the walls, and ultimately lead to formation of stable AQNs in color-superconducting phase of size $R\sim m_{\mathrm{a}}^{-1}$, discussed in detail in [57].

It is precisely this energy, hidden in form of the bound axion domain wall since the moment of formation, that can be released in form of free-propagating axions when the baryon charge from antimatter AQNs are annihilated with matter in the interior of the Earth. The present work is devoted to an analysis of the fate of these, previously hidden, axions. In the following subsection we briefly review the basic features of the AQN model relevant for the present work.

II.2 Observational constraints on $\langle B\rangle$

An AQN is analogous to a gigantic atom, with a central “nucleus” made out of quark condensate in the color-superconducting phase enclosed by a thick domain wall made out of axion and $\eta^{\prime}$ meson. Because the quark core is in general not neutral, it is surrounded by a cloud of leptons in what is known as the “electrosphere”. Once an AQN made of antimatter impacts upon matter annihilation will occur. As a consequence, the AQN will loose mass and emits axions.

The flux of AQNs on Earth is determined by the AQN mass $M_{N}$, or equivalently by the baryon charge $B$, related via $M_{N}\simeq m_{\mathrm{p}}B$. Assuming that AQNs saturate the DM density the AQN flux can be approximated as

[TABLE]

so that direct detection experiments impose lower limits on the value of $\langle B\rangle$ for the baryon-charge distribution of the AQN.

The strongest direct detection limit on AQN is likely set by the non-detection of an AQN flux [60] by the IceCube Observatory, which limits the AQN flux to $\Phi\lesssim 1$ km*-2yr-1*. If we take the local DM mass density to be $\rho\approx 0.3$ GeV/cm3 and assume that AQN make up all DM we may translate the flux constraint obtained from IceCube into a lower limit with $3.5\sigma$ confidence on the mean baryon number of the AQN distribution of

[TABLE]

see Appendix A for details. A similar constraint also comes from the Antarctic Impulsive Transient Antenna (ANITA) [61].

Limits on $\langle B\rangle$ may also be obtained from indirect astrophysical and cosmological observations. Photon production during matter-AQN interactions along a given line of sight implies an observable flux:

[TABLE]

where $R\sim B^{1/3}$ is the typical size of the AQN, which determines the effective cross section of interaction between DM and visible matter, while $n_{\rm DM}\sim B^{-1}$ is the number density of the AQNs. Thus, astrophysical constraints impose a lower bound on $\langle B\rangle$. The annihilation contribution to the galactic spectrum has been analyzed for different frequencies from a radio bands to X-ray and $\gamma$-ray photons. In particular, if $\langle B\rangle\sim 10^{25}$ the AQN model could offer an explanation for the observed 511 keV emission features (including the width, morphology, $3\gamma$ continuum spectrum, etc). In all other cases, where computations have been performed in different frequency bands, the predicted emission is consistent with observations. In fact, in some cases it could explain some, or even the dominant, portion of the observed excess of radiation. All these emissions in different frequency bands are expressed in terms of one and the same integral (10), and therefore, the relative intensities are completely determined by internal structure of the AQNs, which is described by conventional nuclear physics and QED. For further details see the original works [62, 63, 64, 65, 66, 67, 68, 69] with specific computations in different frequency bands in galactic radiation, and a short overview [52].

II.3 Baryon charge distribution of AQN model

Another AQN-connected observation might be related to the so-called “Solar corona heating mystery”. It has been a long-standing puzzle that the corona has a temperature $T\simeq 10^{6}$K that 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. This puzzle may find a resolution within the AQN framework, as recently argued in [54, 70, 55]. The basic idea is that nanoflares, postulated by Parker to explain the Solar corona heating, can be identified with the AQN annihilation events. In [55] it was shown that most of the Solar-incident AQN are are annihilated in the transition region between the photosphere and corona, specifically 2000 km above the photosphere, where drastic changes in temperature and pressure are known to occur. Solar flares are sensitive to the distribution $\mathrm{d}N/\mathrm{d}B$ because the AQN framework allows us to connect the nanoflare-size distribution, $\mathrm{d}N/\mathrm{d}W$, to the baryon-charge distribution. This is possible due to $\mathrm{d}N/\mathrm{d}W$ being modelled via magnetic-hydro-dynamics (MHD) simulations [71, 72] in such a way that the Solar observations match simulations. The connection with AQN is to ensure that the energy distribution of nanoflare events in the Solar corona equates with the baryon charge distribution within AQN framework [54], i.e.

[TABLE]

where $\mathrm{d}N$ is the number of nanoflare events per unit time, with energy between $W$ and $W+\mathrm{d}W$, which occur as a result of complete annihilation of the antimatter AQN carrying the baryon charges between $B$ and $B+\mathrm{d}B$. The parameter $\alpha$ is fixed to match observations [71, 72] of the nanoflare distribution.

In what follows we use the same baryon charge distribution $\mathrm{d}N/\mathrm{d}B$ with the purpose of analyzing the axion emission from AQN annihilated in the interior of the Earth. One should also note that it has been recently argued [73] that the algebraic scaling (11) is a generic feature of the AQN formation mechanism based on percolation theory. The phenomenological parameter $\alpha$ is determined by the properties of the domain wall formation during the QCD transition in the early Universe, but it cannot be theoretically computed in strongly coupled QCD. Instead, it will be constrained based on the observations as discussed below.

If we assume that the baryon charge distribution follows a power-law as suggested in (11) the average baryon number of the distribution is

[TABLE]

where $f(B)$ is normalized and the power-law is taken to hold in the range from $B_{\rm min}$ to $B_{\rm max}$. This range is determined by the complicated formation and evolution of the AQN as discussed in [73] but is not strongly constrained from the theoretical side. Model independent constraints on heavy dark matter candidates requires an average mass above $\sim 55$ g (corresponding to $\langle B\rangle\sim 10^{25}$) while lunar seismology disfavours a significant AQN population in the 10 kg to one ton range ($B\sim 10^{28}-10^{30}$) as reviewed in [52, 74]. However these constraints do not constrain the distribution of AQN within this range.

Rather than introducing a completely new set of distributions we will assume that the baryon number distribution follows the nanoflare distribution as implied by expression (11). Following the flare distribution models considered in [71, 72], we investigate three different choices for the power-law index $\alpha$:

[TABLE]

We also investigate two different choices of $B_{\rm min}$: $10^{23}$ and $3\times 10^{24}$ while we fix $B_{\rm max}=10^{28}$. Therefore, we have a total of 6 different models for $f(B)$ corresponding to different models from [71, 72]. In Table 1 we show the mean baryon charge $\langle B\rangle$ for each of the 6 models.

For the simulations in this work, we will only investigate parameters that give $\langle B\rangle\gtrsim 10^{25}$ because present constraints require $\langle B\rangle>3\cdot 10^{24}$ according to (9). This means excluding two models: that with $B_{\rm min}\sim 10^{23}$ and that with $\alpha=2,5$ and $\alpha=2$.

II.4 Axion emission mechanism

The goal of this subsection is to explain the process of axion emission when an AQN enters impacts upon the Earth and annihilates. We refer the readers to Refs. [1, 2] for the technical details on the computation of the axion spectrum.

First, one should recall that the axion plays a key role in the AQN construction: due to the domain wall pressure the AQN remain stable on cosmological time scales. The same domain walls that provide the stability of AQN also contribute to their energy. The analysis of [59] suggests that the axion domain wall contributes $\sim 1/3$ to its total energy. It is this energy in the form of free, propagating axions that can be released when the antibaryon charge from an AQN is annihilated by surrounding material.

Second, note that an axion domain wall in the equilibrium does not emit any axions as a result of a purely kinematical constraint: the domain wall axions are off-shell axions in equilibrium. Now consider an AQN losing mass due to annihilation processes: the axion portion to the energy is initially unchanged, as it is not being directly annihilated. However, time-dependent perturbations due to annihilation processes obviously change this equilibrium configuration, and the configuration becomes unstable with respect to the axion emission because the total energy of the system is no longer at its minimum when some portion of the baryon charge in the core is annihilated. To retrieve the ground state, an AQN will therefore lower its domain wall contribution to the total energy by radiating axions. To summarize: the emission of axions is an inevitable consequence during the annihilation of antimatter AQNs and is a consequence of AQN minimizing their binding energy.

A mathematically consistent procedure that allows one to compute the emitted axion-velocity spectrum was developed in [2] and can be explained as follows: Consider a general form of a domain wall solution as follows:

[TABLE]

where $R_{0}$ is the radius of the AQN, $\phi_{w}$ is the classical solution of the domain wall, while $\chi$ describes excitations due to the time-dependent perturbation. Note that, $\phi_{w}$ is clearly off-shell time-independent classical solution, while $\chi$ describes on-shell propagating axions. Thus, whenever the domain wall is excited, corresponding to $\chi\neq 0$, freely propagating axions will be produced and emitted by the excitation modes.

Suppose an AQN is travelling in vacuum where no annihilation events take place. We expect the solution stays in its ground state $\phi(R_{0})=\phi_{w}(R_{0})$, which corresponds to the minimum energy state. Since there is no excitation (i.e. $\chi=0$), no free axion can be produced. However, the scenario changes when baryon charge from the AQN is annihilated. As the AQN starts to loose a small amount of mass, its size consequentially shrinks from $R_{0}$ to a slightly smaller radius $R_{\rm new}=R_{0}-\Delta R$. The quantum state $\phi(R_{0})=\phi_{w}(R_{0})$ is then no longer the ground state, because a lower energy state $\phi_{w}(R_{\rm new})$ becomes available. Then, we may write the current state of the domain wall as $\phi(R_{0})=\phi_{w}(R_{\rm new})+\phi_{w}^{\prime}(R_{\rm new})\Delta R$, so the domain wall now has a nonzero exciting mode $\chi=\phi_{w}^{\prime}(R_{\rm new})\Delta R$ and free axions can be produced during oscillations of the domain wall.

The annihilation processes inside the AQN force the surrounding domain wall to oscillate. These oscillations of the domain wall generate excitation modes and ultimately lead to radiation of the propagating axions. The corresponding oscillations have a typical frequency of order the axion mass $\sim m_{\mathrm{a}}$; the size of the AQN is also order of $m_{\mathrm{a}}^{-1}$, such that typical velocity of the oscillations of order $c$. These simple dimensional arguments suggest that the typical velocities of the emitted axions will be of order $c$, i.e. the most of the radiated axions will have relativistic velocities, in contrast with conventional galactic axions. However, the present work is devoted to the low energy portion of this spectrum with $v\lesssim v_{\rm esc}$, which is suppressed as $v^{3}$ at $v/c\ll 1$ [2]. Explicit computations presented in [2] support all these claims of the basic features of the axion emission and the corresponding spectrum.

III Annihilation modeling

In this section we present all equations that will be used in our simulations and calculations of the axion-density distribution around the Earth and the speed spectrum of axions passing through the surface of the Earth. In what follows we adopt the following terminology for the sake of brevity: whenever “per AQN” is referred, we mean “per $\langle B\rangle$ baryon charge”.

III.1 Annihilation of AQNs impacting the Earth

The equations derived in this subsection are similar to equations discussed in [55], which was devoted to annihilation of AQNs in the Solar corona. However, there are two main differences between our work and that of [55]: First, in our work, the motion of AQNs in the interior of the Earth will experience differing strengths of local gravity because we find that the average AQN will only be partially annihilated when passing through Earth. This is in contrast to Ref. [55] where AQNs are almost completely annihilated in the transition region 2000 km above the Solar photosphere and there is an extremely small probability for an AQN to survive until the Solar surface. The second distinct feature of our modeling is related to the different structural properties of the Sun and Earth. The Solar corona is a highly-ioninized plasma, where protons, ions and electrons move freely and experience a long-range Coulomb interaction with AQNs due to the induced electric charge on the AQN at $T\neq 0$, see estimates in Appendix A in ref. [54]. This should be contrasted with the structure of Earth, where most atoms are neutral or weakly ionized. Furthermore, most material on Earth are in the solid or liquid phase, with low mobility of the constituents. As a result, atoms from the interior of the Earth interact with AQNs on as a result of head-on collisions with the cross section $\sigma\sim\pi R^{2}$. This key difference ensures that AQN interactions with Earth have a much smaller effective cross section compared to interactions between AQN and the Solar corona.

The energy loss due to the collision of AQNs with the Earth can be expressed as follows [75]

[TABLE]

, where $E=mv^{2}/2$ is the kinetic energy, $s$ is the path length, $\sigma$ is the effective cross section of the AQN, $\rho$ is the density of the local environment and $v$ is the AQN velocity. Note that, following from the definition, we can express the time derivative of $E$ in the form

[TABLE]

where in the last step, we substitute the conventional rate of mass loss

[TABLE]

Comparing Eqs. (15) and (16), we conclude that the mass loss acts like a friction term in the equation of motion

[TABLE]

In the last step, we use $\sigma\equiv\pi(\epsilon R)^{2}$, where $R$ is the radius of the AQN and $\epsilon$ is a dimensionless parameter that describes deviations from the geometrical cross section $\pi R^{2}$. In the studies of the Solar corona preformed in Ref. [55], $\epsilon_{\rm corona}\gg 1$ due to the high mobility of ions in plasma. In our case we expect $\epsilon\sim 1$, but we leave it as a free parameter to account for our ignorance regarding a number of physical effects that we have neglected but that we believe may be important. The main effects that $\epsilon\neq 1$ could describe are partial ionization of the AQN, the fact that AQN travel at supersonic speeds. Both effects can change the internal-core heating of the AQN and therefore change the cross-section. However, in our case, the rigidity of the surrounding material while the AQN traverses the Earth, implies that $\epsilon$ is likely to remain very close to one. We therefore set $\epsilon=1$ in most simulations by default, but we also investigate a much larger value of $\epsilon=3$ in order to assess how our results are modified when the AQN cross-section is significantly increased.

The complete dynamical equations of motion in vector form are:

[TABLE]

where $G$ is the gravitational constant, and we define

[TABLE]

In Eq. (20a), we follow the parameters adopted in [55], in which $\rho_{\mathrm{n}}=3.5\times 10^{17}$ kg m*-3* is the nuclear density and $m=m_{\mathrm{p}}B$ is the AQN mass. In Eq. (20b), we approximate the local environmental density $\rho(r)$ as discrete step functions due to the discontinuous geological structure of the Earth. The labels $i,j=1,...,5$ correspond to layers, summarized in Table 2. The parameter $r_{i}$ is the radius of the start of the layer as measured from the center of the Earth. $\rho_{i}$ is the average density of the corresponding shell respectively. We make the approximation that the density within each layer is uniform, and therefore take the average value from density at the top/bottom of the shell as the local density. The data in Table 2 is taken from Ref.333Don L. Anderson, Theory of the Earth, Boston, Blackwell Publications (1989). Also see, http://pubs.usgs.gov/gip/interior.

III.2 Axion-emission spectrum in the observer frame

We assume axion emission with non-relativistic velocities from AQN will be spherically symmetric in the AQN frame. However, AQNs move with velocities $v_{\rm AQN}\sim 220\,\mathrm{km}\,\mathrm{s}^{-1}\sim 10^{-3}c$, which greatly exceeds the escape velocity of the axions we are interested in. Therefore, the angular distribution of the emission may be strongly affected by the motion of the AQN in the interior of Earth. This subsection is devoted to analysis of this modification in the angular distribution due to the frame change.

In general, an annihilating AQN will emit axions in a frame moving with respect to an observer on Earth. We introduce the notation $\tilde{K}$ and $K$ for the rest frame of the AQN and the frame of the observer respectively. The axion-emission velocity spectrum in the AQN rest frame is calculated in Ref. [2] and we refer the reader to that paper for the calculation details. The velocity spectrum in the rest frame is, by definition, the derivative of the radiated flux $\Phi_{\rm rad}(\tilde{v}_{\mathrm{a}})$:

[TABLE]

where $\tilde{v}_{\mathrm{a}}$ and $\tilde{E}_{\mathrm{a}}$ are the rest frame axion velocity and energy respectively. The function $H_{l}(\tilde{p},\delta)$ corresponds to partial wave expansion in the approximate solutions, as derived in [2]. In the present work we are interested in non-relativistic AQNs and axions, therefore, only the $l=0$ mode and $\tilde{v}_{a}/c\ll 1$ are considered in Eq. (21). The parameter $\delta\in(0,1)$ is a convenient factor introduced in Ref. [2] as a result of approximations due to absence of simple analytic solutions of the general expressions. Tuning $\delta\in(0,1)$ leads to changes in the velocity spectrum (21) that do not exceed $\sim 20\%$. The velocity spectrum is not known more accurately than this (see Appendix B). The normalization factor $N(\delta)$, depending on the parameter $\delta$, is also known and presented in Appendix B.

In the frame of the observer, an AQN is moving with a velocity $v_{\rm AQN}\lesssim 10^{-3}c$. Thus, we need only to consider a non-relativistic transformation of frames, with relations as follows:

[TABLE]

Working in the manifold of $\mathbf{v}_{\mathrm{a}}$, we know $\rho(v_{\mathrm{a}})$ is normalized within a unit 3-ball $B^{3}$

[TABLE]

where in the second step, we transform our coordinate from $\tilde{v}_{\mathrm{a}}$ to $v_{\mathrm{a}}$ using Eqs. (22). Note that such transformation produces a small error related to a slight shift of spherical center by $\sim 10^{-3}c$. However, the inconsistency is negligible comparing to the uncertainty of approximation in terms of $\delta$. Noting that the last equality in Eq. (23) is completely expressed in terms of $v_{\mathrm{a}}$, the velocity in the frame of the observer. Therefore, we read off the emission spectrum in the frame of the observer

[TABLE]

where $\Omega$ is the solid angle made by $\langle\hat{\mathbf{v}}_{\mathrm{a}},\hat{\mathbf{v}}_{\rm AQN}\rangle$. For a moving AQN with axion emission, it is reasonable to assume azimuthal symmetry, which gives

[TABLE]

where $\tilde{v}_{\mathrm{a}}$, $\tilde{p}$, and $\tilde{E}_{\mathrm{a}}$ are expressed as functions of $v_{\mathrm{a}}$ defined in Eqs. (22). In this work, we are specifically interested in the hierarchy when $v_{\mathrm{a}}\ll v_{\rm AQN}\ll c$, such that one can expand Eq. (25) to arrive at

[TABLE]

and using the fact that $H_{0}(p,\delta)$ is locally quadratic near $p=0$, (Appendix B). These corrections originate from the fact that axions that are eventually trapped are emitted with velocity in the rest frame that greatly exceeds the escape velocity, but they are emitted in the opposite direction to $\mathbf{v}_{\rm AQN}$, so cancellation occurs and these are the axions that end up with a low velocity in the observer frame and therefore end up trapped and contributing to $\rho_{\rm obs}(v_{\mathrm{a}})$. This transformation has a small dipole term, but we have numerically checked that the dipole contribution to the emission is small ($\pm 6\%$ correction for the dipole, while the quadrupole gives $\pm 0.2\%$), and so we neglect it and take axion emission to be isotropic in the observer frame.

To summarize: the emission spectrum the frame of the observer (26) is directly expressed in terms of the spherically symmetric spectrum in the rest frame via equation (21) and computed in [2]. This is expressed in terms $H_{0}(0,\delta)$ given in Appendix B. We emphasize that the axion angular distribution is not spherically symmetric; instead, for each incoming AQN, it is azimuthally symmetric with respect to the direction defined by the velocity of the AQN $\mathbf{v}_{\rm AQN}$. It is the integration over incoming directions in Eq. (24) that leads to the simple relation in equation (26) between the two distributions.

III.3 The axion density distribution

The slowest of the axions that are emitted due to AQN annihilation while travelling through the Earth will be trapped by the gravitational field. As the axions are non-interacting particles, their dynamical equations of motion are simpler than that of the AQNs:

[TABLE]

In other words, the trajectories of axions are entirely determined by the gravitational force from the Earth. Axions are mutually non-interacting, so if they are generated trapped then they remain trapped for all time. The only difference in comparison with standard Kepler problem being that the effective acceleration, $g(r)$, does not scale like $r^{-2}$ in the interior of the Earth. We will use this set of equations and sample large numbers of axion trajectories with initial conditions set by the AQN annihilation events within Earth. This allows us to obtain a distribution function $p_{\mathrm{a}}(r)$ that gives us the probability to find an axion at a given distance $r$ from the center of the Earth. This is then directly related to the axion density profile, which builds up around the Earth over 4.5 Gyrs and is the crucial quantity for the axion search experiments reviewed in the introduction. A full description of the numerical simulations that allow us to compute $p_{\mathrm{a}}(r)$ will be described in Sec. IV, while here we relate $p_{\mathrm{a}}(r)$ to other observables that enter the numerical simulations.

First, we want to relate $p_{\mathrm{a}}(r)$ to the energy density of gravitationally-trapped axions, $\rho_{\mathrm{a}}(r)$:

[TABLE]

where $\langle N_{\rm AQN}^{\rm 4.5Gyr}\rangle$ is the expected total number of antimatter AQNs that impacted the Earth within 4.5 Gyr, and $\langle N_{\mathrm{a}}^{\rm trap}\rangle$ is the mean number of trapped axions emitted per AQN. The first term entering (28) is easy to estimate

[TABLE]

where $\langle\dot{N}\rangle$ is the average flux of AQNs hitting the surface of the Earth, estimated in Appendix C.

The estimation of $\langle N_{\mathrm{a}}^{\rm trap}\rangle$ in (28) is more complicated. We choose to express this in terms of the “heat” distribution function $q(r,v_{\rm AQN})$ defined as the fraction of the mass loss of an AQN moving with local velocity $v_{\rm AQN}$ at a location $r$. We name this distribution “heat” to emphasize that this function describes the annihilation events where locally one should expect about $2m_{\mathrm{p}}c^{2}$ energy release per single event of annihilation. We normalize the heat distribution as follows:

[TABLE]

where $\langle\Delta m_{\rm AQN}\rangle$ is mean total mass loss per AQN as it passes through the Earth. The total number of axions emitted for a given mass loss $\Delta m_{\rm AQN}$ can be estimated as

[TABLE]

where $\langle E_{\mathrm{a}}\rangle\simeq 1.3m_{\mathrm{a}}$ is the average energy of axions computed from the spectrum in equation (21). The coefficient $1/3$ in (31) comes from the fact that approximately $1/3$ of the total AQN energy is stored in the form of axion energy at formation time. This energy is inevitably released as axions during the annihilation.

We are interested in trapped axions with velocities below the local escape velocity $v_{\rm esc}(r)$ in the interior of the Earth, with $v_{\rm esc}(r)$ given by:

[TABLE]

Now we can express $\langle N_{\mathrm{a}}^{\rm trap}\rangle$ entering (28) in terms of the heat distribution $q(r,v_{\rm AQN})$ as follows

[TABLE]

where $\rho_{\rm obs}(v_{\mathrm{a}})$ is the axion-velocity spectrum in the frame of the observer, computed in equation (26). The computation of $q(r,v_{\rm AQN})$ requires Monte Carlo simulations, which we now discuss.

IV Algorithm and simulations

IV.1 Simulating initial conditions of AQNs for the heat emission profile

To simulate the trajectories of AQNs through Earth we start from the velocity distribution function of AQNs in the vicinity of Earth. This distribution is Gaussian in each Cartesian direction with dispersion $\sigma\simeq 110\,\mathrm{km}\,\mathrm{s}^{-1}$ and an additional ‘wind’ component in one direction $\mu\simeq 220\,\mathrm{km}\,\mathrm{s}^{-1}$ due to the motion of the Solar system with respect to the DM halo:

[TABLE]

We consider the accumulation of axions trapped in orbit around Earth over the lifetime of the Solar System ($\sim 4.5\,\mathrm{Gyr}$). Over this time, due to the spin and orbit of the Earth, the orbit of the Solar system around the galactic centre, and perturbations to both of these due to secular effects of gravity, we consider the wind direction to have averaged over all $4\pi$ solid angle, see Fig. 1. Taking the average of Eq. (34) gives us an isotropic distribution for incoming AQN speed that includes the effect of the wind:

[TABLE]

To start our simulations we must use the axion flux distribution through the surface of the Earth. The flux distribution is given by

[TABLE]

We show the derivation of the distributions (35) and (36) in Appendix C. As a consistency check of the isotropized distribution in equation (36) we also carry out calculations in which the wind direction is fixed in space (Appendix D) and obtain consistent results.

To generate AQN initial conditions, we uniform-randomly choose a point on the surface of Earth and then generate an AQN with velocity sampled from distribution (36). We remove AQNs with initial velocity pointing away from the Earth as these correspond to AQNs that would have already experienced an Earth interaction. Following this sampling scheme, we generate $10^{6}$ AQNs and compute $q(r,v_{\rm AQN})$, the heat emission profile of axions at a given location within the Earth, $r$, and originating from an AQN with local velocity $v_{\rm AQN}$ as defined by Eq. (30). As explained in Sec. III.3, this function is a key element to the calculation of the axion density distribution $\rho_{\mathrm{a}}(r)$. Furthermore in the next subsection, we will see that it is $q(r,v_{\rm AQN})$ that solely determines the number of axions we simulate at a given point in the phase space $(r,v_{\rm AQN})$.

For individual AQNs, Eqs. (19) are solved numerically, and a 2D histogram of the $(r,v_{\rm AQN})$ phase space is computed, where each point in a given trajectory is weighted by $\dot{m}(t)$ to give the total amount of annihilation energy going into axions. We refer to Appendix E for further details .

IV.2 Computations of the density distribution of trapped axions

The density distribution of trapped axions $\rho_{\mathrm{a}}(r)$ is directly related to $p_{\mathrm{a}}(r)$, the probability of finding a trapped axion at a given distance $r$, given by Eq. (28). Similar to the calculation of the AQN heat profile $q(r,v_{\rm AQN})$, we obtain $p_{\mathrm{a}}(r)$ by simulation, but the numerical procedure is more sophisticated. First, it is clear that the number of axions to simulate at a given point in the phase space $(r,v_{\rm AQN})$ is proportional to two weighting factors – the rate of axion production [i.e. the heat profile function $q(r,v_{\rm AQN})$], and the rate of producing trapped axions [i.e. the spectrum $\rho_{\rm obs}(r)$ and the local escape speed $v_{\rm esc}$]. Therefore, we sample axion initial conditions from the distribution within the integral

[TABLE]

The interpretation of Eq. (37) is as follows: we first choose the initial conditions $(r,v_{\rm AQN})$ with probability weighted by $q(r,v_{\rm AQN})v_{\rm AQN}/\langle v_{\rm AQN}\rangle$, and then we shoot an axion at the given position $r$ at a random angle,with a speed in the range $[0,v_{\rm esc}(r)]$ weighted by $\langle v_{\rm AQN}\rangle/v_{\mathrm{a}}\rho_{\rm rest}(v_{\mathrm{a}})$.

Monte Carlo sampling can be performed in 2D space due to the azimuthal symmetry. At $\mathbf{r}(t=0)=(0,r)$, the velocity

[TABLE]

is randomly generated using the following distributions:

[TABLE]

Then, we obtain the initial velocity of the sampled axion:

[TABLE]

Note that the sign of $v_{\theta}$ is irrelevant, as $p_{\mathrm{a}}(r)$ preserves symmetry of time reversal.

Unlike for the AQN, the equations of motion for the axion are simple and obey an energy and angular momentum conservation laws, see Eq. (27). We use these conservation laws, together with the fact that the probability, $p_{\mathrm{a}}^{\rm(one)}(r)$, to find a single trapped axion at a distance, $r$ from Earth must be inversely proportional to the radial velocity $v_{r}$ of that axion at that distance:

[TABLE]

where $N$ is a normalization, and $v_{r}$ can be explicitly written in terms of $r$ using the conservation laws for angular momentum and energy:

[TABLE]

We then conclude that

[TABLE]

The full $p_{\mathrm{a}}(r)$ is obtained by averaging over sufficient numbers of $p_{\mathrm{a}}^{\rm(one)}(r)$ using the Monte Carlo technique.

IV.3 Computations of the axion spectrum on the Earth surface

Finally, we describe the method to calculate $F(v_{\mathrm{a}})$, the axion speed spectrum on the Earth surface. The spectrum is important for detection instrument as the de Broglie wavelength of the trapped axion $\lambda_{\mathrm{a}}=h/m_{\mathrm{a}}v_{\mathrm{a}}$ is sensitive to the local speed in the non-relativistic limit.

We follow a similar approach as in the previous subsection. Once the initial conditions are known, it is easy to find the speed of an axion at the surface of the Earth from the equation of motion (42):

[TABLE]

We can then obtain the axion speed profile, $F(v_{\mathrm{a}})$, by computing a histogram of the axions with weighting inversely proportional to the orbital period of each axion (axions with shorter periods pass through the surface more often). The period of a given trajectory is found from (42):

[TABLE]

Therefore, the (half-)period is obtained by numerical integration

[TABLE]

V Axion profiles: density distribution and velocity spectrum

V.1 Heat emission profile

We first discuss the results of the simulation ($10^{6}$ samples) for the heat-emission profile of axions: $q(r,v_{\rm AQN})$. In our simulations we use a number of parameters such as the AQN baryon-charge distribution (parameters $\alpha$ and $B_{\rm min}$), $\epsilon$ (describing the strength of the interaction of AQNs with earth’s material) and the flux distribution of incoming AQNs (characterized by the galactic wind and DM velocity dispersion). Surprisingly, we find that our results are quite insensitive to these parameters. Therefore, in the main body of the paper we only present the case with $\alpha=(1.2,2.5)$ and $B_{\rm min}=3\times 10^{24}$ shown in Fig. 2, while leaving other cases with different parameters to Appendix F.

From Fig. 2, we first observe a Maxwellian-like distribution with mean speed of approximately 300 $\,\mathrm{km}\,\mathrm{s}^{-1}$ and dispersion 200 $\,\mathrm{km}\,\mathrm{s}^{-1}$ in all cases. In addition, we specifically note that $q(r,v_{\rm AQN})$ looks like a stack of multiple Maxwellian-like distributions, such that an abrupt change appears at some specific distance (e.g. $r\sim 0.9$ and 0.5$R_{\oplus}$). One can see, from Table 2, that these jumps correspond to the successive layers of the Earth interior. Consequently, whenever an AQN moves into a new layer with an abrupt change of local density, the mass loss drastically changes. Another observations is that very few AQNs reach the core of the Earth, as the dominant portion of the AQNs propagate at the distances $r\geq 0.5R_{\oplus}$, i.e. far away from the center of the Earth. This is a geometrical effect that stems from the fact that there is more volume in the outer part of the Earth compared to the inner part.

As shown in Appendix F, we reach very similar conclusion for almost all effects we considered. The heat emission profile is insensitive to the baryon-charge distribution of AQN models and to the velocity distribution of incoming DM flux. The effective cross section, parametrized by $\epsilon$, obviously increases the number of annihilation events expressed as $\Delta B/B$ in Table 3. However, the corresponding increase of $\epsilon$ does not lead to significant modification of the axion density on the surface $\rho_{\mathrm{a}}(R_{\oplus})$, which is the main observable of the system. We have verified that this conclusion of insensitivity to variety of physical parameters is a generic feature of the system and this feature will be further elaborated in next subsections.

V.2 Axion density distribution

The key results of our simulations are summarized in Table 3. From the mass loss, we obtain the expected number of trapped axions $\langle N_{\mathrm{a}}^{\rm trap}\rangle$ per AQN by Eq. (33) and substitute this value into Eq. (28), the energy density distribution of axions is therefore

[TABLE]

where $p_{\mathrm{a}}(r)$ is computed using the Monte Carlo simulations discussed in the previous section.

We can parametrize our results for $\rho_{\mathrm{a}}(r)$ using a simple analytical fit as follows

[TABLE]

where $x\equiv r/R_{\oplus}$, with parameters chosen to be $(a_{1},n_{1},k_{1},m_{1})=(3.11919,1.58484,1.11709,2.06927)$ and $(a_{2},n_{2})=(1.000051,1.56429)$ for the best fit. Note that there exists a discontinuous jump at $r=R_{\oplus}$, which is due to the steep change of density near the Earth surface, because axion emission always happens inside the earth. The corresponding plot is presented in Fig. 3. Similar plots obtained with different choice of parameters are shown in Appendix F. It is found that the approximation (48) fits the density profile of all models with $\epsilon\simeq 1$ over a range of radii $r\in[0.5,5]R_{\oplus}$ with a maximum deviation of 10%.

We would like to make following comments:

$\bullet$ We note that the shape of $\rho_{\mathrm{a}}(r)$ is not sensitive to the parameters ($\alpha$, $B_{\rm min}$, and flux distribution) of the AQN mass-distribution function; only its magnitude shows a modest variation. This feature was anticipated before the computations: Indeed, the axions inside the Earth are non-interacting particles (except for the gravity). However, the gravitation becomes weak near the inner core in comparison with the surface of the Earth. Thus, one can treat axions as a free ideal gas with a probability density scaling as $\sim r^{2}$ near $r\lesssim R_{\oplus}$, and as $\sim r^{-2}$ outside the radiation source. We can do this because the axion-density profile is a statistical count of occupation states when the entropy is maximized, and therefore it should not depend on the location of the heat source inside the system (Earth). This argument is generic, model independent and is supported by explicit computations that we have parameterized by simple analytical formula (48).

$\bullet$ The amplitude of the axion-density profile increases with $\epsilon$. However, the observable density $\rho_{\mathrm{a}}(R_{\oplus})$ on the earth’s surface shows very modest changes when $\epsilon$ varies as discussed in Appendix F.

Qualitative explanation of this effect is as follows. In general, a larger $\epsilon$ results in a peak in the density profile closer to the surface of the Earth because more AQNs annihilate closer to the surface of the Earth due to a larger effective cross section. Nevertheless, the observable $\rho_{\mathrm{a}}(R_{\oplus})$ shows only modest variations with respect to these changes as mentioned above. This is not a surprise because $\rho_{\mathrm{a}}(r)$ is predominantly determined by the local annihilation rate of AQNs, and therefore the local density of the Earth. Consequentially, even for the large $\epsilon$ cases, $\rho_{\mathrm{a}}$ becomes more uniform in the inner layers, but drops quickly near the surface due to the reduced rate of local AQN annihilation.

$\bullet$ It is instructive to compare this result with the crude estimate (4) presented in Ref. [2]: there is about 3 orders of magnitude deviation from that naive estimate. This can easily be understood from the work described here: First, [2] assumed that $\Delta B/B\sim 1$ (i.e., most incident AQN are mostly annihilated), however, our simulations show that $\Delta B/B\sim 0.1$ for a typical model with $\epsilon\simeq 1$ (see Table 3). Second, the emitted axions form an extended cloud in space around the Earth, such that the axion density remains high even far away from Earth, at the distances $r\sim 5R_{\oplus}$, for instance see Fig. 3 and corresponding figures in Appendix F. It is the physical size of the extended cloud that is the main source for the discrepancy with the naive estimate in [2], where it was assumed the axion density dropped sharply beyond the surface of the Earth. Finally, there is a suppression factor related to averaging over inclination angles between the AQN velocity and the surface of the Earth that was missing from [2]. The order of magnitude estimate from [2] assumes $\dot{N}\simeq nAv_{\rm AQN}$. The proper treatment of the inclination angle is given in Appendix C and produces an additional suppression factor $\sim\frac{1}{4}$ such that corresponding expression assumes the form $\dot{N}\simeq\frac{1}{4}nAv_{\rm AQN}$.

The results of our Monte Carlo simulations are presented in Table 3. The main result from the point-of-view of axion-search experiments is the axion density on the surface of the Earth $\rho_{\mathrm{a}}(R_{\oplus})$444We do not expect any large enhancement factors due to the gravitational focusing of DM, which may briefly increase the incident AQN flux by factor $10^{6}$, as advocated in [76, 77]. This is because trapped axions are accumulated during 4.5 billion of years and enhancement effects due to gravitational focusing are of short temporal duration and we do not think they will significantly change the mean AQN impact rate over the long accumulation time.

[TABLE]

There is no remaining freedom in the model to significantly modify our final results in eq. (49). In this scenario, the AQN framework is rigid and predictive. This is because dark- and visible-matter densities in the AQN framework must always satisfy the relation (6), irrespective of the axion mass $m_{\mathrm{a}}$. This is in contrast with conventional estimates for the galactic-axion density, which strongly depends on the axion mass and scale as $m_{\mathrm{a}}^{-7/6}$.

The axion velocity distribution has a maximum close to $v_{\mathrm{a}}(R_{\oplus})\simeq 8\,\mathrm{km}\,\mathrm{s}^{-1}$ (see next subsection V.3). The corresponding wavelength $\lambda_{\mathrm{a}}\sim{\hbar}/({m_{\mathrm{a}}v_{\mathrm{a}}})$ is approximately $30$ times greater than for galactic axions, which have a typical velocity of about $\sim 220\,\mathrm{km}\,\mathrm{s}^{-1}$. Therefore, coherent effects can be maintained for a longer time period compared to those for conventional galactic axion searches. One may hope that the feature of having a large coherence length, $\lambda_{\mathrm{a}}\sim v_{\mathrm{a}}^{-1}$, could play a key role in the design of instruments, capable of discovering such gravitationally trapped axions (see concluding Section VI for more comments).

Indeed, the local axion density (49) for the trapped axions is approximately three orders of magnitude smaller than conventional galactic axion density (assuming that the galactic axions with a given mass $m_{a}$ saturate the DM density). If the experimental sensitivity on the axion-photon coupling is inversely proportional to the square root of the local axion density555this case is realized when the design of the axion search experiment is based on coherent accumulation of the signal., one should expect a suppression of order $\sqrt{10^{-3}}\sim 1/30$ or so in the amplitude. However, the enhancement factor due to the longer wave length $(220\,\mathrm{km}\,\mathrm{s}^{-1})/(8\,\mathrm{km}\,\mathrm{s}^{-1})\sim 30$ may cancel the suppression factor in the amplitude mentioned above. It remains to be seen if this potentially strong enhancement factor can be practically realized in real experiments.

V.3 Axion-velocity spectrum features on the surface of the Earth

In this subsection we derive the axion-velocity spectrum $F(v_{\mathrm{a}})$ at the surface of the Earth. Before discussing the simulation results, we note that the axion spectrum $F(v_{\mathrm{a}})$ can be estimated based on a simple toy model. Surprisingly, the estimated profile gives a consistent and precise description which is in good agreement with Monte Carlo numerical simulations. The arguments presented below strongly suggest that $F(v_{\mathrm{a}})$ is not sensitive to any parameters of the incoming AQN flux, and $F(v_{\mathrm{a}})$ is expected to be similar regardless of the distribution functions computed in previous sections. The main reason for this insensitivity to parameters is that the basic physics of low-energy axions is governed by gravity rather than by any specific features related to the AQNs emitting these axions.

We start our argument by rewriting eq. (44) as

[TABLE]

and introduce following notation

[TABLE]

The escape velocity at the surface of the Earth is given by

[TABLE]

The energy in eq. (50) can now be written as

[TABLE]

We now construct a toy model using a number of approximations. First, we assume that the Earth is point mass, which simplifies the gravitation potential:

[TABLE]

This is correct for all distances above the surface of the Earth. As the next step, we use Eqs. (53) and (54) to simplify Eq. (50) and derive

[TABLE]

which should be valid for $r\geq R_{\oplus}$. From this simplified case, one can extract many useful relations that must be valid for all radii above the surface of the Earth. For example, by setting $v_{\mathrm{a}}(r)=0$ in Eq. (55), one can estimate the maximum radial distance $r_{\rm max}$ for a trapped axion as follows

[TABLE]

This relation allows to estimate the orbital period for these axions:

[TABLE]

where $v_{\mathrm{a}}\equiv v_{\mathrm{a}}(R_{\oplus})$ is the speed of a given axion at the surface of the Earth. Also, at $r=R_{\oplus}$, the parameter $u_{0}$, the axion velocity at the center of the Earth, can be expressed as a function of $v_{\mathrm{a}}$ using Eq. (55) as follows:

[TABLE]

We can now estimate the axion spectrum $F^{\rm(est)}(v_{\mathrm{a}})$ at the surface of the Earth:

[TABLE]

where $u_{0}$ is given by Eq. (58). Expression (59) comes from the fact that $F(v_{\mathrm{a}})$ must be proportional to the axion emission spectrum $\rho_{\rm obs}$ and the frequency of the orbital period $1/T$. The last term $(\mathrm{d}u_{0}/\mathrm{d}v_{\mathrm{a}})$ is the Jacobian of the transformation from the basis of $u_{0}$ to $v_{a}$.

Substituting Eq. (58) into Eq. (59), the spectral axion distribution (not normalized) at the surface of the Earth can be written as:

[TABLE]

where $v_{\rm esc}(r_{0})$ is the only unknown parameter affecting the form of $F^{(\rm est)}(v_{\mathrm{a}})$. However, the local escape velocity $v_{\rm esc}(r_{0})$ varies only in a narrow range, from $11\,\mathrm{km}\,\mathrm{s}^{-1}$ (from the surface of the Earth) to $15\,\mathrm{km}\,\mathrm{s}^{-1}$ (from the center of the Earth). Therefore, the spectrum $F^{(\rm est)}(v_{\mathrm{a}})$ given by (60) depends only weakly on $v_{\rm esc}(r_{0})$.

Fig. 4 shows that our simple estimate, Eq. (60), and the simulations are in very good agreement. In Appendix F we show that the agreement holds for different choices of the AQN model parameters. The reason for this insensitivity to large variations in AQN parameters (such as $\alpha$, $B_{\rm min}$, $\epsilon$, incident-flux distribution) is because the velocity spectrum of the gravitationally-bound axions is mostly determined by the gravity of the Earth, rather than by any features of the AQNs. Therefore, the analytical expression (60) remains a valid description of the spectrum, where all specific features related to the AQNs are hidden in $v_{\rm esc}(r_{0})$ which can be assumed to be a constant computed as an average over the entire ensemble of axions emitted from any point in the interior of the Earth by any AQN of any size at any moment. One can therefore use $\langle v_{\rm esc}(r_{0})\rangle\simeq 13.5\,\mathrm{km}\,\mathrm{s}^{-1}$ in Eq. (60), which is in good agreement with numerical simulations.

VI Conclusions and future directions

The main results of our work can be summarized as:

a. We computed the energy-density profile of trapped axions as a function of the distance from the center of the Earth using numerical Monte Carlo simulations. These axions are produced when antimatter AQNs incident on Earth annihilate with the rocky material in the interior. The trapped axions represent the tiny, low-velocity tail of the axions emitted during this annihilation process, a velocity distribution that otherwise peaks at around $0.5c$. However, the number of trapped axions can be appreciable since they build up over the lifetime of the Earth. Results for the density profile are shown in Fig. 3, and in Figs. 8(a), 8(b) in Appendix F for different parameters of the incident AQN population ($\alpha$ and $B_{\rm min}$, governing the AQN-mass distribution; $\epsilon$ governing the AQN-Earth material cross section; incident flux velocity distribution). The final value we calculate for the axion density at the surface of the Earth after 4.5 Gys, $\rho(R_{\oplus})$, is not very sensitive to these AQN parameters, which suggests that one can use the estimate (49) as a baseline for the trapped-axion density that is quite independent of the details of the AQN model. We refer the reader to Table 3 for an overview of the sensitivity of the results to the AQN parameters.

b. We computed the trapped-axion velocity-distribution $F(v_{\mathrm{a}})$, shown on Fig. 4, and Figs. 8(e), 8(f) in Appendix F for different parameters of the incident AQN population. Once again, we find that $F(v_{\mathrm{a}})$ is not sensitive to the incident AQN parameters, implying that the formula (60) can be used as a good approximation for the velocity distribution of the gravitationally bound axions that is quite independent of the details of the AQN model. One can express the important features of the trapped-axion system in terms of the de Broglie wavelength $\lambda_{\mathrm{a}}$:

[TABLE]

which is much greater than the typical wavelength of the conventional galactic axions because $v_{\mathrm{a}}\leq v_{\oplus}\approx 11\,\mathrm{km}\,\mathrm{s}^{-1}$ in formula (61), whereas typical galactic axions have $v_{\mathrm{a}}\simeq 200\,\mathrm{km}\,\mathrm{s}^{-1}$. This key feature of gravitationally-trapped axions may be a decisive element that may drastically increase the discovery potential for these axions. This is because axion-search instruments that are based on idea of coherence could be made more sensitive, and could collect the signal for much longer period of time without losing the coherent features. These experiments may be able to overcome the relatively small axion density we calculate in eq. (49) in comparison with conventional estimates for the galactic axions, due to the possibility that the signal may be collected over a longer period of time.

For example, the conventional cavity quality factor $Q\approx 10^{5}-10^{6}$ is mostly constrained by the finite bandwidth ${\cal{O}}(\beta^{2})\sim 10^{-6}$ for galactic axions with $\beta\equiv v/c\sim 10^{-3}$. Galactic axions can be thought as a highly monochromatic radiation with a coherence time corresponding to $10^{6}$ oscillation periods and coherence length of roughly $10^{3}$ times the axion Compton length $\lambda_{C}={\hbar}/({m_{\mathrm{a}}c})$, see e.g. the analysis [49] of the DM radio experiment and recent review [19]. This should be compared with the gravitationally trapped axions when $\beta_{\oplus}\equiv v_{\oplus}/c\simeq 3.6\times 10^{-5}$, and therefore factor $Q\sim 10^{10}$ could be potentially four orders of magnitude greater than for the galactic axions. This huge enhancement factor is due to the fact that gravitationally trapped axions are much more monochromatic than galactic axions as the coherence length is roughly $\lambda_{\mathrm{a}}/\lambda_{C}\sim 10^{5}$ of the axion Compton length, to be contrasted with factor $10^{3}$ for the galactic axions.

We mentioned in Sec.I a number of instruments existing, upgraded (or planning to be upgraded), designed (but not build yet), or under consideration that in principle are capable of achieving the precision required to detect gravitationally-trapped axions. Along with the list of instruments presented in Introduction we also want to mention another possible design, advocated in [78], which is based on quantum coherence of the propagating axions, and which may also benefit from the large coherence length (61) mentioned above. Finally, we want to mention the idea advocated in [79, 80, 81] that the combination of large number of cavities can drastically enhance the axion signal. The design advocated in [79, 80, 81] may benefit from large wave length (61) of the gravitationally bound axions.

It is the central result of this work that all the questions formulated in the Introduction as items 1-6 have been successfully answered: we have derived the axion density on the surface of the Earth, where detection can be made (as stated in item “a” above) and the velocity-spectral features of these gravitationally trapped axions (as stated in item “b” above). We have also demonstrated that these results are not very sensitive to the incident AQN distribution model, nor to the interaction pattern between the AQNs in the interior of the Earth. Therefore, we think our results are solid and robust predictions within the AQN framework. We emphasize once again that the corresponding properties (49), (60) and (61) are very distinct from conventional galactic axions. In particular, the spectral features (60) and (61) describe the gravitationally trapped axions that are produced in the interior of the Earth; and no other model we know of would generate such spectral properties. The rigidity and predictiveness of the model leads to very limited freedom and flexibility for any modification of (49), (60) and (61).

Why should we take this AQN model seriously? First of all, 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 $\Omega_{\rm DM}\sim\Omega_{\rm visible}$, it may also explain a number of other (naively unrelated) phenomena, such as the excess of galactic emission in different frequency bands as discussed in section II.2. The AQN model may also offer resolutions to some other astrophysical mysteries: the so-called “Primordial Lithium Puzzle” [53], the so-called “The Solar Corona Mystery” [54, 55], the recent EDGES observations of a stronger than anticipated 21 cm absorption features [56], to name just a few.777In fact, there is some conceptual similarity between the “Primordial Lithium Puzzle” and the axions-emission problem which is the topic of the present work. In both cases, a fraction of the constituents is hidden inside the AQNs. The Li nuclei are captured by AQNs after BBN at $T\simeq 20$ keV as discussed in [53], while off-shell axions are hidden in the form of axion domain walls during the formation time at $T\simeq 100$ MeV as discussed in Section I. In first case it leads to the depletion of the visible Li- nuclei as argued in [53]. In the second case the originally hidden axions might become observable propagating on-shell particles as a result of AQN annihilation events similar to those discussed in the present work. These cosmological puzzles could be resolved within AQN framework with the same set of physical parameters advocated in the present work. The only parameter that changes is $\epsilon$ which scales the effective cross section as in expression (18).

In this respect, the proposal advocated in this work: to search for gravitationally trapped axions that are inevitably produced as a result of the annihilation of incident antimatter AQN with the matter of the Earth is a direct manifestation of the AQN model. In fact, the observation of these axions, with their very distinct density and velocity-spectral properties would be the smoking gun supporting the entire AQN framework.

Acknowledgments

AZ is thankful to Konstantin Zioutas for questions which motivated these studies. AZ also thanks Kent Irwin and Yannis Semertzidis for discussions on possible enhancement factors due to the large wave length (61) which drastically increases the coherence time scale. We are also thankful to David Marsh, Alex Millar, Jan Schuette-Engel, Pierre Sikivie and Mike Tobar for comments, correspondence and discussions. This work was supported in part by the National Science and Engineering Research Council of Canada. AM acknowledges support from the Horizon 2020 research and innovation programme of the European Union under the Marie Skłodowska-Curie grant agreement No. 702971.

Appendix A Constraint from IceCube

Here we estimate the constraint of the AQN model from IceCube’s observation. In short, depending on the sensitivity of the detector, we find a mean baryon charge $\langle B\rangle\lesssim 1.6\times 10^{24}$ is excluded up to $3.5\sigma$ level of significance for $\eta=1$. For $\eta=3$ similar constraint reads $\langle B\rangle\lesssim 4.7\times 10^{24}$, see definition for the parameter $\eta$ below. We present the details of such estimation as follows.

To estimate the hit rates of AQNs to the IceCube, we assume majority of AQNs pass through the Earth in the end without annihilation, which is in agreement with the simulation results in the present work.888To be more specific, simulation indicates the trajectories of AQNs are mostly as like free motions and annihilations only slow down the AQNs by a small fraction $\sim 25$-50 $\,\mathrm{km}\,\mathrm{s}^{-1}$. Since slower velocity results in smaller flux, the simplification of zero velocity loss (i.e. free motion) is in fact putting us to a stricter constraint from IceCube. Then, we conclude the flux of AQNs hitting the Icecube must be approximately isotropic in all directions. To simplify the calculation without loss of generality, we further approximate IceCube as a spherical detector with diameter $\sqrt{\eta}\cdot(1\rm km)$. As shown in Fig. 5, we introduce a parameter $\eta$ to account for the effective cross section of the detector. For the modest estimation, we may expect $\eta=1$, while $\eta=3$ is clearly the maximal cross section ever possible. In all cases, the reasonable range of $\eta$ is from 1 to 3, where $\eta=3$ is a strict upper limit of IceCube’s sensitivity. Using average flux density of AQNs computed in Appendix C, we obtain the total hits expected to IceCube in 10 years from Eq. (74):

[TABLE]

Assuming a Poisson distribution, the probability of observing $k$ events within 10 years is

[TABLE]

The probability of seeing zero event in IceCube over 10 years is summarized in Table 4. As an rough estimate, we choose the average value between $\eta=1$ and 3

[TABLE]

as the constraint from IceCube. This is precisely the value presented in the main body of the paper in eq. (9).

Appendix B Spectral properties in the rest frame

In this appendix, we fulfill the technical details in Sec. III.2. In Eq. (21), the $H_{l}(\tilde{p},\delta)$ is the partial wave expansion known as follows [2]:

[TABLE]

where $f(a,b,c;z)$ is defined to be the regularized Gauss hypergeometric function ${}_{2}F_{1}(a,b,c;z)$, i.e. $f(a,b,c;z)\equiv\frac{1}{\Gamma(c)}{}_{2}F_{1}(a,b,c;z)$, see Refs. [82, 83] and recent article [36]. We also note a useful fact in Eq. (65b) that $f(a,b,c;z)$ has a simple quadratic behaviour in the non-relativistic limit $z\rightarrow 0$.

Next, we turn to the normalization factor $N(\delta)$. According to Eq. (65), the normalization factor obviously depends on $\delta$. We refer the detailed calculation to the original work [2], while quoting results as follows

[TABLE]

In this work, we choose the intermediate value $\delta=0.5$ in the simulation.

Appendix C Estimation of AQN flux

We compute the average flux of AQN (in number density). We start with the standard definition of flux

[TABLE]

where $n$ is the number density of AQNs, $A=4\pi b^{2}$ is the cross section of the impact, and $v$ is the velocity of AQN flux. Now we take the average:

[TABLE]

where we take the average number density of an antimatter AQN

[TABLE]

and the incident flux vector is also averaged out as

[TABLE]

where in the last step, we take $\sigma=110\,\mathrm{km}\,\mathrm{s}^{-1}$, $M_{\oplus}=5.972\times 10^{24}\rm kg$, and $R_{\oplus}=6371\rm km$, and the impact parameter $b$ for the Earth is taken into account:

[TABLE]

As an additional note, in Eq. (70) we have taken the angular average over all direction of the galactic wind as

[TABLE]

which gives the distribution (35) in Sec. IV.1. Also we can read off the flux spectrum from Eqs. (68) and (70) as

[TABLE]

which gives the flux distribution (36).

We therefore obtain the hit rate per unit area on Earth surface

[TABLE]

where we take $m_{\mathrm{p}}=938.27$ MeV and the total hit rate is therefore

[TABLE]

Therefore, we can find the expected number of AQNs that have impac the Earth within 4.5 Gyr is

[TABLE]

Appendix D Estimation of AQN flux (fixed wind)

We also present calculation of AQN flux for fixed wind direction.

As shown in Fig. 6, The coordinate system is set up as shown in Fig. 6, which implies the following relation

[TABLE]

Starting with the general velocity distribution

[TABLE]

and noting that

[TABLE]

from the relation (77), we conclude

[TABLE]

and the flux of AQN is

[TABLE]

where in the last step we drop out all numerical factors as normalization will be applied in the simulation, and the $\phi$ component is simply integrated out due to the azimuthal symmetry of setup.

The integral over $\psi$ and $\varphi$ can be performed numerically by Monte Carlo method, as presented in Fig. 7.

Appendix E Heat emission profile simulation

This appendix describes the computation details solving for the heat emission profile $q(r,v_{\rm AQN})$. First, the simulation implements numerical solving of differential equations (19) using Python for individual AQNs, and for each of these AQN trajectories through the Earth the simulation returns the solution in the form of the following array

[TABLE]

where each row corresponds to a time $t$ for the trajectory at which sampling has been done, and denotes the radial position $\mathbf{r}(t)$, the speed $v(t)$, the rate of mass loss $\dot{m}(t)$ and the instantaneous mass $m(t)$, where $t_{i}$ and $t_{f}$ give the AQN’s time of entry into the Earth and the time of exit respectively.

By doing the histogram statistics of the first two columns with weighting factor $\dot{m}$ in the third column, we obtain the unnormalised $q(r,v_{\rm AQN})$. To do so, we used NumPy’s $hist2d$ function to bin the first two columns into a 200$\times$200 pixel grid using weights $\dot{m}(t)$. The simulation for an AQN trajectory was set to have a default sampling rate of $2$ sample points per second. For a given AQN trajectory in the phase space $(r,v_{\rm AQN})$, at sufficiently high velocities this would lead to gaps (empty pixels) being seen in the trajectory on the weighted 2D histogram plot. To tackle this, we use Python’s SciPy package to perform spline interpolation on our points for the radial distances $\mathbf{r}(t)$, instantaneous velocities $v(t)$ and instantaneous masses $m(t)$ on the interval [$t_{i}$,$t_{f}$] and demand a sufficient number of sample points such that for the highest $v(t)$ for the trajectory we would obtain at least one point per pixel. We use these and Eq. (17) to compute the corresponding $\dot{m}(t)$ (directly interpolating $\dot{m}(t)$ is unreliable due to it not being a smooth function). This procedure changes the number of sampling points per second and makes it different for every AQN trajectory, and hence when plotting $q(r,v_{\rm AQN})$ we multiply the weighting factor $\dot{m}$ by 1/(sample points per second) to obtain the new weighting factor $\dot{m}(t)$/(sample points per second) for points corresponding to each AQN. By plotting sufficiently many such weighted trajectories on the 2D histogram, we obtain $q(r,v_{\rm AQN})$. For our simulations, we have chosen to do this for $10^{6}$ AQNs.

To normalize $q(r,v_{\rm AQN})$, we should use the information in the fourth column. We first build up a new array for the total $N$ samples of AQNs

[TABLE]

with each row recorded the mass loss of one sample AQN. Then we can estimate $\langle\Delta m_{\rm AQN}\rangle$, the total mass loss per AQN on average, using the general formula

[TABLE]

We therefore obtain the normalized $q(r,v_{\rm AQN})$ from condition (30).

Appendix F On sensitivity of the main results to the AQN’s parameters

The main goal of this appendix is to argue that the main results of this work are not very sensitive to the parameters of the model, such as size distribution of AQN (represented by parameters $\alpha$ and $B_{\rm min}$), $\epsilon$ (describing the strength of the interaction of AQNs with material), flux distribution of AQNs (depending on the direction of galactic wind ${\bm{\mu}}$ and gravitation attraction from the Solar system $v_{\rm esc}\simeq 42.1\,\mathrm{km}\,\mathrm{s}^{-1}$).

Before presenting the simulation results, we describe the parameters considered in the work as summarized in Table 3. First, we choose the set of $(\alpha,B_{\rm min})$ in Table 1 such that $\langle B\rangle\gtrsim 10^{25}$, subjecting to the constraint (9) from IceCube. Next, we choose different values of $\epsilon=1$ and 3 to study the dependence of the effective cross section of AQN, as discussed in Sec. III.1. Additionally, we consider three different flux distributions of AQNs: isotropic averaging, fixed-wind direction, and acceleration by solar gravitation. Since the first two are fully described in the main text as well as Appendices C and D, we only comment on the third distribution here. The gravitation of the solar system implies an additional escape velocity $\sim 42.1\,\mathrm{km}\,\mathrm{s}^{-1}$ at Earth distance. To take into account of such effect, we add an additional magnitude to the initial speed for each AQN by $\sim 42.1\,\mathrm{km}\,\mathrm{s}^{-1}$ in the simulation based on the isotropic distribuiton (36).

As mentioned earlier, we find the simulation results are largely similar in all cases, with one moderate exception $\epsilon=3$. For the sake of brevity, we only present one common case [$(\alpha,B_{\rm min})=(2.0,3\times 10^{24})$, $\epsilon=1$ and fixed-wind distribution] and the case with $\epsilon=3$ here. As shown in Fig. 8, the heat emission profiles $q(r,v)$ is presented in the first row, following by by the axion density distribution $\rho_{a}(r)$ and the velocity spectrum $F(v_{\mathrm{a}})$ on Earth surface.

We only comment on the sensitivity to the effective cross section $\epsilon$ (column on the right), as the common case (column on the left) presents the same features already discussed in the main text. It is found for a larger $\epsilon$, the heat distribution become distorted in favor of a slower speed and more annihilation events closer to the Earth surface. This behaviour is well anticipated as annihilation becomes more efficient. However, these modifications of the heat emission profile $q(r,v)$ have almost no effect on the axion density distribution $\rho_{\mathrm{a}}(R_{\oplus})$ and the axion spectral velocity distribution $F(v_{\mathrm{a}})$ on the surface of the Earth as we discuss below.

We now move to the axion density distributions (Figs. 8(c) and 8(d)) and velocity spectrums (Figs. 8(e) and 8(f)), which are the main observable. We reach a conclusion that the simulation results are in general insensitive to almost all effects we considered. In fact, it is specifically insensitive to the baryon charge distribution of AQN models or the velocity distribution of DM flux. The effective cross section $\sim\epsilon$ modifies the heat distribution as mentioned above. However, this modification almost has no effect on the axion density $\rho_{\mathrm{a}}(R_{\oplus})$ and spectral properties $F(v_{\mathrm{a}})$ on the surface of the Earth, which are the relevant parameters for the observations.