Axion Quark Nugget Dark Matter: Time Modulations and Amplifications

TL;DR

This paper investigates axion production from Axion Quark Nugget (AQN) dark matter, focusing on unique short-term signal enhancements and their potential detection through time modulation and amplification effects.

Contribution

It introduces a novel short-term amplification mechanism for axion signals in the AQN model and evaluates its implications for detection strategies.

Findings

AQN model predicts short-term burst-like axion signal amplifications up to 10^2-10^3 times.

Annual and daily modulations are confirmed as key signatures for axion detection.

The study identifies a specific mass window for relativistic axions relevant to the AQN model.

Abstract

We study the new mechanism of the axion production suggested recently in [1,2]. This mechanism is 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. We perform numerical simulations to evaluate the axion flux on the Earth's surface. We examine annual and daily modulations, which have been studied previously and are known to occur for any type of dark matter. We also discuss a novel type of short time enhancements which are unique to the AQN model: the statistical fluctuations and burst-like amplification, both of which can drastically amplify the axion signal, up to a factor $\sim10^2-10^3$ for a very short period of time. The present work studies the AQN-induced axions within the mass window $10^{-6}{\rm\,eV}\lesssim m_a\lesssim10^{-3}\rm\,eV$ with typical velocities $\langle v_a\rangle\sim0.6c$. We also comment on the broadband detection strategy to search for such relativistic axions by studying the daily and annual time modulations as well as random burst-like amplifications.

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 : Comparison of potential enhancement: AQN-induced axions vs. WIMPs. Amplifications factors are listed up to order of magnitude estimate.

Potential enhancement	WIMPs	AQN-induced axions
Annual modulation	$1\%-10\%$ [53, 54]	$1\%-10\%$
Daily modulation	$\ll 1\%$ [54]	$1\%-10\%$
Statistical fluctuation	0	$20\%-60\%$
Local flashes	0	${10}^2-{10}^3$
Gravitational lensing	${10}^4-{10}^6$ [55, 56]	$\lesssim 1\%$

Table 3. Table 3 : Summary of statistical fluctuations ( B min = 3 × 10 24 subscript 𝐵 min 3 superscript 10 24 B_{\rm min}=3\times 10^{24} unless specified, 10 3 superscript 10 3 10^{3} trials involved in the two-step Monte Carlo). The uncertainties represent 1 ​ σ 1 𝜎 1\sigma significance. σ m subscript 𝜎 𝑚 \sigma_{m} is the standard deviation of ⟨ Δ ​ m AQN ⟩ delimited-⟨⟩ Δ subscript 𝑚 AQN \langle\Delta m_{\rm AQN}\rangle in the simulation.

$⟨B⟩$	$\alpha$	Other parameters	$⟨\Delta t⟩$ [s]	$⟨\dot{N}⟩⟨\Delta t⟩$	$⟨\Delta m_{\mathrm{AQN}}⟩$ [kg]	${\sigma }_m/⟨\Delta m_{\mathrm{AQN}}⟩$
$8.84\times {10}^{24}$	2.5	–	$38.3\pm 5.6$	$25.8\pm 5.0$	$\left(4.93\pm 1.26\right)\times {10}^{-3}$	25.6%
$8.84\times {10}^{24}$	2.5	$\mu =V_{\odot }-V_{\oplus }$	$41.8\pm 8.5$	$25.2\pm 5.0$	$\left(4.90\pm 1.17\right)\times {10}^{-3}$	23.8%
$8.84\times {10}^{24}$	2.5	$\mu =V_{\odot }+V_{\oplus }$	$35.2\pm 4.5$	$26.1\pm 5.1$	$\left(4.96\pm 1.34\right)\times {10}^{-3}$	27.1%
$8.84\times {10}^{24}$	2.5	Solar gravitation	$32.5\pm 3.8$	$25.0\pm 4.9$	$\left(4.97\pm 1.61\right)\times {10}^{-3}$	32.5%
$2.43\times {10}^{25}$	2.0	–	$37.5\pm 4.8$	$25.2\pm 5.0$	$\left(8.13\pm 4.51\right)\times {10}^{-3}$	55.5%
$4.25\times {10}^{25}$	(1.2, 2.5)	$B_{\min }={10}^{23}$	$62.3\pm 15.2$	$41.9\pm 6.6$	$\left(1.03\pm 0.49\right)\times {10}^{-2}$	47.6%
$1.05\times {10}^{26}$	(1.2, 2.5)	–	$35.4\pm 4.7$	$23.8\pm 4.9$	$\left(2.48\pm 0.98\right)\times {10}^{-2}$	39.5%

Table 4. Table 4 : Estimation of local flashes: the time duration, and the corresponding event rate as a function of amplification factor A 𝐴 A . Here we choose β = 1 𝛽 1 \beta=1 , ⟨ N ˙ ⟩ = 0.672 ​ s − 1 delimited-⟨⟩ ˙ 𝑁 0.672 superscript s 1 \langle\dot{N}\rangle=0.672\rm\,s^{-1} and ⟨ Δ ​ t ⟩ = 40 ​ s delimited-⟨⟩ Δ 𝑡 40 s \langle\Delta t\rangle=40\rm\,s .

$A$	Time Span	Event rate
1	10 s	0.3 ${\min }^{-1}$
$10$	3 s	0.5 ${\mathrm{hr}}^{-1}$
${10}^2$	1 s	0.4 ${\mathrm{day}}^{-1}$
${10}^3$	0.3 s	5 ${\mathrm{yr}}^{-1}$
${10}^4$	0.1 s	0.2 ${\mathrm{yr}}^{-1}$

Table 5. Table 5 : 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 6. Table 6 : Summary of some results ( B min = 3 × 10 24 subscript 𝐵 min 3 superscript 10 24 B_{\rm min}=3\times 10^{24} , ϵ = 1 italic-ϵ 1 \epsilon=1 unless specified, 10 8 superscript 10 8 10^{8} samples), the uncertainties in axion flux ⟨ E a ⟩ ​ Φ a delimited-⟨⟩ subscript 𝐸 𝑎 subscript Φ 𝑎 \langle E_{a}\rangle\Phi_{a} and density ⟨ ρ a ⟩ delimited-⟨⟩ subscript 𝜌 𝑎 \langle\rho_{a}\rangle describe the maximum daily modulation. Note that the AQN-induced energy flux and energy density are independent of axion mass in contrast to the conventional galactic axion models, see discussion at the end of Sec. IV.3 .

$⟨B⟩$	$\alpha$	Other parameters	$⟨\Delta m_{\mathrm{AQN}}⟩$ [kg]	$⟨\Delta B⟩/⟨B⟩$	$⟨E_a⟩{\Phi }_a\left[{10}^{13}\frac{\mathrm{eV}}{{\mathrm{cm}}^{2}s}\right]$	$⟨{\rho }_a⟩\left[{10}^{-6}\frac{\mathrm{GeV}}{{\mathrm{cm}}^3}\right]$
$8.84\times {10}^{24}$	2.5	–	$4.96\times {10}^{-3}$	$33.5\%$	$8.12\pm 0.93$	$4.51\pm 0.52$
$8.84\times {10}^{24}$	2.5	$\mu =V_{\odot }-V_{\oplus }$	$4.93\times {10}^{-3}$	$33.4\%$	$6.99\pm 0.72$	$3.88\pm 0.40$
$8.84\times {10}^{24}$	2.5	$\mu =V_{\odot }+V_{\oplus }$	$5.00\times {10}^{-3}$	$33.8\%$	$9.30\pm 1.17$	$5.17\pm 0.65$
$8.84\times {10}^{24}$	2.5	Solar gravitation777We 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 D for more details.	$4.94\times {10}^{-3}$	$33.4\%$	$9.67\pm 1.05$	$5.37\pm 0.58$
$2.43\times {10}^{25}$	2.0	–	$8.16\times {10}^{-3}$	$20.1\%$	$4.91\pm 0.43$	$2.73\pm 0.24$
$4.25\times {10}^{25}$	(1.2, 2.5)	$B_{\min }={10}^{23}$	$1.04\times {10}^{-2}$	$14.7\%$	$3.60\pm 0.20$	$2.00\pm 0.11$
$1.05\times {10}^{26}$	(1.2, 2.5)	–	$2.48\times {10}^{-2}$	$14.1\%$	$3.46\pm 0.16$	$1.92\pm 0.09$

Equations

$${\rho }_a^{\mathrm{AQN}}\sim 10^{-4}\left(\frac{\Delta B}{B}\right)\frac{\mathrm{GeV}}{\mathrm{c}{\mathrm{m}}^3}[\mathrm{AQN}-\mathrm{induced}],$$

$$⟨E_a⟩{\Phi }_a^{\mathrm{AQN}}(t)\simeq 10^{14}A(t)\left[\frac{\mathrm{eV}}{\mathrm{c}{\mathrm{m}}^2\mathrm{s}}\right],⟨E_a⟩\simeq 1.3m_a,$$

$$A_{(\mathrm{a})}(t)\equiv [1+{\kappa }_{(\mathrm{a})}\cos {\Omega }_a(t-t_0)],$$

$$A_{(\mathrm{d})}(t)\equiv [1+{\kappa }_{(\mathrm{d})}\cos ({\Omega }_{d}t-{\varphi }_0)],$$

$$\frac{⟨\dot{N}⟩}{4\pi R_{\oplus }^2}=\frac{0.4}{\mathrm{k}{\mathrm{m}}^2\mathrm{yr}}\left(\frac{10^{24}}{⟨B⟩}\right)\left(\frac{{\rho }_{\mathrm{DM}}}{0.3\frac{\mathrm{GeV}}{\mathrm{c}{\mathrm{m}}^3}}\right)\left(\frac{⟨v_{\mathrm{AQN}}⟩}{220\mathrm{km}{\mathrm{s}}^{-1}}\right).$$

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

$$⟨B⟩={\int }_{B_{\min }}^{10^{28}}\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}$$

$$\mu (t)\simeq V_{\odot }+b{V}_{\oplus }\cos {\omega }_{\mathrm{a}}(t-t_0)$$

$$m_a{\Phi }_a\simeq \frac{v_a}{c}\frac{⟨\dot{N}⟩⟨\Delta m_{\mathrm{AQN}}⟩}{16\pi R_{\oplus }^2},$$

$$⟨E_a⟩{\Phi }_a\sim 10^{14}\left[\frac{\mathrm{eV}}{\mathrm{c}{\mathrm{m}}^2\mathrm{s}}\right],$$

$$\mathrm{d}{m}_{\mathrm{AQN}}=-\sigma \rho v\mathrm{d}t=-\sigma \rho \mathrm{d}s,$$

$$\frac{{\sigma }_{\mathrm{entry}}}{{\sigma }_{\mathrm{exit}}}\simeq \frac{{\sigma }_{\mathrm{entry}}}{({\sigma }_{\mathrm{entry}}-\Delta \sigma )}\simeq \left(1+\frac{\Delta \sigma }{\sigma }\right)\simeq 1+\frac{2}{3}\frac{\Delta B}{B},$$

$$\mathrm{daily}\mathrm{modulations}\equiv {\kappa }_{(\mathrm{d})}\simeq \frac{1}{3}\frac{⟨\Delta B⟩}{⟨B⟩},$$

$$\mathrm{Prob}(N)\sim \frac{{\lambda }^N}{N!}{e}^{-\lambda },\lambda =⟨\dot{N}⟩⟨\Delta t⟩,$$

$$\frac{\mathrm{d}}{\mathrm{d}S}\Delta N_a$$

$$\simeq \frac{\beta }{2\pi d}\frac{1}{v_{\mathrm{AQN}}}\frac{{\dot{m}}_{\mathrm{AQN}}}{4m_a}$$

$$\mathrm{d}{N}_a\simeq \frac{1}{4m_a}\mathrm{d}{m}_{\mathrm{AQN}}=\frac{1}{4m_a}\frac{{\dot{m}}_{\mathrm{AQN}}}{v_{\mathrm{AQN}}}\mathrm{d}z,$$

$$\Delta {\Phi }_a(d)\simeq \frac{v_a}{\tau c}\frac{\mathrm{d}}{\mathrm{d}S}\Delta N_a=\beta \frac{v_a}{m_{a}c}\frac{{\dot{m}}_{\mathrm{AQN}}}{16\pi d^2},$$

$$A(d)\equiv \frac{\Delta {\Phi }_a(d)}{{\Phi }_a}\simeq \frac{\beta }{⟨\dot{N}⟩⟨\Delta t⟩}{\left(\frac{R_{\oplus }}{d}\right)}^2,$$

$$\tau \simeq \frac{2d}{v_{\mathrm{AQN}}}\simeq {\left(\frac{⟨\Delta t⟩}{⟨\dot{N}⟩}\right)}^{1/2}{A}^{-1/2},$$

$$\mathrm{Prob}(z\le d)\simeq {\left(\frac{d}{R_{\oplus }}\right)}^2\simeq \frac{A^{-1}}{⟨\dot{N}⟩⟨\Delta t⟩},.$$

$$\mathrm{Event}\mathrm{rate}$$

$$\simeq \frac{A^{-3/2}}{\sqrt{⟨\dot{N}⟩⟨\Delta t{⟩}^3}},.$$

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

$$\sigma \simeq \pi R^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,$$

$$m\frac{\mathrm{d}v}{\mathrm{d}t}=-\frac{1}{2}\sigma \rho v^2\simeq -\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),$$

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, 20, 21]. In this model, conventional dark-matter (DM) axions are produced either by the misalignment mechanism [22, 23, 24], when the cosmological field $\theta(t)$ oscillates and emits cold axions before it settles at a minimum, or via the decay of topological objects [25, 26, 27, 28, 29, 30, 31]. There are some uncertainties in the estimation of the axion abundance from these two channels and we refer the reader to the original papers for discussions111According to recent computations [31] 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 [30] on this matter. One should also emphasize that the computations [25, 26, 27, 28, 29, 30, 31] have been performed with assumption that Peccei-Quinn symmetry was broken after inflation.. In both production mechanisms, the axions are radiated as non-relativistic particles with typical galactic velocities $v_{\rm axion}/c\sim 10^{-3}$, and their contribution to the cosmological DM density scales as $\Omega_{\rm axion}\sim m_{a}^{-7/6}$. This scaling implies that the axion mass must be tuned to $m_{a}\simeq 10^{-5}$ eV in order to saturate the observed cosmological DM density today. Higher axion masses will contribute very little to $\Omega_{\rm DM}$ and lower axion masses will over-close the Universe which would be in strong conflict with cosmological data [32]. Cavity-type experiments have the potential to discover these non-relativistic axions.

Note that the type of axions being discussed here are conventional Quantum Chromo-Dynamics (QCD) axions with a mass range of ($10^{-6}{\rm eV}\lesssim m_{a}\lesssim 10^{-3}{\rm eV}$). It should be contrasted with Ultra Light Axions (ULAs) which were suggested as another Dark Matter candidate, also called fuzzy dark matter, but it is not the kind of axions considered in this study. ULAs do not solve the strong CP problem, and their primary motivation is not rooted in QCD.

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

Recent work [1] has suggested a fundamentally novel mechanism for axion production in planets and stars, with a mechanism rooted in the so-called axion quark nugget (AQN) dark matter model [35]. The AQN construction in many respects is similar to the original quark nugget model suggested by Witten [36] long ago (refer to [37] for a review). This type of DM model is “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 [36, 37]. First, there is an additional stabilization factor for the nuggets provided by the axion domain walls which are produced in great quantity during the QCD transition which help to alleviate a number of problems with the original nugget model222In particular, a first-order phase transition is not required as the axion domain wall plays the role of the squeezer. Another problem with [36, 37] 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 [36, 37], which must be stable at zero external pressure.. Another feature of AQNs 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 [35] because the DM density, $\Omega_{\rm DM}$, and the baryonic matter density, $\Omega_{\rm visible}$, automatically assume the same order of magnitude $\Omega_{\rm DM}\sim\Omega_{\rm visible}$ without any fine tuning. This is because they have the same QCD origin and are both proportional to the same fundamental dimensional parameter $\Lambda_{\rm QCD}$ which ensures that the relation $\Omega_{\rm DM}\sim\Omega_{\rm visible}$ always holds.

The existence of both AQN species explains the observed asymmetry between matter and antimatter as a result of separation of the baryon charge when some portion of the baryonic charge is hidden in the form of AQNs. Both AQNs made of matter and antimatter serve as dark matter. It should be contrasted with the conventional baryogenesis paradigm when extra baryons (1 part in $10^{10}$) must be produced during the early stages of the evolution of the Universe. 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. The AQNs may also offer a resolution to the so-called “Primordial Lithium Puzzle” [38] and to the “The Solar Corona Mystery” [39, 40] and may also explain the recent EDGES observation [41], which is in some tension with the standard cosmological model. It may also explain the DAMA/LIBRA annual modulation as argued in [42]. It is the same set of physical parameters of the model which were used in aforementioned phenomena that will be adopted for the present studies.

We refer to the original papers [43, 44, 45, 46] on the formation mechanism when the separation of baryon charge phenomenon plays the key role. The AQN framework aims at resolving two fundamental problems at once: the nature of dark matter and the asymmetry between matter and antimatter. This is irrespective of the specific details of the model, such as the axion mass $m_{a}$ or misalignment angle $\theta_{0}$.

AQNs are composite objects made out of axion field and quarks and gluons in the color superconducting (CS) phase, squeezed by a domain wall (DW) shell. It represents a cosmologicaly stable configuration state assuming the lowest energy state for a given baryon charge. An important point, in the context of the present study, is that the axion portion of the energy contributes to about 1/3 of the total mass in form of the axion DW. This DW is a stable time-independent configuration which kinematically cannot convert its energy to freely propagating (time-dependent) axions. However, any time-dependent perturbation results in the emission of real propagating axions via DW oscillations. The axion emission happens whenever annihilation events occur between antimatter AQNs with surrounding material.

In particular when AQNs cross the Earth the corresponding axion energy density is estimated [1] as:

[TABLE]

where $\Delta B/{B}$ is the portion of the baryon charge being annihilated during the passage of the AQN through the Earth. This portion $\Delta B/{B}$ is estimated on the level $(10\%-30\%)$ depending on the size distribution of the AQNs [2]. If the conventional galactic axions saturate the DM density $\rho_{\rm DM}\simeq 0.3\,\rm GeV\,cm^{-3}$ then the estimate (1) suggests that the AQN-induced axion density is about four orders of magnitude smaller than the conventional galactic axions density.

The key distinct feature between the conventional galactic axions and AQN-induced axions is the spectrum. In the former case $v_{a}\sim 10^{-3}c$, while for the latter, the spectrum is very broad and demonstrates a considerable variation in the entire interval $v_{a}\in(0,c)$ with average velocity $\langle v_{a}\rangle\simeq 0.6c$ [47]. This crucial difference in spectrum requires a different type of instruments and drastically different search strategies. A possible broadband detection strategy for the relativistic AQN-induced axions (1) has been discussed recently in an accompanying paper [48], where it is suggested to analyse the annual and daily modulations and different options are proposed in order to discriminate the true signals from spurious signals by using the global network. Appendix A reviews some of the ideas advocated in [48].

The main goal of the present work is the computation of the intensity and time-dependence of all the effects discussed in [48], assuming that the broadband detection strategy will be available in the future. It is clear that the AQN-induced axion density (1) is very small and the signal will be strongly suppressed in conventional cavity axion search experiments such as ADMX, ADMX-HF [49] or HAYSTAC [50] which explicitly depend on the axion density. However, in other type of experiments such as CASPEr [51] or QUAX [52] when the observables are proportional to the axion velocity $\mathbf{v_{a}}\sim\bm{\nabla}a$ the AQN-induced relativistic axions will clearly enhance the signal, as shown by Eq. (46) in Appendix A. In that case, it is the axion flux $\Phi^{\rm AQN}_{a}$, proportional to axion velocity $\mathbf{v_{a}}$, rather than the axion density becomes the relevant observable. In any case, axion density and flux are closely related: $v_{a}\rho_{a}^{\rm AQN}\approx m_{a}\Phi^{\rm AQN}_{a}$, and we present our basic numerical results in the main body of the text in terms of both.

We would like to emphasize that the main result of this work is the computation of the annual, daily modulations of the induced axions, as well as rare bursts-like amplifications. While the annual modulations have been discussed in the literature [53, 54], the daily modulations are normally ignored in DM literature because in WIMPs based models the effect is negligible. This is no longer the case with the AQN model when both effects are large, of the order of $10\%$. More specifically, we computed the AQN-induced axion flux on the Earth’s surface, which can be conveniently presented as follows

[TABLE]

where $A(t)$ is the modulation/amplification time dependent factor. We normalize the flux such that $\langle A(t)\rangle=1$ if averaged over very long period of time, much longer than few years.

The effect of the annual modulation has been known since [53, 54]. For the AQN-model we computed the annual modulation parameter $\kappa_{\rm(a)}$ defined as follows:

[TABLE]

where $\Omega_{a}=2\pi\,\rm yr^{-1}$ is the angular frequency of the annual modulation and label $``a"$ in $\Omega_{a}$ stands for annual. The $\Omega_{a}t_{0}$ is the phase shift corresponding to the maximum on June 1 and minimum on December 1 for the standard galactic DM distribution.

We have also computed the daily modulations with parameters defined as follows:

[TABLE]

where $\Omega_{d}=2\pi\,\rm day^{-1}$ is the angular frequency of the daily modulation, while $\phi_{0}$ is the phase shift similar to $\Omega_{a}t_{0}$ in (3). It can be assumed to be constant on the scale of days. However, it actually slowly changes with time due to the variation of the direction of DM wind with respect to the Earth. In addition to annual and daily modulations, we also show that the factor $A(t)$ can be numerically large for rare bursts-like events, the so-called “local flashes”. These short bursts resulting from the interaction of the AQN hitting the Earth in a close vicinity of a detector.

It is important to emphasize that, in the present work, along with the accompanying paper [48], we are studying the relativistic axions, $v_{a}\sim 0.6c$, which represent a direct manifestation of the AQN model. It should be contrasted with indirect manifestations of the AQN model mentioned above. An observation of axions with very distinct spectral properties compared to those predicted from conventional galactic axions with $v_{a}\sim 10^{-3}c$ would be a smoking gun for the AQN framework and this may then answer a fundamental question on the nature of DM.

The presentation is organized as follows. We start with a brief overview of the axion emission mechanism in the AQN framework in Sec. II. In Sec. III we compute the relevant parameters describing the modulations and amplifications as announced in this Introduction. The corresponding derivation of analytical equations and algorithm of numerical simulation are left to Sec. IV and V, with detailed results presented in VI. Finally, we conclude with some thoughts on possible future developments in Sec. VII.

II The AQN model and its axion emission mechanism

The AQNs hitting the Earth surface is given by [2]:

[TABLE]

Eq. (5) shows that conventional DM detectors are too small to detect AQNs directly. However, axions will be emitted when the AQN crosses the Earth interior, due to the annihilation processes that will lead to time-dependent perturbations of the axion DW. The time-dependent perturbations due to annihilation processes will change the equilibrium configuration of the axion DW shell, and axions will be emitted 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. The resulting emitted axions can be detected by conventional haloscope axion search experiments.

The emitted axion velocity spectrum was calculated in [47] using the following approach: Consider a general form of a domain wall:

[TABLE]

where $R_{0}$ is the radius of the AQN, $\phi_{w}$ is the classical solution of the domain wall, and $\chi$ describes excitation due to the time-dependent perturbation. $\phi_{w}$ is the DW time-independent classical solution, while $\chi$ describes on-shell propagating axions. Suppose an AQN is travelling in the vacuum where no annihilation events are taking place. The DW solution will remain in its ground (minimum energy) state $\phi(R_{0})=\phi_{w}(R_{0})$. Since there is no excitation (i.e. $\chi=0$), no free axion can be produced. When some baryon charge of the AQN is annihilated, the AQN starts loosing mass, its size decreases 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. The state of the domain wall then becomes $\phi(R_{0})=\phi_{w}(R_{\rm new})+\phi_{w}^{\prime}(R_{\rm new})\Delta R$, and the domain wall acquires a nonzero exciting mode $\chi=\phi_{w}^{\prime}(R_{\rm new})\Delta R$ leading to the production of free axions. Thus, whenever the domain wall is excited, corresponding to $\chi\neq 0$, freely propagating axions will be produced and emitted by the excited modes. The emission of axions is therefore an inevitable consequence of the annihilation of antimatter AQNs and of AQNs minimizing their binding energy. We refer the readers to Refs. [47] for the technical details of the emission mechanism and the calculation of its axion spectrum. For convenience, we also review the important results from [47] in Appendix B.

Let $\mathrm{d}N/\mathrm{d}B$ be the number of AQNs which carry the baryon charge [$B$, $B+dB$]. We shall use the same models for $\mathrm{d}N/\mathrm{d}B$ as in [2], and we refer to that paper for a description of the baryon charge distributions still allowed. Following [2], the mean value of the baryon charge $\langle B\rangle$ is given by

[TABLE]

where $f(B)$ is properly normalized distribution and the $\alpha$ is power-law index which assumes the following values:

[TABLE]

One should note that that the algebraic scaling (7) is a generic feature of the AQN formation mechanism based on percolation theory [46]. The 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, the parametrization (8) is based on fitting the observations of the Extreme UV emission from the solar corona as discussed in [40, 2].

Another parameter that defines the distribution function (7) is the minimum baryonic charge $B_{\rm min}$. We use the same models as discussed in [40, 2] and take $B_{\rm min}=10^{23}$ and $B_{\rm min}=3\times 10^{24}$. Therefore, we have a total of 6 different models for $f(B)$. In Table 1 we show the mean baryon charge $\langle B\rangle$ for each of the 6 models.

For simulations in this work, we will only investigate parameters that give $\langle B\rangle\gtrsim 10^{25}$ to be consistent with IceCube and Antarctic Impulsive Transient Antenna (ANITA) experiments as discussed in [2]. Therefore, in our numerical studies we exclude two models corresponding $B_{\rm min}\sim 10^{23}$ with power-law index $\alpha=2.5$ and $\alpha=2.0$.

III modulation and amplifications: WIMPs and the AQN-induced axions

In this section we quantify the time modulations and amplifications of the axion flux that can potentially be realized in the AQN model, and compare it to the conventional DM candidates such as weakly interacting massive particles (WIMPs). We examine time modulations and possible enhancements in descending order of importance, from annual to daily modulation effects [53, 54] specific to the AQN model as well as other new phenomena.

Our results are presented in Table 2, and the resulting comparison between AQN-induced axions and conventional DM enhancements can be summarized as follow:

The annual modulation discussed in Sec. III.1 has a similar amplitude in both cases, but the daily modulation presented in Sec. III.2 is much stronger for AQN-induced axions than for conventional WIMP-like models. Sec. III.3 and III.4 introduce two new, time-dependent phenomena, that are unique to the AQN model, and do not exist in conventional DM: the statistical fluctuations and the “local flashes” respectively. The latter effect becomes operational when the axion detector happens to be in the vicinity of the point where AQN enters or exits the Earth surface. We consider this novel phenomenon as the most promising effect which may drastically enhance the discovery potential for the axion search experiments and provide a decisive test of the AQN model. Gravitational lensing is another modulation effect, originally discussed in [55, 56]. We revisit this effect in subsection III.5 assuming a conventional galactic DM distribution (in coordinate and momentum spaces) which explains our result of a negligible effect in comparison with huge enhancement reported in [55, 56] where nonconventional streams of DM were considered.

III.1 Annual modulation

Conventional DM candidates such as galactic axions and WIMPs are affected by annual modulation with change estimated to be $\mathcal{O}(1\%-10\%)$ based on standard halo model (SHM) of the galactic halo [53, 54]. To understand this result, we first note the local speed of DM stream is not constant and subject to an annual modulation due to the motion of the Earth [53, 54]:

[TABLE]

where $V_{\odot}=220\,\mathrm{km}\,\mathrm{s}^{-1}$ is the oribital speed of the Sun around the galactic center, $V_{\oplus}=29.8\,\mathrm{km}\,\mathrm{s}^{-1}$ is the orbital speed of the Earth around the Sun, $\omega_{\rm a}=2\pi\,\rm yr^{-1}$ is the angular frequency of the annual modulation, and $|b|\leq 1$ is a geometrical factor associated with the direction of local velocity $\bm{\mu}$ relative to the orbital plane of Earth. Hence, it is natural to expect the amplification must be of order $\mathcal{O}(V_{\oplus}/V_{\odot})\sim 10\%$, as the incoming flux of particles depends on the incident speed $\mu$.

Similar arguments apply to the AQN-induced axions. We first note the flux of AQN-induced axions (derived in Sec. IV.3) is:

[TABLE]

where $v_{a}\sim 0.6c$ is speed of the emitted axions, $\langle\dot{N}\rangle$ is the expected hit rate of AQNs on Earth, and $\langle\Delta m_{\rm AQN}\rangle$ is the total average mass loss333For sake of brevity, in this work we adopt the following terminology: whenever “per AQN” is referred, we mean “per $\langle B\rangle$ baryon charge”. per single AQN with baryon charge $B$ such that $\Delta m_{\rm AQN}\approx m_{p}\Delta B$. Numerically, the corresponding flux of the AQN-induced axions is

[TABLE]

as the presented in Eq. (2) with $\langle A(t)\rangle=1$.

We should mention here that the exact formula to be derived in Sec. IV.3 contains some additional features (e.g. angular dependence), but it is sufficient to use Eq. (10) for qualitative discussion in this section. The linear relation (10) simply states that the the output axion flux rate $m_{a}\Phi_{a}$ is proportional to the amount of AQN flux supplied $\langle\dot{N}\rangle$ and its mass loss $\langle\Delta m_{\rm AQN}\rangle$. Clearly, the hit rate $\langle\dot{N}\rangle$ is, by definition, linearly proportional to the magnitude of incident speed $\mu$, resulting in a modulation up to order $\mathcal{O}(V_{\oplus}/V_{\odot})\sim 10\%$. On the other hand, we note the mass loss $\langle\Delta m_{\rm AQN}\rangle$ is a speed independent quantity by its conventional definition:

[TABLE]

where $\sigma$ is the effective cross section of the AQN, $\rho$ is local density of the Earth, $v$ is the speed of the AQN, $t$ and $s$ is propagation time and the path length of the AQN respectively. Thus, the farther an AQN travels underground the more axions emitted due to mass loss, but this is insensitive to the AQN speed.

To illustrate this, we plot the annual modulation fraction, defined as the size of the modulation amplitude relative to its average flux density, as a function of time in Fig. 1. In this work, we choose two extreme cases $\mu=V_{\odot}\pm V_{\oplus}$ to demonstrate the existence of annual modulation in our numerical simulations.

III.2 Daily modulation

The effect of daily modulation is rarely considered important in conventional discussion of DM candidates because the additional correction (at most $0.5\,\mathrm{km}\,\mathrm{s}^{-1}$ rotational velocity on the surface of the Earth) is much smaller than the Earth’s orbital velocity $V_{\oplus}$. However, the scenario is drastically different in case of axion emission from AQNs. The basic reason is that the AQNs are the composite microscopical objects444This should be contrasted with conventional DM candidates which are microscopic fundamental particles characterized by a cross section that is independent from orientation and position of the Earth with respect to Sun and galactic centre. which lose a portion of their material while crossing the Earth. As a result the size of an AQN decreases while it crosses the Earth.

The basic mechanism of the daily modulation of the AQN-induced axions can be describe as follow: Eq. (12) implies that the flux of AQN-induced axions is sensitive to the local mass loss of AQNs. Furthermore, the same equation suggests that the local mass loss is not constant during the passage of the AQN through Earth as the cross section $\sigma$ decreases along the trajectory. Consequently, the annihilation rate, and therefore the axion production rate, also decreases along the trajectory. Since the AQN cross section decreases as they traverse Earth and lose mass, there is then more axions emitted facing the wind compared to facing away from it. Fig. 2 illustrates this mechanism: the AQN flux impacts the Earth at a fixed angle $63\degree$ due to the alignment of the celestial equator (the plane of Earth’s equator) relative to the Galactic plane. Consider the downward motion as shown in Fig. 2. For this case we expect more heat (and therefore more axions) to be emitted in the upper hemisphere (the half sphere facing the AQN wind) than the lower one (the half sphere opposing the AQN wind). The difference can be as large as $\sim 10\%$ as we will soon estimate. Then, as the Earth rotates daily, we expect to see a daily modulation of order $10\%$ especially for detectors built in lower latitude regions with respect to the north pole.

In order to estimate the amplitude of this effect, we consider the ratio between the annihilation cross section when an AQN enters and when AQN exits the Earth (assuming the same orientation for the wind $\bm{\mu}$):

[TABLE]

where we use the relation $\sigma\propto R^{2}\propto B^{2/3}$. Therefore, the total fluctuation of daily modulation deviated from the mean value is half of the above estimation (13):

[TABLE]

which can be trusted as long as the factor on the right hand side is numerically small and expansion (13) is justified. From Monte Carlo simulation, we calculated that the typical fraction of total mass loss $\langle\Delta B\rangle/\langle B\rangle$ is about 30%. Therefore, we expect an amplitude modulation of order 10% from the mean value in mass loss. The estimate (14) is consistent with numerical simulations, later demonstrated in Sec. VI, which supports our interpretation in terms of the daily modulation. The corresponding daily modulation can be described with the parameter $A_{\rm(d)}(t)$ as defined in Eq. (4). $A_{\rm(d)}(t)$ can be assumed to be constant on the scale of days. However, it actually slowly changes with time due to the variation of the direction of DM wind with respect to the Earth’s position and orientation.

Lastly, we would like to point out another interesting consequence: there also exists a “spatial” modulation: axion flux is slightly more intense (the same $10\%$) in the Northern Hemisphere compared to the south, because the DM wind points to the northern portion of the Earth, see Fig. 2.

III.3 Statistical fluctuation

In this subsection we study a new type of time-dependent DM signal which has not been discussed before in the context of conventional DM: AQNs are macroscopically large objects with a number density approximately 23 orders of magnitude smaller in comparison to WIMPs with mass $\sim 10^{2}$ GeV. Therefore Poisson fluctuations, which are completely irrelevant for conventional DM, could potentially be important for the AQN model. The tiny flux (5) of AQNs is an explicit manifestation of this unique feature. The details of the AQN statistical fluctuation depend on the AQN size (and mass) distribution, trajectories and velocities.

In order to study the statistical fluctuation of AQN-induced axions, we consider the flux formula (10) and conjecture that the two quantities, $\langle\dot{N}\rangle$ and $\langle\Delta m_{\rm AQN}\rangle$, may have large statistical fluctuation due to the low number statistics. The average hit rate of AQNs hitting the Earth is about $1\rm\,s^{-1}$ [2], while the crossing time of an AQN inside the Earth interior is of order $\Delta t\sim R_{\oplus}/v_{\rm DM}\sim 30\rm\,s$. Hence, we estimate that there are only about 30 AQNs in the Earth interior at any given moment. Statistical fluctuation can therefore be important for AQNs, unlike conventional DM.

We numerically simulate the statistical fluctuations using a two-step Monte Carlo. First, we simulate $N$, the number of AQNs in the Earth interior at a given moment, by Poisson distribution:

[TABLE]

where $\langle\dot{N}\rangle=0.672{\rm~{}s}^{-1}$ is the average hit rate of AQNs [2], and $\langle\Delta t\rangle\sim 30$ s is the average time duration of an AQN inside the Earth555A subtle point is, $\langle\dot{N}\rangle$ is also proportional to $\mu$, the mean speed of AQN. Therefore for some model parameters (e.g. annual modulation) that modify $\mu$, $\langle\dot{N}\rangle$ should be also modified correspondingly.. And then, we obtain $\Delta m_{\rm AQN}$, the average mass loss per AQN from this sample size $N$. By repeating the above process, we obtain the standard deviation of $\langle\Delta m_{\rm AQN}\rangle$ where fluctuation in $\langle\dot{N}\rangle$ is already accounted.

The results are shown in the last columns of Table 3. We find statistical fluctuations of the order the ${\cal O}(20\%-60\%)$, depending on the size distribution model (7), (8). Numerically, the effect is greater than the enhancements assessed in the previous subsections for annual and daily modulations.

It is interesting to observe that smaller values of power-law index $\alpha$ corresponds to a larger fluctuation of mass loss $\langle\Delta m_{\rm AQN}\rangle$. Qualitatively this can be explained by the fact that a smaller value of $\alpha$ corresponds to higher AQNS average mass and consequently the rate of AQNs hitting Earth is reduced. This implies a larger dispersion, in agreement with explicit numerical computations presented in Table 3. In addition, we note that the statistical fluctuations of $\langle\Delta m_{\rm AQN}\rangle$ are almost insensitive to the AQN mean speed $\mu$ (DM wind).

III.4 Local flashes

Now we turn to the most interesting and most promising enhancement effect. Sharing a similar origin to the statistical fluctuation, we note that the detection signal of axion flux emitted by AQNs may be greatly amplified via a “local flash” on rare occasions when an AQN hits (or exits) the Earth surface in the vicinity of an axion search detector.

To understand the local flash, we first note that the intensity of axion flux is inversely proportional to distance square from the source. As estimated from the preceding section, there are only about 30 AQNs inside the Earth at any moment. Most of these AQNs do not loose much momentum and they travel a distance of order $R_{\oplus}$ along straight path without changing their original directions. Now, consider a case when an AQN is moving from a distance $d$ close enough to the axion detector. The intensity of the axion flux will be greatly enhanced by factor of $(R_{\oplus}/d)^{2}$ for a short period of time. This short enhancement in intensity is called a “local flash”. In the following, we will estimate the signal amplification factor of a local flash compared to the average axion flux induced by AQNs, and we will derive the instrumental requirements for a possible detection.

As shown in Fig. 3, suppose an AQN is moving close to the detector with a minimum distance $d$. The intensity of the emitted axion flux is maximized within $z\lesssim d$. Assuming the time is sufficiently short, the mass loss rate $\dot{m}_{\rm AQN}$ and velocity $\dot{z}\simeq v_{\rm AQN}$ can be treated as constants. We obtain the total number (per surface area) of emitted axions $\Delta N_{a}$ within $z\in[-z_{\rm cut},z_{\rm cut}]$:

[TABLE]

where $\beta$ is the angle related to $z_{\rm cut}$ as shown in Fig. 3, and we have implicitly used the relation

[TABLE]

when the integration $\mathrm{d}z$ is replaced by integration $\mathrm{d}m_{\rm AQN}$, see Sec. IV.3 for details. Finally, the axion flux density emitted by this AQN is

[TABLE]

where $\tau\simeq 2d/v_{\rm AQN}$ is the approximate travel time for an AQN inside the interval $[-z_{\rm cut},z_{\rm cut}]$.

The expression (18) represents the axion flux emitted by a single AQN travelling at distance $d\ll R_{\oplus}$ from detector. We want to compare this “local flash” with the average flux (10) by introducing the amplification factor $A$ defined as the ratio:

[TABLE]

where we approximated $\dot{m}_{\rm AQN}\simeq\langle\Delta m_{\rm AQN}\rangle/\langle\Delta t\rangle$. Using $\langle\dot{N}\rangle\langle\Delta t\rangle\sim 30$ (see Table 3), and for $z_{\rm cut}\sim d$, corresponding to $\beta\sim 1$, the typical amplification becomes significant if $d\lesssim 0.1R_{\oplus}$. The time duration $\tau$ of a local flash as a function of amplification $A$ is given by:

[TABLE]

where we assumed $v_{\rm AQN}\simeq 2R_{\oplus}/\langle\Delta t\rangle$ and $\beta\sim 1$ for simplicity. Table 4 shows various values for the time duration $\tau$ as a function of amplification factor $A$. We see that a standard detection signal without amplification ($A\sim 1$) lasts for 10 seconds, while a strongly enhanced signal amplified by $A\sim 10^{4}$ flashes for 0.1 second. However, it is a very rare event as it happens only once in every 5 years. More realistic amplification factors are somewhere between these two limiting cases as presented in Table 4.

Let us now estimate the event rate of a local flash for a given amplification $A$. The probability of observing an AQN for $z\leq d$ is given by:

[TABLE]

where we use Eq. (19) to express $d$ in terms of $A$. The event rate can be expressed in terms of amplification parameter $A$,

[TABLE]

where averages $\langle\dot{N}\rangle$ and $\Delta t$ have been numerically computed for different size distribution models, and can be found in Table 3.

Table 4 shows the event rate calculated for a few values of the amplification factor $A$. Specifically, we conclude there is about one event every two days for local flash amplified by $\sim 10^{2}$ if the detector has a time resolution of 1 second.

In conclusion, the amplification by local flashes is a unique feature of the AQN-induced axions. Moreover, these relativistic axions, with $v_{a}\sim 0.6c$, have very different spectral properties in comparison to the conventional DM candidates when $v_{a}\sim 10^{-3}c$. Therefore, even with an AQN-induced density (1) smaller than galactic DM axions the amplification could be sufficiently large because the AQN-induced flux could produce stronger signal when observables of the experiments is proportional to the axion velocity, see Eq. (46) in Appendix A and corresponding discussion in the Introduction. In addition, studying the time correlation of the local flashes can effectively distinguish the true signals from background noise as they are uncorrelated, while such approach does not apply to the conventional axion experiment and most cold DM searches because the distribution of cold DM halo is always uncorrelated and featureless in time. We refer the detailed study of correlation and broadband detection to Ref. [48].

III.5 Gravitational lensing

This subsection is partly motivated by Refs. [55, 56] where it has been claimed that the Sun or Jupiter can focus the flux of DM particles with speed $\sim(10^{-3}-10^{-1})c$ by an amplification factor up to $\sim 10^{6}$. The claim was based on two key assumptions:

The deflection angle $\gamma$ due to gravitational focusing is small, namely $\gamma\ll 1$; 2. 2.

The DM flux is colinear.

Assumption 1 requires that the bending angle caused by stars and planets is always very small. Assumption 2 strongly enforces the gravitational focusing as a result of (assumed) high level of coherency of the DM flux when all particles move in a highly colinear way with the same direction. To strengthen the focusing effect, Ref. [55] additionally assumed the DM particles are non-interacting and can pass through opaque objects such as the Sun and planets, because a transparent Sun has a shorter focal length by one order of magnitude compared to the case of opaque Sun (i.e. for interacting particles) and this results in a correspondingly stronger gravitational focusing correspondingly. Under these assumptions, calculations are greatly simplified and can be carried out analytically.

It is important to point out that the two assumptions are not well justified in SHM due to the misalignment of the ecliptic plane in the Milky Way. Therefore to satisfy the requirements of perfect focusing alignment and high coherency in propagating direction, Refs. [55, 56] considered special streams of slow-moving particles originating from distant point-like sources, such as stars, distant galaxies, and cluster of galaxies, etc..

One can imagine that these two assumptions also apply to the case of AQN-induced axions as conjectured in the work [1]: Assuming the existence of special emitting source of AQNs for ideal gravitational focusing, one could suspect that the stream of AQNs gains enhancement (up to $\sim 10^{6}$) by gravitational lensing666In Ref. [1], it is discussed the amplification factor can even go up to $\sim 10^{11}$ given the Sun as gravitational lens for source at a very specific distance., and consequently the flux of emitted axions is enhanced by the same amplification factor when the AQN stream impacts the Earth. However, the existence of such special emitting source is not of interest in the present work because a SHM is already assumed and implemented. Further investigations carried out in this work shows that in SHM the two assumptions are invalid for AQNs, and therefore the AQN-induced axions. Similarly, the two assumptions are also strongly violated for conventional WIMPs. Therefore, our conclusion is that amplification by gravitational focusing generally do not apply to DM particles in SHM, but the ideas advocated in Refs. [55, 56] are not excluded because nonconventional streams of DM were considered.

In what follows, we clarify the reasons that the two assumptions are violated for SHM. First, assumption 1 is invalid because of the inclination of the ecliptic plane relative to the DM wind: there is a $60\degree$ angle between the ecliptic plane and the Galactic plane, facing the DM wind direction. Consequently, for DM particles to be gravitationally focused to Earth, the deflection angle will have to be of the same order of magnitude. This is illustrated on Fig. 4. Such large deflection angle are impossible to obtain with lenses like normal stars and planets for the large majority of incoming DM particles, consequently, calculation in Refs. [55, 56] is no longer applied as the gravitational lensing is strong. Next, assumption 2 is violated because AQNs (and in fact most DM candidates) have a large velocity dispersion $\sim 110\,\mathrm{km}\,\mathrm{s}^{-1}$ that is comparable to their mean galactic velocity. This large velocity dispersion is requirement of the Virial theorem. Thus, the enhancement of DM flux by gravitational lensing as advocated in Refs. [55, 56] has a very narrow window for applicability, even for conventional DM candidates like WIMPs, for which both assumptions becomes invalid in the present framework of SHM.

The above arguments have been verified by analytical and numerical calculations, and it is discussed in Appendix C. Fig. 5 shows the resulting enhancement factor $A$ caused by gravitational lensing, where $A$ is plotted as a function of the velocity spectrum of AQN flux. Noting that $A$ implicitly depends on $\theta$ and $\sigma$: the angle between its incident direction and the plane of the Solar system, and the dispersion in velocity respectively.

In the realistic case $(\theta=60\degree,\sigma=110\,\mathrm{km}\,\mathrm{s}^{-1})$ such that both assumptions are violated, the amplification factor never exceeds $10^{-2}$, so the actual enhancement by gravitational lensing is small and can be neglected in the current study. In addition, we observe the existence of a resonance in the parameter space, namely when $\log_{10}(v/c)\sim-2.7$, the amplification is maximized. This is because the deflection angle is very large and only velocity within a specific range will be deflected within such narrow window $\gamma\simeq 60\degree$.

Note that if we considered an ideal case when the assumption 1 is externally enforced $(\theta=0\degree,\sigma=110\,\mathrm{km}\,\mathrm{s}^{-1})$ we indeed see some amplification. This happens because the DM wind-solar system alignment is perfect and the amplification indeed becomes much stronger (and grows as a power law of $v$), see orange line on Fig. 5. Furthermore, if assumption 2 is also enforced $(\theta=0\degree,\sigma=1\,\mathrm{km}\,\mathrm{s}^{-1})$ which is precisely the case considered in [55, 56], the AQN flux is perfectly aligned and highly colinear, see green line on Fig. 5. In this case we indeed find a strong amplification $\sim 10^{4}$, in agreement with computations of [55, 56]. However, it is very hard to imagine how such conditions can be justified, and how a colinear stream of DM could be formed in realistic (not ideal) galactic environment at least in SHM.

The results of the present subsection III.5 were discussed in [57], referring to the amplification $\sim 10^{11}$ mentioned in [1]. This factor is well above the enhancement $10^{4}$ derived in the present work. There is no contradiction: [1] did not derive the factor $10^{11}$, but it was mentioned as a possible strong enhancement motivated by the analogy with monochromatic coherent electromagnetic emission with $\lambda\simeq 1\,\mu\rm m$ from a distance source. This case of a monochromatic EM radiation emitted by a distance source has drastically different physics from the smooth distribution of ordinary DM assumed in present work. Furthermore, the authors of [56, 57] specifically distinguish “invisible DM” from “ordinary DM”, the latter being the subject of the present work. This difference is reflected in the assumptions 1 and 2 formulated above.

To conclude this subsection: we do not expect any substantial amplifications due to the gravitational lensing, at least within SHM. It is quite possible that future studies accounting for a number of other effects which are ignored in the present SHM (such as accounting for planets) may modify our conclusion. However, this topic is well beyond of the present studies.

IV Annihilation modeling

In this section we present all the equations that will be used in our simulations and calculations of the axion-density distribution around the Earth and the flux spectrum of axions passing through the interior of the Earth.

IV.1 Annihilation of AQNs impacting the Earth

The equations of motion describing the AQNs impacting the Earth are derived in Ref. [2]. We give a brief overview here.

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

[TABLE]

where $E=mv^{2}/2$ is the kinetic energy, $s$ is the path length, $\rho$ is the density of the local environment, $v$ is the AQN velocity, and $\sigma$ is the effective cross section of the AQN. In case of AQN annihilation within the Earth, the surrounding material is very rigid. Therefore, the cross section is the geometrical cross section,

[TABLE]

where $R$ is the radius of the AQN. Equation (24) implies that all nuclei which are on the AQN’s path will get annihilated while the AQN traverses the Earth. We refer readers to Ref. [2] for technical details and justification of this argument.

To proceed, following from the definition, we can express the time derivative of $E$ in the form

[TABLE]

where in the last step, we utilize the conventional rate of mass loss (12). Comparing Eqs. (23) and (25), we conclude that the mass loss acts like a friction term in the equation of motion

[TABLE]

The complete dynamical equations of motion in vector form are:

[TABLE]

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

[TABLE]

In Eq. (28a), we adopt the parameters in [40, 2], in which $\rho_{\mathrm{n}}=3.5\times 10^{17}$ kg m*-3* is the nuclear density and $m=m_{p}B$ is the AQN mass. In Eq. (28b), 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 5. 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. We make the approximation that the density within each layer is uniform, and therefore take the mean value from density at the top/bottom of the shell as the local density. The data in Table 5 is taken from Ref. [59].

IV.2 Axion-emission spectrum in the observer frame

We assume that the emission of AQN-induced axions is relativistic with velocity $v_{a}\sim 0.6c$ and dominantly spherically symmetric in the AQN frame [1, 47]. In fact, spherical symmetry is well preserved even in the observer frame due to the relatively slow speed of AQNs $v_{\rm AQN}\sim 10^{-3}c\ll v_{a}$. To demonstrate this, we analyze the 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. [47] 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}_{a})$:

[TABLE]

where $\tilde{v}_{a}$ and $\tilde{E}_{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 [47]. As claimed in the beginning, annihilation of AQN is assumed to preserve spherical symmetry in the AQN rest frame, therefore only the $l=0$ mode is considered in Eq. (29). The parameter $\delta\in(0,1)$ is a convenient factor introduced in Ref. [47] 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 (29) that do not exceed $\sim 20\%$. The velocity spectrum is not known to better precision (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}_{a}$, we know $\rho(v_{a})$ is normalized within a unit 3-ball $B^{3}$

[TABLE]

where in the second step, we transform our coordinate from $\tilde{v}_{a}$ to $v_{a}$ using Eqs. (30). Note that such transformation produces a small error related to a slight shift of spherical center by $\sim 10^{-3}c$. However, this inconsistency is negligible comparing to the uncertainty of approximation in terms of $\delta$. Noting that the last equality in Eq. (31) is completely expressed in terms of $v_{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}}_{a},\hat{\mathbf{v}}_{\rm AQN}\rangle$. From Eq. (30b), we expand Eq. (32) as

[TABLE]

where ${\cal O}(\frac{v_{\rm AQN}}{v_{a}})\sim 0.1\%$, the correction is negligible (comparing to the uncertainty of $\delta$) as claimed at the beginning of this subsection.

To summarize: the emission spectrum in frame of the observer (33) has an identical form, up to negligible correction, to the spherically symmetric spectrum in the rest frame (29) and computed in [47]. This is expressed in terms $H_{0}(p,\delta)$ given in Appendix B. One important fact to emphasize is, although spherical symmetry is well preserved in velocity spectrum of axion emission, the symmetry does not extend to the flux profile on Earth’s surface due to asymmetric effects such as the daily modulation (see Sec. III.2). In the next subsection, we introduce azimuthally symmetric function $P_{a}(\theta)$ to account for such effect.

IV.3 The axion flux density on Earth’s surface

The trajectory of an axion can be approximated as free motion because the gravity of the Earth is too weak to modify the relativistic axions: $v_{a}\sim 0.6c\gg v_{\rm esc}=11\,\mathrm{km}\,\mathrm{s}^{-1}$. The axion flux density on the surface of Earth is

[TABLE]

where $\langle v_{a}\rangle\simeq 0.6c$ is the average speed of emitted axion flux (see Appendix B for numerical estimation), $\langle\dot{N}\rangle$ is the total hit rate of AQNs to surface of Earth, $\langle N_{a}\rangle$ is the total number of axions emitted per AQN, and $P_{a}(\theta)$ is the azimuthal distribution of the axion flux due to daily modulation (see Sec. III.2). The first quantity is already estimated in [2]:

[TABLE]

The second quantity is proportional to the average mass loss per AQN $\langle\Delta m_{\rm AQN}\rangle$:

[TABLE]

where $\langle E_{a}\rangle\simeq 1.3m_{a}c^{2}$ is the average energy of axion emitted given by the spectrum (33), also see Ref. [47]. The last quantity assumes an azimuthal symmetry of system and is defined by the normalization condition:

[TABLE]

Lastly, the energy flux on Earth’s surface is

[TABLE]

Note that both the energy flux and energy density of the AQN-induced axions is independent of axion mass $m_{a}$, as the $m_{a}$ dependence in $\langle E_{a}\rangle$ and $\langle N_{a}\rangle$ cancels. This is a distinct feature comparing to conventional galactic axions which are $m_{a}$-sensitive.

V algorithm and simulation

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

The flux of relativistic axions is emitted instantaneously once an AQN looses mass through baryonic annihilation with matter on Earth. As shown in Fig. 6, to simulate the trajectories of AQNs through Earth we use the flux distribution function of AQNs in the vicinity of Earth with wind coming at fixed direction [2]:

[TABLE]

where $C$ is a normalization constant, $\mu$ is the mean speed of AQN wind as defined in (9), the notation of angles ($\theta$, $\psi$, and $\varphi$) is presented in Fig. 6.

To generate AQN initial conditions, we pick a point on the Earth’s surface determined by the flux distribution function - the integrand of (39). We sample initial conditions $v$, $\theta$, $\psi$ and $\varphi$ from the function defined on a 4D space of $v$, $\theta$, $\psi$ and $\varphi$. The function has no dependence on $\phi$ due to our approximation of the Earth as a sphere, so we set it to 0 for our simulation. We remove AQNs with initial velocities pointing away from the Earth as these correspond to AQNs that would have already experienced an Earth interaction. Following this sampling scheme, we generate $2\times 10^{5}$ AQNs and compute $q(r,\theta)$, the heat emission profile of axions at a given location within the Earth $(r,\theta)$, subject to the following normalization condition

[TABLE]

$\theta$ is the angle between $\mathbf{r}$ and the (opposite-direction) galactic wind $\bm{\mu}$, see Fig. 6. Note that the normalization condition of $q(r,\theta)$ is not defined in a standard way where the solid angle term $\sin\theta$ is implicitly absorbed into $q(r,\theta)$ for simplicity.

We also present a few AQN trajectories in Fig. 7. In the plot the AQNs are seen to be entering the Earth from one point on the Earth’s surface for plotting purposes. However, the initial entering velocities of the AQNs, and their positions on the Earth’s surface are determined by the flux distribution function and have been accounted for in our simulations. Their trajectories are determined by Eqs. (27). As expected, the vast majority of AQNs penetrate and escape from the Earth in linear motions due to the weakness of gravitation compared to the kinetic energy of AQNs ($v_{\rm esc}/v_{\rm AQN}\sim 0.1$). We also present a few trajectories of trapped AQNs (up to negligible portion) for illustrative purpose.

V.2 Computations of the surface axion flux density

As shown in Eq. (38) the magnitude of axion flux on the surface of Earth is determined by the mass loss $\langle m_{\rm AQN}\rangle$ simulated in the preceding subsection. However, the azimuthal component is directly related to $P_{a}(\theta)$, the surface probability of finding an axion at a given azimuthal angle $\theta$, given by Eq. (37). To obtain $P_{a}(\theta)$, we use Monte Carlo simulations based on the heat profile $q(r,\theta)$. First, it is clear that the number of axions to simulate at a given point at position $(r,\theta)$ is proportional to the rate of axion production [i.e. the heat emission profile function $q(r,\theta)$]. Therefore, we sample axion initial conditions from the distribution within the integral

[TABLE]

As shown in Fig. 8, we choose the initial position vector of emitted axion to be on the $y$-$z$ plane by azimuthal symmetry of $q(r,\theta)$. In addition, we know gravitation of Earth is negligible because axion has relativistic speed $v_{a}\sim 0.6c$ and largely exceeds the escape velocity $\sim 10^{-5}c$. Therefore given an emitted axion starting at $(r_{0},\theta_{0})$, its trajectory is simply a straight line of uniform motion

[TABLE]

where the initial velocity vector is sampling as follows:

[TABLE]

Here ‘Unif’ stands for uniform distribution. Note that the equation of motion (42) gives the intercept point(s) at the Earth surface $r=R_{\oplus}$ and hence the angle $\theta$ in the azimuthal distribution $P_{a}(\theta)$. Therefore solving for

[TABLE]

we find the angle $\theta$ on surface of the Earth:

[TABLE]

Note that the solution gives two roots in Eq. (44): one positive and one negative solution of $t$ in Eq. (42). By simulating sufficient numbers of axions based on initial conditions (42) and (43), we obtain the the azimuthal distribution $P_{a}(\theta)$.

As an additional note, the simulation is not efficient by picking out the positive-time solution only each time from the sampling array. In practice, we do not distinguish the two solutions and put both equally into statistics. In other words, we launch a pair of axions, instead of one only, with opposite initial velocity at position $(r_{0},\theta_{0})$ in the simulation code.

VI Discussion of results

VI.1 Heat emission profile

We first discuss the results of the simulation ($2\times 10^{5}$ samples) for the heat-emission profiles of axions $q(r,\theta)$. In our simulations we use a number of parameters such as the AQN baryon-charge distribution (parameters $\alpha$ and $B_{\rm min}$) and flux distribution of incoming AQNs (e.g. the annual/daily modulation and solar gravitation). 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,B_{\rm min})=(2.5,3\times 10^{24})$ shown in Fig. 9, while leaving other cases with different parameters to Appendix D. From Fig. 9, we observe the profile is in general as like a typical collection set of linearly moving trajectories in polar coordinate. This is as expected because the gravitational force is weak compared to the kinetic energy of AQNs.

The heat emission is greater near the upper pole $\theta=0$ with respect to the wind direction, despite of the fact that every annihilation event is assumed to be spherically symmetric. This is the cause of daily modulation explained in Sec. III.2: the AQN emits more heat when it enters Earth compared to when it leaves due to a larger initial cross section.

In addition, we observe abrupt changes of heat emission at some specific distance (e.g. $r\sim 0.9R_{\oplus}$ and $0.5R_{\oplus}$). One can see, from Table 5, that these jumps precisely correspond to the successive layers of the Earth interior. Whenever an AQN moves into a new layer with an abrupt change of local density, the mass loss drastically changes.

Another observation 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}$. This is a geometrical effect that stems from the fact that very few AQNs hit the surface with zero incident angle.

In Appendix D we study the sensitivity of our results to the parameters of the system. We observe that almost all the effects discussed in the paper are insensitive to specific choices for the parameters of the model such as the baryon-charge distribution of AQN models and the velocity distribution of incoming DM flux. Numerical simulations confirm this behavior.

VI.2 Axion flux density

The key results of our simulations are summarized in Table 6. The axion flux $\langle E_{a}\rangle\Phi_{a}$ is obtianed from Eq. (38) and the energy density $\langle\rho_{a}\rangle$ is the axion flux divided by $\langle v_{a}\rangle\simeq 0.6c$. We also present the azimuthal distribution $P_{a}(\theta)$ as defined in Eq. (37) for $(\alpha,B_{\rm min})=(2.5,3\times 10^{24})$. Similar to the conclusion in the preceding subsection, we find that the results are not sensitive to the modeling parameters: the axion flux $\langle E_{a}\rangle\Phi_{a}$ in general constrains within $(10^{13}-10^{14})\,\rm eV\,cm^{-2}s^{-1}$ depending on the size distribution, similar result applies for the density $\langle\rho_{a}\rangle\sim(10^{-6}-10^{-5})\,\rm GeV\,cm^{-3}$. While the main effect of modulations and enhancements on the axion flux have been discussed in Section III, here we want to make few additional comments on these results (we refer the reader to Appendix D for further technical details on the sensitivity to the parameters of the system).

Regarding the annual modulation, we explored two extreme cases of mean galactic velocity $\mu$ with $\mu=V_{\odot}\pm V_{\oplus}$ [see Eq. (9) and Table 6]. We find the flux is modulated by ${\cal O}(10\%)$ as expected from discussion in Sec. III.1. Regarding the daily modulations, the results are shown in Table 6 as the fluctuations $\pm$ in the last column. We find that the magnitude of daily modulation is of order ${\cal O}(10\%)$ and proportional to the ratio of average mass loss $\langle\Delta B\rangle/\langle B\rangle$, which is consistent with our simplified estimate (14).

We also studied the related effect on the azimuthal distribution $P_{a}(\theta)$ plotted on Fig. 10. It shows nearly linear dependence on $\cos\theta$ which is consistent with our interpretation because the mass loss is indeed proportional to the path length $s$ and therefore to $\cos\theta$, according to Eq. (12). We also note there is a noticeable jump near the upside pole ($\cos\theta=1$). This is interpreted as a sharp impact caused by the initial stage of the annihilation process.

Finally, it is instructive to compare the result from Table 6 to the order of magnitude estimate (1) presented in Ref. [1]. There is approximately two orders of magnitude deviation in numerical factor from that naive estimate. The reason for the discrepancy can be understood as follow: First, the estimate in Ref. [1] assumed that $\Delta B/B\sim 1$ (i.e., most incident AQN are completely annihilated). However, our simulations show that $\Delta B/B\sim 0.1$ in most cases, see Table 6. Similarly there is a numerical factor of 3/5 that is neglected in Ref. [1] as only the AQNs made out of antiquarks will be annihilated underground. Finally, there is a geometrical suppression factor related to averaging over inclination angles between the AQN velocity and the surface of the Earth that was ignored in [1].

VII Conclusions and future directions

The goal of the present work was to perform the Monte Carlo simulations describing the distributions of the AQN trajectories and annihilation events in the Earth’s interior in order to calculate the AQN-induced axions with $\langle v_{a}\rangle\simeq 0.6c$. This allowed us to compute a number of time dependent modulation and amplification effects which have been heavily used in [48]. Our goal was not to discuss the detection of these AQN-induced axions with a specific instrumentation, we refer to [48] where the detection issues were addressed. The main results of the present work can be summarized as follow:

a. We computed the axion flux, the spectrum and the angular distribution of the axions generated by antimatter AQNs crossing the Earth using Monte Carlo simulations. These axions are produced when the antimatter AQNs annihilate with the rocky material in the Earth interior. The main results are presented in Table 6.

b. We computed the time-dependent the annual and daily modulations in the axion intensity, see Table 2. Some of the predicted effects are unique to the AQN framework and constitute a decisive test of the model. The corresponding results have been heavily used in accompanying paper [48] where the broadband detection strategy has been proposed.

c. We also computed the time-dependent burst like enhancement, which we coined as the “local flash”. This amplification could be enormous $(10^{2}-10^{4})$, see Table 4. This strong enhancement occurs when AQN hits the earth close enough from the position of the axion search detector. This effect has been heavily used in accompanying paper [48] where a powerful test to discriminate the true axion signal from a spurious background noise, has been suggested.

The effects discussed in items a-c depend very weakly on the parameters of the model such as the nuggets size-distribution, consequently, there is little to no flexibility in the predictions presented in Tables 2, 4 and 6.

Why should we consider the AQN model seriously? From an observational cosmology angle, this model is consistent with all available cosmological, astrophysical, satellite and ground based constraints, where AQNs could leave a detectable electromagnetic signature. It provides a simple explanation for the observed relation $\Omega_{\rm DM}\sim\Omega_{\rm visible}$ and the baryon asymmetry without the need for fundamental changes or extensions of the standard model. The baryogenesis is replaced by “charge separation” effect. It was also showed that the AQNs can form and survive the very early epochs of the evolution of the Universe, such that they may serve as the DM candidates today. The same AQN framework may also explain a number of other (naively unrelated) phenomena, such as the excess of galactic emission in different frequency bands; it may also offer resolutions to some other astrophysical mysteries such as “Primordial Lithium Puzzle”, so-called “The Solar Corona Mystery”, the DAMA/LIBRA observed annual modulation, as well as the recent EDGES observations of a stronger than anticipated 21 cm absorption features, see Introduction for references and details.

Interestingly, the results of the present work could also be used for different purposes, not directly related to the axion searches, such as analysis of the annual modulation in the AQN-induced neutrino intensity. Such a study could be a key element in explanation [42] of the 20 years old DAMA/LIBRA puzzling observation of the annual modulation.

Acknowledgments

This work was supported in part by the National Science and Engineering Research Council of Canada. AZ thanks Yannis Semertzidis for explaining the role of different time scales (cavity storage time, axion coherence time etc) in axion search experiments. AZ also thanks many participants of the workshop “Axion Experiments in Germany”, August 19-22, 2019 and the “IBS Conference on Dark World”, Daejeon, Korea, November 4-7, 2019 where this work has been presented, for discussions and large number of good questions. 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. We thank the authors of ref. [57] for correspondence which resulted in clarification of our claims in Section III.5.

Appendix A Few comments on the broadband detection strategy.

The goal here is to highlight few ideas on broadband detection strategy as suggested in [48]. First of all, we remark that the axion field $a(\mathbf{r},t)$ can be treated as a classical field because the number of the AQN-induced axions (10), (11) accommodated by a single de-Broglie volume is very large in spite of the fact that the de-Broglie wavelength $\lambda$ for relativistic AQN-induced axions is much shorter than for conventional galactic axions,

[TABLE]

We start our overview with QUAX [52], CASPEr [51] ideas when the basic coupling is the interaction between the gradient of the axion field and the spin,

[TABLE]

In formula (46) we use, for generality, a single parameter $g_{a}$ which may assume the value $g_{\rm a}\equiv g_{\rm aee}$ for electrons in case of QUAX or $g_{\rm a}\equiv g_{\rm aNN}$ for nucleons in case of CASPEr. The coupling is proportional to the axion velocity $\mathbf{v}_{a}\propto\bm{\nabla}a(\mathbf{r},t)$ such that it is 3 orders of magnitude enhancement for the AQN axions with $v_{a}\sim 0.6c$ in comparison with galactic axions with $v_{a}\sim 10^{-3}c$. For the coupling (46) one can use a broadband detection technique if a single photon detectors for GHz range become available, which is claimed to be the case [60, 61].

There is another type of instruments such as ABRACADABRA [62], LC Circuit [63], DM Radio [64] and others which can operate in the resonance as well as in broadband regimes. A similar design [65] was also suggested as an instrument capable to measure the intensity of axion field ${a}(\mathbf{r},t)$ by using the Topological Casimir Effect. All these designs are based on the idea that the time-dependent axion field generates an additional electric current $\mathbf{j}_{a}\propto\dot{a}(\mathbf{r},t)\,\mathbf{B}$ in the presence of the background magnetic field $\mathbf{B}$. All the instruments which are mentioned above are very different in sizes and designs. However, their common feature is that these instruments are capable to measure the induced electric current in pickup loop [62] which is obviously represents a broadband capability not shared by the cavity type instruments.

We reiterate the main essence of the proposal. The conventional cavity detection technique search assumes a scanning to be sensitive to very narrow resonance line. The broadband strategy [48] can be used to probe the AQN emitted axions using the daily or annual modulation computed in the present work to select some specific frequency bin where modulation is observed. This broadband strategy allows to effectively remove the noise by fitting the data to the expected time dependent modulation (3) or (4). As the next step one can use a resonance based cavity experiment to scan a single frequency bin where modulation is observed to pinpoint a precise value for the axion mass.

Appendix B Spectral properties in the rest frame

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

[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. [66, 67] and recent article [68]. We also note a useful fact in Eq. (47b) 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. (47), the normalization factor obviously depends on $\delta$. We refer the detailed calculation to the work [47], while quoting results as follows

[TABLE]

In this work, we choose the intermediate value $\delta=0.5$ in the simulation, and we find the corresponding mean velocity of emitted axions is $\langle v_{a}\rangle\simeq 0.6c$.

Appendix C Detailed calculation of strong gravitational lensing

It is shown in Refs. [55, 56] the signal of dark matter (DM) detection can be amplified as large as $\sim 10^{6}$ by weak gravitational lensing for colinear stream of slow-moving particles. However, as shown in Sec. III.5 assumptions in Refs. [55, 56] do not apply to the case of AQNs in SHM. Unlike Refs. [55, 56] where the lensing is weak, we need to calculate strong gravitational deflection because there is a fixed $60\degree$ incident angle as shown in Fig. 4. Therefore, we will derive the amplification from first-principle Newtonian gravity. In what follows, we present a pedagogical approach from simple case to realistic situation.

C.1 Simple case: weak gravitational lensing and coherent velocity

We first analyze the simplest situation: a directed aligned flux with one and the same velocity $\mathbf{v}$, see Fig. 11. Because of the good alignment, we expect the lensing is weak and the amplification should reach a similar agreement with Refs. [55, 56].

Following the notation in Fig. 11. we first write down the equations of motion from first-principle Newtonian gravity:

[TABLE]

The set of differential equations describes a DM flux comes from a far distance $d_{\rm LS}$ with speed $v$ and impact parameter $b$. The Earth is approximated as 2D disk with radius $R_{\oplus}$ as we only concern the AQNs that have impact to the Earth. Solving for the differential equations, we have

[TABLE]

where the sign $\pm$ depends on the time before/after sthe distance of closest approach $r_{\rm min}$ defined as:

[TABLE]

Now, denoting a set of dimensionless parameters for convenience:

[TABLE]

Integrating out Eq. (50b), we obtain an explicit equation of $\theta$ as a function of $r$:

[TABLE]

As a quick check, note the deflection angle $\gamma$:

[TABLE]

which recovers to classic formula of (Newtonian) deflection angle in the weak lensing limit $b\gg\frac{GM_{\odot}}{v^{2}}$.

Now, we want to determine condition such that the flux will hit the Earth, that is

[TABLE]

Perturbatively solving the inequality (55), we conclude to hit the Earth the impact parameter $b$ must be constraint within a small ring $(b\pm\delta b)$:

[TABLE]

From Eqs. (56), the amplification factor $A_{0}(v)$ is just the ratio of the incident ring to the cross section without gravity:

[TABLE]

Therefore, we conclude the amplification (weak lensing and coherent velocity) is of order $\sim 10^{3}-10^{4}$ for a typical galactic speed $220\,\mathrm{km}\,\mathrm{s}^{-1}$. We conclude we reach a similar, but not the same, magnification as calculated in Refs. [55, 56]. But we warn the readers such comparison for illustrative purpose only as the assumptions in this work are different from Refs. [55, 56], as explained in the Sec. III.5.

Lastly, we comment on one subtle point regarding to the validity of equations (49). That is, the initial position of DM flux $d_{\rm LS}$ to be “sufficiently” far distance such that we can safely assume the initial velocity $v$ is not affected by gravity. The condition of validity can be expressed as:

[TABLE]

alternatively we can define a dimensionless parameter $\varepsilon$ as:

[TABLE]

Although it sounds to be a trivial point, the parameter $\varepsilon$ turns out to be an important parameter in the next subsection. As we will see, the amplification derived in the realist case will be expressed in terms of $\varepsilon$.

C.2 Realistic case: strong gravitational deflection and dispersive velocity

In reality the AQN flux is obviously dispersive and the gravitational deflection is large because of the fixed $60\degree$ incident angle, see Fig. 4. To justify this realistic case, we first consider the configuration as shown in Fig. 12. We assume initially there is a small tilted angle $\delta\psi$. We ask how big $\delta\psi$ is allowed? To answer this question, we trace the particle trajectory after its launch for some distance. Due to the existence the gravitation, the particle goes back to the parallel aligned case at $(z_{0},s_{0})$. Therefore, the only requirement is

[TABLE]

such that the particle can still hit the Earth.

Now we want to find $(z_{0},s_{0})$ as a function of $\delta\psi$. We modify the equation of motion in Eqs. (49) as

[TABLE]

In the last step, we use the fact $|z|\simeq d_{\rm LS}\gg s$. We present only the radial component $s$, which is of interest

[TABLE]

Solving this equation, we obtain

[TABLE]

where we define $\Delta s\equiv s_{0}-b$ and utilize the fact $v_{s,0}=0$ and $v_{s}=v\sin\delta\psi\simeq v\delta\psi$ for small $\delta\psi$ in the second step, and $s_{0}=b(v)+\mathcal{O}(\Delta s)$ in the last receptively. To ensure the particle still hit the Earth, we require $\Delta s\leq\delta b$. Therefore using Eqs. (56) and (63), we have the following inequality:

[TABLE]

where $\varepsilon\ll 1$ is a convenient parameter defined in Eq. (59).

From the flux distribution (39), we now define the angular flux spectrum as follows:

[TABLE]

Then a more realistic estimate of enhancement is therefore

[TABLE]

Substituting inequality (64) insto Eq. (66), we arrive

[TABLE]

The condition is clearly sensitive to the value of $\varepsilon$. Subjecting to the constraint (59), $\varepsilon$ has to be small. In Fig. 5 we choose the marginal value $\varepsilon=\frac{1}{3}$, namely the gravitational effect does not contribute to the initial total energy more than 25%. We note that for $\varepsilon$ greater than this value, the flux spectrum (65) may fail because strong gravitation becomes comparable to initial kinetic energy. Finally, we also comment that our estimation agrees with Ref. [69] although different assumptions and models are adopted.

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

The main goal of this appendix is argue that thee main resultss of this work are not very sensitive to the parameters of the model, such as size distribution of AQN (paramters $\alpha$ and $B_{\rm min}$) and flux distribution of incoming AQNs (e.g. the annual/daily modulation and solar gravitation).

Before presenting the simulation results, we describe the parameters considered in the work as summarized in Table 6. First, we choose the set of $(\alpha,B_{\rm min})$ in Table 1 such that $\langle B\rangle\gtrsim 10^{25}$, subjecting to the observational constraints by IceCube and and ANITA as discussed in Sec. II. Next, we consider two extreme cases of annual modulation by modifying the mean galactic velocity $\mu$ by adding (subtracting) the orbital speed of the Earth. In addition, we also consider the gravitation of solar system because it 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 flux distribution (39).

As mentioned earlier, we find the simulation results are largely similar in all cases. For the sake of brevity, we only present one case of annual modulation [$(\alpha,B_{\rm min})=(2.0,3\times 10^{24})$, and $\mu=V_{\odot}+V_{\oplus}$], and two cases with different choices of $(\alpha,B_{\rm min})$: $(2.0,3\times 10^{24})$ and $({\rm piesewise},10^{23})$ as shown in Figs. 13. Apparently, the heat emission profiles $q(r,\theta)$ and the azimuthal distributions $P_{a}(\theta)$ share similar features as presented in Figs. 9 and 10. We conclude the main results are not senstive to the parameters of model.