# Gravitational focusing effects on streaming dark matter as a new detection
concept

Abaz Kryemadhia, Marios Maroudasb, Andreas Mastronikolisc, and Konstantin
Zioutasb a Dept. Computing,Math & Physics, Messiah University, Mechanicsburg
PA 17055, USA
b Physics Department, University of Patras, GR 26504,Patras-Rio, Greece
c Dept. Physics & Astronomy, University of Manchester, Manchester M13 9PL, UK

###### Abstract

This work suggests how to search for streaming constituents from the dark
sector. Following cosmological reasoning, a very large number of dark fine
grained streams are expected in the Galaxy. Within the solar system, planetary
gravitational focusing is effective for the slow DM particles, with possible
strong flux enhancements only for streaming dark constituents. The present
simulation incorporates gravitational effects preceding their encounter with
Earth bound detectors. The gravitational self-focusing by the inner Earth can
be large, and therefore must also be taken into account. Spatiotemporal
transient density enhancements are also derived, which are of interest for DM
experiments. Due to the built-in gravitational flux enhancements, advantages
are expected when compared to the detection schemes designed for the smooth
Standard Halo Model (SHM). In addition, a network of sensors would have a much
better discovery potential as long as the hypothetical streams of dark
constituents, along their rest mass and coupling strength to normal matter are
unknown. This work is relevant for all ongoing DM searches, including state-
of-the-art detection schemes like GNOME, EDM experiments, and the ECHO idea.
Also astrophysical searches for the dark sector could benefit more when
targeting planetary and solar/stellar atmospheres. The present simulation
focuses on fine grained streams but it can be applied to any configuration of
invisible matter propagating inside the solar system.

##  I Introduction

The dark matter constitutes about 80% of the matter in the universe and plays
an important role in structure formation. In the standard halo model (SHM)
dark matter is described as a cold collision-less fluid with a smooth
Maxwellian velocity distribution ($v\approx 0.001c$), which is an
approximation. At the last scattering as the universe expands and the DM fluid
cools it would have had a very small velocity dispersion Vogelsberger:2020 ;
Stucker:2020 .

A lot of progress has been made in studying the small scale structure of DM in
the galaxy particularly the velocity distribution in the halo despite the
N-body computational challenges related to finer resolution requirement
Springel:2021 ; Vogelsberger:2020 . The study of small-scale structure of DM
spanning the size of solar system is important as it can shed light on local
DM density distribution and eventually on the flux of particles in Earth bound
experiments Vogelsberger:2009 ; Vogelsberger:2011 .

Vogelsberger et. al. Vogelsberger:2009 ; Vogelsberger:2011  used a geodesic
deviation equation and the N-body equations of motion to study extensively the
evolution of phase space of the halos and concluded that many caustics and
fine grained streams could be present at the position of the solar system. The
simulations carried out by Vogelsberger et al. Vogelsberger:2011  for fine
grained streams indicate that up to $\sim 10^{12}$ such streams could be
present in the solar system, and about $10^{6}$ streams contain half the
average local DM density. Consequently, the DM distribution at a point in the
halo is described as the superposition of many fine-grained streams, each of
which has a small velocity dispersion. This can be an advantage for some
experimental searches. Fine grained streams would appear in detectors with
very good energy or velocity resolution as narrow peaks at specific velocities
and directions. The primordial velocity dispersions of the streams are $\sim
10^{-10}$c for WIMPs and $\sim 10^{-17}$c for axions Vogelsberger:2009 ;
Vogelsberger:2008 ; Sikivie:1995 ; SikivieVelo:1999 ; natarajan:2008 .
Existing detectors might not have sufficient resolution to distinguish the
discrete nature following the iso-thermal model. Though in future, detectors
designed with appropriate sensitivity could unravel the fine grained
configuration. In order to probe the fine-grained nature one would need high
sensitivity detectors with extremely good resolution. This work suggests to
take into account for the first time gravitational focusing effects by the
intervening solar system bodies before a stream encounters a detector at a
near earth location. Hence, the flux enhancement can be by several orders of
magnitude and may allow to observe the so far overlooked DM. Here the
quantitative estimates of the expected flux enhancements are given. We wish to
stress that DM experiments can benefit from, as this might only require to
modify the analysis code or adopt the detector design in order to be sensitive
to transients due to streaming DM.

DM streams propagating inside the solar system experience gravitational
effects by the Sun, Earth, Moon and the other planets of the solar system.
Their small dispersion is advantageous to reach large density enhancements by
focusing effects of solar bodies as they can focus non-relativistic DM
constituents on Earth Hoffmann:2003 ; alenazi:2006 ; Patla:2014 ; Prezeau:2015
; sofue:2020 . Moreover, the small dispersion is attractive for axion
haloscopes because of the coherence length and time increase with smaller
velocity dispersion. This is described as a quality factor benefit and it is
attributed to the fact that axion to photon conversions in the cavity will
grow coherently over large coherence time. (See Refs. slee:2020 ; admx:2018 ;
CAST:2017 )

Gravitational impact on DM by solar system bodies have been studied before. We
mention here briefly a number of earlier relevant publications, e.g.:
Hoffmann, Jacoby, and Zioutas Hoffmann:2003 , and Patla et al. Patla:2014
studied gravitational effect of streaming non-relativistic particles by the
Sun and other solar bodies. Alenazi and Condolo, studied phase-space
distribution of unbound dark matter near the Sun producing sky maps with
density enhancement due to gravitational focusing by the Sun alenazi:2006 .
Lee at al. lee:2014  studied the effect of gravitational focusing on annual
modulation in dark-matter direct-detection experiments. Prezeau Prezeau:2015
used the geodesic equation to show strong enhancement by Earth, Jupiter and
other compact bodies, also discussed feasibility for the number of probes
needed to detect these large amplification regions. Prezeau came to the
conclusions that a huge number of probes were needed making it hard to detect
but encouraged further studies. Prezeau used particles with speeds of 220 km/s
in the calculations quoting mostly a far away focus distance from the Earth.
The coverage area needed scales as distance2, so aside from the challenge of
sending probes so far away in space, a very large number of them are needed to
cover such a large area. Prezeau dubbed regions in space with large density
enhancements as DM hairs Prezeau:2015 . Sofue discussed self-focusing by the
inner Earth on DM particles, which is of particular consequence for direct
detection of DM. Sofue demonstrated that it’s possible to get high
amplifications for dispersion-less velocity particles when one explores
planetary self-focusing effects sofue:2020 . We highlight the fact that the
difference of the smooth SHM case caused by gravitational focusing (GF) is a
flux change of some % and at certain times of the year, while here the
enhancements can be up to several orders of magnitude (see also Refs.
Hoffmann:2003 ; Patla:2014 ; Prezeau:2015 ; sofue:2020 ).

Prompted by the prospect for large density enhancements with low velocity
streams which could be within reach of space and ground based experiments, in
this work we developed simulation tools to study at more extensive level the
gravitational effects. The Earth’s inclined ecliptic plane with respect to
galactic plane imposes kinematical constrain for streaming DM particles to
reach the Earth. Moreover, the combined GF by the Sun and its planets
including the Moon might open a new window of opportunity in DM research. And,
this is one novelty of the present work to also include GF by the Sun and
Earth together; the Sun modifies the phase-space of these particles falling on
its gravitational well at 1 A.U. after which the particles experience
additional self-focusing by the inner Earth. This is important for Earth bound
DM searches. Furthermore, also the dispersion effects on density enhancements
are taken into account. We pivoted the study on fine grained streams of axions
and WIMPs, however the tools could be broadly used to evolve the phase-space
of DM near the Sun and the Earth, in particular.

Furthermore, we refer to gravitational (G) produced (P) density (D)
enhancement (E) as GPDE regions. Also whenever we discuss streams we imply
fine grained streams as compared to for instance tidal streams that are caused
by halo disruptions Gaia:2017 . Therefore, streams of cosmological origin with
much smaller velocity dispersion seem to be more advantageous for direct DM
detection. The spatial extent of streams is considered to be at least the size
of solar system and the spatial distribution of particles in a stream is
assumed to be uniform in the absence of any additional gravitational effects
throughout this work.

Self-focusing plays an important role for DM searches at the Earth location.
The reason is that it is always at work.

This paper is organized as follows: section II describes the simulation
approach used, followed by the results for the Sun Earth system in section
III. We then discuss frequency and prospect for detection of these GPDE
regions of different densities in IV and the summary with conclusions is
outlined in section V.

##  II Simulation Framework

The first step is to compute the combined gravitational focusing (GF) effects
by the Sun and the Earth on the phase-space density of fine grained streams at
the lab frame given any known but arbitrary chosen input phase-space far away
from the solar system in the galactic rest frame. This allows to explore the
potential of DM experiments after DM particles fall in the gravitational well
of the Sun, and then eventually focus at the Earth’s position. The combined
gravitational focusing effects by the Sun and the Earth yield regions of high-
density near the Earth. Figure 1 illustrates how DM particles from one stream
are curved by the Sun and the Earth and a density enhancement results on the
opposite side of the Earth relative to the incoming direction of the stream.

Figure 1: Top: Schematics of DM particles being first gravitationally
deflected by the Sun and then self-focusing occurs by the inner Earth. Bottom:
A map in equatorial coordinates of simulated DM particles in the SHM. The
hurtling of solar system towards Cygnus causes the typical bias towards more
particles coming from its direction confirming that our simulation is setup
correctly, as it produces a widely accepted effect.

We recall that the velocity distribution in the SHM is modeled with a Gaussian
function where DM is on average at rest with respect to the galaxy and the
standard deviation is taken to be
$\sigma_{\textrm{v}}\approx{\textrm{2}30km/s}/\sqrt{2}\approx
160~{}{\textrm{k}m/s}$. While the SHM is the most simplified commonly used
model for DM distribution, there are also models where a fraction of DM is in
a co-rotating dark disk to mimic the baryonic counterpart. For instance,
following Ref. purcel:2009 , up to 32% of the average local DM could be in the
co-rotating dark disk (dd). For generality, in our simulation we construct the
total velocity distribution function as a combination of an SHM component and
a co-rotating fraction as it is given by:

$f_{\textrm{t}ot}(\textbf{v}_{g})=(\frac{1-\rho_{\textrm{d}d}}{\rho_{0}})f_{\textrm{S}HM}(\bf
v_{g})+\frac{\rho_{\textrm{d}d}}{\rho_{0}}f_{\textrm{d}d}(\bf v_{g}),$ (1)

where $\rho_{\textrm{d}d}$ is the DM density of the dark disc content and
$\rho_{0}$ is the mean nominal DM density of
$\rho_{\text{0}}=0.45\frac{\textrm{GeV}}{\textrm{cm}^{3}}$. We have assumed
that the dark disc fraction is $\frac{\rho_{\textrm{d}d}}{\rho_{0}}=0.32$
purcel:2009  as an upper limit aiming to probe a wide range of the phase
space. We used multivariate Gaussian for both $f_{\textrm{S}HM}(\bf v_{g})$
and $f_{\textrm{d}d}(\bf v_{g})$. Table 1 lists velocities and dispersion used
in our simulations of the SHM and the co-rotating dark disk parameters where
dd is 32% (dd32) of the nominal DM density purcel:2009 .

Table 1: Input Simulation Parameters in galactic frame. Parameter | SHM | dd32
---|---|---
$\sigma_{r}$ | 160 km/s | 133 km/s
$\sigma_{\phi}$ | 160 km/s | 78 km/s
$\sigma_{z}$ | 160 km/s | 63 km/s
$v_{\phi}$ | 0 km/s |  $v_{\textrm{L}SR}$+11.7 km/s
$v_{r}$ | 0 km/s | 0 km/s
$v_{z}$ | 0 km/s | 0 km/s

We now turn our attention to the velocity distribution of fine grained
streams. According to Refs. Sikivie:1995 ; Vogelsberger:2008  a fraction of
cold DM would have thermalized due to undergoing many orbits in the Milky Way
and acquired Maxwellian velocity distribution. Cold DM with fewer orbits would
have sharper peaks with accordingly smaller velocity dispersion. It is not
known yet, however what fraction of cold DM has undergone fewer orbits. In
addition the high density streams might be more likely to have performed fewer
orbits compared to low density ones and therefore have different velocity
distribution from low density ones. In the absence of knowledge about the
fraction of DM in streams or whether velocity distribution is shaped
differently for high density streams compared to low density ones, in our
simulation we sample streams with discrete velocities following the total
velocity distribution function such that

$f_{\textrm{t}ot}(\bf v_{g})=\sum_{i}N_{si}\delta^{3}(\bf v_{g}-\bf v_{si}),$
(2)

where $\delta^{3}(\bf v_{g}-\bf v_{si})$ is the Dirac-delta function, $\bf
v_{si}$ is the velocity of i-th stream, and $N_{si}$ is the number of
particles in that particular stream. The Dirac-delta function is implemented
as a narrow Gaussian with a width representing the dispersion of the stream.
The velocity of a DM particle unperturbed by the solar system coming from far
away with respect to the Sun is given by

${\textbf{v}_{\infty}}={\textbf{v}_{g}-\textbf{v}_{\odot}},$ (3)

where ${\textbf{v}_{g}}$ represent the velocity of the dark matter with
respect to the galaxy and

${\textbf{v}_{\odot}=\textbf{v}_{pec}+\textbf{v}_{\textrm{L}SR}},$ (4)

represent the velocity of the Sun with respect to the Galaxy, with peculiar
velocity ${\bf v_{\textrm{p}ec}}=(11,12,7)$ km/s, and velocity of standard at
rest ${\bf v_{\textrm{L}SR}}=(0,230,0)$ km/s chare:2014 . Then, applying
conservation of energy and Laplace-Runge-Lenz vector to account for the
gravitational effects by the Sun alenazi:2006 . This relates
${\textbf{v}_{\infty}}$ with velocity of DM at Earth’s location v in a Sun-
centered frame of reference

${\textbf{v}_{\infty}}=\frac{\text{v}_{\infty}^{2}{\textbf{v}}+{\text{v}_{\infty}}(GM/r){\hat{\bf{r}}}-{\text{v}_{\infty}}{\textbf{v}}\left({\textbf{v}}\cdot{\hat{\textbf{r}}}\right)}{{\text{v}_{\infty}^{2}}+\left(GM/r\right)-{\text{v}_{\infty}}\left(\textbf{v}\cdot{\hat{\textbf{r}}}\right)},$
(5)

where

$v_{\infty}=\sqrt{v^{2}-\frac{2GM}{r}}$ (6)

relates the speed at far away with the speed at the position of the Earth with
respect to the Sun. We follow a free fall trajectory towards the Sun until the
Earth’s position. Finally to find the velocity of dark matter in the
laboratory frame we use:

${\textbf{v}_{lab}=\textbf{v}-\textbf{V}_{\oplus}},$ (7)

where ${\textbf{V}_{\oplus}}$ is the velocity vector of the Earth with respect
to the Sun. We used scipy.optimize.fsolve to solve for v from Eq. 5 given
${\textbf{v}_{\infty}}$ from Eq. 3. We show in Fig. 1 a map of the arrival
direction of simulated particles during a period of one year (2022) in
equatorial coordinates where Cygnus constellation is shown for reference. For
positioning the Earth with respect to the Sun in the galaxy, we used a code
developed by C. O’Hare chare:2014 . We observe the typical pattern as
discussed by other authors chare:2014 ; chare:2018 ; chare:2020 ; mayet:2016
where a higher DM density is expected to arrive from the direction of Cygnus
due to the boost from lab velocity around the galaxy (The solar system is
moving towards the Cygnus constellation). Figure 2 shows a distribution of
velocity of particles in galactic frame and laboratory frame for the SHM and
SHM + dd32.111When we refer to SHM and SHM + dd32, we imply discreet large
number of streams ($10^{12}$) whose velocities are sampled from Eq. 1. The
mean value of the densities of the various overlapping streams reproduces the
local DM density $\rho_{0}$. These models in that respect are different from
smooth isothermal SHM.  As expected, the mean velocity is higher in the lab
compared to the galactic frame, and the SHM + dd32 has higher percentage of
particles in the low velocity tail compared to the SHM.

Figure 2: Galactic and laboratory velocity distributions for two different
simulation models. Top is only the SHM and bottom is the SHM+dd32. A shift
towards higher speed is noticed between galactic and lab frame as expected.
For comparison, the SHM+dd32 has a slightly more pronounced low speed trail
than SHM. The dd32 notation stands for a dark disk contribution of 32%.

After the free fall in Sun’s gravitational field, all the way until 1 A.U.
from the Sun, the streams reach the Earth where they also experience
gravitational self-focusing by the interior mass distribution of the Earth as
shown in Fig. 3. Using Newtons second law the trajectories are reconstructed
using python scipy.odeint package.

${\bf v_{i}}=\frac{d\bf r_{i}}{dt}$ (8) $\frac{d{\bf
v_{i}}}{dt}=\begin{cases}-G\frac{M_{E}}{r_{i}^{3}}{\bf
r_{i}}&~{}~{}~{}~{}r_{i}\geq R_{E}\\\\[5.0pt]
-G\frac{4\pi\eta_{E}(r_{i})}{3}{\bf r_{i}}&~{}~{}~{}~{}r_{i}<R_{E},\\\
\end{cases}$ (9)

where ${\bf r_{i}}$ represents the position of the $i^{th}$ DM particle and
$\eta_{E}(r_{i})$ is the local Earth density. The Earth is located at the
origin of the reference frame. For Earth’s inner density an approximate PREM
model prem:1981  is used, with densities $12000~{}kg/m^{3}$,
$7000~{}kg/m^{3}$, and $5000~{}kg/m^{3}$, for inner distances from 0-3400 km,
3400-5700 km, and from 5700-6340 km, respectively. Throughout this work a
possible self-interaction of DM particles is ignored.

In our approach so far we used Eq. 5 to solve the velocity distribution of
particles affected by the gravity of the Sun. Particles from the same stream
are taken to have the same direction. Particles then from same stream at some
distance ( we chose 30 $R_{E}$ as sufficiently far from the Earth) away from
the Earth are injected in the Earth gravitational field using Eq. 9 and
density features are studied afterwards (see sections III and IV).

As a cross-check we also simulate the trajectory of DM particles in the
gravitational field of the Sun and the Earth by solving namely the Eq. 10
which describes DM particles experiencing both gravitational fields. As
expected the enhancement factors are not significantly different from when
using the combined approach of Eq. 5 and 9. The simulations using Eq. 10 took
a long computer time since the majority of simulated particles ended up
outside of the Earth due to its’ small aperture compared to solar system size.
The results we quote in the following are based on combination of computation
based on Eq. 5 and simulation based on Eq. 9. The Newton’s second law equation
guiding DM particles trajectories in the gravitational field of the Sun and
the Earth is given by the following:

${\bf\frac{d{\bf
v_{i}}}{dt}=\left\\{\begin{array}[]{ll}-\frac{GM_{\odot}}{(r_{i}-r_{S})^{3}}(\bf
r_{i}-r_{S})-\frac{GM_{E}}{(r_{i}-r_{E})^{3}}(\bf r_{i}-r_{E})&\mbox{outside
Sun and Earth},\\\\[5.0pt] -\frac{4\pi}{3}G\eta_{S}(\bf
r_{i}-r_{S})-\frac{GM_{E}}{(r_{i}-r_{E})^{3}}(\bf r_{i}-r_{E})&\mbox{inside
Sun, outside Earth},\\\\[5.0pt] -\frac{4\pi}{3}G\eta_{E}(\bf
r_{i}-r_{E})-\frac{GM_{\odot}}{(r_{i}-r_{S})^{3}}(\bf
r_{i}-r_{S})&\mbox{inside Earth, outside Sun},\\\ \end{array}\right\\}}$ (10)

where $r_{E}$, $r_{S}$ are position for the Earth and the Sun, and $\eta_{S}$,
$\eta_{E}$ are the densities of the Sun and the Earth.

##  III Gravitational Focusing

After carrying out the Sun’s gravitational effects using Eq. 5, some DM
particles move towards the Earth; their trajectories are derived using Eq. 9.
Figure 3 shows the XZ and YZ view of the trajectories for some particles from
one stream as they experience the gravitational self-focusing by the inner
Earth mass distribution prem:1981 . A caustic-like GPDE region is formed at
the focus.

Figure 3: A stream experiences gravitational effects by the Sun and then a
gravitational self-focusing by the inner Earth; for a stream velocity of 20
km/s in laboratory frame incident from negative x. The self-focusing brings
the particles some distance away from the Earth. The top plot shows the XY
projection while the bottom plot shows the YZ projection at the focus. The
number of particles used in top and bottom plot is different.

The flow rate of particles $\rho Sv$ is constant, assuming that the velocity
(v) does not change appreciably while propagating from one side of the Earth
to the opposite one, the amplification factor is given by

$A\approx\frac{S_{0}}{S_{f}}=\frac{\pi r_{0}^{2}}{\pi r_{\text{f}}^{2}},$ (11)

where $S_{0}$, $S_{f}$ are the areas subtended by particles before and after
gravitational self-focusing at the focal point, respectively.

Figure 4 shows the amplification factor as a function of the distance from
geo-center to the ”root” of the self-focused stream (”hair”) by the Earth for
different stream velocities. The six streams in Fig. 4 are all assumed to be
dispersionless, while the amplifications peak at about $10^{9}$ and are nearly
independent on the velocity. However, the root of the focus distance is
velocity dependent (as in Refs. Prezeau:2015 ; sofue:2020 ).

Figure 4: Amplification as a function of geo-center distance in Earth radii
for different velocities. The streams in this case are all assumed
dispersionless. The amplification peaks at about $10^{9}$ with a peak height
nearly independent from velocity, the location of the focus point (”root”) is
however strongly dependent on velocity.

The DM density within a GPDE region varies with the largest value being at the
focal point and decreases away from the focal point as shown in Fig. 3 and 4.

For simplicity reasons, we define GPDE regions by an average length, an
average width, and an average amplification. A schematic of definitions used
is apparent from Fig. 5 marking an average width, length and how far away
these GPDE regions start from the geo-center. In our discussion we choose
amplification averages from 10 to $10^{8}$, with the implication that a region
with small average amplification is larger in size and a region with high
amplification is smaller in size due to dilution of densities away from focus
region. Table 2 summarizes the sizes of GPDE regions when choosing different
amplification averages.

Figure 5: Schematics of GPDE regions for two different streams at different incident velocity and incident directions. The schematics illustrates how the arrival direction of particles from one stream forms the GPDE region on the opposite side of the Earth. Table 2: Derived parameters for GPDE regions for different amplifications. Amplification Average | 10 | 100 | 1000 | 1e4 | 1e6 | 1e7 | 1e8
---|---|---|---|---|---|---|---
Length (L) ($R_{E}$) | 20.3 | 9.3 | 3.0 | 1.0 | 0.76 | 0.4 | 0.2
Width Radius ($r_{ave}$) (km) | 47 | 46 | 32 | 13 | 9.5 | 5 | 2.5

Figure 6 describes the amplification as a function of the geo-distance for two
stream velocities (50 km/s and 200 km/s). Particles in this case are injected
with different velocity dispersion values sampled from Gaussian functions. The
dispersion values chosen are $\delta v=10^{-17}$c, and $\delta v=10^{-10}$c
corresponding to primordial fine grained streams of axions and WIMPs,
respectively. As a cross-check we simulated also the case of $\delta v\sim
0.00053$ c corresponding to Maxwellian velocity dispersion, which is the
velocity dispersion for the smooth SHM. Dispersion is negligible for
primordial axion streams with amplification factor reaching up to $10^{9}$
just as in the dispersionless case. The WIMP streams have amplifications
factors of about a factor of 10 lower than axion streams due to their slightly
larger velocity dispersion. For the quite wide Maxwellian velocity profile as
expected the gravitational focusing effects become negligible.

Figure 6: The amplification for different distance from the Earth is shown for
particles with mean velocity of 50 km/s (top) and for particles with mean
velocity of 200 km/s (bottom). Three values of dispersion are chosen for each
plot to include primordial axions and WIMPs, and Maxwellian velocity
dispersion. Note that for the Maxwellian (green) there are cases where
particles are diverging at that particular point (Amplification less than 1).

At a particular location a large number of streams contribute to the local
nominal DM density, and when discussing GF effects, each one of the multiple
streams is taken into account at that location. The density of DM before
gravitational focusing is defined as

$\rho_{0}=\sum_{i}\rho_{si},$ (12)

where $\rho_{si}$ represents the density contribution to overall dark matter
from the i-th stream. Once gravitational focusing effects are calculated, the
total density is represented by

$\rho=\sum_{i}A_{i}\rho_{si},$ (13)

where $A_{i}$ represents the amplification factor which could be smaller or
larger than 1. After GF, to calculate the density from the many streams at
some arbitrary location, all these amplification coefficients are obtained
from simulation of multiple overlapping streams at that region. In most
locations where there is no dominant term from Eq. 13, the contribution of so
many streams should average to $\rho\sim\rho_{0}$. If there is a dominant term
however, Eq. 13 becomes

$\rho=A_{j}\rho_{sj}+\sum_{i\neq j}A_{i}\rho_{si},$ (14)

where the j-th stream has been singled out from the rest. If there are regions
where the amplification from the j-th stream is much larger than
amplifications from all the other streams $A_{j}\gg A_{i}$ and the density
after amplification from one stream dominates over the others
$A_{j}\rho_{sj}\gg\rho_{0}$ than instead the second term in 14 is smaller than
the first and Eq. 14 becomes

$\rho\sim A_{j}\rho_{sj}$ (15)

In these regions the simple expression given by Eq. 15 is sufficient even in a
multi stream environment.

##  IV Results & Discussion

For DM detection it is important to know how many streams are expected on
Earth’s location and how many of them form GPDE regions close to the Earth,
how large the expected density enhancements are compared to the average local
DM density, and how long they last. Firstly, the number of streams with a
particular density ratio $\rho_{s}/\rho_{0}$ at Earth’s position is taken from
Ref. Vogelsberger:2011 . Table 3 gives a summary of the abundance of fine
grained streams at the solar system and their probability. The number of
streams is calculated following the approximate relation:

$N_{s}\cdot\frac{\rho_{s}}{\rho_{0}}\sim F(>\rho_{s})N_{ts},$ (16)

where $F(>\rho_{s})$ is the fraction of streams with densities above
$\rho_{s}$, and $N_{ts}$ is the total number of streams Vogelsberger:2011 .
For instance, following Table 3 a stream of density $\rho_{s}=0.01\rho_{0}$
has a 20% chance of being present in our vicinity, and the stream number
becomes $2\cdot 10^{4}$ and $2\cdot 10^{6}$ for a density
$\rho_{s}=10^{-5}\rho_{0}$, and $\rho_{s}=10^{-7}\rho_{0}$, respectively.

Table 3: Streams with different densities and their probability in solar neighbourhood. Density ratio ($\rho_{s}/\rho_{0}$) | Number of streams | Probability (%)
---|---|---
1 | 1 | 0.002
0.1 | 1 | 0.2
0.01 | 1 | 20
10-3 | 10 | 100
10-4 | 500 | 100
10-5 | $2\cdot 10^{4}$ | 100
10-6 | $4\cdot 10^{5}$ | 100
10-7 | $2\cdot 10^{6}$ | 100

It is also of interest to know whether the GPDE regions produced by streams
are isotropic or have preferred directions. Figure 7 shows the equatorial
coordinates map for GPDE regions, which form on the opposite side to the
incoming stream directions. Figure 7 which includes all possible stream
velocities shows more GPDE regions in the direction opposite Cygnus (on the
right below the ecliptic). This is due to the kinematic boost where larger
number of particles arrive from Cygnus direction, consequently more streams
are sampled to arrive from this direction, which leads to more GPDE regions in
the opposite side. Figure 8 represents the same maps but limiting the
velocities to less than 50km/s, implying streams whose GPDE regions are near
Earth’s surface. The distributions in this case are approximately isotropic.
The Sun causes a large bending effect on low velocity streams strongly
deviating their directions compared to the unperturbed directions. For near
the Earth experiments we conclude that GPDE regions are approximately
isotropic for both models (for velocities $\leq$ 50 km/s).

Figure 7: Equatorial coordinates map distribution for GPDEs of streams of all
possible velocities, for SHM and SHM + dd32. This is contrary to Figure 1
which shows the incoming direction of streams. The GPDE regions are on the
opposite side of the incoming direction as illustrated in Figure 5. A larger
density of GPDE regions is observed in the direction opposite Cygnus which is
consistent with more streams arriving from Cygnus direction. For reference,
the Cygnus is depicted with red stars while the ecliplic with orange circles.
Figure 8:  Equatorial coordinates map distribution for GPDEs of streams with
velocities less than 50 km/s, for the SHM and the SHM + dd32. The GPDEs are on
the opposite side of the incoming direction of streams. The Cygnus is depicted
with red stars while the ecliplic with orange circles. The distribution in
both models is approximately isotropic for low velocities.

Next question is how many streams have their GPDE regions near the Earth’s
surface as being relevant for Earth bound detectors; therefore two distances
from geo-center are chosen ($\sim 1R_{E}$) or in orbit ($\sim 1.5R_{E}$).
Table 2 refers to different GPDE sizes for amplification range between 10 and
$10^{8}$, a small amplification corresponds to a large GPDE size and vice
versa. We define a cutoff velocity as the upper value of a stream velocity
such that its particles fall just inside the GPDE region. Larger size regions
are expected to have larger value for the cutoff velocity. Table 4 gives the
cutoff velocities for two near Earth observation distances from geo-center
($\sim 1R_{E}$ and $\sim 1.5R_{E}$ ) for different amplifications. The cutoff
velocities vary from about $v_{\text{cut}}\leq 75$ km/s to $v_{\text{cut}}\leq
10$ km/s for amplification factors from 10 to $10^{8}$ given the observation
region is at ($1R_{E}$). The cutoff velocities for the observation region
($1.5R_{E}$) vary from $v_{\text{cut}}\leq 78$ km/s to $v_{\text{cut}}\leq 14$
km/s.

Table 4: Cutoff velocities for GPDE regions for different amplifications and different observation distances. Amplification Average | 10 | 100 | 1000 | 1e4 | 1e6 | 1e7 | 1e8
---|---|---|---|---|---|---|---
$v_{\text{cut}}\leq$ (km/s) | 75 | 51 | 29 | 18 | 16 | 13 | 10
($d_{\text{obs}}\simeq 1R_{E}$ ) |  |  |  |  |  |  |
$v_{\text{cut}}\leq$ (km/s) | 78 | 52 | 32 | 22 | 19 | 17 | 14
($d_{\text{obs}}\simeq 1.5R_{E}$ ) |  |  |  |  |  |  |

Figure 9 and 10 give the number of streams whose GPDE regions are at $\sim
1R_{E}$ and $\sim 1.5R_{E}$ from geo-center. Thus, each line represents a
fixed value for the density taken from Table 3. The x-axis represents the
density ratios of the GPDE region density as derived by equation 15. Density
ratio for WIMPs and axions is given by the x-axis (top) and (bottom)
respectively. 222We use the expression from equation 15 to calculate the
density ratio. We emphasise that this expression works in a multi-stream
environment as long as $A_{j}\gg A_{i}$, and $A_{j}\rho_{sj}\gg\rho_{0}$ which
would be fine for density ratios of greater than about 10. Above 10 or so the
graphs in Figs. 9, 10 are a better representation of the multi-stream
environment.  Note that for the same incoming stream density we expect
different GPDE density ratios for the fact that we have different average
amplification assumptions as in the following Tables 2 and 4. Also the minimum
(10) and the maximum ($10^{8}$) amplification factors, limit the range of
density ratios differently for each incoming stream density assumption, as
observed in ranges of each line in Figs. 9 and 10.

Figure 9: The number of streams expected at the Earth’s location, where their
GPDE regions overlap at the Earth’s surface ($\sim 1R_{E}$). The top plot
refers to SHM while the bottom plot is for the SHM+dd32 model. The different
lines represent different assumptions for stream densities taken from Table 3.
The variations within one line are due to different amplification values (
from 10 to $10^{8}$) (Table 4). The x-axis represents the ratio of the GPDE
density over the nominal mean DM density which is given by
$\rho_{s}A/\rho_{0}=\rho/\rho_{0}$. The amplification factor applied to axions
is used in the bottom x-axis and for WIMPs at the top x-axis. As an example,
out of $4\cdot 10^{5}$ axion streams with density of
$\rho_{s}=10^{-6}\rho_{0}$ before GF, only 1.4 of them result to a GPDE region
with density ratio $\rho/\rho_{0}=100$ with an observation distance of $\sim
1R_{E}$ in SHM. The afore-mentioned number of streams becomes 6 for the
SHM+dd32 Figure 10: Number of streams expected on Earth where the GPDE regions
are at $\sim 1.5R_{E}$. The top plot refers to SHM while the bottom plot is
for the SHM+dd32 model. The different lines represent different assumptions
for the densities of incoming streams taken from table 3. The changes within
one line are due to different amplification assumptions as described in Table
4. The x-axis represents the ratio of the GPDE density over nominal mean DM
density. The amplification factor applied to axions is used in the bottom
x-axis and for WIMPs at the top x-axis. As an example, out of $4\cdot 10^{5}$
axion streams with density $\rho_{s}=10^{-6}\rho_{0}$ before GF, only about 7
of them will have a GPDE region with density ratio $100$ with an observation
distance of $\sim 1.5R_{E}$ in SHM. This value for similar conditions is about
20 for the SHM+dd32. Such plots are relevant for detectors in space.

As an example, we analyze $\rho_{s}=10^{-6}\rho_{0}$ line on Fig. 9
corresponding to the same incoming stream density ($\rho_{s}=10^{-6}\rho_{0}$)
in SHM. In this case, there are about $4\cdot 10^{5}$ streams passing through
the Earth continuously following Table 3. Interestingly, out of all these
streams only about 1.4 on average have GPDE regions at an observation distance
of $\sim 1R_{E}$, with an average amplification of $10^{8}$. This number
increases to 2 requiring an average amplification of $10^{7}$. The density
ratios of these regions are then given by Eq. 15, and, they are $100$ and $10$
for axions, and, $10$ and $1$ for WIMPs, respectively. For the SHM+dd32 shown
in Fig. 9 these stream numbers increase and become 6 and 10, respectively. For
the larger observation distance ($\sim 1.5R_{E}$), 7 GPDE regions in SHM, and
about 20 in SHM+dd32, have axion density of $100\rho_{0}$ (see Fig. 10).

Figures 9 and 10 show all possible GPDE regions with particular stream
densities for a full sky coverage. Since, no experiment alone can have full
sky coverage, while the experiment’s field-of-view moves as the Earth rotates
and evolves. Experiments need to know the likelihood of detecting a GPDE event
per day given streams with a particular density, and how long these events
last. The probability of encountering one GPDE region per day is given by

$P_{GPDE}/day\sim N_{\text{GPDE}}N_{\text{exp}}P_{\text{loc}},$ (17)

where $N_{\text{GPDE}}$ is the number of GPDE regions being present at the
observation distances which are obtained from Figs. 9 and 10. $N_{\text{exp}}$
is the number of experiments, and $P_{\text{loc}}$ is the distribution
probability of GPDE regions at a particular location on Earth which could be
position dependent on experiment’s geographical location and time dependence
due to Earth’s orbital motion around the Sun and its rotation. From Fig. 8 it
follows that, for the case of small velocities (less than 50 km/s), the
expected distribution of GPDE regions during the entire year (2022) is
approximately isotropic. By assuming a uniform distribution map (as in Fig.
8), the location probability is given as

$P_{\text{loc}}\sim\frac{2\pi d_{\text{obs}}cos(l)(2r_{ave})}{4\pi
d^{2}_{\text{obs}}},$ (18)

where $r_{\text{ave}}$ is the average radius of the cylinder shaped GPDE
region for a particular amplification factor as given by Table 2, and $l$ is
the latitude of the detector. The numerator represents the area covered by the
GPDE per day, and the denominator represents the corresponding area for full
coverage. Figures 11 and 12 give the probability to encounter a GPDE region
per day, for one single experiment being in two different observation
distances. These plots assume a mid latitude observer location at about $l\sim
40^{o}$

Figure 11: The probability to encounter a GPDE region on Earth per day for one
single experiment at an observation distance of $\sim 1R_{E}$ as a function of
relative density over the nominal halo for different density enhancement
values.  Figure 12: The probability to encounter a GPDE region per day on
Earth, per one single experiment at an observation distance of $\sim 1.5R_{E}$
as a function of relative density over the nominal halo for different density
enhancement values. Such plots are relevant for detectors in space.

The probability to encounter one GPDE region per day, for one experiment at
$1R_{E}$, and a density ratio $10$ for axions or $1$ for WIMPs is $\sim$0.005
for SHM and $\sim$0.02 for SHM+dd32. The same probability for one experiment
at $\sim 1R_{E}$, and a density ratio $100$ for axions or $10$ for WIMPs is
$\sim 4\cdot 10^{-4}$ for SHM, and $\sim 2\cdot 10^{-3}$ for SHM+dd32. The
corresponding values for one experiment at $\sim 1.5R_{E}$, and a density
ratio $100$ for axions or $10$ for WIMPs, becomes $1.5\cdot 10^{-3}$ for SHM
and $4\cdot 10^{-3}$ for SHM+dd32.

As the Earth rotates so does also the detector, and therefore the GPDE regions
can give rise to transient signals. A detector at the surface of the Earth
moves at a speed of $\sim 0.5{\text{km/s}}\cdot cos(l)$. The transient time
for an experiment that is moving through a GPDE region at an observation
distance $d_{\text{obs}}$ is given by

$t_{enc}\sim\frac{2\cdot r_{ave}}{(d_{\text{obs}}/R_{E})\cdot 0.5cos(l)}$ (19)

For an $r_{\text{ave}}=2.5\text{km}$, the transient time is $\sim$13 seconds
for $1R_{E}$, and $\sim$9 seconds for $1.5R_{E}$.

Apparently, a network of many experiments seems promising, in particular as
long as the mass of the DM particles and their interaction strength or their
stream configuration are not known. For instance the GNOME network gnome:2021
, the electric dipole moments (EDM), storage ring experiments yannis:2018 ;
yannis:2021 , and the ECHO idea ariel:2019 , etc., fit the reasoning of the
present work.

Another possible implication of the remarkable gravitational self-focusing
effect of DM constituents by the inner Earth’s mass distribution, results to
possible overlap with, say, conventional cosmic ray studies. Because, the
appearance of caustics-like shapes downstream of an incident low speed stream
might correlate with upward escaping cosmic-ray like events from the Earth’s
surface. In this context it is intriguing to note the observed anomalous
events by the Antarctic Impulse Transient Antenna (ANITA) Anita:2016 ;
Anita:2018 ; Anita:2020  that appear to be like energetic cosmic showers
emerging from the Earth’s surface. Interestingly, recently Liang and
Zhitnitsky advocate as a possible explanation for these apparently anomalous
events within the axion quark nugget (AQN) model Ariel_Zh:2022 . Of note, the
AQN model was invented in 2003 Ariel_Zh:2003 , and solves the DM problem and
the charge asymmetry in the Universe.

Finally, we also wish to stress that the present first simulation results with
a streaming DM component applies in principle to DM searches for any other DM
constituent to be invented. Also, should a DM particle be identified, one can
use the same study presented here to eventually derive the fine structure of
the DM distribution in our galaxy going eventually back to cosmological times.
Our results show that transient GPDE regions of high density around the Earth
(dubbed “hairs”) would form only if the streaming nature of DM is true,
compared to the smooth SHM configuration. The simulation tools developed in
this work can be adapted to any kind of dark matter particle, since they all
couple gravitationally. 333The simulation code can be provided to interested
researchers upon request.

###  IV.1 Numerical Examples

In this subsection we mention a couple figures of merit from this section that
one can use to initiate a promising search for GPDE regions.

At the surface of the Earth, there are present between 1 and 6 GPDE regions
with radius of $\sim 2.5$ km and length of $0.2R_{E}$, and with densities
$100\rho_{0}$ for axions, or $10\rho_{0}$ for WIMPs. A single experiment at
the surface of the Earth has a probability per day of anywhere between $4\cdot
10^{-4}$ to $2\cdot 10^{-3}$ to run into such a region as the Earth rotates
and the experiment too. The transit time is determined from the size of the
GPDE region and the Earth’s rotation velocity and it is about 13 seconds.

At a distance of $\sim 1.5R_{E}$ there are anywhere between 7 and 20 GPDE
regions with radius of $\sim 2.5$ km and length of $0.2R_{E}$, and with
densities $100\rho_{0}$ for axions, or $10\rho_{0}$ for WIMPs. A single
experiment at this location has a probability per day of anywhere between
$1.5\cdot 10^{-3}$ to $4\cdot 10^{-3}$ to enter such a region as the Earth
rotates. The transit time is about 9 seconds. 444More precisely, the transit
time for this case depends on relative speed of satellite with respect to the
Earth.

##  V Summary & Conclusion

This simulation work addresses the question as to whether the invisible
Universe being also in form of streams offers a better DM detection
sensitivity than the direct DM searches, which have targeted so far the DM
within the isotropic SHM. The performed simulations show certain detection
advantages for the streaming DM configuration compared with the SHM. For
example:

  * •

There are enhanced DM fluxes following planetary gravitational lensing
effects, resulting to caustics-like configurations in particular when the
gravitational self-focusing effects by the inner Earth are at work.

  * •

The velocity selection of streaming DM via the focal plane position. E.g., at
the Earth’s surface the speed is about 12 km/s, 10 km/s, and 14 km/s at
$1R_{E}$, $0.9R_{E}$, and $1.5R_{E}$, respectively.

  * •

Transient events due to Earth’s rotation. The expected signatures last about
10 seconds and could become an additional signature in the search of streaming
DM.

  * •

A potential new method to detect and identify DM tidal or cosmological streams
with Earth bound experiments like the ongoing or future DM axion searches.

  * •

Search for similar planetary relationships in exo-solar systems, but also
correlations between our solar system and other nearby exo-solar systems
perrymankonstantin:2021 . This indicates a far going perspective, once data of
long time series from near exo-planetary systems becomes available.

The main advantage to expect from a dark sector being at least partially in
form of streams comes from the fact that with such a large number of streams
in our vicinity (up to $10^{12}$), some of them experience occasionally strong
gravitational flux enhancement by the various solar system bodies including in
particular the self-focusing by the inner Earth. When some planet or the Sun
intervene along a stream incident towards the Earth, it enhances temporally
downstream the DM flux, improving the detection sensitivity. Of note, if some
stream(s) last relatively long, say few months or more, following the
aforementioned self-focusing effects by the inner Earth, an experiment can
accumulate many temporally enhanced fluxes, accumulating very long exposure
times unthinkable to reach with conventional searches targeting DM within SHM.
In addition, the expected transient events can become the signatures to
unravel the existence of invisible streams and their structure.

All these are new possibilities and an important by-product from the
undertaken simulation in this work.

###### Acknowledgements.

We sincerely thank Mark Vogelsberger for the numerous feedback about setting
up simulation inputs for the fine grained streams. Yannis Semertzidis is
acknowledged for encouraging discussions about DM searches. AK acknowledges
Messiah University Scholarship Programs for support of the work. KZ thanks
Antonis Gardikiotis for his help.

## References

  * (1) M. Vogelsberger, F. Marinacci, P. Torrey and E. Puchwein, _Cosmological simulations of galaxy formation_ , _Nature Reviews Physics_ 2 (2020) 42.
  * (2) J. Stucker, O. Hahn, R.E. Angulo and S.D. White, _Simulating the Complexity of the Dark Matter Sheet I Numerical Algorithms_ , _MNRAS_ 295 (2020) 4943.
  * (3) V. Springel, R. Pakmor, O. Zier and M. Reinecke, _Simulating cosmic structure formation with the gadget-4 code_ , _MNRAS_ 506 (2021) 2871.
  * (4) M. Vogelsberger, A. Helmi, V. Springel, S.D.M. White, J. Wang, C.S. Frenk et al., _Phase-space structure in the local dark matter distribution and its signature in direct detection experiments_ , _MNRAS_ 395 (2009) 797.
  * (5) M. Vogelsberger and S. White, _Streams and caustics: the fine-grained structure of $\Lambda$ cold dark matter haloes_, _MNRAS_ 413 (2011) 1419.
  * (6) M. Vogelsberger, S. White, A. Helmi and V. Springel, _The fine-grained phase-space structure of cold dark matter haloes_ , _MNRAS_ 385 (2008) 236.
  * (7) P. Sikivie, I. Tkachev and Y. Wang, _Velocity peaks in the cold dark matter spectrum on earth_ , _Phys. Rev. Lett._ 75 (1995) 2911.
  * (8) P. Sikivie, _Velocity peaks and caustic rings_ , _Nuclear Physics B-Proceedings Supplements_ 72 (1999) 110.
  * (9) A. Natarajan and P. Sikivie, _A further look at particle annihilation in dark matter caustics_ , _Phys. Rev. D_ 77 (2008) 043531.
  * (10) D. Hoffmann, J. Jacoby and K. Zioutas, _Gravitational lensing by the Sun of non-relativistic penetrating particles_ , _Astropart. Phys._ 20 (2003) 73.
  * (11) M.S. Alenazi and P. Gondolo, _Phase-space distribution of unbound dark matter near the Sun_ , _Phys. Rev. D_ 74 (2006) 083518.
  * (12) B.R. Patla, R.J. Nemiroff, D.H.H. Hoffmann and K. Zioutas, _Flux Enhancement of Slow-moving Particles by Sun or Jupiter: Can they be Detected on Earth?_ , _ApJ_ 780 (2014) 158.
  * (13) G. Prezeau, _Dense Dark Matter Hairs Spreading Out from Earth, Jupiter and Other Compact Bodies_ , _ApJ_ 814 (2015) 122.
  * (14) Y. Sofue, _Gravitational Focusing of Low-Velocity Dark Matter on the Earth’s Surface_ , _Galaxies_ 8 (2020) 42.
  * (15) S. Lee, S. Ahn, J. Choi, B. Ko and Y. Semertzidis, _Axion Dark Matter Search around 6.7 $\mu$eV_, _Phys. Rev. Lett._ 124 (2020) 101802.
  * (16) ADMX Collaboration, N. Du, N. Force, R. Khatiwada, E. Lentz, R. Ottens et al., _A Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 120 (2018) 151301.
  * (17) CAST Collaboration, V. Anastassopoulos, S. Aune, K. Barth, A. Belov, H. Bräuninger et al., _New CAST limit on the axion–photon interaction_ , _Nature Physics_ 13 (2017) 584.
  * (18) S.K. Lee, M. Lisanti, A.H.G. Peter and B.R. Safdi, _Effect of Gravitational Focusing on Annual Modulation in Dark-Matter Direct-Detection Experiments_ , _Phys. Rev. Lett._ 112 (2014) 011301.
  * (19) C. Mateu, A.P. Cooper, A.S. Font, L. Aguilar, C. Frenk, S. Cole et al., _Predictions for the detection of tidal streams with Gaia using great-circle methods_ , _MNRAS_ 469 (2017) 721.
  * (20) C.W. Purcell, J.S. Bullock and M. Kaplinghat, _The dark disk of the Milky Way_ , _ApJ_ 703 (2009) 2275.
  * (21) C. O’Hare and A. Green, _Directional detection of dark matter streams_ , _Phys. Rev. D_ 90 (2014) 123511.
  * (22) C. O’Hare, C. McCabe, N. Evans, G. Myeong and V. Belokurov, _Dark matter hurricane: Measuring the S1 stream with dark matter detectors_ , _Phys. Rev. D_ 98 (2018) 103006.
  * (23) C. McCabe, C. O’Hare, W. Evans, G. Myeong and V. Belokurov, _Velocity substructure from Gaia and direct searches for dark matter: Dark Shards_ , _Phys. Rev. D_ 101 (2020) 023006.
  * (24) F. Mayet, A.M. Green, J.B.R. Battat, J. Billard, N. Bozorgnia, G.B. Gelmini et al., _A review of the discovery reach of directional Dark Matter detection_ , _Phys. Rept._ 627 (2016) 1.
  * (25) A.M. Dziewonskia and D.L. Anderson, _Preliminary reference Earth model_ , _Physics of the Earth and Planetary Interiors_ 25 (1981) 297.
  * (26) GNOME Collaboration, S. Afach, B.C. Buchler, D. Budker, C. Dailey, A. Derevianko et al., _Search for topological defect dark matter with a global network of optical magnetometers_ , _Nature Physics_ 17 (2021) 1396.
  * (27) S. Haciomeroglu and Y.K. Semertzidis, _A hybrid ring design in the storage-ring proton electric dipole moment experiment_ , _Phys. Rev. Accel. Beams_ 22 (2019) 034001.
  * (28) P.W. Graham, S. Haclu, D.E. Kaplan, Z. Omarov, S. Rajendran and Y.K. Semertzidis, _Storage ring probes of dark matter and dark energy_ , _Phys. Rev. D_ 103 (2021) 055010.
  * (29) A. Arza and P. Sikivie, _Production and detection of an axion dark matter echo_ , _Phys. Rev. Lett._ 123 (2019) 131804 [1902.00114].
  * (30) ANITA Collaboration, P.W. Gorham, J. Nam, A. Romero-Wolf, S. Hoover, P. Allison et al., _Characteristics of four upward-pointing cosmic-ray-like events observed with anita_ , _Phys. Rev. Lett._ 117 (2016) 071101.
  * (31) ANITA Collaboration, P.W. Gorham, B. Rotter, P. Allison, O. Banerjee, L. Batten et al., _Observation of an unusual upward-going cosmic-ray-like event in the third flight of anita_ , _Phys. Rev. Lett._ 121 (2018) 161102.
  * (32) ANITA Collaboration, P.W. Gorham, A. Ludwig, C. Deaconu, P. Cao, P. Allison et al., _Unusual near-horizon cosmic-ray-like events observed by anita-iv_ , _Phys. Rev. Lett._ 126 (2021) 071103.
  * (33) X. Liang and A. Zhitnitsky, _Anita anomalous events and axion quark nuggets_ , _Phys. Rev. D_ 106 (2022) 063022.
  * (34) A. Zhitnitsky, _‘Nonbaryonic’ dark matter as baryonic colour superconductor_ , _JCAP_ 10 (2003) 010.
  * (35) M. Perryman and K. Zioutas, _Gaia, Fundamental Physics, and Dark Matter_ , 2106.15408.