# Multi-Modal Clustering Events observed by Horizon-10T and Axion Quark
Nuggets

Ariel Zhitnitsky  Department of Physics and Astronomy, University of British
Columbia, Vancouver, V6T 1Z1, BC, Canada

###### Abstract

The Horizon-10T collaboration [1, 2, 3, 4, 5] have reported observation of
Multi-Modal Events (MME) containing multiple peaks suggesting their clustering
origin. These events are proven to be hard to explain in terms of conventional
cosmic rays (CR). We propose that these MMEs might be result of the dark
matter annihilation events within the so-called axion quark nugget (AQN) dark
matter model, which was originally invented for completely different purpose
to explain the observed similarity between the dark and the visible components
in the Universe, i.e. $\Omega_{\rm DM}\sim\Omega_{\rm visible}$ without any
fitting parameters. We support this proposal by demonstrating that the
observations [1, 2, 3, 4, 5], including the frequency of appearance,
intensity, the spatial distribution, the time duration, the clustering
features, and many other properties nicely match the emission characteristics
of the AQN annihilation events in atmosphere. We list a number of features of
the AQN events which are very distinct from conventional CR air showers. The
observation (non-observation) of these features may substantiate (refute) our
proposal.

##  I Introduction

In this work we discuss two naively unrelated stories. The first one is the
study of a specific dark matter (DM) model, the so-called axion quark nugget
(AQN) model [6], see a brief overview of this model below. The second one
deals with the recent puzzling observations [1, 2, 4, 3, 5] by the Horizon-10T
(H10T) collaboration of the Multi-Modal Events (MME). We overview the
corresponding events below in details. We also highlight some difficulties in
interpretation of these events in terms of the standard CR air showers. The
unusual features of MME include:

1. “clustering puzzle”: Two or more peaks separated by $\sim 10^{2}$ ns are present in several detection points, while entire event may last $\sim 10^{3}$ ns. It can be viewed as many fronts separated by $\sim(10^{2}-10^{3})$ ns, instead of a single front;

2. “particle density puzzle”: The number density of particles recorded at different detection points apparently weakly dependent on distance from Extensive Air Showers (EAS) axis;

3. “pulse width puzzle”: The width of each individual pulse is around $(20-35)$ ns and apparently does not depend on distance from EAS axis;

4. “intensity puzzle”: The observed intensity of the events (measured in units of a number of particles per unit area) is of order $\rho\sim(100-300)\rm m^{-2}$ when measured at distances $(300-800)$ m from EAS axis. Such intensity would correspond to the CR energy of the primary particle on the level $E_{p}\gtrsim 10^{19}$ eV which would have dramatically different event rate in comparison with observed rate recorded by H10T.

Before we proceed with our explanation of the proposal to view the MMEs as the
AQN events one should mention that similar unusual features of the CR air
showers have been noticed long ago for the first time by Jelley and Whitehouse
[7] in 1953. Later, EAS exhibiting the unusual time structures were studied by
several independent experiments, see e.g. [8]. We refer to ref. [2] for
overview and references of the older observations and studies of such unusual
events.

The considerable recent advances in rising the resolution (on the level of a
few ns) has allowed the H10T collaboration dramatically improve the collection
and analysis of the MMEs. In particular, during $\sim 3500$ hours of operation
the H10T collaboration collected more than $10^{3}$ MMEs.

In what follows we overview the observation [1, 2, 3, 4, 5] by emphasizing the
puzzling nature of these events if interpreted as conventional EAS events. At
the same time the same observations can be explained in very natural way if
interpreted in terms of the AQN annihilation events as we shall argue in this
work. One should mention here that a similar conclusion has been also reached
for a different type of unusual CR-like events. First, it has been argued in
[9, 10] that the Telescope Array (TA) “mysterious bursts” [11, 12] can be
naturally interpreted as the AQN events. Secondly, it has been also argued in
[13] that the Antarctic Impulse Transient Antenna (ANITA) observation [14, 15,
16] of two anomalous events with noninverted polarity might be also a
consequence of the AQN annihilation events. Important comment here that in all
those cases the basic parameters of the AQN model remain the same as they have
been fixed long ago from dramatically different observations in a very
different context.

Our presentation is organized as follows. In next Section II we explain why
the observation [1, 2, 3, 4, 5] listed as 1-4 puzzles are very mysterious
events if interpreted as conventional EAS events. In Section III we give a
brief overview of the AQN framework with emphasize on the key elements
relevant for the present studies. In Section IV we formulate our proposal on
identification of the unusual Multi-Modal Events with the upward moving AQN
events. In Section V we estimate the event rate. Our main Section is VI where
we estimate a variety of relevant time scales and the intensity of the events.
In the same section we also confront our proposal with observations and argue
that puzzles 1-4 can be naturally understood within the AQN framework. Section
VII is our conclusion where we suggest several tests which may substantiate or
refute our proposal on identification of the Multi-Modal Cluster Events with
the upward moving AQNs.

##  II Conventional CR picture confronts the MME observations [1, 2, 3, 4]

Figure 1: top: Aerial view and geometry of H10T instrument with location of 10
detectors, adopted from [3]; bottom: A typical MME event recorded on April
6.2018 by H10T instrument at point $\\#9$, adopted from [3]. All pulses are
recorded at a single detection point. Delay times, the width of each peak
$\tau_{0.8}$ in ns, and the particle density $\rho_{0.8}$ per $\rm m^{-2}$
within $\tau_{0.8}$ are also shown in the table.

.

Figure 2:  Solid lines: the particle density distribution $\rho(R)$ in
simulated EAS disk versus distance from axis for different energies shown by
different colours, depending on energy of the CR. Blue dots: cumulative
particle density for each bimodal pulse vs. distance to EAS axis, adopted from
[3].

.

Figure 3:  Solid lines: simulated EAS disk width versus distance from axis for
different energies shown by different colours. Pulse width $\tau_{0.8}$ in ns
for the first pulse (in blue) and second pulse (in red) for the bimodal pulses
versus distance to the EAS axis, adopted from [3].

.

This section is devoted to the first part of our story where we describe the
MME observations and argue why the observed events are inconsistent with
conventional CR interpretation, while the next section III is devoted to the
second part of our story, the AQN dark matter model.

We start by reviewing the conventional picture of the EAS. It is normally
assumed that EAS can be thought as a disk -pancake with well defined EAS axis.
It is also assumed that EAS represents an uniform, without any breaks
structure. The standard picture also assumes that the particle density drops
smoothly with distance when moving away from the core, while the thickness of
the EAS pancake increases with the distance from the core. Now we want to see
why this conventional picture is in dramatic contradiction with observations
of the MME events.

1. Indeed, a typical MME is shown on Fig. 1 where a complicated temporal features are explicitly seen. Several peaks separated by $\sim 10^{2}$ ns in a single detector represent the “clustering puzzle”, listed above. In conventional EAS picture one should see a single pulse in each given detector with the amplitude which depends on the distance from the EAS axis. It is not what actually observed by H10T.

2. The manifestation of the “particle density puzzle” is as follows. On Fig. 2 we show the particle density distribution $\rho(R)$ in simulated EAS disk versus distance from axis for different energies shown by solid lines with different colours, depending on energy of the CR. In particular, for energy of the primary particle on the level of $10^{17}$ eV one should expect a strong suppression $\sim 10^{3}$ when distance $R$ to the EAS axis changes from $R\simeq 100$ m to $R\simeq 700$ m. It is not what actually observed by H10T: the density of the particles $\rho(R)$ is not very sensitive to the distance to the EAS axis, and remains essentially flat for the entire region of observations. Furthermore, the magnitude of the density $\rho(R)$ is much higher than it is normally expected for CR energies $\sim(10^{17}-10^{18})$ eV.

3. The manifestation of the “pulse width puzzle” is as follows. On Fig. 3 we show the simulated EAS disk width versus distance from axis for different energies shown by solid lines in different colours. As explained above the thickness of the EAS pancake increases with the distance from the core. Therefore, the pulse width must also increase correspondingly as shown by sold lines on Fig. 3 for different energies. It is not what actually observed by H10T: all MME events show similar duration of the pulse width on the level of $(20-35)$ ns irrespective to the distance to the AES axis. This observation is in dramatic conflict with conventional picture as outlined above.

4. The manifestation of the “intensity puzzle” is as follows. The Fig. 2 suggests that the charged particle density $\rho(R)$ varies between $(30-300)$ particles per $m^{2}$ at the distances $(500-700)$ m from the EAS axis. This is at least factor $10^{2}$ above the expected $\rho(R)$ for the primary particle with energy in the interval $(10^{17}-10^{18})$ eV. Only the particles with energies well above $10^{19}$ eV could generate such enormous particle density as shown on Fig. 2. However the frequency of appearance for such highly energetic particles is only once every few years. This represents the “intensity puzzle” when the intensity of the events estimated from $\rho(R)$ for MME is several orders of magnitude higher than the energy estimated by the event rate. This shows a dramatic inconsistency between the measured intensity and observed event rate carried out by one and the same detector.

We conclude this section with the following comment. From the puzzles listed
above it must be obvious that the events detected by H10T are not the
conventional CR air showers as they demonstrate enormous inconsistency with
standard CR interpretation. What it could be? Before we put forward our
proposal we would like to briefly overview in next section III the second part
of our story- the AQN dark matter model.

##  III The AQN Dark matter model

We start with few historical remarks and motivation of the AQN model in
subsection III.1, while in subsection III.2 we overview recent observations of
CR-like events (such as puzzling bursts observed by the TA experiment and
ANITA’s two anomalous events with non-inverted polarity) of some mysterious
events which could be explained by the AQN events hitting the Earth. Finally,
in section III.3 we briefly 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
MMEs as the AQN events.

###  III.1 The basics

The AQN dark matter model [6] was invented long ago with a single motivation
to explain in natural way the observed similarity between the dark matter and
the visible densities in the Universe, i.e. $\Omega_{\rm DM}\sim\Omega_{\rm
visible}$ without any fitting parameters. We refer to recent brief review [17]
on the AQN model. Here we want to mention few key elements which are important
for this work.

The AQN construction in many respects is similar to the Witten’s quark
nuggets, see [18, 19, 20], and review [21]. 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 [18, 19, 20, 21]. First new element is the presence of the axion domain
walls which are copiously produced during the QCD transition111The 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
[22, 23, 24, 25, 26, 27, 28], and recent reviews [29, 30, 31, 32, 33, 34, 35,
36, 37].. The axion field $\theta(x)$ plays a dual role in this framework:
first it serves as an additional stabilization factor for the nuggets, which
helps to alleviate a number of problems with the original nugget construction
[18, 19, 20, 21]. 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 matter as well as 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 antimatter over another species made of matter. The
preference is determined by the initial sign of the $\theta$ field when the
formation of the AQN nuggets 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 as they both proportional
to one and the same dimensional parameter $\Lambda_{\rm QCD}$. We refer to the
original papers [38, 39, 40, 41] 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, see also brief review [17] covering these topics.

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. The strongest direct detection limit222Non-detection of etching
tracks in ancient mica gives another indirect constraint on the flux of dark
matter nuggets with mass $M<55$g [42]. 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 [43]:

$\displaystyle\langle B\rangle>3\cdot 10^{24}~{}~{}~{}[{\rm
direct~{}(non)detection~{}constraint]},$ (1)

where we formulated the constraints in term of the AQN’s baryon charge $B$
rather than in terms of its mass $M\approx m_{p}B$. Similar limits are also
obtainable from the ANITA and from geothermal constraints which are also
consistent with (1) as estimated in [44].

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
dark matter 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 CR are much better suited for our studies of the
AQNs as we discuss in next section III.2.

###  III.2 When the AQNs hit the 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 deposition 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. In this short subsection we want to mention
several observed phenomena which could be related to the AQN annihilation
events when the nuggets propagate in the Earth’s atmosphere.

We start with the Telescope Array (TA) experiment which has recorded [11, 12]
several short time bursts of air shower like events. These events are very
unusual and cannot be interpreted in terms of conventional 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 [11, 12]. It has been argued in [9, 10] that
these bursts represent the direct manifestation of the AQN annihilation
events.

This feature, in many aspects, is very similar to the intensity puzzle listed
in previous Section II, where analogous dramatic inconsistency (between the
observed intensity and measured event rate) is recorded by H10T instrument.

Our next example is the observed seasonal variations of the X ray background
in the near-Earth environment. To be more precise, the XMM-Newton at
$11\sigma$ confidence level [45] had recorded the seasonal variations in the
2-6 keV energy range at distances $r\gtrsim 8R_{\oplus}$ where the
measurements are effectively had been performed. The authors [45] argue that
conventional astrophysical sources have been ruled out. Furthermore, the XMM-
Newton’s operations exclude pointing at the Sun and at the Earth directly,
diminishing a possible direct X ray background from the Earth and the Sun. It
has been argued in [46] that this seasonal variation may be naturally
explained by the AQNs exiting the Earth. The AQNs continue the emission of the
X rays at $r\gtrsim 8R_{\oplus}$, long after the nuggets exit the Earth’s
atmosphere. It has been also shown that the spectrum and the intensity
computed in [46] almost identically match the observed spectrum [45].

The final example we want to mention here is the observation by ANITA [14, 15,
16] of two anomalous events with noninverted polarity. These two events had
been identified as the Earth emergent upward going CR-like events with exit
angles of $-27\degree$ and $-35\degree$ relative to the horizon. These
anomalous events are in dramatic tension with the standard model because
neutrinos are exceedingly unlikely to traverse through Earth at a distance of
$\gtrsim 5\times 10^{3}\,$km with such ultrahigh energy, even accounting for
the $\nu_{\tau}$ regeneration [14]. The analysis [47] 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
level. We advocated in [13] that these two anomalous events with noninverted
polarity observed by ANITA could be interpreted as the upward going AQNs.

The applications [46] to the X ray emission in the near-Earth environment and
the ANITA anomalous events [13] are especially relevant for the present work
because both these phenomena are originated from the AQN upward going (Earth
emergent) events when the AQNs traversed through the Earth interior and exit
the Earth surface. As we shall argue in this work the MME events could be also
the consequence of the same upward moving AQNs.

However, there is an additional technical challenging problem to be addressed
for the application to MME analysis. The point is that the ANITA anomalous
events are generated by AQNs which just crossed the surface and entered the
Earth’s atmosphere (upward-going events). Exactly at this instant large number
of electrons instantaneously released into atmosphere and generated the radio
pulse measured by ANITA. In contrast, the X rays are observed by XMM-Newton at
very large distances $r\gtrsim 8R_{\oplus}$ from Earth. In this case the
impact of the atmospheric material at $r\lesssim 60$ km can be ignored, and
the cooling of the AQN is dominated by nugget’s propagation in an empty space.
The application to MME, which is the topic of the present work, deals with an
intermediate stage between these two cases: it is not the first instant when
the AQN enters the atmosphere emerging from the interior, and it is not the
asymptotically far away when the AQN propagates in empty space. In other
words, we need to know the emission pattern and the energy deposition rate for
the AQNs when they propagate in the atmosphere, and the impact of the
atmospheric surrounding material cannot be ignored. We formulate our proposal
in Section IV where we argue that emission of the electrons by the AQN in form
of the well-separated “bunches” may explain the unusual features of the MME
observations.

However, before we present our arguments supporting this interpretation of the
MME unusual events (as highlighted in section II) we should overview the basic
characteristics of the AQNs traversing the Earth, which represents the topic
of the next subsection III.3.

###  III.3 Internal structure of the AQN

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 [48] 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}}$ and atmosphere with
$n_{\rm air}\sim 10^{21}~{}{\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 [48] and
it is given by

$F_{\text{tot}}\approx\frac{16}{3}\frac{T^{4}\alpha^{5/2}}{\pi}\sqrt[4]{\frac{T}{m}}\,,$
(2)

where $\alpha\approx 1/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 black body radiation, see [48],
see also Appendix A 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. (2) must balance the flux of energy onto the nugget

$\displaystyle F_{\text{tot}}(4\pi R^{2})\approx\kappa\cdot(\pi
R^{2})\cdot(2~{}{\rm GeV})\cdot n\cdot v_{\rm AQN},$ (3)

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 (2) 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_{p}c^{2}\approx 2~{}{\rm GeV}$ energy. In Eq. (3) 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, neutrinos333In a neutral dilute galactic environment
considered previously [48] the value of $\kappa$ was estimated as
$\kappa\approx 0.1$.. 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 [49] 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 (3). 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
[46] explaining the seasonal variations of the X rays observed the XMM-Newton
at $11\sigma$ confidence level [45] at distances $r\sim(6-10)R_{\oplus}$ from
the Earth surface.

The crucial element for the study [46] was the emission rate by AQNs at very
high temperatures $T\approx(100-200)$ keV. In this case the emission is not
determined by simple Bremsstrahlung radiation given by (2) which is only valid
for relatively low temperatures. In the high temperature regime a number of
many-body effects in the electrosphere, that were previously ignored, become
important. In particular, it includes: the generation of the plasma frequency
in electrosphere, which suppresses the low frequency emission. Also, the the
AQNs get ionized which dramatically decreases the density of positrons in
electrosphere, and consequently suppresses the emission, among many others
effects. All these complications have been carefully considered in [46] in
applications to the seasonal variations of the X rays as observed the XMM-
Newton, and we refer to that paper for the details. Here we quote the result
of these studies. The rate of emission can be effectively represented as
follows:

$\displaystyle\frac{dE_{\rm emiss}}{dt}\sim(4\pi
R^{2})\eta(T,R)F_{\text{tot}}(T)$ (4)

where factor $\eta(T,R)$ is a result of strong ionization of the
electrosphere, which leads to the corresponding suppression of the emission.
The details of this suppression are explained in Appendix A and given by (36).
This suppression is a direct consequence of high internal temperature $T$ when
a large number of weakly bound positrons are expanded over much larger
distances order of $R$ rather than distributed over much shorter distances of
order $m^{-1}$ around the nugget’s core. This basically determines the
suppression factor $\eta\sim(mR)^{-1}\sim 10^{-6}$.

In both previously considered applications (to ANITA [13] and to XMM-Newton
[46] experiments) when the upward moving AQNs with high $T\approx(100-200)$
keV play the key role in the explanations of the observations, we ignored an
additional444We use the term“additional” to emphasize that the AQN had
accumulated a huge amount of energy during the transpassing the Earth interior
as the capacity of the quark core nugget is very large [46]. Precisely this
accumulated energy will be emitted in form of the X rays when the AQN
propagates in empty space at distances $r\sim(6-10)R_{\oplus}$ from the Earth
surface. annihilation events with atmospheric material as we already mentioned
at the end of the previous subsection III.2. Now we want to include the
corresponding physics in our analysis.

First of all we want to estimate the number of direct head-on collisions of
the atmospheric molecules with AQN per unit time. It can be estimated as
follows:

$\displaystyle\frac{dN}{dt}\simeq(\pi R^{2})N_{\rm m}v_{\rm AQN}\simeq
0.8\cdot 10^{18}\left(\frac{n_{\rm air}}{10^{21}~{}{\rm cm^{-3}}}\right)\rm
s^{-1},~{}~{}~{}$ (5)

where $N_{m}\simeq 2.7\cdot 10^{19}~{}{\rm cm^{-3}}$ is the molecular density
in atmosphere when each molecule contains approximately 30 baryons such that
the typical density of surrounding baryons in air is $n_{\rm air}\simeq
30\cdot N_{m}\simeq 10^{21}~{}{\rm cm^{-3}}$. The dominant portion of these
collisions are the elastic scattering processes rather than successful
annihilation events suppressed by parameters $\kappa$ as discussed above. The
energy being deposited to the AQN per unit time as a result of annihilating
processes can be estimated as follows

$\displaystyle\frac{dE_{\rm deposit}}{dt}\approx(2~{}{\rm
GeV})\cdot\kappa\cdot\frac{dN}{dt}$ $\displaystyle\approx$ $\displaystyle
1.6\cdot 10^{18}\cdot\kappa\cdot\left(\frac{n_{\rm air}}{10^{21}~{}{\rm
cm^{-3}}}\right)\rm\frac{GeV}{s},$

One should emphasize that this energy is being deposited into the AQN which is
already at full capacity with very high temperature $T\simeq(100-200)$ keV
accumulated during a long journey in much denser Earth’s interior.

What is the mechanism to release this extra energy? In other words: In what
form the energy (III.3) could be emitted into the surrounding atmosphere as it
cannot be easily transferred to the AQN’s quark core (as it is already
saturated, see footnote 4)?

There are two qualitatively different regimes which can be normally realized
in such circumstances. If a typical time scale $\tau_{deposit}$ to deposit a
specific amount of energy (let us say, 1 GeV) determined by (III.3) is longer
than the time scale $\tau_{\rm cool}$ to release the same amount of energy,
i.e. $\tau_{deposit}\gtrsim\tau_{\rm cool}$ than a continuous cooling process
takes place and the temperature slowly decreases. This is a thermodynamically
equilibrium behaviour which can be analyzed using conventional technical
tools555For example, the AQNs moving in empty space and slowly cooling by
emitting the X rays as computed in [46] belongs to this class..

The second option is realized when the time scale to deposit the energy
$\tau_{deposit}$ is much shorter than the time scale of cooling processes,
i.e. $\tau_{deposit}\ll\tau_{\rm cool}$, in which case the burst-like
(explosion like, eruption-like blasts) non- equilibrium processes must occur.
These eruption-like events will release the extra accumulated energy in form
of very short pulses, which cannot be described as conventional cooling
equilibrium processes. These short pulses must alternate with much longer
periods of accumulation energy which inevitably end with subsequent eruption-
like events. In other words, in this case, the non-equilibrium fast
equilibration process can be thought as the eruption-like event666An analogy
for such eruption-like event is the lightning flash under thunderclouds. The
clouds accumulate the electric charge in form of the ionized molecule very
efficiently. The corresponding time scales plays the role of $\tau_{deposit}$
in our system. The neutralization of these ions is less efficient process,
which is analogous to our $\tau_{\rm cool}$. If some conditions are met (the
so-called runaway breakdown conditions are satisfied, see [50, 51] for review)
the discharge occurs in form of the eruption which is the lightning strike in
our analogy. The system is getting neutralized in form of the non- equilibrium
lightning event (eruption) on the time scales which are much shorter than any
time scales of the problem. This analogy is in fact, quite deep, and will be
used in section VI when we discuss different time scales of MMEs as observed
by H10T detector..

Now we give a numerical estimation for the cooling rate determined by (4) to
find out what regime is realized in given circumstances. The order of
magnitude estimate can be expressed as follows:

$\displaystyle\frac{dE_{\rm emiss}}{dt}\sim 10^{16}\left(\frac{T}{100~{}\rm
keV}\right)^{\frac{17}{4}}\rm\frac{GeV}{s}.$ (7)

Comparing this estimate for cooling rate with estimate for the deposition rate
(III.3) one concludes that the rate of energy deposition is much shorter than
the rate at which the system is capable to release the same amount of energy
in form of the photon’s emission. As a result, we arrive to the conclusion
that second option when $\tau_{deposit}\ll\tau_{\rm cool}$ is realized. This
conclusion inevitably implies the eruption like events must occur, and these
very short explosions must be alternated with much longer periods of time when
the energy is being accumulated by the AQNs. These longer periods of the
energy accumulations also must end with consequent eruption-like blasts. This
alternating pattern continues as long as following condition holds:

$\displaystyle\tau_{deposit}\ll\tau_{\rm
cool}~{}~{}\Rightarrow~{}~{}{\rm(short~{}eruptions~{}occur)}.$ (8)

With these preliminary comments with estimates from the previous studies on
the AQN features we are now in position to formulate the proposal interpreting
the MME events observed by X10T collaboration in terms of the upward moving
AQNs, which is the topic of the next section IV.

##  IV MME as the AQN event

In this Section we formulate the basic idea of our proposal on identification
of the unusual Multi-Modal Events with the upward moving AQN events, while
supporting estimates on the event rate, the intensity, and the variety of time
scales will be presented in following sections V and VI.

The AQNs which propagate in the Earth’s interior are very hot. Their
temperature could be as hot as $T\approx(100-200)$ keV at the moment of exit
on the Earth’s surface as we discussed in [46] in the application to studies
of the seasonal variations of the X rays observed by XMM-Newton at very large
distances $r\sim(6-10)R_{\oplus}$ from the Earth. At such large distances the
AQN’s cooling can be described by conventional thermodynamically equilibrium
processes with well defined spectrum, see footnote 4. It should be contrasted
with our application to ANITA anomalous events [13] when the same upward
moving AQNs just crossed the surface and entered the Earth’s atmosphere. In
this case the dominant cooling mechanism is realized in form of short burst-
like events according to (8). The outcome of these short burst-like events is
emission of numerous bunches (clumps) of relativistic electrons. These
eruption like events are analogous to lightning flashes from footnote 6.

In our studies [13] we focused on an estimation of the number of emitted
electrons $N$ and their energy $\langle E\rangle\sim 10~{}{\rm MeV}$ which,
according to the proposal [13], essentially determine the intensity and the
spectral features of the ANITA anomalous events. These electrons will be
always accompanied by the emission of much more numerous number of photons
with typical energy of order $\omega\sim T\approx(100-200)$ keV, which in fact
represent the dominant fast non- equilibrium mechanism of cooling of the
system in these circumstances (8).

Precisely these multiple explosion-like emissions of photons and electrons
will be identified with Multi-Modal Events with their unusual features
reviewed in Section II. The mean free path of the photons with $\omega\sim T$
is relatively short, measured in meters. Therefore, they cannot be directly
observed. In contrast, the mean free path of the energetic electrons at sea
level with $\langle E\rangle\sim 10~{}{\rm MeV}$ is around several kilometres
in atmosphere (and even longer at the elevation of 3346 m where H10T is
located). The propagation of these ultra-relativistic electrons takes place in
the background of the geomagnetic field ${\cal{B}}\sim 0.5~{}\rm G$ which
determines the instantaneous curvature $\rho\sim 3$ km of the electron’s
trajectories with such energies, see section V with relevant estimates. Due to
the large component of magnetic field parallel to the Earth’s surface the
electrons may move in downward direction even if they are initially emitted in
upward direction. Our proposal is that precisely these downward moving
energetic electrons could mimic the CR and could be responsible for the Multi-
Modal pulses as detected by H10T instrument.

Any precise computation during a short period of time of emission is very hard
problem of non-equilibrium dynamics, which is beyond the scope of the present
work. Fortunately, the observable intensity, the event rate, the time delays
between the pulses are determined in the AQN framework by the basic
characteristics of the model and are not very sensitive to the details of this
non-equilibrium mechanism of production and its corresponding time scales.
This is because all these observables depend on several parameters such as
typical energy $\langle E\rangle\sim 10~{}{\rm MeV}$ of the emitted electrons,
typical number of emitted electrons $N$, the deposition rate (III.3) and the
cooling rate (7) which had been previously estimated for very different
purposes in different circumstances. We shall use the same parameters in the
present work.

In particular, the number of electrons $N\sim(10^{8}-10^{9})$ emitted by AQNs
is very hard to compute from the first principles. However it can be fixed by
the observed intensity and spectral features of the anomalous ANITA radio
pulses assuming, of course, that these observed radio pulses with inverted
polarities are due to the upward moving AQNs [13]. The typical energy $\langle
E\rangle\sim 10~{}{\rm MeV}$ of these electrons is also consistent with
duration of the observed pulses [13]. The next section V is devoted to the
estimation of the event rate of MME, while section VI is devoted to
explanations of the puzzles which demonstrate a dramatic deviation from
conventional CR picture as formulated in section II. We shall argue that the
observed MME features are consistent with our AQN-based interpretation.

##  V Event Rate of MME

This section is devoted to estimate of the event rate of MMEs within the AQN
framework. We anticipate that this evaluation is expected to be very
qualitative estimation due to large uncertainties in parameters and rare
occurrence of the observed MMEs. In particular, it is known that key
parameter, the DM density locally may dramatically deviate from the well
established average global value $\rho_{\rm DM}\approx
0.3{\rm\,GeV\,cm^{-3}}$. Nevertheless we would like to present such estimate
to demonstrate that our interpretation of MMEs as an outcome of upward-going
AQNs is at least a self-consistent proposal. In what follows we use the same
formulae we used for the estimates of the frequency of appearance of the ANITA
Anomalous Events [13] such that the numerous uncertainties related to the
local DM density or/and the AQN size distribution would not dramatically
affect our estimates below if they are normalized to the ANITA event rate.

We start with the same formula from [13] for the expected number ${\cal N}$ of
the MMEs assuming that they are induced by the AQNs:

${\cal N}\approx{\cal{A}}_{\rm eff}{\cal
T}\Delta\Omega\frac{\mathrm{d}\Phi}{\mathrm{d}A\mathrm{d}\Omega}\,,$ (9)

where $\Delta\Omega\approx 2\pi$ for the isotropic AQN flux, and the
expression for the local rate of upward-going AQNs per unit area is given by:

$\displaystyle\frac{\mathrm{d}\Phi}{\mathrm{d}A\mathrm{d}\Omega}\approx\frac{\Phi}{4\pi
R_{\oplus}^{2}}\approx 4\cdot 10^{-2}\left(\frac{10^{25}}{\langle
B\rangle}\right)\rm\frac{events}{yr\cdot km^{2}},$ (10)
$\displaystyle\Phi\approx\frac{2\cdot 10^{7}}{\rm yr}\left(\frac{\rho_{\rm
DM}}{0.3{\rm\,GeV\,cm^{-3}}}\right)\left(\frac{v_{\rm
AQN}}{220\,\mathrm{km}\,\mathrm{s}^{-1}}\right)\left(\frac{10^{25}}{\langle
B\rangle}\right),$

where $\rho_{\rm DM}$ is the local density of DM and $R_{\oplus}=6371\,$km is
the radius of the Earth and $\Phi$ is the total hit rate of AQNs on Earth
[43]. In formula (9) the ${\cal T}$ is the time of operation while
${\cal{A}}_{\rm eff}$ is the effective area to be estimated below.

The estimation of the effective area ${\cal{A}}_{\rm eff}$ is a complicated
task as it is not simply determined by the area of the detectors similar to
standard CR analysis. Instead, it is determined by the area along the AQN’s
path where it can emit the bunches of electrons in form of short pulses which
can mimic the CR as outlined in previous section IV. The area ${\cal{A}}_{\rm
eff}\sim L_{\rm AQN}\cdot D_{\rm AQN}$ can be thought as a strip of length
$L_{\rm AQN}$ and width $D_{\rm AQN}$.

To estimate these parameters let us consider an AQN moving in upward direction
with angle $\theta_{\rm AQN}$ with respect to the Earth’s surface such that
upward velocity component is $(v_{\rm AQN}\cdot\sin\theta_{\rm AQN})$. This
vertical velocity component determines the maximal height $h$ of the
atmosphere where the atmospheric density is still sufficiently high such that
condition (8) holds. We estimate $h\approx(30-50)$ km where the atmospheric
density falls by two order of magnitude such that conventional X-ray emission
(which had been used in computations [46] to explain the seasonal variations
observed by XMM Newton) becomes the dominant cooling process at higher
altitudes. This condition determines the maximal length $L_{\rm AQN}$ and time
scale $\tau_{\rm AQN}$ when the AQNs can emit the bunches of electrons in form
of pulses which we identify with MMEs and which can mimic the CR air showers
as outlined in previous section IV. Numerically these parameters are estimated
as follows:

$\displaystyle\tau_{\rm AQN}$ $\displaystyle\sim$
$\displaystyle\frac{h}{(v_{\rm AQN}\cdot\sin\theta_{\rm AQN})}\sim 0.2\rm s,$
(11) $\displaystyle L_{\rm AQN}$ $\displaystyle\sim$ $\displaystyle
h\cot\theta_{\rm AQN}\sim(50-90)\rm km,$

where we used $\theta_{\rm AQN}\simeq 30^{0}$ for an order of magnitude
numerical estimates.

To estimate the width of the strip $D_{\rm AQN}$ we recall that the
geomagnetic field ${\cal B}$ in location of the H10T instrument can be
characterized by two numbers: it has strong component in up to down direction
${\cal B}_{\rm down}\simeq 0.49$ G, and strong component from south to north
direction ${\cal B}_{\rm north}\simeq 0.24$ G, see e.g. [52].

Assuming that the energetic electrons with $\langle E\rangle\sim 10~{}{\rm
MeV}$ will be emitted along the ${\mathbf{v}}_{\rm AQN}$-direction one can
estimate instantaneous radius of curvature $\rho$, see e.g. Jackson [53]:

$\rho\approx\frac{\gamma mc}{e{\cal B}\sin\theta_{\cal B}}\approx
2.8{\rm\,km}\left(\frac{\gamma}{20}\right),~{}~{}~{}~{}~{}~{}\gamma\equiv\frac{\langle
E\rangle}{m}$ (12)

where ${\cal B}\approx{\cal B}_{\rm north}\approx 0.24$ G is the local
magnetic field strength, which changes the direction of the electrons from
upward to downward moving such that they can mimic the CR air showers. The
angle $\theta_{\cal B}$ is the angle between the particle electron velocity
$\mathbf{v}$ and magnetic direction. We choose $\theta_{\cal B}\approx
30\degree$ in (12) for the numerical estimates.

Now we estimate the relevant scale $\lambda$ which determines the survival
pattern of the electron’s bunches. It is mostly determined by the Coulomb
elastic scattering with cross section $\sigma_{\rm Coul}$

$\displaystyle\sigma_{\rm
Coul}\approx\frac{\alpha^{2}}{E^{2}\theta^{4}}\approx 3\cdot
10^{-27}\left(\frac{1/2}{\theta}\right)^{4}\left(\frac{20}{\gamma}\right)^{2}\rm
cm^{2}.$ (13)

The electrons with $\theta\gtrsim 1/2$ may strongly deviate from their main
paths and cease to stay with majority of particles forming the bunch (which
eventually becomes the pulse being interpreted as MME). The corresponding
length scale $\lambda(h)$ at altitude $h$ (accounting for the air-density
$n_{\rm air}(h)$ variation) is estimated as follows

$\displaystyle\lambda(h)\sim\frac{1}{\sigma_{\rm Coul}n_{\rm air}(h)}\sim
3\cdot\exp\left({\frac{h}{8~{}\rm km}}\right)~{}\rm km.$ (14)

The physical meaning of $\lambda$ is the length distance particles propagate
at which the majority of the particles in the bunch remain within the bunch
before the dispersing to much larger distances when they cease to be a part of
the pulse (burst). Important comment here is that $\lambda\gtrsim\rho$ such
that the majority of electrons within the bunch do not strongly re-scatter,
and therefore survive the change in orientation due to the geomagnetic field
${\cal B}$ (from upward to downward direction).

The width of the strip $D_{\rm AQN}$ now can be roughly estimated as follows:

$\displaystyle D_{\rm AQN}\sim\lambda\sim 3\cdot\exp\left({\frac{h}{8~{}\rm
km}}\right)\cdot\left(\frac{\gamma}{20}\right){\rm\,km}.$ (15)

Combining all the estimates above one arrives to the following order of
magnitude estimation for the expected number of MME events:

${\cal N}\approx 0.5\cdot 10^{2}\rm events\left(\frac{{\cal{A}}_{\rm
eff}}{500~{}\rm km^{2}}\right)\left(\frac{{\cal{T}}}{0.4~{}\rm yr}\right).$
(16)

This order of magnitude estimate should be compared with observed frequency of
appearance of the MMEs. According to [2] the observations from start until
August 2016 during ${\cal T}\approx\,3,500\rm h\approx 0.4\rm yr$ of operation
the H10T instrument had detected ${\cal N}_{\rm obs}\simeq 10^{3}$ MMEs.
According to [3] the observations from Feb 15, 2018 to May 12, 2018 (which
corresponds to ${\cal{T}}\approx{0.25~{}\rm yr}$) the H10T instrument had
detected ${\cal N}_{\rm obs}=217$ MMEs.

The AQN based estimations (16) suggest that number of events ${\cal N}$ is
almost one order of magnitude lower than ${\cal N}_{\rm obs}$ events observed
by H10T. Nevertheless, we consider this order of magnitude estimation (16)
being consistent with our proposal due to many uncertainties which enter this
estimate. There are several reasons for this optimistic view.

First, as we already mentioned the parameters entering (10) are in fact not
precisely known. The fundamental parameters such as $\rho_{\rm DM}$ and
$\langle v_{\rm AQN}\rangle$ only have accuracy up to order one as the local
flux distribution of DM and size distribution of AQN remain unknown to date as
the DM density locally may dramatically deviate from the well established
average global value $\rho_{\rm DM}\approx 0.3{\rm\,GeV\,cm^{-3}}$, see
introduction in [54] for the references and details. Furthermore, the size
distribution factor entering (10) in form $\langle B\rangle^{-1}$ had been
fixed from dramatically different physics (including solar corona heating
puzzle) and can easily deviate by large factor. The estimation of the
effective area ${\cal{A}}_{\rm eff}$ also suffers enormous uncertainty.

Furthermore, the numerical suppression factor $\sim 0.1$ (between computed
${\cal N}$ and observed ${\cal N}_{\rm obs}$ values) which appears in our
order of magnitude estimates (16) is very similar to suppression factor $\sim
0.1$ which occurred in analogous estimates for the mysterious TA bursts [9]
and the ANITA anomalous events [13]. This similarity hints on a common origin
for all these phenomena. Therefore, if one normalize ${\cal N}$ to the
observed ANITA anomalous events (assuming that all three phenomena originated
from the same AQN based physics) one arrives to correct order of magnitude
estimate as the suppression factor $\sim 0.1$ is approximately the same for
all three phenomena estimated for the mysterious TA bursts in [9], the ANITA
anomalous events in [13], and Multi-Modal Events estimated in (16).

Our main arguments, however, leading to the identification of the MMEs with
the AQN events are not based on an order of magnitude estimate of the
frequency of appearance (16) which is presented in this work exclusively for
illustrative purposes. Rather, our main arguments are based on specific
qualitative features (formulated as puzzles) which have been listed in Section
I, and which cannot be explained in terms of the conventional CR air showers
as highlighted in Section II. We consider a qualitative agreement between the
observations and our theoretical estimates (to be discussed in next section
VI) as the strong arguments supporting our identification.

##  VI AQN proposal confronts the MME observations [1, 2, 3, 4, 5]

We start this section by estimating the particle number density $\rho(R)$
which appears in formulation of the puzzles in Section II. The particle number
density in the AQN framework is determined by the number of particles which
could be observed at the detector site. Assuming that the total number of
electrons being emitted in form of a pulse (as a result of erupted burst) is
similar to the number of particles which was used in our analysis [13] of the
ANITA anomalous events $N\simeq(10^{8}-10^{9})$ we arrive to the following
estimate for $\rho(R)$:

$\displaystyle\rho(R)\sim\frac{N}{(\lambda\Delta\theta)^{2}}\sim\frac{(10^{8}-10^{9})}{(2.5~{}\rm
km)^{2}}\sim\frac{(20-200)}{\rm m^{2}}$ (17)

where the area which is hit by the bunch of particles from a single burst is
estimated as $(\lambda\Delta\theta)^{2}$ with $\lambda(h)$ given by (14) at
the detector site with $h\simeq 4$ km. It is assumed that the particles are
propagating within the cone with angle $\Delta\theta\lesssim 1/2$.

Our next task is the estimation of the time delay between the pulses within
the same cluster. The corresponding time delay can be estimated from the
condition that the energy deposited into the AQN during time $\tau_{\rm
delay}$ must be released in form of the eruption when $N\simeq(10^{8}-10^{9})$
electrons being emitted in form of a single pulse. As we mentioned previously,
the dominant portion of the energy in such eruption is released in form of the
photons with energy $\omega\sim T$ which always accompany the electron’s
emission as a result of these bursts. The relative ratio $r$ between these two
components is estimated in Appendix B, see eq. (38). Numerically, the estimate
can be presented as follows:

$\displaystyle r\equiv\frac{\langle E_{e^{+}e^{-}}\rangle}{\langle
E_{\gamma}\rangle}\sim\frac{\alpha}{2\pi}\exp\left(-\frac{2m}{T}\right)\sim
0.7\cdot 10^{-5},$ (18)

where for numerical estimates we use $T\simeq 200$ keV. We must emphasize that
(18) is really an order of magnitude estimate as mechanism of emission is
determined by a very complex physics as highlighted in Appendix B.
Furthermore, the corresponding eruption event is not a thermodynamically
equilibrium process as we already mentioned previously. The dominant
contribution in form of the photon’s emission during a single pulse can be
estimated in terms of these parameters as follows

$\displaystyle\langle E_{\gamma}\rangle\sim N\cdot(2m)\cdot r^{-1}\sim
1.4\cdot(10^{10}-10^{11})~{}\rm GeV,~{}~{}~{}~{}~{}~{}$ (19)

where we use the numerical value for $N\in(10^{8}-10^{9})$ extracted from
analysis of the ANITA anomalous events [13] and parameter $r$ is given by
(18). In this estimate we assumed that the energy of the pair $\langle
E_{e^{+}e^{-}}\rangle\sim 2mN$ computed at the moment of the pair creation.
The released energy during a single pulse (19) determines the time delay
$\tau_{\rm delay}$ between the pulses within the same cluster:

$\displaystyle\tau^{-1}_{\rm delay}\sim\frac{{dE_{\rm deposit}}/{dt}}{\langle
E_{\gamma}\rangle}\sim\kappa\frac{10^{18}\cdot\rm{GeV}s^{-1}}{(10^{10}-10^{11})~{}\rm
GeV},$ (20)

where we used numerical value for $({dE_{\rm deposit}}/{dt})$ from (III.3)
with density of atmosphere at $h\simeq 4$ km where the detector H10T is
located. Our final, the order of magnitude estimate for the time delay
$\tau_{\rm delay}$ between the pulses can be represented as follows

$\displaystyle\tau_{\rm
delay}\sim\left(\frac{0.1}{\kappa}\right)\cdot\left(10^{2}-10^{3}\right)\rm
ns,$ (21)

where we use $\kappa\simeq 0.1$ in (21) for numerical estimates which was
previously extracted from studies of the emission in the dilute galactic
environment, see footnote 3.

While this is very crude, an order of magnitude estimate, the important point
here is that the time scale $\tau_{\rm delay}$ may dramatically vary from
event to event as it strongly depends on many parameters entering the problem.
In particular, it obviously depends on intensity (number of emitted particles
$N$) of a previous pulse. Essentially the parameter (21) should be interpreted
as a preparation time the system requires for the next eruption.

The strong variation of the $\tau_{\rm delay}$ should be contrasted with
another parameter which is the time duration of a single pulse $\tau_{\rm
pulse}$. In the AQN framework this parameter must not vary much from one event
to another as it is entirely determined by internal dynamics of the AQN during
the eruption irrespective to the intensity of the previous events, and
irrespective to the prehistory of the AQN propagation as long as condition (8)
is met777This approximately constant parameter $\tau_{\rm pulse}$ can be
understood using the analogy with lightnings mentioned in footnote 6. Indeed,
the time scale between lightning flashes may dramatically vary (measured in
minutes) during the same thunderstorm while the lightning strikes themselves
are much shorter and characterized approximately by the same duration
(measured in ms).. In different words, our proposal suggests that

$\displaystyle\tau_{\rm pulse}\approx\rm constant,$ (22)

where the constant cannot be computed from first principles as it is entirely
determined by the non-equilibrium dynamics of the AQN at the instant of
eruption888The analogy with lightnings mentioned in footnotes 6 and 7 can be
useful here: the instability in form of the runaway breakdown mechanism in
lightnings which is responsible for the flashes is similar to non-equilibrium
dynamics of the AQN at the instant of eruption. Furthermore, in both cases
there must be some trigger which initiates the eruption. In case of lightning
events the trigger is thought to be related to the cosmic rays, though this
element remains a part of controversy, see [50, 51] for review. In case of the
AQN eruption events such triggers could be several consequent successful
annihilation events with atmospheric material..

Now we are prepared to confront the basic consequences of our proposal
(identifying the MMEs with the AQN annihilation events) with observations.

1. “clustering puzzle”: The multiple number of events is a very generic feature of the system as explained at the very end of subsection III.3 as long as condition (8) is met. The time delays between the pulses $\tau_{\rm delay}$ dramatically fluctuate between different events. The window for these variations is huge according to (21). In fact, it could be even well outside of this window. The time delays $\tau_{\rm delay}$ may considerably vary even for different events within the same cluster at the same detection point, which is consistent with observed cluster shown on Fig. 1. Furthermore, the AQN itself remains almost at the same location as the displacement $\Delta l_{\rm AQN}$ during entire cluster of events is very tiny

$\displaystyle\Delta l_{\rm AQN}\sim v_{\rm AQN}\cdot\tau_{\rm delay}\sim
20~{}\rm cm,$ (23)

which implies that all individual bunches making the cluster are likely to be
emitted by the AQN along the same direction, and can be recorded and
classified as MME by H10T. Each event can be viewed as an approximately
uniform front as mere notion of the “EAS axis” does not exist in this
framework, see also next item. However, each individual event may appear to
arrive from slightly different direction due to the inherent spread of the
emitted electrons at the moment of eruption.

2. “particle density puzzle”: Particle density distribution $\rho(R)$ in the AQN framework is estimated by (17) and it shows strong fluctuation from one event to another event. These variations are mostly related to the intensity of the individual bursts being expressed by the number of electrons in the bunch $N$. However, the distinct feature of the distribution in the AQN framework is that it does not depend on $R$, as the notion of the “EAS axis” does not exist in this framework as we already mentioned. All these generic features of the AQN framework are perfectly consistent with H10T observations as presented on Fig. 2. However, these observed features are in dramatic contradiction with conventional EAS prediction shown by solid lines on Fig. 2 with different colours, depending on energy of the CR.

Furthermore, the magnitude of the density $\rho(R)$ in the AQN framework
(which is mostly determined by parameter $N$) is much higher than one normally
expects for CR energies $\sim(10^{17}-10^{18})$ eV. The corresponding
parameter $N$ representing the number of electrons in the bunch was not fitted
for the present studies to match the observations. Instead, it was extracted
from different experiment in dramatically different circumstances (in proposal
[13] to explain the ANITA anomalous events as the AQN events).

3. “pulse width puzzle”: In the AQN framework the width of the pulse (22) cannot vary much from one event to another. It is a fundamental feature of the framework because the duration of the pulse is entirely determined by internal dynamics of the AQN during the blast as long as condition (8) is met, see footnote 7 with a comment. This feature is in perfect agreement with observations [3]

$\displaystyle\tau_{\rm pulse}\approx(20-35)\rm
ns~{}~{}~{}\Leftarrow~{}~{}~{}[\rm observations]$ (24)

for all recorded MMEs. At the same time, this feature is in dramatic conflict
with conventional picture when the duration of the pulse must depend on the
distance to the EAS axis as shown by sold lines on Fig. 3. This basic
prediction of the conventional CR analysis is due to increase of the thickness
of the EAS pancake with the distance from the EAS axis. As we already
mentioned the mere notions such as the “EAS axis” and the “thickness of the
EAS pancake” do not exist in the AQN framework.

4. “intensity puzzle”: Particle density distribution $\rho(R)$ in the AQN framework is estimated by (17). The corresponding event to event fluctuations do not depend on the distance to the EAS axis as we already mentioned. Such intensity of the events as given by (17) in the AQN framework is consistent with observations shown on Fig. 2. However, the observations are in dramatic conflict with conventional CR analysis when such huge intensity could be generated by a primary particle with energy well above $10^{19}$ eV with dramatically lower event rate on the level of once every few years. The frequency of appearance in the AQN framework is estimated in section V, and it is consistent with observed event rate.

We conclude this section with the following comment. All our formulae
presented in this section are the order of magnitude estimations at the very
best, as they include many inherent uncertainties which are inevitable
features of any composite system (such as AQN) propagating in a complex
environment (such as Earth atmosphere) with very large Mach number $M\equiv
v_{\rm AQN}/c_{s}\gg 1$ where $c_{s}$ is speed of sound.

Nevertheless, the emergent picture suggests that all the puzzles formulated in
Sections I and II can be naturally understood within the AQN framework as
explained above in items 1-4. Needless to say that all phenomenological
parameters used in the estimates had been fixed long ago for dramatically
different observations in different circumstances for different environments
as overviewed in Section III.

##  VII Conclusion and future development

Our basic results can be summarized as follows. We have argued that all the
puzzles formulated in Sections I and II can be naturally understood within the
AQN framework as explained in items 1-4 in previous Section, and we do not
need to repeat these arguments again in this Conclusion. Instead, we want to
discuss the drastic differences between the events induced by conventional CR
showers and the AQNs. These dramatic distinct features can be tested by the
future experiments, such that our proposal can be discriminated from any other
proposals and suggestions. We list below the following typical features of the
AQN events and contrast them with any other possible mechanisms which could be
responsible for MMEs.

1\. The events which are generated by the bunches of electrons as a result of
eruption of the propagating AQN in the Earth atmosphere suggests an enormous
number of possibilities to generate different clusters when each event within
a given cluster may have very different intensity from a previous and
consequent events with very different time delays between the events. In other
words, the AQN proposal suggests that there should be large variety of shapes
and delays between the events with very different patterns due to the
complexity of the AQN system. It should be contrasted, for example, with
hypothesis of “delayed particles” (which was originally suggested to explain
the MMEs) in which case all clusters must be the same as they should be
determined by a specific pattern of decaying fundamental particle of unknown
nature.

2\. A “rule of thumb” suggests that a typical 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 the energy of the primary
particle $E_{\rm CR}\approx(10^{17}-10^{18})$ eV. This estimate suggests that
any detector which is designed to study the EAS with energies $E_{\rm
CR}\approx(10^{17}-10^{18})$ eV are, in principle, capable to study MMEs if
the resolution of the detectors is in $\sim 10$-ns level, similar to H10T, see
also item 4 below as an alternative option to properly select and discriminate
the MMEs.

3\. In particular, we expect that the extension of the H10T detector would
produce more multiple pulses (at each given detector) instead of simple
bimodal pulses. We also expect that more detectors in the area will be
recording MMEs because the area covered by each individual pulse is relatively
large (few kilometres) according to (17), which is well above the present size
of H10T instrument.

4\. A large number of charged particles $N\approx(10^{8}-10^{9})$ in the
background of the geomagnetic field ${\cal{B}}\sim 0.5$ gauss will produce the
radio pulse in both cases: the CR-induced radio pulse [55, 56] as well as AQN-
induced radio pulse [13]. However, these pulses can be easily discriminated
from each other as argued in [13].

The main reason for the dramatic differences between these two radio pulses is
that the AQN event could be viewed as an (approximately) uniform front of size
$\lambda\Delta\theta\sim$km with a constant width, while EAS is characterized
by central axis. In different words, the number of particles per unit area
$\rho(R)$ in the AQN case does not depend on the distance from the central
axis, in huge contrast with conventional CR air showers when $\rho(R)$
strongly depends the distance from the central axis. The width of the
“pancake” in CR air shower also strongly depends on $R$. As a result, the
effective number of coherent particles contributing to the radio pulse is
highly sensitive to the width of the “pancake” when it becomes close to the
wavelength of the radio pulse. These distinct features lead to very different
spectral properties of the radio pulses in these two cases, which can be
viewed as an independent characteristic of MMEs. In fact, this unique feature
can be used in future studies for purpose of discrimination and proper
selection of the Multi-Modal Clustering Events.

If future studies and tests (including the detecting of the synchronized radio
pulses with MMEs as suggested above) indeed substantiate our proposal it would
be a strong argument supporting the AQN nature of the MMEs.

We conclude this work with the following comment. We estimated the event rates
for three dramatically different puzzling CR-like events: for mysterious TA
bursts [9], for the ANITA anomalous events [13], and finally for the Multi-
Modal Events in this work as estimated in (16). All three puzzling phenomena
are proportional to one and the same AQN flux (10). The self-consistency
between all three estimates hints on a common nature for these puzzling CR-
like events. We interpret this self-consistency in the event rates as an
additional indirect argument supporting the AQN nature for all three
mysterious phenomena, while our direct arguments are presented in section VI
and listed as item 1-4. We finish on this optimistic note.

## Acknowledgements

I am very thankful to Vladimir Shiltsev of FNAL who brought to my attention
the original references [1, 2, 3, 4, 5] and asked me many questions on
possible nature of MMEs, which essentially motivated and initiated this work.
This research was supported in part by the Natural Sciences and Engineering
Research Council of Canada.

##  Appendix A On suppression of emissivity due to the ionization

The main goal of this section is to produce simple estimates of the
suppression parameter $\eta(T,R)$ which enters the cooling rate (4) in case of
high temperature and large ionization. It was introduced in [46] as a
suppression parameter accounting for the ionization of the AQN after it
travelled through the Earth. It was also shown that $\eta(T,R)\sim 10^{-6}$
produces very reasonable results consistent with intensity of the radiation
observed by XMM-Newton [45].

The original computations of the Bremsstrahlung radiation were performed in
[48] in case of low temperature and low ionization as given by (2). The key
element of the computations was the electrosphere density

$n(z)=\frac{T}{2\pi\alpha}\frac{1}{(z+\bar{z})^{2}},$ (25)

with

$\bar{z}^{-1}=\sqrt{2\pi\alpha}\cdot
m\cdot\left(\frac{T}{m}\right)^{1/4},~{}~{}~{}~{}n(z=0)\simeq\left(mT\right)^{3/2},$

where $n(z)$ is the positron number density of the electrosphere, $z=0$
corresponds to the surface of the nuggets, and $n(z=0)$ reproduces an
approximate formula for the plasma density in the Boltzmann regime at the
temperature $T$ as computed in [48]. The 1d approximation formulated in terms
of distance $z$ from the surface was more than sufficient sufficient to study
the low temperature behaviour.

However, after an AQN crosses the Earth, it acquires a very high internal
temperature of the order $T\simeq(100-200)$ keV, and approximation (25) is not
sufficient. This is because temperature’s rise will cause the electrosphere to
expand well beyond the thin layer $\sim\bar{z}$ surrounding the nugget
surface. Some positrons will leave the system but the majority will stay in
close vicinity of the moving, negatively charged nugget core. Consequently,
the nugget will acquire a negative charge of approximately $-|e|Q$ with the
number of positrons $Q$ estimated as:

$\displaystyle Q\simeq 4\pi R^{2}\int^{\infty}_{0}n(z)dz\sim\frac{4\pi
R^{2}}{\sqrt{2\pi\alpha}}\cdot\left(mT\right)\cdot\left(\frac{T}{m}\right)^{1/4}.~{}~{}~{}~{}~{}~{}$
(26)

The distance $\rho$ at which the positrons remain attached to the nugget is
given by the capture radius $R_{\rm cap}(T)$, determined by the Coulomb
attraction:

$\displaystyle\frac{\alpha Q(\rho)}{\rho}>\frac{mv^{2}}{2}\approx
T~{}~{}~{}{\rm for}~{}~{}~{}\rho\lesssim R_{\rm cap}(T).$ (27)

In order to estimate $\eta$ entering (4) we start our analysis with an
estimation of the positron density $n(\rho,T)$ when the AQNs enter the Earth
atmosphere moving in upward direction. Formula (27) shows that the capture
radius $R_{\rm cap}(T)$ could be much larger than $R$, i.e. we have
$R\lesssim\rho\lesssim R_{\rm cap}$. The expression (25), which is valid when
$z=|\rho-R|\ll R$, in close vicinity to the core, breaks down for $\rho\gtrsim
R$. In that case the curvature of the nugget surface cannot be neglected and
one should use a truly 3-dimensional formula for $n(\rho,T)$ instead of the
1-dimensional approximation (25). We will assume that the density $n(\rho,T)$
behaves as a power law for $\rho\gtrsim R$ with exponent $p$:

$\displaystyle n(\rho,T)\simeq
n_{0}(T)\left(\frac{R}{\rho}\right)^{p},~{}~{}~{}\rho\gtrsim R,$ (28)

where $n_{0}(T)\equiv n(\rho=R,T)$ is a normalization factor and $p$ a free
parameter. Eq. (28) is consistent with our previous numerical studies [57] of
the electrosphere with $p\simeq 6$. It is also consistent with the Big Bang
Nucleosyntheses (BBN) study discussed in [58] and with the conventional
Thomas-Fermi model999In [59], the dimensionless function $\chi(x)$ behaves as
$\chi\sim x^{-3}$ at large $x$. The potential $\phi=\chi(x)/x$ behaves as
$\phi\sim x^{-4}$. The density of electrons in Thomas-Fermi model scales as
$n\sim\phi^{3/2}\sim x^{-6}$ at large $x$. at $T=0$ [59]. While keeping $p$ as
a free parameter, we will show below that our main claim is not very sensitive
to the precise value of $p$. The normalization $n_{0}(T)$ can be estimated
from the condition that the finite portion of the positrons satisfying Eq.
(27) is such that the total number of positrons surrounding the nugget is
approximately equal to the ionization charge $Q$ determined by (26).
Therefore, we arrive to the following expression for the normalization factor
$n_{0}(T)$:

$\displaystyle 4\pi\int_{R}^{\infty}\rho^{2}d\rho~{}n(\rho,T)\simeq Q(T),$
(29)

Note that the integrant $\rho^{2}n(\rho,T)$ is mostly dominated by the inner
shells $\rho\sim R$ such that the intgration can be extended to infinity with
very high accuracy, instead of cutting off at $R_{\rm cap}$. The resulting
estimate for $n_{0}(T)$ assumes the form:

$\displaystyle
n_{0}(T)\simeq\frac{\left(mT\right)}{\sqrt{2\pi\alpha}}\cdot\frac{(p-3)}{R}\cdot\left(\frac{T}{m}\right)^{1/4},$
(30)

which, as expected, is smaller than $n(z=0)\simeq\left(mT\right)^{3/2}$
because the positrons now are distributed over a distance of order $R$ from
the nugget core surface, rather than over a distance of order $\bar{z}$. This
suppression is convenient to represent in terms of the dimensionless ratio:

$\displaystyle\frac{n_{0}(T)}{\left(mT\right)^{3/2}}\simeq\frac{(p-3)}{\sqrt{2\pi\alpha}}\cdot\frac{1}{(mR)}\cdot\left(\frac{m}{T}\right)^{1/4}\ll
1.$ (31)

The next step is the calculation of the Bremsstrahlung emissivity by the
nuggets due to the strong suppression of the plasma density (31) as a result
of the electrosphere’s expansion. The spectral surface emissivity is denoted
as $dF/d\omega=dE/dtdAd\omega$, representing the energy emitted by a single
nugget per unit time, per unit area of nugget surface and per unit frequency.
For low temperature, when approximation (25) for the flat geometry is
justified, the corresponding expression assumes the form [48]:

$\frac{dF}{d\omega}(\omega,T)=\frac{1}{2}\int_{0}^{\infty}dz~{}n^{2}(z)~{}{\cal{K}}(\omega,T)$
(32)

where $n(z)$ is the local density of positrons at distance $z$ from quark
nugget surface and function ${\cal{K}}(\omega)$ does not depend on $z$ and
describes the spectral dependence of the system:

${\cal{K}}(\omega,T)=\frac{4\alpha}{15}\left(\frac{\alpha}{m}\right)^{2}2\sqrt{\frac{2T}{m\pi}}\left(1+\frac{\omega}{T}\right){\rm
e}^{-\omega/T}h\left(\frac{\omega}{T}\right).$ (33)

The dimensionless function $h\left(\frac{\omega}{T}\right)$ in (33) is a
slowly varying logarithmic function of ${\omega}/{T}$ which was explicitly
computed in [48]. In order to calculate the emissivity with the positron
density (28), when the nugget core is ionized, we have to replace the integral
$\int dz$ in (32) by the spherical integration over $\rho$:

$\displaystyle 4\pi
R^{2}\int_{0}^{\infty}dz~{}n^{2}(z)~{}~{}\Rightarrow~{}4\pi
n^{2}_{0}(T)\int_{R}^{\infty}\rho^{2}d\rho\left(\frac{R}{\rho}\right)^{2p}~{}~{}~{}$
(34)

where the normalization $n_{0}(T)$ is determined by (30). Using the positrons
distributed according to Eq. (28) the nugget surface emissivity can be
calculated as follows:

$\displaystyle\frac{dF}{d\omega}(\omega,T)=\frac{{\cal{K}}(\omega,T)}{2R}\frac{(mT)^{2}}{2\pi\alpha}\left[\frac{(p-3)^{2}}{2p-3}\right]\left(\frac{T}{m}\right)^{\frac{1}{2}}.~{}~{}~{}~{}~{}~{}~{}$
(35)

The key result here is the strong suppression factor $R^{-1}$ which was not
present in the no-ionization case (32) when the positrons are localized in a
thin layer. It is instructive to compare the emissivity given by eq. (35)
accounting for ionization of the nugget core with the original formula (32).
The dimensionless ratio is given by

$\displaystyle\eta\equiv\frac{\frac{dF}{d\omega}^{(\rm
ion)}}{{\frac{dF}{d\omega}^{(\rm
no~{}ion)}}}=\frac{\left[\frac{(p-3)^{2}}{2p-3}\right]}{3(mR)}\frac{1}{\sqrt{2\pi\alpha}}\left(\frac{m}{T}\right)^{\frac{1}{4}}\sim
10^{-6}.~{}~{}~{}~{}~{}~{}~{}~{}$ (36)

The total emissivity integrated over all frequencies in case of low
temperature without ionization has been computed in [48] and it is given by
(2). In case of strong ionization, Eq. (36) implies:

$\displaystyle F^{(\rm ion)}_{\rm tot}(T)=\eta(T,R)\cdot F^{(\rm
no~{}ion)}_{\rm tot}(T).$ (37)

The critical element leading to the suppression factor (36) is the very small
quantity $(mR)^{-1}\sim 10^{-6}$. It is a direct result of the expanding
electrosphere and high level of ionization at very high temperature $T$.

The expression (37) with numerical suppression (36) is precisely the formula
(4) which we used in main body of the text to proceed with our key arguments.

##  Appendix B On $e^{+}e^{-}$ emission at high temperature in high density
QCD phases.

The main goal of this Appendix is to overview the photon emission and the
$e^{+}e^{-}$ production at high temperature $T\gtrsim 10^{2}$ keV from high
density QCD phases. The corresponding studies [60, 61, 62, 63, 64, 65, 66]
have been carried out in the past in context of the quark stars. In context of
the present work all the key ingredients relevant for $e^{+}e^{-}$ production
are also present in the system. Indeed, the AQN is characterized by very high
temperature $T\gtrsim 10^{2}$ keV, the quark core is assumed to be in CS dense
phase, and a strong internal electric field is also present in the system.
However, we cannot literally use the results from the previous studies
obtained in context of the quark stars because the size of the AQN is much
smaller than the relevant mean free paths for all elementary processes as
discussed in details in [13]. As a result the thermal equilibration cannot be
achieved in the AQN system, and entire physics is determined by non-
equilibrium dynamics in high temperature regime. It should be contrasted with
large size quark stars where the thermal equilibrium is maintained.

Nevertheless, it is very instructive to review the relevant results from the
previous studies [60, 61, 62, 63, 64, 65, 66] on quark stars due to the
following reasons. First, it explicitly shows the role of the main ingredients
of the system, such as temperature and the large electric field. Secondly, it
demonstrates the complexity of the problem when even a much simpler case of
the bare quark star remains to be a matter of debates.

The idea on possibility of the $e^{+}e^{-}$ emission at high temperature from
quark stars was originally suggested in [60, 61]. The temperature range
considered in [60, 61] includes the typical temperatures $T\gtrsim 10^{2}$ keV
which is expected to occur in our case when the AQN exits the Earth’s surface
as mentioned in Section III.3. In refs. [62, 63] the authors argue that
bremsstrahlung radiation from the electrosphere could be much more important
than $e^{+}e^{-}$ emission. It has been also argued that a number of effects
such as the boundary effects, inhomogeneity of the electric field, and the
Landau-Pomeranchuk-Migdal (LPM) suppression may dramatically modify the
emission rate. In [64] it has been argued that the Pauli blocking will
strongly suppress the bremsstrahlung emission. Finally, in refs. [65, 66] it
has been argued that the so-called mean field bremsstrahlung could be the
dominant mechanism.

It is not the goal of the present work to critically analyze all these
suggested mechanisms of the emission. Rather, our goal is to demonstrate that
even a relatively simple system of the bare quark star remains to be a matter
of debates. Our case of the emission from AQN at high temperature is even more
complicated as it is determined by non-equilibrium dynamics.

In the present work we do not even attempt to solve this very ambitious
problem of estimating the absolute intensity rate of the $e^{+}e^{-}$ and
$\gamma$ emissions. Rather, the absolute number $N\simeq(10^{8}-10^{9})$ of
the $e^{+}e^{-}$ pairs in a single pulse was extracted from ANITA anomalous
events as explained in Section VI. In particular, this number $N$ enters the
estimate (17) for the density of the observed particles by H10T experiment [1,
2, 3, 4]. Our goal here is much less ambitious as we try to estimate the
relative ratio $r$ between the $e^{+}e^{-}$ pair production and the $\gamma$
radiation assuming that the emission of both components is mostly originated
from the region in electrosphere where the Boltzmann regime is justified, and
the Pauli blocking suppression is not very dramatic and can be ignored.

In this case, the relative ratio $r$ between these two components contributing
to the emission can be estimated as follows. It is assumed that the dominant
contribution to the $\gamma$ emission comes from the Bremsstrahlung radiation
($e^{+}e^{+}\rightarrow e^{+}e^{+}+\gamma$) resembling the low temperature
case considered in [48]. The $e^{+}e^{-}$ pair is produced by a similar
mechanism through the virtual photon ($e^{+}e^{+}\rightarrow
e^{+}e^{+}+\gamma^{*}$) which consequently decays to the pair:
$\gamma^{*}\rightarrow e^{+}e^{-}$. If this is the dominant mechanism of
radiation the pair production is expected to be suppressed by a factor
${\alpha}/{2\pi}$ as a result of conversion of the virtual photon
$\gamma^{*}\rightarrow e^{+}e^{-}$ to pair101010The cross section for pair
production as a result of collisions of two particles is known exactly[67]. It
contains, of course factor ${\alpha}/{2\pi}$ entering (38) along with many
other numerical factors reflecting a complex kinematics of the process.. It
must be also suppressed by a factor $\exp(-{2m}/{T})$ describing the
suppression in the density distribution because the virtual photons must have
sufficient energy $\gtrsim 2m$ to produce $e^{+}e^{-}$ pair. This
oversimplified estimate leads to the following expression

$\displaystyle r\equiv\frac{\langle E_{e^{+}e^{-}}\rangle}{\langle
E_{\gamma}\rangle}\sim\frac{\alpha}{2\pi}\exp\left(-\frac{2m}{T}\right)\sim
0.7\cdot 10^{-5},$ (38)

where for numerical estimates we use $T\simeq 200$ keV, and ignored a possible
complicated function which could depend on parameter $T/m\sim 0.4$, which is
order of unity for our case. We must emphasize that (38) is really an order of
magnitude estimate as mechanism of emission is not a thermodynamically
equilibrium process as we already previously mentioned. Formula (38) is used
in the main body of the text in (18) for our estimate (21) for $\tau_{\rm
delay}$ which is measured by H10T experiment [1, 2, 3, 4] as it describes a
typical time delay between subsequent pulses representing the Multi-Modal
Clustering Events.

## References

  * Beznosko  _et al._ [2020] D. Beznosko _et al._ , J. Phys. Conf. Ser. 1342, 012007 (2020), arXiv:1605.05179 [physics.ins-det] .
  * Beznosko  _et al._ [2017] D. Beznosko, R. Beisembaev, K. Baigarin, E. Beisembaeva, O. Dalkarov, V. Ryabov, T. Sadykov, S. Shaulov, A. Stepanov, M. Vildanova, N. Vildanov, and V. Zhukov, in  _European Physical Journal Web of Conferences_, European Physical Journal Web of Conferences, Vol. 145 (2017) p. 14001.
  * Beznosko  _et al._ [2019] D. Beznosko, R. Beisembaev, E. Beisembaeva, O. D. Dalkarov, V. Mossunov, V. Ryabov, S. Shaulov, M. Vildanova, V. Zhukov, K. Baigarin, and T. Sadykov, PoS ICRC2019, 195 (2019).
  * Beisembaev  _et al._ [2019] R. Beisembaev, D. Beznosko, K. Baigarin, A. Batyrkhanov, E. Beisembaeva, O. Dalkarov, A. Iakovlev, V. Ryabov, T. Sadykov, S. Shaulov, M. Vildanova, T. Uakhitov, and V. Zhukov, in  _European Physical Journal Web of Conferences_, European Physical Journal Web of Conferences, Vol. 208 (2019) p. 06002.
  * Beisembaev  _et al._ [2019] R. U. Beisembaev _et al._ , Phys. Atom. Nucl. 82, 330 (2019).
  * Zhitnitsky [2003] A. R. Zhitnitsky, JCAP 10, 010 (2003), hep-ph/0202161 .
  * Jelley and Whitehouse [1953] J. V. Jelley and W. J. Whitehouse, Proceedings of the Physical Society. Section A 66, 454 (1953).
  * Linsley and Scarsi [1962] J. Linsley and L. Scarsi, Phys. Rev. 128, 2384 (1962).
  * Zhitnitsky [2020] A. Zhitnitsky, J.Phys.G:Nucl.Part.Phys. (2020), 10.1088/1361-6471/abd457, arXiv:2008.04325 [hep-ph] .
  * Liang and Zhitnitsky [2021a] X. Liang and A. Zhitnitsky, (2021a), arXiv:2101.01722 [hep-ph] .
  * Abbasi  _et al._ [2017] R. Abbasi _et al._ (Telescope Array Project), Phys. Lett. A 381, 2565 (2017).
  * Okuda [2019] T. Okuda, Journal of Physics: Conference Series 1181, 012067 (2019).
  * Liang and Zhitnitsky [2021b] X. Liang and A. Zhitnitsky, (2021b), arXiv:2105.01668 [hep-ph] .
  * Gorham  _et al._ [2016] P. W. Gorham _et al._ (ANITA), Phys. Rev. Lett. 117, 071101 (2016), arXiv:1603.05218 [astro-ph.HE] .
  * Gorham  _et al._ [2018] P. W. Gorham _et al._ (ANITA), Phys. Rev. Lett. 121, 161102 (2018), arXiv:1803.05088 [astro-ph.HE] .
  * Gorham  _et al._ [2021] P. W. Gorham _et al._ (ANITA), Phys. Rev. Lett. 126, 071103 (2021), arXiv:2008.05690 [astro-ph.HE] .
  * Zhitnitsky [2021] A. Zhitnitsky, Mod. Phys. Lett. A 36, 2130017 (2021), arXiv:2105.08719 [hep-ph] .
  * Witten [1984] E. Witten, Phys. Rev. D 30, 272 (1984).
  * Farhi and Jaffe [1984] E. Farhi and R. L. Jaffe, Phys. Rev. D 30, 2379 (1984).
  * De Rujula and Glashow [1984] A. De Rujula and S. L. Glashow, Nature (London) 312, 734 (1984).
  * Madsen [1999] J. Madsen, in  _Hadrons in Dense Matter and Hadrosynthesis_, Lecture Notes in Physics, Berlin Springer Verlag, Vol. 516, edited by J. Cleymans, H. B. Geyer, and F. G. Scholtz (1999) p. 162, astro-ph/9809032 .
  * Peccei and Quinn [1977] R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977).
  * Weinberg [1978] S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek [1978] F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Kim [1979] J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman  _et al._ [1980] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B 166, 493 (1980).
  * Dine  _et al._ [1981] M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B 104, 199 (1981).
  * Zhitnitsky [1980] A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980), [Yad. Fiz.31,497(1980)].
  * Van Bibber and Rosenberg [2006] K. Van Bibber and L. J. Rosenberg, Physics Today 59, 30 (2006).
  * Asztalos  _et al._ [2006] S. J. Asztalos, L. J. Rosenberg, K. van Bibber, P. Sikivie, and K. Zioutas, Annual Review of Nuclear and Particle Science 56, 293 (2006).
  * Sikivie [2008] P. Sikivie, in  _Axions_ , Lecture Notes in Physics, Berlin Springer Verlag, Vol. 741, edited by M. Kuster, G. Raffelt, and B. Beltrán (2008) p. 19, astro-ph/0610440 .
  * Raffelt [2008] G. G. Raffelt, in  _Axions_ , Lecture Notes in Physics, Berlin Springer Verlag, Vol. 741, edited by M. Kuster, G. Raffelt, and B. Beltrán (2008) p. 51, hep-ph/0611350 .
  * Sikivie [2010] P. Sikivie, International Journal of Modern Physics A 25, 554 (2010), arXiv:0909.0949 [hep-ph] .
  * Rosenberg [2015] L. J. Rosenberg, Proceedings of the National Academy of Science 112, 12278 (2015).
  * Marsh [2016] D. J. E. Marsh, Phys. Rep. 643, 1 (2016), arXiv:1510.07633 .
  * Graham  _et al._ [2015] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, Annual Review of Nuclear and Particle Science 65, 485 (2015), arXiv:1602.00039 [hep-ex] .
  * Irastorza and Redondo [2018] I. G. Irastorza and J. Redondo, Prog. Part. Nucl. Phys. 102, 89 (2018), arXiv:1801.08127 [hep-ph] .
  * Liang and Zhitnitsky [2016] X. Liang and A. Zhitnitsky, Phys. Rev. D 94, 083502 (2016), arXiv:1606.00435 [hep-ph] .
  * Ge  _et al._ [2017] S. Ge, X. Liang, and A. Zhitnitsky, Phys. Rev. D 96, 063514 (2017), arXiv:1702.04354 [hep-ph] .
  * Ge  _et al._ [2018] S. Ge, X. Liang, and A. Zhitnitsky, Phys. Rev. D 97, 043008 (2018), arXiv:1711.06271 [hep-ph] .
  * Ge  _et al._ [2019] S. Ge, K. Lawson, and A. Zhitnitsky, Phys. Rev. D 99, 116017 (2019), arXiv:1903.05090 [hep-ph] .
  * Jacobs  _et al._ [2015] D. M. Jacobs, G. D. Starkman, and B. W. Lynn, Mon. Not. Roy. Astron. Soc. 450, 3418 (2015), arXiv:1410.2236 [astro-ph.CO] .
  * Lawson  _et al._ [2019] K. Lawson, X. Liang, A. Mead, M. S. R. Siddiqui, L. Van Waerbeke, and A. Zhitnitsky, Phys. Rev. D 100, 043531 (2019), arXiv:1905.00022 [astro-ph.CO] .
  * Gorham [2012] P. W. Gorham, Phys. Rev. D 86, 123005 (2012), arXiv:1208.3697 [astro-ph.CO] .
  * Fraser  _et al._ [2014] G. W. Fraser, A. M. Read, S. Sembay, J. A. Carter, and E. Schyns, Mon. Not. Roy. Astron. Soc. 445, 2146 (2014), arXiv:1403.2436 [astro-ph.HE] .
  * Ge  _et al._ [2020] S. Ge, H. Rachmat, M. S. R. Siddiqui, L. Van Waerbeke, and A. Zhitnitsky, (2020), arXiv:2004.00632 [astro-ph.HE] .
  * Fox  _et al._ [2018] D. B. Fox, S. Sigurdsson, S. Shandera, P. Mészáros, K. Murase, M. Mostafá, and S. Coutu, (2018), arXiv:1809.09615 [astro-ph.HE] .
  * Forbes and Zhitnitsky [2008] M. M. Forbes and A. R. Zhitnitsky, Phys. Rev. D 78, 083505 (2008), arXiv:0802.3830 .
  * Raza  _et al._ [2018] N. Raza, L. van Waerbeke, and A. Zhitnitsky, Phys. Rev. D 98, 103527 (2018), arXiv:1805.01897 [astro-ph.SR] .
  * Gurevich and Zybin [2001] A. V. Gurevich and K. P. Zybin, Physics-Uspekhi 44, 1119 (2001).
  * Dwyer and Uman [2014] J. R. Dwyer and M. A. Uman, Phys. Rep. 534, 147 (2014), the Physics of Lightning.
  * [52] NOAA, National Centers for Environmental Information,  _Magnetic Field calculator_ ,  .
  * Jackson [1999] J. D. Jackson,  _Classical electrodynamics_, 3rd ed. (Wiley, New York, NY, 1999).
  * Liang  _et al._ [2021] X. Liang, E. Peshkov, L. Van Waerbeke, and A. Zhitnitsky, Phys. Rev. D 103, 096001 (2021), arXiv:2012.00765 [hep-ph] .
  * Huege and Falcke [2003] T. Huege and H. Falcke, Astron. Astrophys. 412, 19 (2003), arXiv:astro-ph/0309622 .
  * Huege and Falcke [2005] T. Huege and H. Falcke, Astroparticle Physics 24, 116 (2005).
  * Forbes  _et al._ [2010] M. M. Forbes, K. Lawson, and A. R. Zhitnitsky, Phys. Rev. D 82, 083510 (2010), arXiv:0910.4541 .
  * Flambaum and Zhitnitsky [2019] V. V. Flambaum and A. R. Zhitnitsky, Phys. Rev. D 99, 023517 (2019), arXiv:1811.01965 [hep-ph] .
  * Landau and Lifshitz [1975] L. D. Landau and E. M. Lifshitz,  _The classical Theory of Fields, 4-th edition_ (Butterworth-Heinemann, Oxford, 1975).
  * Usov [1998] V. V. Usov, Phys. Rev. Lett. 80, 230 (1998), arXiv:astro-ph/9712304 .
  * Usov [2001] V. V. Usov, Astrophys. J. Lett. 550, L179 (2001), arXiv:astro-ph/0103361 .
  * Jaikumar  _et al._ [2004] P. Jaikumar, C. Gale, D. Page, and M. Prakash, Phys. Rev. D 70, 023004 (2004).
  * Harko and Cheng [2005] T. Harko and K. S. Cheng, The Astrophysical Journal 622, 1033 (2005).
  * Caron and Zhitnitsky [2009] J.-F. Caron and A. R. Zhitnitsky, Phys. Rev. D 80, 123006 (2009), arXiv:0907.4715 [astro-ph.HE] .
  * Zakharov [2011] B. G. Zakharov, J. Exp. Theor. Phys. 112, 63 (2011), arXiv:1007.4697 [hep-ph] .
  * Zakharov [2010] B. G. Zakharov, Phys. Lett. B 690, 250 (2010), arXiv:1003.5779 [hep-ph] .
  * Berestetskii  _et al._ [1982] V. B. Berestetskii, L. P. Pitaevskii, and E. M. Lifshitz,  _Quantum Electrodynamics, 2-nd edition_ (Butterworth-Heinemann, Oxford, 1982).