\section{Introduction}\label{sec:introduction}
 The  ANITA   collaboration have reported \cite{Gorham:2016zah,Gorham:2018ydl,Gorham:2020zne} observation of two anomalous events   that appear to be energetic cosmic showers emerging from the Earth with large exit angles.  We overview the corresponding events  below in details.
  We also highlight some difficulties in interpretation of these events in terms of the standard model (SM) physics and in terms of beyond-the-SM (BSM) physics. Specifically, some BSM explanations were proposed exclusively with a single goal to explain the observed ANITA anomalous events (AAEs).

 In the present work we advocate an alternative idea  that the ANITA anomalous events   could be a direct consequence of the axion quark nugget (AQN) model which was  invented long ago without any relation  to the ANITA observations. Rather, it was invented to
 to   explain the observed  similarity between the DM and the visible densities in the Universe, i.e. $\Omega_{\rm DM}\sim \Omega_{\rm visible}$ without any fitting parameters.
 Nevertheless, we will argue in this work that  the  observations  \cite{Gorham:2016zah,Gorham:2018ydl,Gorham:2020zne}, including the    frequency, intensity and  time duration  of the radio pulses    can be explained in terms     of the upward (Earth-emergent) AQN events. It is important to emphasize that the relevant parameters of the model have been fixed long ago for completely different purposes in a very different context in dramatically different systems, without any references to the ANITA anomalous events. It should be   contrasted with many recent proposals  when specific  BSM fields and interactions  were specifically introduced   to match the observations.

To date, the ANITA experiment has completed four flights \cite{Gorham:2016zah,Gorham:2018ydl,Gorham:2020zne} and reported two  anomalously steeply upward-going, radio-detected  CR-like air shower events that is compatible with a $\nu_\tau$ neutrino interpretation of energy $\sim\,$EeV at exit angles of $-27\degree$ and $-35\degree$ relative to horizontal in the first \cite{Gorham:2016zah} and third \cite{Gorham:2018ydl} flights respectively. The radio pulse of an ANITA anomalous event observed by the ANITA balloon payload at an altitude of $\sim{35}\,$km is of order $(0.1-1)\,$mV/m in electric field strength and  $(1-10)\,$ns in time duration. The observed frequency spectrum of the signal is in range $(40-800)\,$MHz, and attenuates sharply beyond the critical frequency near $800\,$MHz.

Now we overview some  suggestions \cite{Motloch:2016yic,Fox:2018syq,Romero-Wolf:2018zxt,Chauhan:2018lnq,Dasgupta:2018dzp,Heurtier:2019git,deVries:2019gzs,Hooper:2019ytr,Cline:2019snp,Shoemaker:2019xlt,Anchordoqui:2019utb,Smith:2020ecb,Anchordoqui:2018ucj,Dudas:2018npp,Huang:2018als,Connolly:2018ewv,Cherry:2018rxj,Anchordoqui:2018ssd,Collins:2018jpg,Chauhan:2018lnq,Heurtier:2019git,Hooper:2019ytr,Heurtier:2019rkz,Chipman:2019vjm,Abdullah:2019ofw,Anchordoqui:2019utb,Borah:2019ciw,Esmaili:2019pcy,Esteban:2019hcm,Altmannshofer:2020axr,Aartsen:2020vir}   to explain the ANITA anomalous events. We also  mention some of the difficulties which occur with many of these  proposals.

The AAEs are in critical tension with the standard model because neutrinos is exceedingly unlikely to traverse through Earth at a distance of $\gtrsim5\times10^3\,$km with such ultrahigh energy, even accounting for the $\nu_\tau$ regeneration \cite{Gorham:2016zah}. The analysis \cite{Fox:2018syq} reviewed the high-energy neutrino events from the IceCube Neutrino Observatory and inferred  that the $\nu_\tau$ interpretation is excluded by at least $5\sigma$ confidence. Similar study in Ref. \cite{Romero-Wolf:2018zxt} estimated ANITA acceptance to a $\nu_\tau$ flux, and concluded an at least two orders of magnitude above the upper limit from Pierre Auger Observatory and IceCube. More recently, IceCube also published severe constraints of astrophysical explanation for the AAEs under SM assumptions \cite{Aartsen:2020vir}.

Alternative   explanations such as transition radiation \cite{Motloch:2016yic,deVries:2019gzs} remains unconfirmed, and some are either disfavoured \cite{Dasgupta:2018dzp} or largely excluded by ANITA \cite{Shoemaker:2019xlt,Smith:2020ecb}. Notably, several BSM explanations are proposed \cite{Anchordoqui:2018ucj,Dudas:2018npp,Huang:2018als,Connolly:2018ewv,Cherry:2018rxj,Anchordoqui:2018ssd,Collins:2018jpg,Chauhan:2018lnq,Heurtier:2019git,Hooper:2019ytr,Heurtier:2019rkz,Chipman:2019vjm,Abdullah:2019ofw,Anchordoqui:2019utb,Borah:2019ciw,Esmaili:2019pcy,Esteban:2019hcm,Altmannshofer:2020axr}. In most cases, it suggests an origin of a (or a group of) massive hypothetical particle(s), which is strongly constrained by the IceCube and Auger bounds \cite{Cline:2019snp}. The other common problem is that the models are largely fine-tuned to match the observation of the AAEs and lack of a  natural motivation.

In the present work we put forward a proposal that the   AAEs  could be a direct consequence of the AQN model \cite{Zhitnitsky:2002qa}. The model   was originally invented to explain the observed  similarity between the DM and the visible densities in the Universe,  see Sec. \ref{AQN}  with more details on the model. The important feature of our proposal is that the AQN model involves no fine-tuning of parameters because this model is not devised to match the observation of AAEs, but rather its properties are well studied and constrained by numerous unrelated phenomena in the previous studies. Our presentation is organized as follows. In
subsection \ref{basics}  we overview the basics of the AQN model, while in subsections \ref{earth} and \ref{upward} we overview the features of the AQN model which will be   relevant for the present work.  We formulate our proposal on identification of AAE with AQN upward-going events in Sec. \ref{proposal}, while in Sec.  \ref{sec:Expected number of events}  we estimate the corresponding event rate
and in Sec \ref{sec:Radio signals induced by AQN} we
carry out explicit computations of the spectral properties of the radio pulse and estimate the intensity, which support our interpretation of the AAEs as Earth-emergent  AQN events.  We  conclude with Sec. \ref{conclusion}
where we explicitly formulate some dramatic differences between AQN   events and conventional CR   events. We also suggest possible tests   which may support or refute our proposal.

 \section{The AQN   DM model }\label{AQN}
 We start with few historical remarks and motivation of the AQN model in subsection \ref{basics},  while  in subsection
 \ref{earth} we overview recent observations (such as puzzling  bursts observed by the Telescope Array experiment) of some mysterious events which could be explained by the AQN events hitting the Earth. Finally, in section \ref{upward} we overview some  specific features of the AQNs  traversing the Earth (such as internal temperature, level  of ionization, etc).  These characteristics  will be important  for the present study  interpreting the ANITA anomalous events   as the Earth-emergent AQN events.

 \subsection{The basics}\label{basics}
   The   AQN DM  model  \cite{Zhitnitsky:2002qa} was  invented long ago with a single motivation
 to naturally explain the observed  similarity between the DM and the visible densities in the Universe, i.e. $\Omega_{\rm DM}\sim \Omega_{\rm visible}$ without any fitting parameters.
 The AQN construction in many respects is
similar to the Witten's quark nuggets, see  \cite{Witten:1984rs,Farhi:1984qu,DeRujula:1984axn},  and    review \cite{Madsen:1998uh}. This type of DM  is ``cosmologically dark'' not because of the weakness of the AQN 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 the  AQN model compared to the original models \cite{Witten:1984rs,Farhi:1984qu,DeRujula:1984axn,Madsen:1998uh}. First new element is the presence of the  axion domain walls which   are copiously produced  during the  QCD  transition\footnote{The axion field had been introduced into the theory to resolve the so-called the strong ${\cal CP}$ problem which  is related to the fundamental initial parameter $\theta_0\neq 0$.
  This source of ${\cal CP}$ violation  is no longer
available at the present time  as a result of the axion's  dynamics in early Universe.
One should mention that the axion  remains the most compelling resolution of the strong ${\cal CP}$ problem, see original papers
 on the axion \cite{axion1,axion2,axion3,KSVZ1,KSVZ2,DFSZ1,DFSZ2}, and   recent reviews \cite{vanBibber:2006rb, Asztalos:2006kz,Sikivie:2008,Raffelt:2006cw,Sikivie:2009fv,Rosenberg:2015kxa,Marsh:2015xka,Graham:2015ouw, Irastorza:2018dyq}.}. This domain wall plays a dual role: first  it serves  as  an   additional stabilization factor for the nuggets,     which helps to alleviate a number of  problems with the original nugget construction  \cite{Witten:1984rs,Farhi:1984qu,DeRujula:1984axn,Madsen:1998uh}.  Secondly, the same axion field $\theta (x)$ generates the strong and coherent $\cal{CP}$ violation in the entire visible Universe.

 This is because the $\theta (x)$ axion field before the QCD epoch could be thought as   classical $\cal{CP}$ violating field correlated on the scale of the entire Universe.
The axion field starts to oscillate at the QCD transition by emitting the propagating axions. However, these oscillations remain coherent on the scale     of the entire Universe. Therefore, the $\cal{CP}$ violating phase remains coherent on the same enormous scale.

  Another feature of the  AQN model which plays
absolutely crucial role for  the present work  is that nuggets can be made of {\it matter} as well as {\it antimatter} during the QCD transition.
Precisely the coherence of the $\cal{CP}$ violating field on large scale mentioned above provides a preferential production of one species of nuggets made of {\it antimatter}   over another
species made of {\it matter}. The preference    is determined by the initial sign of the $\theta$ field when the formation of the AQN   starts.
The direct consequence of this feature along with coherent $\cal{CP}$ violation in entire Universe  is that  the DM density, $\Omega_{\rm DM}$, and the visible    density, $\Omega_{\rm visible}$, will automatically assume the  same order of magnitude densities  $\Omega_{\rm DM}\sim \Omega_{\rm visible}$ without any fine tuning. We refer to the original papers   \cite{Liang:2016tqc,Ge:2017ttc,Ge:2017idw,Ge:2019voa} devoted to the specific questions  related to the nugget's formation, generation of the baryon asymmetry, and
survival   pattern of the nuggets during the evolution in  early Universe with its unfriendly environment.

One should emphasize that AQNs are absolutely stable configurations on cosmological scales. Furthermore, the antimatter which is hidden  in form of the very dense nuggets is unavailable for annihilation unless the AQNs hit the stars or the planets.

However, when the AQNs hit  the stars or the planets it may lead to observable phenomena. In particular,  the injection of the energy due to the AQNs hitting  the Sun
  may explain\footnote{In fact, to resolve this problem  Parker conjectured  long ago  \cite{Parker} that ``nanoflares"   are   identified with the  annihilation events in the AQN framework.
    The luminosity  of the Extreme UV (EUV) radiation from corona due to these annihilation events is unambiguously determined by the DM density. It is very nontrivial consistency check that  the computed  luminosity from the corona nicely matches with observed EUV radiation.     The same   events  of annihilation are  also   manifested themselves as
     the radio impulsive events in quiet solar corona as  recently recorded    by Murchison Widefield Array Observatory \cite{Mondal-2020}, which also represents a very nontrivial consistency check of the proposal  \cite{Zhitnitsky:2017rop,Raza:2018gpb} on the ``Solar Corona Mystery" resolution and accompanying radio impulsive events in quiet solar corona \cite{Ge:2020xvf}.}  the ``Solar Corona heating problem" as advocated in \cite{Zhitnitsky:2017rop,Raza:2018gpb,Ge:2020xvf}.
 There are also very rare events of annihilation in the center of the galaxy, which, in fact, may explain some observed galactic excess emissions in different frequency bands, including famous 511 keV line.

 The strongest direct detection limit\footnote{Non-detection of etching tracks in ancient mica gives another indirect constraint on the flux of   DM nuggets with mass $M< 55$g   \cite{Jacobs:2014yca}. This constraint is based on assumption that all nuggets have the same mass, which is not the case  as we discuss below.
The nuggets with small masses represent a tiny portion of all nuggets in this model, such that this constraint is easily  satisfied with any reasonable nugget's size distribution.} is  set by the IceCube Observatory's,  see Appendix A in \cite{Lawson:2019cvy}:
\begin{eqnarray}
\label{direct}
\langle B \rangle > 3\cdot 10^{24} ~~~[{\rm direct ~ (non)detection ~constraint]}.
\end{eqnarray}
Similar limits are   also obtainable
from the   ANITA
  and from  geothermal constraints which are also consistent with (\ref{direct}) as estimated in \cite{Gorham:2012hy}. It has been also argued in \cite{Gorham:2015rfa} that that AQNs producing a significant neutrino flux
in the 20-50 MeV range cannot account for more than 20$\%$ of the DM
density. However, the estimates \cite{Gorham:2015rfa} were based on assumption that the neutrino spectrum is similar to  the one which is observed in
conventional baryon-antibaryon annihilation events which typically produce a large number of pions and  muons and thus generate a
significant number of neutrinos and antineutrinos in the 20-50 MeV range where
SuperK has a high sensitivity. However, the critical difference in the case of AQNs  is that the annihilation proceeds within the
colour superconducting (CS) phase where the energetics are drastically  different \cite{Lawson:2015cla}. The main point is that, in most CS phases, the lightest pseudo Goldstone mesons
(the pions and kaons) have masses in the  20 MeV range,
rather than 140 MeV in hadronic phase. This dramatically changes
entire spectrum such that the  main assumption of \cite{Gorham:2015rfa} on similarity of the neutrino's spectrum in both phases is incorrect. The  resulting flux computed in \cite{Lawson:2015cla} is perfectly consistent with observations. Furthermore, precisely these low energy ($\lesssim 20~ \rm MeV$) AQN-induced neutrinos produced in the Earth's interior might be responsible  for explanation of the long standing puzzle of the DAMA/LIBRA  observation  of the annual modulation at $9.5\sigma$ confidence level as argued in \cite{Zhitnitsky:2019tbh}.

The authors of Ref. \cite{SinghSidhu:2020cxw} considered a generic constraint for the nuggets made of antimatter (ignoring all essential  specifics of the AQN model such as quark matter  CS phase of the nugget's core). Our constraints (\ref{direct}) are consistent with their findings including the Cosmic Microwave Background (CMB), Big Bang Nucleosynthesis (BBN), and others, except the constraints derived from    the so-called ``Human Detectors".
As explained in \cite{Ge:2020xvf}
  the corresponding estimates of Ref. \cite{SinghSidhu:2020cxw} are   oversimplified   and do not have the same status as those derived from CMB or BBN constraints\footnote{In particular, the rate of energy deposition was estimated in \cite{SinghSidhu:2020cxw} assuming that the annihilation processes between antimatter nuggets and baryons are similar to $p\bar{p}$ annihilation process. It is known that it cannot be the case because the annihilating objects have drastically different internal structures (hadronic phase versus CS phase). It has been also assumed in \cite{SinghSidhu:2020cxw} that  a typical X ray energy  is around  1 keV, which is  much lower than direct computations in the AQN model would  suggest \cite{Budker:2020mqk}. Higher energy X rays have much longer mean free path, which implies that the dominant portion of the energy will be deposited outside the human body. Finally, the authors of Ref.  \cite{SinghSidhu:2020cxw} assume that an antimatter nugget will result in ``injury   similar to a gunshot". It is obviously a wrong picture as the size of a typical nugget is only $R\sim 10^{-5} {\rm cm}$ while the most of the energy is deposited in form of the X rays on centimeter  scales \cite{Budker:2020mqk} without making a large hole similar to  bullet as assumed in \cite{SinghSidhu:2020cxw}. In this case a human's  death may occur as a result of  a large dose of radiation with a long time delay, which would make  it hard to identify the cause of the death. This argument (about the time delay and difficulties with possible death identification) should be contrasted   with the main assumption of \cite{SinghSidhu:2020cxw} that all such cases would be unambiguously and quickly identified. While in the journal version of \cite{SinghSidhu:2020cxw} the constraint had been weaken in comparison with the arXiv version, we still think it is not adequately reflect the complex physics as outlined  above.}.

 While ground based direct searches
offer the most unambiguous channel
for the detection of quark nuggets
the flux of nuggets   is inversely proportional to the nugget's mass   and
consequently even the largest available conventional DM detectors are incapable  to exclude    the entire potential mass range of the nuggets. Instead, the large area detectors which are normally designed for analysing     the high energy cosmic rays are much better suited for our studies of the AQNs as we discuss in next section \ref{earth}.

 \subsection{When the AQNs hitting the Earth...}\label{earth}
    For our present work, however,  the  most relevant studies  are related to the effects which  may occur when the AQNs made of antimatter hit the Earth  and continue to  propagate in deep underground in very dense environment.
   In this case the  most of the energy injection will   occur  in the Earth's interior.    The corresponding signals are very hard to detect  as the photons, electrons and positrons will be quickly absorbed by surrounding dense  material deep underground, while the emissions of the very weakly interacting neutrinos and axions are hard to recover. Nevertheless, as we already mentioned,   the AQN-induced neutrinos produced in the Earth's interior might be responsible  for explanation of the long standing puzzle of the DAMA/LIBRA  observation  of the annual modulation  \cite{Zhitnitsky:2019tbh}. The AQN-induced axions from deep interior can be recovered by analyzing    the daily and annual modulations  as suggested in \cite{Budker:2019zka} and elaborated in \cite{Liang:2020mnz}.      The AQN annihilation events in the Earth's atmosphere could produce  infrasound and seismic acoustic waves     as discussed in   \cite{Budker:2020mqk,Figueroa:2021bab}  when the infrasound and seismic acoustic waves indeed have been recorded  by  dedicated instruments\footnote{A single observed event properly recorded by  the Elginfield Infrasound Array  (ELFO) which was  accompanied by correlated seismic waves is dramatically different from conventional meteor-like  events. In particular, while the event was very intense it has not been detected by a synchronized  all-sky camera network (visible frequency bands)  which ruled out a meteor source. At the same time this event is consistent  with   interpretation of the    AQN-induced    event because  the visible frequency bands must be strongly suppressed when  AQN propagates in atmosphere \cite{Budker:2020mqk} }.  Furthermore,   the AQN annihilation events  my explain  the mysterious bursts  (representing a cluster of 3+ individual  cosmic ray-like events) observed  by Telescope Array (TA) experiment  \cite{Abbasi:2017rvx,Okuda_2019} in terms of the AQN annihilation events under the thunderstorm\footnote{These events are very unusual and cannot be interpreted  in
terms of  conventional cosmic ray (CR) single showers. In particular, if one tries to fit the observed bursts (cluster events) with conventional code, the energy for CR events should be in $10^{13}$ eV energy range to match the frequency of appearance, while the observed bursts correspond to $(10^{18}-10^{19})$ eV energy range as estimated by signal amplitude and distribution. Therefore, the estimated energy from individual events within the bursts is five to six orders of magnitude higher than the energy estimated by event rate  \cite{Abbasi:2017rvx,Okuda_2019}. Furthermore, the bursts events do not have sharp edges in waveforms in comparison with  conventional CR events. Also, all 10 recorded bursts occur under thunderstorm, which is very hard to interpret in terms of the conventional cosmic rays events which should not be modified by low-energy physics of thunderclouds.} as argued in \cite{Zhitnitsky:2020shd,Liang:2021wjx}.

Finally, the seasonal variations of the X ray   background in the near-Earth environment in the 2-6 keV energy range as observed by the  XMM-Newton at  $ 11\sigma$ confidence level  \citep{Fraser:2014wja} may be also naturally explained within the same AQN framework as argued \cite{Ge:2020cho}. This  application to the X ray emission in the near-Earth environment  is especially relevant for the present work because the AQN-induced  X rays   according to the proposal \cite{Ge:2020cho} are  originated from the AQN upward (Earth emergent) events when the AQNs traversed  through the Earth interior and exit the Earth surface. Such events could, in principle, be responsible for the ANITA mysterious events \cite{Gorham:2016zah,Gorham:2018ydl,Gorham:2020zne}   with exit angles  of $ -27\degree$ and $ -35\degree$   relative to the horizon    as advocated in present work. Before we present our arguments in next section we have to highlight the basic characteristics of the AQNs traversing the Earth, which is the topic of the next subsection.

     \subsection{Upward (Earth-emergent) events}\label{upward}

      The goal  here is to explain the basic features of the AQNs when they enter the dense regions of the surrounding material and annihilation processes start.
       The related  computations originally have been carried out in \cite{Forbes:2008uf}
 in application to the galactic environment with a typical density of surrounding visible baryons of order $n_{\rm galaxy}\sim 300 ~{\rm cm^{-3}}$ in the galactic center, in dramatic contrast with dense region in the Earth's interior  when $n_{\rm rock}\sim 10^{24} ~{\rm cm^{-3}}$. We review  these computations with few additional elements which must be implemented in case of propagation in the  Earth's atmosphere  and interior when the density of the environment is much greater than in the galactic  environment.

The total surface
emissivity from electrosphere has been computed in \cite{Forbes:2008uf} and it is given by
\begin{equation}
  \label{eq:P_t}
  F_{\text{tot}} \approx
  \frac{16}{3}
  \frac{T^4\alpha^{5/2}}{\pi}\sqrt[4]{\frac{T}{m}}\,,
\end{equation}
 where $\alpha\approx1/137$ is the fine structure constant, $m=511{\rm\,keV}$ is the mass of electron, and $T$ is the internal temperature of the AQN.
 One should emphasize that the emission from the electrosphere is not thermal, and the spectrum is dramatically different from blackbody radiation, see \cite{Forbes:2008uf},  see also Appendix \ref{details} with more details.

A typical internal temperature   of  the  AQNs for very dilute galactic environment can be estimated from the condition that
 the radiative output of Eq. (\ref{eq:P_t}) must balance the flux of energy onto the
nugget
\begin{eqnarray}
\label{eq:rad_balance}
    F_{\text{tot}}  (4\pi R^2)
\approx \kappa\cdot  (\pi R^2) \cdot (2~ {\rm GeV})\cdot n \cdot v_{\rm AQN},
\end{eqnarray}
where $n$ represents the number density of the environment.  The left hand side accounts for the total energy radiation from the  AQN's surface per unit time as given by (\ref{eq:P_t})  while
 the right hand side  accounts for the rate of annihilation events when each successful annihilation event of a single baryon charge produces $\sim 2m_pc^2\approx 2~{\rm GeV}$ energy. In Eq.\,(\ref{eq:rad_balance}) we assume that  the nugget is characterized by the geometrical cross section $\pi R^2$ when it propagates
in environment with local density $n$ with velocity $v_{\rm AQN}\sim 10^{-3}c$.

The factor $\kappa$ is introduced to account for the fact that not all matter striking the  AQN will
annihilate and not all of the energy released by an annihilation will be thermalized in the  AQNs by changing the internal temperature $T$. In particular,  some portion of the energy will be released in form of the axions, neutrinos and electron-positron pairs by the mechanism  discussed below.
  In a neutral dilute environment considered  previously \cite{Forbes:2008uf}  the value of $\kappa$ cannot exceed $\kappa \lesssim 1$ which would
correspond to the total annihilation of all impacting matter into thermal photons. The high probability
of reflection at the sharp quark matter surface lowers the value of $\kappa$. The propagation of an ionized (negatively charged) nugget in a  highly ionized plasma (such as solar corona)   will
 increase  the effective cross section. As a consequence,   the value of  $\kappa$ could be very large as discussed in \cite{Raza:2018gpb} in application to the solar corona heating problem.

 The  internal AQN temperature  had been estimated previously for a number of  cases.  It  may assume dramatically different values, mostly due to the huge difference in number density $n$ entering (\ref{eq:rad_balance}). In particular, for the galactic environment $T_{\rm galaxy}\approx 1$\,eV, while in
 deep Earth's interior it could be as high as  $T_{\rm rock}\approx  (100-200) $\,keV. Precisely this value of  $T$  had been used as initial temperature
of the nuggets in  the proposal \cite{Ge:2020cho} explaining the seasonal variations of the X rays observed the  XMM-Newton at  $ 11\sigma$ confidence level  \citep{Fraser:2014wja} at distances $r\sim (6-10) R_{\oplus}$ from the Earth surface. For our estimates in the present  work we shall use the same
 $T_{\rm rock}\approx  (100-200) $\,keV for explaining the ANITA anomalous events\footnote{\label{cooling}It is very unlikely for the temperature to reach  much higher  values    because a different cooling mechanism ($e^+e^-$ pair production) becomes very  efficient as $T$ reaches the region which is   relatively close to $m=511$ keV   such that the anticipated  suppression $\exp(-m/T)$   becomes less dramatic for $T\gtrsim 10^2$ keV
 as we discuss in Appendix \ref{details}.}. The difference with previous studies  \cite{Ge:2020cho} is that
 in this work we are interested in the instant  when  the very hot AQNs   just cross the Earth surface and enter the atmosphere (upward events).
 As we argue below the radio pulse is emitted precisely at this moment.  It should be contrasted with  the   applications considered in \cite{Ge:2020cho} when the  nuggets had travelled $\sim 10^4$ km in empty space before   emitted X rays  could be  detected by XMM-Newton observatory.

  One more feature   we want to mention here which is relevant for our studies is the ionization  property. Ionization, as usual, may occur in the system   as a result of  the high internal temperature $T$ discussed above. To be more precise,    the high  internal temperature $T$
  excites   a large number of positrons $\sim Q$ from electrosphere. These positrons  are    weakly bound particles which   can easily leave the system. The corresponding parameter $Q$ can be estimated as follows:
\begin{eqnarray}
  \label{Q}
Q\approx 4\pi R^2 \int^{\infty}_{0}  n(z, T)\rmd z\sim \frac{4\pi R^2}{\sqrt{2\pi\alpha}}  \left(m T\right)   \left(\frac{T}{m}\right)^{\frac{1}{4}} , ~~
  \end{eqnarray}
where $n(z, T)$ is the local density of positrons at distance $z$ from the nugget's surface, which has been computed in the mean field approximation  in \cite{Forbes:2008uf,Forbes:2009wg} and has the following form
\begin{equation}
\label{eq:nz0}
n(z, T)=\frac{T}{2\pi\alpha}\frac{1}{(z+\bar{z})^2},
\end{equation}
where $\bar{z}$ is the integration constant is chosen to match the Boltzmann regime at sufficiently  large $z\gg \bar{z}$. Numerical
studies ~ \cite{Forbes:2009wg}  support the approximate analytical  expression (\ref{eq:nz0}).
\begin{eqnarray}
\label{eq:zbar}
    \bar{z}^{-1}\approx \sqrt{2\pi\alpha}\cdot m  \cdot \left(\frac{T}{m}\right)^{\frac{1}{4}}, ~ n(z=0)\approx \left(mT\right)^{\frac{3}{2}}. ~~~
\end{eqnarray}
  In the equilibrium  with small  annihilation rate the positrons will normally occupy very thin layer  of order $\bar{z}$ around the  AQN's quark core as computed in \cite{Forbes:2008uf,Forbes:2009wg}. However, in our case when the AQN enters the Earth's atmosphere and further the interior   a large number of non-equilibrium processes (such as generation of the shock wave as a result  of  large Mach number $M$) are expected to occur. Furthermore,    the positron's cloud  is expected to expand well beyond the  thin layer  around the core's nugget
 as a result of  the  direct collisions with atmospheric molecules, in which case some positrons will be kicked off and   leave the system.

In what follows  we   assume that, to first order, the  finite portion of positrons  $\sim Q$   leave the system as a result of these complicated processes, in which case    the AQN as a system   acquires a negative  electric  charge $\sim -|e|Q$  and  get   partially  ionized,  see also Appendix \ref{details} with more details.

\section{ANITA anomalous  events as AQN Earth-emergent    events}\label{proposal}
With these preliminary comments  from previous section on the AQN features we  are  now in position to formulate the proposal  interpreting  the ANITA anomalous  events in terms of  the   upward moving   AQNs. We refer to some technical details to the Appendix \ref{details}, while here we highlight the basic idea which goes as follows.

The AQNs which propagate in the Earth's interior are very hot. Their temperature could be as hot as $T\approx (100-200)$\,keV as we discussed in
 \cite{Ge:2020cho} in the application to studies of the seasonal variations of the X rays observed   by XMM-Newton at very large distances from Earth, see footnote \ref{cooling}. At such large distances the AQNs (which traversed the Earth and continued to travel in empty space) are already sufficiently cold and emit mostly X rays with well defined spectrum. It should be contrasted with our present studies when we are interested
 in the same AQNs at the instant when the AQNs just crossed the surface and entered the Earth's atmosphere (upward-going events).  It is expected that these AQNs are still very hot with $T\approx (100-200)$\,keV, and their main cooling mechanism at the instant of surface's crossing  remains the same ($e^+e^-$ pair production).

 Exactly at this instant the electrons experiencing the repulsion force due to AQNs unscreened (negative) electric charge (\ref{Q}) which represents inevitable feature of the AQNs with high internal temperature $T$. The numerical value of the charge can be estimated as  (\ref{Q_eff}) where electric field is very large. This electric field may accelerate the produced electrons to very high energy   $\langle E\rangle \sim 10~ {\rm  MeV}$ as given by (\ref{E}). The same electric field has opposite effect for positrons when    the produced positrons will experience  the attractive  force and will assume a  location close to the nugget's surface.  One should emphasize that the  positrons produced by such mechanism have dramatically different properties from  the positrons from  electrosphere (\ref{density}) with much smaller bound energies and which could be localized far away from the AQN's surface\footnote{Precisely those weakly bound positrons could be kicked off as a result of  the elastic collisions with atmospheric molecule, which was  the topic of the proposal \cite{Zhitnitsky:2020shd} interpreting the mysterious burst events  detected by the Telescope Array Experiment  as the AQN events under a thunderstorm.}. One should mention here that the instant when AQN crosses  the boundary is very unique in a sense that all the characteristics of the  AQN (the internal temperature, the pressure, the Mach number $M=v_{\rm AQN}/c_s\gg 1$, where $c_s$ is the speed of sound in a given environment) remain the same as the AQN had before the crossing in a course  of   propagation in dense Earth's  interior with the density  $n_{\rm rock}\sim 10^{24} ~{\rm cm^{-3}}$. At the moment of crossing    the AQN suddenly enters
 a new environment with dramatically different density $n_{\rm air}\sim 10^{21} ~{\rm cm^{-3}}$.

 We expect that the result of this instant characterized by the  dramatic perturbation will be the  emission of energetic electrons with typical energy $\langle E\rangle \sim 10~ {\rm  MeV}$ as given by (\ref{E}). One should emphasize that such radiation  of photons (and accompanied emission of the $e^+e^-$ pairs)   is expected to be the dominant  cooling mechanism for the AQN propagating in the deep interior before the crossing. There are two  new elements (in comparison with a case when the AQNs propagate in deep interior with approximately constant density) which occur at the instant when the  AQNs cross  the hard surface and enter the Earth's atmosphere: 1. the emitted energetic electrons can propagate for several kilometres in atmosphere being the source  of the geosynchrotron radiation, i.e. they can produce  observable effects to be measured, in contrast with propagation in deep underground when the photons and the electrons will be quickly absorbed; 2.  the electrons are  mostly emitted   in the direction of the AQN velocity
 at the  instant of very short imbalanced  pressure at the moment of crossing.

Any precise computation of this instant of crossing is very hard problem  of non-equilibrium dynamics, which is beyond the scope of the present work. Fortunately, the observable radio  signal  (which is a  result  of the emitted electrons) is not very sensitive to the details of this non-equilibrium mechanism and corresponding time scales as it depends on several basic parameters such as typical energy $\langle E\rangle \sim 10~ {\rm  MeV}$ of the emitted electrons as estimated in (\ref{E}). Another parameter which enters all our formulae below is the number of emitted electrons $N$ which generates    a coherent radio signal. The parameter $N$ determines the intensity of the radiation, and cannot be computed from the first principles due to the very large uncertainties  of the  complicated non-equilibrium dynamics such as turbulence, shock waves, strong ionization as mentioned above and further  elaborated  in Appendix
 \ref{details}. We treat this parameter as a phenomenological unknown parameter which must satisfy constraint $N\ll N_{\rm max}$ where  $N_{\rm max}$ is the maximal number of potentially available electrons   (\ref{N}) which could be, in principle, liberated from the AQN at the instant of crossing the Earth's surface.

  The number of electrons $N$ which are emitted by a conventional cosmic ray showers of energy   $E_{\rm CR}\sim (10^{17}-10^{18})~ \rm eV$ is of order $N\sim (10^8-10^{9})$. We anticipate   a similar magnitude for the number of electrons $N\sim (10^8-10^{9})$ emitted by AQNs  as the observed intensity of the field strength  for the anomalous ANITA events of order $\rm mV/m$ which agrees with the value of $N$ assuming the  shower energy is $E_{\rm CR}\sim (10^{17}-10^{18})~ \rm eV$.

  In context of our work the  value  of $N$ is treated as the phenomenological unknown parameter, as we already mentioned. It  must be    much smaller than $N_{\rm max}$  which can be estimated (\ref{N}).  Important arguments supporting our proposal that the AQN induced events could mimic the anomalous radio signals  observed by ANITA are  based on very specific qualitative characteristics such as the  spectrum and pulse time duration, rather than on a precise estimation of parameter $N$. The  event rate  of such anomalous events is also shown to be  consistent  with our AQN-based interpretation, see next section \ref{sec:Expected number of events}.

 Furthermore, the average  electron's  energy in the energetic CR events with $E_{\rm CR}\sim (10^{17}-10^{18})~ \rm eV$ is around 30 MeV which is in the same energy range of the AQN-induced electrons as estimated in (\ref{E}). Therefore, the electrons which are released as a result of AQN crossing the boundary in upward (Earth emergent) event could mimic the radio signal of the conventional CR shower events as detected by ANITA. The estimation of the specific properties of the AQN-induced radio signal such as the spectrum, time pulse duration and the intensity    is the topic of      Sec.  \ref{sec:Radio signals induced by AQN}.

\section{Event Rate of AAE }
\label{sec:Expected number of events}
 This section is devoted to estimate the event rate of AAEs within the AQN framework. It is expected to be a qualitative estimate up to an order-of-magnitude check due to large uncertainties in parameters and rare occurrence of the observed AAEs. Nevertheless we would like to present such estimate to demonstrate that our interpretation of AAEs as a consequence of upward-going AQNs is at least a self-consistent proposal.

The expected number of the AAEs assuming that they are induced by the AQNs   can be  estimated as follows:
\begin{equation}
\label{eq:cal N}
{\cal N}
\approx {\cal{A}}_{\rm eff} {\cal T}\Delta\Omega
\frac{\rmd\Phi}{\rmd A\rmd\Omega}\,,
\end{equation}
where  ${\cal{A}}_{\rm eff}\approx4{\rm\,km}^2$ is the effective area of ANITA \cite{Gorham:2018ydl}, ${\cal T}\approx48.75\,$days is the combined exposure time of ANITA\footnote{\label{ft:exposure time}The effective exposure time was 17.25 days for ANITA-I \cite{Schoorlemmer:2015afa}, 7 days for ANITA-III \cite{Gorham:2018ydl}, and 24.5 days for ANITA-IV \cite{Gorham:2020zne}, where we exclude ANITA-II (28.5 days) as it is not sensitive to upward-going air showers.},  $\Delta\Omega\approx2\pi$ for isotropic flux of AQNs, $\Phi$ is the total hit rate of AQNs on Earth \cite{Lawson:2019cvy}:
\begin{equation}
\label{Phi}
\begin{aligned}
\Phi
&\approx2.12\times10^7{\rm\,yr}^{-1}  \\
&\quad\times\left(\frac{\rho_{\rm DM}}{0.3{\rm\,GeV\,cm^{-3}}}\right)
\left(\frac{v_{\rm AQN}}{220\kmps}\right)
\left(\frac{10^{25}}{\langle B\rangle}\right)\,,
\end{aligned}
\end{equation}
where $\rho_{\rm DM}$ is the local density of DM. The local rate of upward-going AQNs per unit area depends on the flux distribution of the AQN:
\begin{equation}
\label{Phi1}
\begin{aligned}
\frac{\rmd\Phi}{\rmd A\rmd\Omega}
=\frac{\eta}{4\pi R_\oplus^2}\Phi  =  4\cdot 10^{-2}\left(\frac{10^{25}}{\langle B\rangle}\right)\rm \frac{\eta ~events}{yr\cdot  km^2}
\end{aligned}
\end{equation}
where $R_\oplus=6371\,$km is the radius of the Earth, and $\eta$ is a parameter that characterizes the local flux distribution of AQN. Specifically, $\eta\approx1$  for isotropic distribution, and $\eta\approx2$ for the so-called fixed-wind distribution based on standard halo model so that more AQNs enter the Earth from the northern hemisphere and exit in the southern hemisphere  \cite{Liang:2019lya}. Note that the survival rate of an AQN traverses through Earth is also taken into account in $\eta$ implicitly, but it is an exceedingly minor effect comparing to the flux distribution because AQNs can penetrate through the Earth easily based on Monte Carlo simulation in Ref. \cite{Lawson:2019cvy}.

Combining the estimates above, one should have
\begin{equation}
\label{eq:cal N ::num}
{\cal N}
\approx0.28
\left(\frac{\eta}{2}\right)
\left(\frac{\rho_{\rm DM}}{0.3{\rm\,GeV\,cm^{-3}}}\right)
\left(\frac{v_{\rm AQN}}{220\kmps}\right)
\left(\frac{10^{25}}{\langle B\rangle}\right)\,.
\end{equation}
The expected number of events ${\cal N}\approx0.3$ is    almost one order of magnitude lower than  ${\cal N}_{\rm obs}=2$ events  observed by ANITA. Nevertheless, we consider this order of magnitude estimation  (\ref{eq:cal N ::num}) being  consistent with our proposal due to many uncertainties which
enter this estimate.

First, the parameters in Eq. \eqref{eq:cal N ::num} are in fact not precisely known. Essential parameters such as $\eta$, $\rho_{\rm DM}$, $\langle v_{\rm AQN}\rangle$, and $\langle B\rangle$ only have accuracy up to order one as the local flux distribution of DM and size distribution of AQN remain unknown to date.    In fact, there are numerous hints suggesting that   $\rho_{\rm DM}$ locally in solar system could be much larger from its canonical value, see a short comment on this with the references in the last paragraph   of this section.
Similarly, the effective area of detection ${\cal{A}}_{\rm eff}$ may double depending on the exit angle (see estimate in e.g. Ref. \cite{Huang:2018als}), and the effective exposure time is potentially longer than our conservative estimate if we take into account the ANITA-II flight in the estimate (see footnote \ref{ft:exposure time}).
In addition, the total number of observed AAEs could be a statistical fluctuation  due to its rare occurrence  (only 2).

We   consider this order of magnitude estimate (\ref{eq:cal N ::num}) as a highly nontrivial consistency check of our proposal  as the basic numerical factors entering  (\ref{eq:cal N ::num}) had been fixed from dramatically different physics (including solar corona heating puzzle) and  can easily deviate by large factor. More importantly, our main arguments leading to the  identification of the AAE with the AQN induced radio pulses
are based on specific qualitative features such as frequency dependence and duration of the pulse which are not sensitive   to these huge uncertainties  in the normalization factor (\ref{eq:cal N ::num}).  We consider the agreement between the observations and our theoretical estimates (to be discussed in next section \ref{sec:Radio signals induced by AQN})  for these specific  characteristics  as the  strong  arguments supporting our identification.

It is also interesting to note that the extra numerical factor $0.1$ (between computed and observed values)  which appears in our  order of magnitude estimates  (\ref{eq:cal N ::num}) is very similar to extra factor $0.1$ which occurred in analogous  computations \cite{Zhitnitsky:2020shd} of a  number of mysterious bursts  observed by the Telescope Array. This similarity hints on a common   origin of both phenomena though the  physics  for  these two phenomena  are dramatically different.  This is because  the main normalization factor representing  the   DM flux in the form of the AQN induced events   (\ref{Phi1})    is identically the same for both estimates, for AAE   (\ref{eq:cal N ::num}) and for Telescope Array mysterious event count \cite{Zhitnitsky:2020shd}.  If future studies support our identification of the AAE with the AQN induced radio pulses this numerical suppression factor might be a hint that the Standard Halo Model (which is used in  estimations for the DM flux)     underestimates the local DM density in solar system.  The true DM density locally  may dramatically deviate  from  average global value $\rho_{\rm DM}\approx 0.3{\rm\,GeV\,cm^{-3}}$, see Introduction in \cite{Liang:2020mnz} for the references and details.

\section{Radio signals induced by AQN}
\label{sec:Radio signals induced by AQN}
 It is well known that the frequency spectrum of geosychrotron radiation by an ultrarelativistic charged particle is equivalent to that emitted by a particle moving instantaneously at constant speed on an appropriate circular path with instantaneous radius of curvature $\rho$, see  e.g. Jackson \cite{jackson}:
\begin{equation}
\label{eq:rho}
\rho
\approx\frac{\gamma mc}{e{\cal B}\sin\theta_{\cal B}}
\approx 0.8{\rm\,km}
\left(\frac{\gamma}{20}\right)
\,,
\end{equation}
where ${\cal B}\approx0.5{\rm\,gauss}$ is the local magnetic field strength, and $\theta_{\cal B}$ is the angle between the particle velocity $\mathbf{v}$ and magnetic direction. We choose $\theta_{\cal B}\approx60\degree$ in this work, as the magnetic field direction in Antarctica is approximately vertical and the exit angle of the AAEs are typically of order $30\degree$.

The geometry follows from Fig. \ref{fig:geometry}, the segment of trajectories lies in the $x$-$y$ plane. The $\theta$ is the observation angle  between $\mathbf{v}$ and the direction of the observer $\mathbf{n}$. An ultrarelativisitc particle with Lorentz factor $\gamma\gg1$ has a narrow emission angle $\theta\lesssim\gamma^{-1}$, beyond which the intensity of radiation is exponentially suppressed. The $\boldsymbol{\epsilon}_\parallel$ and $\boldsymbol{\epsilon}_\perp$ are the two directions of polarization as shown on Fig. \ref{fig:geometry}.
\begin{figure}[h]
	\centering
	\captionsetup{justification=raggedright}
	\includegraphics[width=0.8\linewidth]{geometry}
	\caption{}.
	\label{fig:geometry}
\end{figure}

For an observer with distance ${\cal R}$   from $N$ coherent charged particles, the spectral component of the electric field $\mathbf{E}(\omega)$ as a function of  frequency  $\omega$ is given by \cite{jackson}
\begin{equation}
\label{eqs:abs E(omega) etc.}
\begin{aligned}
&|\mathbf{E}(\omega)|
=N\left(\frac{4\pi}{c}\right)^{1/2}\frac{1}{{\cal R}}|\mathbf{A}(\omega)|\,, \\
&\mathbf{A}(\omega)
=\frac{-ie\omega}{\sqrt{8c}\,\pi}\left[
-\boldsymbol{\epsilon}_\parallel A_\parallel(\omega)
+\kappa\,\boldsymbol{\epsilon}_\perp A_\perp(\omega)
\right]\,,
\end{aligned}
\end{equation}
where $A_\parallel(\omega)$ and $A_\perp(\omega)$ corresponds to the amplitudes of two polarization directions in terms of modified Bessel functions:
\begin{subequations}
\label{eqs:A_parallel etc.}
\begin{equation}
A_\parallel(\omega)
=i\frac{2\rho}{\sqrt{3}\,c}
\left(\frac{1}{\gamma^2}+\theta^2\right)K_{2/3}(\xi)\,,
\end{equation}
\begin{equation}
A_\perp(\omega)
=\theta\frac{2\rho}{\sqrt{3}\,c}
\left(\frac{1}{\gamma^2}+\theta^2\right)^{1/2}K_{1/3}(\xi)\,,
\end{equation}
\end{subequations}
with
\begin{equation}
\label{eq:xi}
\xi
=\frac{\rho\omega}{3c}\left(\frac{1}{\gamma^2}+\theta^2\right)^{3/2}\,.
\end{equation}
Note that Eqs. \eqref{eqs:A_parallel etc.} are only well-defined for $\omega>0$, the values in the negative domain is defined by   $\mathbf{A}(-\omega)=\mathbf{A}^*(\omega)$.
The parameter $\kappa$ in (\ref{eqs:abs E(omega) etc.}) is introduced to characterize the screening effect of $e^+e^-$ pair, where the $\boldsymbol{\epsilon}_\perp$-component is effectively cancelled ($\kappa\approx0$) when $e^+e^-$ pairs are predominantly formed in conventional CR events \cite{Huege:2003up}. In the opposite limit, we choose $\kappa\approx1$ in case of AQN-induced signal as it is primarily initiated by electrons such that the screening effect is diminished. The spectrum of electric field is therefore
\begin{equation}
\label{eq:abs E(omega)}
\begin{aligned}
|\mathbf{E(\omega)}|
&\approx\sqrt{\frac{2}{3\pi}}\frac{Ne\rho\omega}{c^2{\cal R}}
\left(\frac{1}{\gamma^2}+\theta^2\right)K_{2/3}(\xi)  \\
&\quad\times\sqrt{
1+\frac{\gamma^2\theta^2}{1+\gamma^2\theta^2}
\left(\frac{K_{1/3}(\xi)}{K_{2/3}(\xi)}\right)^2
}\,.
\end{aligned}
\end{equation}
The observational distance from ANITA balloon payload is of order 35 km and the size of the detector is of order 10 m, therefore the effective observation angle is tiny comparing to the emission angle $\theta\sim10^{-4}\ll\gamma^{-1}$.

Choosing $\theta=0$, we plot the spectrum \eqref{eq:abs E(omega)}
in Fig. \ref{fig:E_omega}  with different values of $\gamma$ for a specific value of $N=5\times 10^8$. One can see that the spectrum is very flat: the absolute value $|\mathbf{E(\omega)}|$ changes by a factor of 3 or so when the frequency $\nu$ varies by 2 orders of magnitude. Furthermore, the total strength  of the electric field integrated over entire frequency band agrees with the observed value on the level of  $|\mathbf{E}|\sim \rm mV/m$, see Refs. \cite{Gorham:2016zah,Gorham:2018ydl,Gorham:2020zne}.

Another  generic feature of synchrotron radiation is the exponential  suppression of the emission  beyond the critical frequency \cite{jackson}
\begin{equation}
\label{eq:nu_c}
\nu_{\rm c}
\equiv \frac{3\gamma^3c}{4\pi\rho}
\approx 0.7 {\rm\,GHz}
\left(\frac{\gamma}{20}\right)^2\,.
\end{equation}
This qualitative consequence of our  proposal is also consistent with
ANITA
observations\cite{Gorham:2016zah,Gorham:2018ydl,Gorham:2020zne}.

\begin{figure}[h]
	\centering
	\captionsetup{justification=raggedright}
	\includegraphics[width=1\linewidth]{E_omega_N58}
		\caption{Spectrum of electric field $|E(2\pi\nu)|$ from Eq. \eqref{eq:abs E(omega)}, with $\theta=0$, ${\cal R}=35\,$km, and $N=5\cdot 10^8$. The Lorentz factor is chosen to be $\gamma=10$ (blue), 20 (orange), and 60 (green).}
\label{fig:E_omega}
\end{figure}

The time-dependent radio pulse can be reconstructed from the frequency spectrum by an inverse Fourier transform:
\begin{equation}
\label{eq:bold E(t)}
\begin{aligned}
\mathbf{E}(t)
&=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^\infty b(\omega)\mathbf{E}(\omega)e^{-i\omega t}\rmd\omega  \\
&\approx -\boldsymbol{\epsilon}_\parallel
\frac{2Ne\rho}{\sqrt{3}\pi c^2\gamma^2 {\cal R}}
{\rm\,Re}\left[
\int_0^\infty b(\omega)\,\omega {\rm K}_{2/3}(\xi)e^{-i\omega t}\rmd\omega
\right]
\end{aligned}
\end{equation}
where $b(\omega)$ is the  filter   characterizing   the receiver,  similar to analysis in \cite{Huege:2003up}. For illustrative purpose, we plot the time-dependent electric field in Fig. \ref{fig:E_t3} by assuming an idealized rectangle filter spanning (40-80)\,MHz and (200-600)\,MHz. The pulse has an amplitude $|\mathbf{E}_0|\sim{\rm\,mV/m}$ and the time duration $\tau\sim{\rm ns}$, which is consistent with observed features of the anomalous pulses observed by ANITA \cite{Gorham:2016zah,Gorham:2018ydl,Gorham:2020zne}. Because $|\mathbf{E}(\omega)|$ is approximately flat in frequency range below $\nu_{\rm c}$, $\mathbf{E}(t)$ is essentially determined by the    inverse Fourier transform of $b(\omega)$. In case of a rectangle-like filter with frequency bandwidth $\Delta\nu$, the time duration of the pulse  is determined by $\Delta \nu$ as follows:
\begin{equation}
\label{eq:tau Delta nu}
\tau
\approx\frac{1}{\Delta\nu}
\approx 2 {\rm\,ns}\left(\frac{600{\rm\,MHz}}{\Delta\nu}\right)\,.
\end{equation}
The numerical value for  the time scale $\tau$ is not   very sensitive to the parameters of the model as explained above due to the flatness of the spectrum  below $\nu_{\rm c}$. In contrast to  the time duration   $\tau$      the absolute value of the electrical field  $|\mathbf{E}_0|\sim{\rm\,mV/m}$  is sensitive to $\gamma$ as shown on Fig. \ref{fig:E_t3}. It assumes the values  which
are also consistent with observations   \cite{Gorham:2016zah,Gorham:2018ydl,Gorham:2020zne}.

\begin{figure}[h]
	\centering
	\captionsetup{justification=raggedright}
	\includegraphics[width=1\linewidth]{E_t3}
		\caption{ Time-dependent electric field from Eq. \eqref{eq:bold E(t)} with $\theta=0$, ${\cal R}=35\,$km, and $N=5\cdot 10^8$, using an idealized rectangle filter. The Lorentz factor is chosen to be $\gamma=10$ (blue) and 20 (orange);
		filter: (40-80)\,MHz and (200-600)\,MHz.    }
	\label{fig:E_t3}
\end{figure}

\begin{figure}[h]
	\centering
	\captionsetup{justification=raggedright}
	\includegraphics[width=1\linewidth]{E_t}
		\caption{ Time-dependent electric field from Eq. \eqref{eq:bold E(t)} with $\theta=0$, ${\cal R}=35\,$km, and $N=5\cdot 10^8$, using an idealized rectangle filter. The Lorentz factor is chosen to be $\gamma=10$ (blue) and 20 (orange). Solid: filter spanning (40-80)\,MHz and (200-800)\,MHz; dashed: (200-800)\,MHz.  }
	\label{fig:E_t}
\end{figure}

 It is instructive to understand the temporal  features of the electric field $\mathbf{E}(t)$ by rewriting the integral (\ref{eq:bold E(t)}) in terms of the dimensionless variable $\xi$ defined by (\ref{eq:xi}) as follows:
\begin{equation}
\label{eq:bold E(t)_temp}
\begin{aligned}
\mathbf{E}(t)
=-\boldsymbol{\epsilon}_\parallel\frac{18}{\sqrt{3}\pi \rho}
\frac{\gamma^4Ne}{{\cal R}}{\rm\,Re}\left[
\int_{\xi_{\rm min}}^{\xi_{\rm max}}\xi e^{-ia\xi}{\rm K}_{2/3}(\xi)\rmd \xi
\right]
\end{aligned}
\end{equation}
where  we assume $b(\omega)$ as an idealized rectangle filter,  $a\equiv3\gamma^3c t/\rho$ and  $(\xi_{\rm min}, \xi_{\rm max})$ are determined by corresponding values of   $\omega_{\rm min}$ and $\omega_{\rm max}$ characterizing the filter $b(\omega)$. From   (\ref{eq:bold E(t)_temp}) one can explicitly see that the typical time duration is determined by parameter
\begin{equation}
\label{eq:tau}
 (a\xi)\approx 2\pi~~\Rightarrow ~~\tau\approx (2-4) {\rm ns},
\end{equation}
and the combination $(a\xi)$ which is the phase entering (\ref{eq:bold E(t)_temp}) is indeed $\gamma$-independent for $\theta\ll \gamma^{-1}$.

The same formula (\ref{eq:bold E(t)_temp}) also shows that the absolute value of the field $\mathbf{E}(t)$ and time duration is mostly determined by the region of the largest values of $\xi$ close to $\xi_{\rm max}$, while the low energy portion of the spectrum does not play a role.   We explicitly checked this feature by  plotting on Fig. \ref{fig:E_t} the electric field  $\mathbf{E}(t)$ with two different models for the filter. Solid line includes both filters spanning the low   and high frequency modes, while the dashed line corresponds to a single filter describing exclusively  high frequency modes. The difference between the two curves is negligible as claimed. Another property worth to be mentioned is that the absolute value of the field  $\mathbf{E}(t)$ in the peak increases with extending the upper value for $\omega_{\rm max}$ to higher value (800 MHz on Fig. \ref{fig:E_t} versus 600 MHz on Fig  \ref{fig:E_t3}). This is also expected behaviour as the integral  (\ref{eq:bold E(t)_temp}) is saturated
by the region of the largest values of $\xi$ close to $\xi_{\rm max}$, as already mentioned.
   \begin{figure}[h]
	\centering
	\captionsetup{justification=raggedright}
	\includegraphics[width=1\linewidth]{P_omega}
	\caption{ Spectrum of power density from Eq. \eqref{eq:d2 P d omega dA}, with $\tau=4{\rm\,ns}$, $\theta=0$, ${\cal R}=35\,$km, and $N=5\cdot10^8$. The Lorentz factor is chosen to be $\gamma=10$ (blue), 20 (orange), and 60 (green).}.
	\label{fig:P_omega}
\end{figure}

\exclude{
\begin{figure}[h]
	\centering
	\captionsetup{justification=raggedright}
	\includegraphics[width=1\linewidth]{E_t1}
		\caption{\Red{Electric field with $\theta=0$, $R=35\,$km, and $N=5\cdot 10^8$.}}
	\label{fig:E_t1}
\end{figure}

\exclude{
\begin{figure}[h]
	\centering
	\captionsetup{justification=raggedright}
	\includegraphics[width=1\linewidth]{E_t2}
		\caption{\Red{Blue: $\gamma=10$, Orange: $\gamma=20$;
		Dashed: (200-800)MHz,
		Solid:(40-80)MHz + (200-800)MHz. (Dashed).}}
	\label{fig:E_t2}
\end{figure}
}

\begin{figure}[h]
	\centering
	\captionsetup{justification=raggedright}
	\includegraphics[width=1\linewidth]{E_t4}
		\caption{\Red{Blue: $\gamma=10$, Orange: $\gamma=20$;
		Dashed: (200-1500)MHz,
		Solid:(40-80)MHz + (200-1500)MHz.}}
	\label{fig:E_t4}
\end{figure}
}

\exclude{
\Red{Explanation as follows:

First, the effective phase condition is $\nu_{\rm c}\xi t\approx\frac{1}{2}$. Then, by definition of $\xi$ and $\nu_c$ [Eqs. \eqref{eq:xi} and \eqref{eq:nu_c}], we have $\nu_c\xi=\frac{\omega}{4\pi}=\frac{\nu}{2}$ by setting $\theta=0$, so no $\gamma$-dependent. The effective time is then proportional to $t\approx \nu^{-1}$. If the filter matches the dominant range of the frequency spectrum, $t$ is solely determined by the upper frequency cutoff of the filter, which is $\gamma$-insensitive. To support this argument, we have an additional plot Fig. \ref{fig:E_t4} with a (200-1500) MHz filter. Obviously the time width is determined by the upper frequency cutoff of the filter. }

\Red{I attempted to write up the above argument formally below Eq. \eqref{eq:bold E(t)}. Also, I would suggest to remove the last line (integral in terms of $\xi$) and express the integral in terms of $\omega$. It is a lot clearer.  }}

The AQN-induced signal is in many aspects similar to the conventional CR shower,
with crucial and dramatic  difference  being the {\it non-inverted polarity}. It well matches the observation of the anomalous events in ANITA experiment.  We list a number of distinct  features (between AQN-induced events and conventional  CR air showers events) in concluding Sec.  \ref{conclusion}.

For completeness, we also derived the spectrum of power emission similar to the estimate in Ref. \cite{Huege:2003up}. The power is related to the electric field by Poynting vector:
\begin{equation}
\mathbf{S}(t)
=\frac{c}{4\pi}\mathbf{E}(t)\times\mathbf{B}(t)
\end{equation}
where we use the Gaussian units. The power density is given by
\begin{equation}
\frac{\rmd P(t)}{\rmd A}
=|\mathbf{S}(t)|
=\frac{c}{4\pi}|\mathbf{E}(t)|^2\,.
\end{equation}
Averaging over the time duration of the pulse $\tau$, it gives
\begin{equation}
\begin{aligned}
\left<\frac{\rmd P(t)}{\rmd A}\right>_\tau
=\frac{c}{4\pi\tau}\int |\mathbf{E}(t)|^2\rmd t
\approx\frac{c}{4\pi\tau}\int|\mathbf{E}(\omega)|^2\rmd\omega\,.
\end{aligned}
\end{equation}
where the last step follows from Parseval’s theorem.   The spectrum of power density is therefore
\begin{equation}
\label{eq:d2 P d omega dA}
\frac{\rmd^2P}{\rmd \omega\rmd A}
\approx\frac{\rmd}{\rmd\omega}\left<\frac{\rmd P(t)}{\rmd A}\right>_\tau
\approx\frac{c}{4\pi\tau}|\mathbf{E(\omega)}|^2\,.
\end{equation}

The spectrum for the power of the emission shares similar properties as the one of electric field as it is expressed in terms of the electric field $|\mathbf{E(\omega)}|$.  The results for the power density is presented on Fig. \ref{fig:P_omega}, where  we use $\tau\approx 4$\,ns for numerical estimates  taken from (\ref{eq:tau}), Fig  \ref{fig:E_t3}  and  Fig. \ref{fig:E_t}. This value     is also consistent with ANITA observations shown on Fig. 2 of Ref. \cite{Gorham:2018ydl}. The power density      of order $(0.2-0.3)\rm\,pW\,m^{-2}\,MHz^{-1}$ for $\gamma=20$ and frequencies $\nu\gtrsim 100$ MHz, then it falls sharply beyond the critical frequency \eqref{eq:nu_c}. This behaviour is consistent with ANITA's results presented on Fig. 4 of Ref. \cite{Gorham:2018ydl}.
\exclude{
\begin{figure}[h]
	\centering
	\captionsetup{justification=raggedright}
	\includegraphics[width=1\linewidth]{P_omega}
	\caption{Spectrum of power density from Eq. \eqref{eq:d2 P d omega dA}, with $\tau=1{\rm\,ns}$, $\theta=0$, $R=35\,$km, and $N=10^8$. The Lorentz factor is chosen to be $\gamma=10$ (blue), 20 (orange), and 60 (green).}.

	\end{figure}
	}

We finish   this section with few comments on accuracy of our estimates.  It should be emphasized that the computations carried out    in this section are  oversimplified as realistic signals can be severely modified by numerous factors such as geometry of the beam, relative position of the observer, and frequency characteristic of the receiver's filters  \cite{Huege:2003up,HUEGE2005116}. Nonetheless most factors, such as lateral structure and inclined axis of the beam, are near-field effects and they are eliminated since the observation angle is small. There are many other factors, such as the electron energy distribution of the beam as a function of $\gamma$, may  modify our predictions. We cannot  predict the corresponding behaviour as it is determined by very complex non-equilibrium dynamics as discussed in Appendix \ref{details}. However, the basic features of the AQN-induced radio  emission remain the same. In other words, the  qualitative picture presented in this work cannot be dramatically modified as most AQN parameters used in the present estimates had been fixed by previous studies of numerous different and unrelated systems as reviewed in Sec.    \ref{AQN}.

\exclude{
The most important far-field effect in our case is the longitudinal distribution of the AQN-induced electrons. For an electron beam with length $\sim c\tau$, the spectrum of the electric field involves an additional attenuation factor $S(\omega)$ \cite{Huege:2003up}:
\begin{equation}
S(\omega)
\sim\frac{1}{\sqrt{1+\omega^2\tau^2}}
\end{equation}
where we assumed the the beam follows an asymptotic $\Gamma$-distribution similar to conventional CR showers. \Red{I attached Fig. \ref{fig:E_omega_All} [including $S(\omega)$] in comparison to Fig. \ref{fig:E_omega}, which one is better?} Evidently this attenuation does not affect our conclusion as it only applies to moderate modification to the higher frequency regime $\nu\gtrsim10^2{\rm\,MHz}$.

\begin{figure}[h]
	\centering
	\captionsetup{justification=raggedright}
	\includegraphics[width=1\linewidth]{E_omega_All}
	\caption{\Red{$\tau=1{\rm\,ns}$}}.
	\label{fig:E_omega_All}
\end{figure}
}

\section{Conclusion and future development}\label{conclusion}
Our basic results can be summarized as follows.
We argued that the two anomalous events    observed by  the ANITA collaboration  \cite{Gorham:2016zah,Gorham:2018ydl,Gorham:2020zne}    can be interpreted as the  AQN-induced radio pulses which  result from    the AQN traversing the Earth and going in the upward direction.
The basic qualitative characteristics (such as the frequency of emission, the strength of the electric field and the duration of the radio pulses as presented in previous section \ref{sec:Radio signals induced by AQN})  of the observed AAE are consistent with our proposal and this  interpretation.

We also want to compare the results presented in this work with  the  studies of the conventional CR-induced radio emission \cite{Huege:2003up,HUEGE2005116}. This comparison shows that the AQN-induced radio pulses can be easily discriminated from conventional sources, such as CR air showers.

But first, we start with similarities between CR-induced and the AQN-induced radio pulses.
A ``rule of thumb" suggests that the maximal number of charged particles (mostly electrons and positrons) in CR air shower is $E_{\rm CR}/{\rm GeV}$, which implies that $N\approx (10^8-10^9)$ for a $E_{\rm CR}\approx (10^{17}-10^{18})$\,eV shower, see e.g. \cite{Huege:2003up}. This number of electrons from CR shower is close  to the number of electrons being emitted by the AQN according to (\ref{N}). Furthermore, the typical average energy in CR shower is 30 MeV which is also very similar to our estimates for the AQN-induced spectrum of the electrons with $E\in (1-10^2)$ MeV with the peak around 10 MeV according to (\ref{N}).
Therefore, it should not be a surprise that the radio emission intensity and the strength of   electric field is very similar in both cases as the geomagnetic field ${\cal B}\approx0.5{\rm\,gauss}$ (which represents the source of the acceleration and consequent   radio emission)   is obviously the same in the same location.

Now we want to discuss the drastic differences between the pulses induced by conventional CR showers and the AQNs. These dramatic distinct features can be tested in future experiments, such that our proposal  can be discriminated from any other suggestions. We list below   the following typical spectral features of the CR-induced radio pulses and contrast them with the AQN induced radio pulses:

1. The generic spectral feature of the CR-induced radio emission is the presence of oscillations which normally start around 100 MHz (depending on the distance from the shower axis), see e.g. Fig. 1 in  \cite{HUEGE2005116}. These oscillations are due to the coherence diminishing as the wavelength becomes shorter (in comparison to the ``pancake" size in CR shower).  While it is obviously affected by the detector's filter, this feature is a physical effect due to  changing  number of coherent particles with different wavelengths.
Such   picture being typical for the CR-induced radio emission  is not expected to occur for the AQN-induced radio signal as the notion of a ``pancake" does not exist in our case, see also item 3 below with an argument that the notion of a ``central axis"  also does not exist for the AQN-induced events.

2. Another typical feature of the CR-induced radio emission is that the most of the power is emitted at frequencies around $(20-30)$ MHz for $E_{\rm CR}\approx 10^{17} $\,eV shower, see Fig. 2 in  \cite{HUEGE2005116}. It is a  result of very strong cutoff frequency $\nu_0\lesssim 50$ MHz which strongly depends  on features of the shower, see Eq. (2) and (12) in  \cite{HUEGE2005116}.
It should be contrasted with our case when the cutoff frequency $\nu_{\rm c}\sim 0.7$ GHz  is determined by dramatically different physics as expressed by (\ref{eq:nu_c}).

3. Final and the most important difference between these two cases is that cutoff frequency $\nu_0 $
in CR air showers strongly depends on many parameters of the shower such as distance from the central axis when number of particles per unit area strongly depends on this parameter. It must be contrasted with our case of the AQN-induced radio signal when all electrons emitted from the same point at the same instant  are moving along the same direction.

This picture suggests that the event  could be viewed as a uniform front of size $\rho\sim 0.8$ km defined by (\ref{eq:rho})  rather than a CR air shower with a well defined central axis.  In different words, the number of particles per unit area in the  bunch of electrons does not depend on the distance from the central axis, in huge contrast  with conventional CR air showers. The notion ``central axis" simply does not exist  for the AQN induced electrons as the number of particles per unit area is approximately the same for all electrons generating the radio pulse.

  These  features of the AQN-induced radio events
  are very distinct from conventional CR induced radio pulses, and it should be easily discriminated by future analysis with more quality data. It can be achieved, for example, by placing two or more independent but synchronized antennas at a distance to study the
  same events from different locations. If future studies indeed support and substantiate    our proposal it would be a strong argument supporting the AQN nature of the AAEs.

  We conclude with few comments on the possibility to test the AQN model with other instruments (in addition to  the TA experiment mentioned previously in Section \ref{earth}) such as the Pierre Auger Observatory and several projects of the Joint Experiment Missions for Extreme Universe Space Observatory (JEM-EUSO), including  the currently operating Mini-EUSO detector, the planned EUSO on a Super Pressure Balloon II Mission (EUSO-SPB2), and the future Probe Of Extreme Multi-Messenger Astrophysics (POEMMA). As advocated recently, these instruments can potentially detect the UV photons from direct interaction between atmospheric molecules and macroscopic DMs (such as the AQN) \cite{Sidhu:2018auv,Anchordoqui:2021xhu} or ultrahigh energy DMs \cite{Xu:2019sgg}.

  Our original comments related to these studies are as follow. First, detection of macroscopically large quark nuggets, such as the AQNs,   would require proper adjustment of bin time of detectors as argued in  Ref. \cite{Sidhu:2018auv}. Furthermore, the AQN  mostly emits photons in X ray bands, such that a signal cannot  be observed by a fluorescence detector
  which is designed to detect the visible  and UV  light\footnote{In fact, it was precisely the argument presented in \cite{Budker:2020mqk} why the all sky camera  has not observed any signal from a powerful meteor-like event detected by the infra-sound  dedicated instrument.}. Finally, the X ray network suggested in \cite{Liang:2020mnz} is capable to detect the X rays from the AQNs, and even determine the directionality and velocity distribution of the DM. It gives a real opportunity to test   the Standard Halo Model (SHM) locally, which is known may dramatically deviate from conventional picture.
  This is because the standard  parameters of the SHM  are  properly fixed only on the global (not local) scales, see Introduction in   \cite{Liang:2020mnz} for references and the details.    We conclude with this optimistic note.

     \section*{Acknowledgements}

This research was supported in part by the Natural Sciences and Engineering
Research Council of Canada, and X.L. also by the UBC four year doctoral fellowship.