A Proposed Network to Detect Axion Quark Nugget Dark Matter

A network of synchronized detectors can increase the likelihood of discovering the QCD axion, within the Axion Quark Nugget (AQN) dark matter model. A similar network can also discriminate the X-rays emitted by the AQNs from the background signal. These networks can provide information on the directionality of the dark matter flux (if any), as well as its velocity distribution, and can therefore test the Standard Halo Model. We show that the optimal configuration to detect AQN-induced axions is a triangular network of stations 100 km apart. For X-rays, the optimal network is an array of tetrahedral units.


INTRODUCTION
The Standard Halo Model (SHM) is locally characterized by an average dark matter (DM) density ρ DM ≈ 0.3 GeVcm −3 and a velocity distribution insensitive to the environment, e.g. to the dynamics of the planets and the Sun (see [1]). Significant modifications of these local DM features would dramatically impact the entire field of direct DM searches. The SHM is consistent with the direct measurement of the average DM mass density from the latest GAIA data release, which finds ρ DM = 0.6 ± 0.4 GeVcm −3 [2]. However, this direct measurement of ρ DM remains imprecise, even with the best kinematics measurement of the nearby stars available to date. Moreover, this measurement cannot rule out large variations of the local DM density at the scale of the Solar System [3,4]. Furthermore, it has been suggested that the Sun and planets may play a key role in the local DM distribution [5][6][7][8] 1 . In addition, the Milky-Way merging history might locally alter the DM distribution, and even at the galactic scale, in nearby galaxies, some observed characteristics of DM are not satisfactorily explained by N-body simulations [9].
In this work, we propose a framework that can test some fundamental assumptions of the SHM, such as the magnitude of the DM density and the DM velocity distribution in the local environment. These assumptions cannot be tested with a single experiment, even a very large one. It requires a network of spatially distributed instruments. These detectors must be synchronized in order to study the directions of the DM particles hitting the Earth, and their velocity distribution.
In a recent proposal [10,11] it was suggested to use a broadband detection strategy to search for axions within * xunyul@phas.ubc.ca † e.peshkov@alumni.ubc.ca ‡ waerbeke@phas.ubc.ca § arz@phas.ubc.ca 1 In particular, [5] found strong correlations between the position of the planets and solar flare activities, and proposed that an "invisible matter stream" is responsible for them. Furthermore, [6] claimed that stratospheric temperature anomalies cannot be explained by conventional atmospheric phenomena and are correlated with planetary positions.
the context of the Axion Quark Nugget (AQN) dark matter model. The main problem with broadband detection is the difficulty to distinguish between the AQN signal and large random noise and spurious events, but the dominant random background noise can be eliminated with the following observational strategy: 1-separate the DM signal from the larger noise by explicitly looking for annual and daily modulations, which are specific to the former; 2-use synchronized detectors to discriminate between the DM signal and the larger noise, by looking at the time delays recorded by two or more nearby detectors, as these delays are unambiguously fixed by the distances between the detectors. Developing such a detection strategy is relevant to the search of any DM particles that leave a detectable trace, being neutrinos [12], infrasonic/ seismic 2 events [13], or X-rays. We must emphasize that we will discuss correlated events that are synchronized but not coherent, for example distinct X-ray photons emitted in different directions at slightly different instants. Being emitted by the same AQN, they can be thought as clustering events which do not satisfy a conventional Poisson distribution. In contrast, noise spurious events are random and follow a conventional Poisson distribution, which allows to discriminate the AQN-induced signal from much larger noise events, provided one collects sufficient statistics of true signals.
Other networks searching for axions, such as The Global Network of Optical Magnetometers for Exotic physics searches (GNOME) [14,15] differ from our proposal. First, their search window is up to the kHz range, while the preferred value for the AQN axion mass corresponds to a frequency of 24 GHz. Second, they search for rare coherent single events, while we study correlated clustering events on a much smaller scale, with a relatively high occurrence rate.
The topic of the present study is two new elements, which have not been discussed previously in the literature: 1. Three (or better four) neighbouring synchronized detectors allow to reconstruct the trajectory of an AQN, the source for the signals detected by all the stations. This allows to study the DM flux distribution, including its velocity distribution and directionality. 2. When an AQN propagates and annihilates in the Earth's atmosphere, it emits weakly coupled axions and neutrinos with very large mean free paths (as discussed in [11] and [12] respectively). The AQN also emits Xrays, with a mean free path between 5 and 50 meters at sea level, depending on the energy, which interact strongly with the Earth environment. The study of crosscorrelated X-ray signals can be more effective than current axion detectors 3 , because X-ray stations are easy to build and to operate, in comparison to more expensive and demanding axion detectors. In order to discriminate AQN-induced X-rays from background X-rays, one should arrange several synchronized detectors separated by less than 50 meters, and study the correlated X-ray signals emitted by one and the same passing AQN. 4 In the next Section, we overview the basic ideas of the AQN framework, section III is an overview of the known spectral features of the axion and X-ray emissions which play a key role in the present work. Our proposed network design is presented in Section IV and concluding remarks are discussed in Section V.

II. THE AQN FRAMEWORK
First, we start with an overview of a novel mechanism for axion production, within the AQN framework, as recently suggested in [10,[16][17][18]. This new mechanism always accompanies the two conventional and well established axion productions : by a misalignment mechanism ( when the cosmological field θ(t) oscillates and emits cold axions), or via the decay of topological objects 5 .
3 provided the detectors we advocate can be built and be sensitive to relevant frequency bands. This is not a trivial assumption : the axion haloscopes known to be working, at present, are based on cavity type detection, rather than on the broadband technology required by our study, see Section III below for details. 4 A portion of the AQN-induced X-rays will be transferred to sound waves, similar to conventional blasts resulting from bombs or large meteors, as recently argued in [13] : the material vaporizes, rapidly expands, which eventually contributes to a small scale shock-wave (of order 10 −3 cm). In our case, similar processes occur but at larger scales (of order 5 meters or so [13]). AQN-induced sound waves can propagate over large distances (of order 100 kilometers [13]), to contrast with X-rays with an attenuation length measured in meters. It means that one can study AQN-induced infrasounds and seismic correlated signals using synchronized stations. Such a network allows to study the DM distribution, its directionality and velocity distribution, but this is beyond the scope of this paper. 5 The AQN framework assumes the conventional Peccei-Quinn resolution of the strong CP problem with the axion mass assumed to be in the classical window ma ∈ (10 −6 −10 −3 ) eV, see original papers [19][20][21][22][23][24][25] and recent reviews [26][27][28][29][30][31][32][33][34]. The AQN framework also assumes that inflation occurs during or after the PQ phase The AQN construction is in several aspects similar to the original quark nugget model suggested by Witten [35] (see [36] for a review). It is a type of "cosmologically dark" DM, but not because nuggets are weakly interacting. AQN is a strongly-interacting DM candidate with a small cross-section-to-mass ratio, that reduces many observable consequences.
There are two additional elements in the AQN model compared to the old proposal [35,36]. First, there is an additional stabilization factor for the nuggets, provided by the axion domain walls (with the QCD substructure) which are copiously produced during the QCD transition. This extra pressure alleviates a number of problems with the original [35,36] nugget model. Another feature of AQNs is that during the QCD transition, nuggets can be made of matter as well as antimatter . The direct consequence is that the DM density, Ω DM , and the baryonic matter density, Ω visible , will automatically assume the same order of magnitude ( Ω DM ∼ Ω visible ) , irrespective of the parameters of the model such as the axion mass m a or misalignment angle θ 0 .
We refer to the original papers [37][38][39][40] for questions related to AQN formation, generation of the baryon asymmetry, and survival pattern of AQNs in the hostile environment of the early Universe. AQNs are characterized by a baryon charge B ∈ (10 23 − 10 28 ) and a mass M N m p |B|. We emphasize that the AQN model is consistent with all presently available cosmological, astrophysical, satellite and ground-based constraints. This model offers a manifold of explanations for seemingly unrelated modern puzzles. It alleviates tensions with observations, such as the observed excesses in galactic emission in various frequency bands. It may resolve some longstanding cosmological puzzles, such as the "Primordial Lithium Puzzle" [41] and the "The Solar Corona Mystery" [42,43]. It could explain "impulsive radio events", recently observed in the quiet solar corona by the Murchison Widefield Array ( [44]). It may also explain three unsolved puzzles: the annual modulation recorded by DAMA/LIBRA (at a 9.5σ confidence level) [12]; the Xray seasonal variations recorded by the XMM-Newton observatory (at a 11σ confidence level) [45] ; the mysterious bursts observed by the Telescope Array [46]. Finally, it could explain mysterious explosions, known as the "sky-quakes" [13], already mentioned above.
We must emphasize that all the above studies are based on one and the same set of AQN parameters, even though they propose explanations for a wide variety of phenomena, in dramatically different environments where parameter differ by many orders of magnitude. In this work, we use again exactly the same set of parameters.
For our estimates in the present work, we need to quote the AQN hitting rate, assuming a conventional dark mattransition, such that one and the same misalignment angle θ 0 occupies entire visible Universe. where B 10 25 is the average baryon charge for a given size distribution determined by the function f (B) ∝ B −α , with α ≈ (2 − 2.5) ( see [12] for details and references). The averaging over all types of AQNtrajectories, with different masses M N m p |B|, different incident angles, different initial velocities and a variety of models for size distribution, does not modify this estimate noticeably. The result (1) suggests that AQNs hit the Earth's surface with a frequency of approximately one a day per (100km) 2 area, which represent a typical "small" event. The hitting rate for large sized objects with B B is suppressed by the distribution function f (B) ∝ B −α , while small sized objects with B B represent more frequent but less intense events, according to the same distribution function f (B).
In present work we study typical events with B 10 25 . In the section IV, we will discuss how a network of synchronized detectors should be able to record the crosscorrelations from such a single trans-passing nugget, and could allow to analyze the directionality and velocity distributions of the AQNs. But first, in the next section III, we will overview the spectral features of the axion and related X-ray emissions.

III. SPECTRAL PROPERTIES OF THE AXION AND X-RAY EMISSIONS
In this section, we overview the spectral properties of the axion and X-ray emissions during the passage of an antimatter AQN through the Earth, as these properties constrain the design of the instruments and the synchronized network to be discussed in section IV.

A. Axion velocity distribution
We start our review with AQN-induced axions. We represent the time-dependent axion flux as follows [10]: where the amplitude A(t) is normalized to unity when averaged over long period of time, i.e. A(t) = 1. The background flux in Eq. (2) is determined for axions emitted from deep Earth's underground using Monte Carlo simulations (as performed in [10,18]) , for different size AQNs and different trajectories. The atmospheric portion of the AQN's path was ignored in these simulations, because the density of the Earth's atmosphere is three orders of magnitude lower than the density of the planet itself. This should be contrasted with the next subsection III B where the atmospheric portion of the path plays the dominant role. Here, the mean free path for X-rays in solids does not exceed 1 cm, while it is three orders of magnitude higher in the atmosphere. There is a number of time-dependent effects which influence A(t). First, there is a well-established annual modulation [1,47] due to the changing position and velocity of the Earth relative to the Sun. It can expressed as : where ω = 2π yr −1 is the angular frequency of the annual modulation. The phase shift ωt 0 corresponds to a maximum on July 1 and a minimum on December 1, assuming that the DM distribution is homogeneous and isotropic in the solar system environment. Second, there is a daily modulation, unique to the AQN framework, that can be described as [10] A where ω = 2π day −1 is the angular frequency of the daily modulation. The phase shift φ 0 is similar to ωt 0 in (3). It changes slowly with time due to the variations of the direction of the DM wind with respect to the Earth's position and velocity. It can be assumed to be constant on the scale of days. In both cases, A does not deviate from its average value by more than (10 − 20)%. These modulations do not constitute the topic of the present study.
The main topic of our work corresponds to cases where A(t) 1 for a short period of time τ . These rare burstslike events are called "local flashes" in Ref. [10]. Typically, for A 10 2 , τ 1 second (see Table I). They result from AQNs hitting the Earth and play a key role in our study. Let's consider a network of detectors recording the cross-correlated signals emitted by a trans-passing AQN. Strong signals should be observed by nearby detectors with a short time delays ∆t ∼ d/v , where d is the distance between detectors ( d ∼ 100 km, see [11]) and v the velocity of the trans-passing AQN, v ∼ 10 −3 c, which must be distinguished from the velocity of the AQN-induced axion, v a .
AQN-induced axions must be distinguished from galactic axions, who have smaller velocities ( v gal a ∼ 10 −3 c) : AQN-induced axions are relativistic ( v a ∼ 0.6c) while galactic axions are not ( see Fig. 1) . The two types of axions have drastically different spectral features and these properties are important for reconstructing the trajectory of the trans-passing AQN.
While the width of the spectrum shown on Fig. 1 does not influence the time delays ∆t, this feature drastically modifies the energy of the photons induced by the axion to photon conversion (a → γ) in the detectors. Indeed, the axion's energy, for the dominant portion of the  (2). The corresponding time duration τ and event rate depend on A, which itself is determined by the shortest distance from the nugget's trajectory to the detector. The table is adopted from [10]: The normalized velocity distribution. The plot is adopted from [17]. The four different curves reflect the dependence on model parameters.
axions in Fig. 1, may vary from E a m a to approximately E a 2m a 6 . It implies that the photons resulting from the a → γ conversion will be distributed with a bandwidth ∆ν ν. This bandwidth should be contrasted with the narrow line (with ∆ν/ν ∼ 10 −6 ) in the search for galactic axions, in cavity type detectors. This is an essential difference between the two types of axions : for AQN-induced axions, the EM signal is very broad. Therefore, one should use a broadband detection strategy instead of the conventional cavity type detectors.
There is a number of proposals to design and build broadband instruments : QUAX [48], CASPEr [49], ABRACADABRA [50], LC Circuit [51], DM Radio [52], along with more recent ideas [53][54][55][56]. All of these instruments differ from cavity type instruments. Some suggest to use different observables, such as induced currents in pickup loop (as advocated in [50][51][52][53][54][55]). While there are presently no broadband experiments operating in the window of interest (10 −6 eV m a 10 −3 eV), we do not foresee any fundamental obstacles to designing and 6 We use the natural units = c = k B = 1 throughout this paper.
building such an instrument in future. In what follows, we assume that broadband axion detectors sensitive to a ∆ν ν distribution can be built and assembled in a network (see section IV).

B. The X-ray emission and its spectral features
Here, we overview the spectral characteristics of Xrays emitted by the annihilation of an AQN entering the Earth's atmosphere. In the previous section III A , the dominant portion of the axions was produced deep underground, while the atmospheric portion of the AQN's path was neglected. In the present subsection, on the contrary, the dominant portion of X-rays is emitted when the AQN propagates in the atmosphere, while the underground path can be ignored.
The spectrum of AQNs at low temperatures was analyzed in [57] and was found to be 4 45 where Q(ω, z) describes the emissivity per unit volume from the AQN's electrosphere. It is characterized by a density n(z, T ), where z measures the distance from the quark core of the nugget. Function h(x) in Eq. (5) is a slow varying logarithmic function (explicitly computed in [57]). Integrating over ω gives the total surface emissivity: where the temperature of the electrosphere 7 , for an AQN propagating in the atmosphere, is (10 − 40) keV, depending on the altitude and model-dependent parameters described in Appendix A (see [13] for details). The corresponding spectrum is presented on Fig. 2. The spectrum in Fig. 2 exhibits three important features. First, it is almost flat in the region ω T , which is a direct manifestation of the soft Bremsstrahlung radiation (the emission of a photon with frequency ω is proportional to dω/ω). Second, for large ω T , the exponential suppression term, exp(−ω/T ), becomes the dominant component of the spectrum. Third, due to the function h(x), the emission is strongly suppressed at very small ω ω p T , where ω p is the plasma frequency The drastic intensity drop at small ω T implies that visible light emitted by the AQN with ω ∼ (1 − 10) eV corresponds to T = 50 keV. Adopted from [13].
is strongly suppressed in comparison to the X-ray emission. It implies that AQNs cannot be observed by optical monitoring. This low frequency suppression plays a key role in the interpretation of sky-quakes. For example, [13] could explain why the sky-quake event observed by ELFO was not observed by all-sky cameras, and therefore was not a conventional meteoroid. Looking at sky-quakes, one might ask if AQNs can be detected by seismic and infrasound instruments. Skyquakes are identified with annihilations of AQNs with a very large baryon charge (e.g. B 10 27 for the quake observed by ELFO). Such events only occur once every 10 years. According to (1), annihilations of AQNs with B 10 25 occur approximately once a day in area ∼ (100 km) 2 . Conventional seismic or infrasound instruments are not sufficiently sensitive to record such numerous but weak events. For those, we propose to use a network of X-ray detectors, described in the next section.

IV. TRAJECTORY AND VELOCITY RECONSTRUCTION.
When an event is detected at time t, the instrument records a burst-like signal with bandwidth τ , as shown in Fig. 3. In a synchronized network of detectors, the time delay (t i − t j ) between two nearby stations located at positions R i and R j is Here, the trajectory of an AQN is assumed to be linear [18] with a constant velocity v over the short duration to be compared to the typical time (∼ 10 s) required by an AQN to traverse the Earth. Here, d is the closest distance from a detector to the detected AQN. When detecting axions, the optimal distance is d 100 km, otherwise the benefit of local flashes diminishes [11]. When detecting X-rays, d 50 m is no more than its attenuation length at the sea level.
Depending on the radiation source, instruments may record more features than just t i and τ i , e.g. the incident direction of the flux. This is beneficial but not essential, as we will see t i and τ i are sufficient to determine the trajectory of an AQN. In this work, we assume that t i and τ i are the only data accessible, while additional detected information are devoted to improving the precision in t i and τ i .

A. Reconstruction from a synchronized network
In a synchronized network, it is convenient to choose one station (labeled as "station 1") as the reference location and as the starting time: Hence, the trajectory of an AQN is where d 1 is the closest vector distance from the R 1 to r(t).
In the trajectory (11), there are 6 scalar unknowns: the components of v and d 1 . At least six independent scalar equations are therefore needed to reconstruct the trajectory (11). For a network of three synchronized stations (see Fig. 4), one can show that: FIG. 4: Coordinate system used in Eqs. (12). An AQN passes nearby three stations (grey circles) located at R i (i = 1, 2, 3), along a linear trajectory r(t) (dashed). The distance from each station to r(t) is denoted by d i (i = 1, 2, 3) respectively.
and substituting into Eqs. (12c) and (12d) respectively. These six equations Eqs. (12) have no unique solution because they consist in an implicit Z 2 × Z 2 symmetry, related to reflection and (modified) time reversal (see Appendix A for details). Consequently, for a given observation Eqs. (12) yields four degenerate solutions, one of which being the true trajectory (11). To eliminate this degeneracy, one can either introduce additional physical constraints or add more detectors, as discussed in subsection IV C.

B. Optimal separation distance between stations
In a network of synchronized detectors, the separation distance ∆R between two neighboring stations should not be greater than an effective distance d: where d 100 km for axions and d 50 m for X-rays, otherwise the event rate of correlated signals diminishes due to the sharp reduction in their shared cross section [11]. On the other hand, ∆R cannot be too small because it is inversely proportional to the relative precision in the measurement of the time delay, t i . Hence, optimal performance is achieved for a narrow range of ∆R, which we will discuss below.
To estimate how the trajectory (11) is affected by the input data t i and τ i , we perturb Eqs. (12). We find a qualitative relation where, according to (10), δt 1 = 0. Only the uncertainty in v matters, because d 1 [of order 100 km (resp. 50 m) for axion (resp. X-rays)] is much smaller than the size of the Earth. A small perturbation δd 1 corresponds to a negligible shift of the trajectories at the Earth surface.
Presumably, |δv| can be expressed in the empirical form of independent noise: where C is a dimensionless factor of order one. Monte Carlo simulation suggests C 1.2 for 95% of samples (see Appendix B 1). From Fig. 3, we expect δt to be proportional to the bandwidth τ i : where ε ∈ (0, 1) relates to the efficiency of measurement. For a moderately large ε 0.1, it corresponds to a low signal-to-noise ratio. This is particularly true for axion detection. In this case, the δt i /t i terms dominate inequality (16): where definitions (8) and (9) were substituted. ∆R is defined as the average separation distance between two nearest stations, and the root mean square (R ·v) −1 12 from Monte Carlo simulation for 95% of samples (see Appendix B 1). Precision of measurement can be enhanced by installing more than one detectors in a station. Assuming there is a total of N detectors in the synchronized network, to ensure that |δv|/v is sufficiently small, the constraint on nearest separation distance is N must be at least 4 (or 3 with additional physical constraints) to obtain a unique solution for Eqs. (12), as discussed in the previous subsection. Hence, even for a moderately low signal-to-noise ratio (ε 0.1), three or four stations are sufficient to measure trajectories, provided their closest separation distance is of order d.
To summarize, constraints (14) and (19) imply that the optimal value for ∆R is of order d, such that both the event rate and the precision of measurement are maximized.

C. Optimal configurations of a synchronized network
The optimal configuration of a synchronized network depends on the type of radiation. For axions, a triangular network is more cost effectiveness because it has the highest ratio event rate per station, especially when operating on an extensible global network. If uniqueness of solution is a priority, a square network is better for axion detection. Another option is to construct a dense network of X-ray detectors consisting of tetrahedral units. X-ray technologies are well developed and overall a lot less expensive.
A summary of configurations is presented in Table II, where the event rate is estimated as similar to Ref. [10], with a more sophisticated estimation based on the Monte Carlo evaluation of cross section.
Here η ∈ (0, 1) is the efficiency of the detection and is assumed to be of order 1. σ and σ 0 are the effective cross section of the network and a single station, respectively (see Appendix B 2 for details). R ⊕ = 6371 km is the radius of Earth.
In what follows, we investigate the optimal configurations for axions and X-rays.

Axion: network of broadband stations
As discussed in subsection IV A, the trajectory (11) can be identified by synchronized measurements from three stations, based on set of Eqs. (12). In Appendix A, it is proven that Eqs. (12) yields four degenerate solutions. To correctly identify the true trajectory out of the four, one should either introduce additional physical constraints or add stations to the network.
When antimatter AQNs annihilate with normal matter inside Earth, there is an additional constraint: axion emission only takes place underground. This implies three additional constraints for axions: Simulation indicates that (21) is sufficient to ensure 94% of the trajectories are uniquely determined by Eqs. (12) (see Appendix B 3 for details). The optimal configuration is a three-station network: an equilateral triangle with d 100 km (Fig. 5a). The event rate is 1.7η day −1 . An alternative way to eliminate the degeneracy is to construct a four-station network, with no additional constrains (see Appendix B 3). Fig. 5b illustrates the optimal configuration of such a network: each station is located at the vertices of a square, with d 100 km. The trade off is the limited improvement in event rate (2.0η day −1 ) compared to the triangular network ( 1.7η day −1 ). If uniqueness is not a major concern, it is more efficient to enhance the event rate by extending a triangular network than using a square one. 2. X-rays: array of ground-based detectors X-rays propagate a shorter distance in air (d 50 m) than axions. It is therefore possible to construct a regular tetrahedral network ( Fig. 6). A tetrahedral network is optimal as it has both a high event rate (due to its triangular structure) and a unique solution (it consists of four detectors). It can serve as a base unit to form a large array of detectors. For increased performance, the array should be lifted 50 meters above ground so that the detection covers all directions in space.
According to simulation in Appendix B 2, the event rate for a single tetrahedral unit is 4.2 × 10 −7 η day −1 (ground based) and 1.2 × 10 −6 η day −1 (50 m above ground). For an array covering a surface area by (100 km) 2 , the event rate is estimated to be 2.9η day −1 .

V. CONCLUSION
The main goal of this work is to develop a new search strategy which allows to study the directionality and ve- II: Synchronized networks. The incoming AQN flux is assumed to be isotropic. For axions (resp. X-rays), the nearest distance between detectors is chosen to be 100 km (resp. 50 m). The typical time delay ∆t and bandwidth τ are of order 0.1 seconds (resp. 0.1 milliseconds). The efficiency parameter η is chosen to be 0.2. ground-based X-ray Regular tetrahedron 9.48 × 10 −5 1.15 × 10 −4 0.35 2.5 × 10 −7 50 m above ground X-ray X-ray arrays (10 km) 2 1.14 × 10 −4 1.37 × 10 −4 0.64 5.8 × 10 −3 10 5 X-ray detectors X-ray X-ray arrays (100 km) 2 1.14 × 10 −4 1.37 × 10 −4 0.64 0.58 10 7 X-ray detectors d FIG. 6: A base unit of the X-ray array of detectors. It consists of four detectors (grey circles) located at the vertices of a regular tetrahedral, separated by d = 50 m. locity distribution of DM particles. To reach this goal, we suggest using a network of synchronized detectors and we apply this idea to two specific cases: 1. a broadband detection strategy to search for axions, within the context of the Axion Quark Nugget (AQN) dark matter model; 2. a network of X-ray detectors that allows to reconstruct the velocity vector of incident AQN particles. We argue that by looking at the time delays recorded by two or more nearby detectors, such a network allows to discriminate between the DM signal and larger noise.

Radiation
In case 1, we focused on typical "small" events, i.e. AQN with a baryon charge B ∼ 10 25 , propagating through the Earth and its atmosphere once a day per (100 km) 2 . We showed that a network of four synchronized detectors is necessary to uniquely determine the directionality of the incident flux and the velocity distribution of the DM flux. On the other hand, we showed that an optimal configuration for a direct detection of AQN-induced axions, when AQNs penetrate the Earth, is a synchronized network of three stations placed at the vertices of an equilateral triangle, 100 km apart. For an array covering a surface area of (100 km) 2 , the event rate is estimated to be 1.7η day −1 . Such a network uniquely determine the directionality of the incident DM particle in 94% of cases.
AQNs propagating through the Earth atmosphere are accompanied by signature X-rays, which therefore constitute an indirect detection of AQNs. X-ray stations are easier to build and operate then axion detectors. In case 2, we showed that the optimal configuration to detect this emission is an array of tetrahedral units. Each unit consist in four synchronised detectors positioned at the vertices of a regular tetrahedron, 50 m apart. Such a network allows to uniquely determine the directionality of the emission. For an array covering a surface area of (100 km) 2 , the event rate is estimated to be 2.9η day −1 . One should emphasize that such a network can eliminate the dominant portion of spurious signals. Therefore, recording synchronized events would be a very strong argument supporting their unconventional nature 8 . Detecting DM particles with those networks can have fundamental consequences since reconstructing the directionality and velocity distribution of DM particles leads to testing the Standard Halo Model.

ACKNOWLEDGMENTS
AZ is thankful to D. Budker and V. Flambaum for infinitely long, never ending discussions about possible networks, which motivated the present studies. This work was supported in part by the National Science and Engineering Research Council of Canada and X.L. by the UBC four year doctoral fellowship. We are grateful to Benoite Pfeiffer for a careful reading of the manuscript.
Appendix A: Degenerate solutions in Eqs. (12) To demonstrate the existence of degenerate solutions, we first solve for Eqs. (12) perturbatively in the following limit of mirror symmetry (i.e., 2 ↔ 3 with similar values): R = R(x cos γ −ŷ sin γ), R 3 = R(x cos γ +ŷ sin γ) , (A1b) where ε t , ε τ 1, and 2γ is the angle between R 2 and R 3 . One of the solution is where S and T are dimensionless parameters determined by t i , τ i and γ: The three remaining solutions are obtained by applying one of, or a combination of, the following operators to solution (A2): The reflection R is an obvious symmetry in Eqs. (12), while the modified time reversal T is less apparent. Note that in Eqs. (12), the time reversal cannot be obtained by simply flipping the sign of time t → −t because one has to also reverse the labeling of each station (1, 2, 3) → (3, 2, 1) at the same time. Therefore, to reverse the time without changing the labeling, we will need the modified "time" T as defined in Eqs. (A3b), and then we take the time reversal of T . We make an intuitive plot of the four degenerate solutions (Z 2 × Z 2 ) in Fig. 7. Qualitatively, solving Eqs. (12) is equivalent to looking for trajectories simultaneously tangent to three sphere of given radii d i (i = 1, 2, 3), where d i is the distance of the AQN trajectory from the i-th detector and is determined by the input parameter (τ i , t i ). Apparently, there are always four tangent lines that satisfy such criteria. In practice, T {d i } are moderately different from the original {d i } because the transformation accounts for a modified time scale.  (12). Allowed trajectories are simultaneously tangent to three spheres of given radii d i (i = 1, 2, 3), where d i is the distance of the AQN trajectory from the i-th detector and is determined by the input parameter (τ i , t i ). In practice, the radii change moderately upon the modified time reversal T transformation to account for the rescaled (v, d 1 ).

Appendix B: Details of numerical simulation
The details of all numerical simulations based on Monte Carlo method are presented in this appendix. The distribution of the AQN trajectories follows an isotropic distribution as derived in Ref. [18], namely cos θ ∼ uniform(−1, 1), φ ∼ uniform(0, 2π) , (B1c) where µ = 220 km s −1 is mean galactic velocity, σ = 110 km s −1 is the galactic dispersion, the effective distance of detection d is 100 km (resp. 50 m) in case of axions (resp. X-rays). Following the simulation conditions (B1), we define two angleŝ n ≡ (sin θ cos φ, sin θ sin φ, cos θ) , v ≡ (sin ψ cos ϕ, sin ψ sin ϕ, cos ψ) , so that a trajectory (11) is generated with given v and Sec. IV C suggests the optimal configurations of a synchronized network are equilateral triangle, square, and regular tetrahedron with separation distance of order d.
For illustrative purpose, we set the distance between two neighboring stations to be ∆R = d (and d/ √ 2 in the square case) in the following simulations.
One more implicit constraint is that all stations should have the same range of effective distance of detection, namely for each i-th station. From the definition (13), the time input parameters are related to the simulated v and d 1 in what follows: (B5) For simulations in this appendix, the trajectories are generated following the above recipe described.

Sensitivity to time measurement
In addition to a simulated trajectory (11) parameterized by (v, d 1 ), a random perturbation (δv, δd 1 ) should be generated: where ε cut is a cutoff of perturbation and we choose its upper limit to be 0.5 in the numerical simulation. Now we definê n ≡ (sin θ cos φ , sin θ sin φ , cos θ ) , v ≡ (sin ψ cos ϕ , sin ψ sin ϕ , cos ψ ) , (B7a) δv = δvv , δd 1 = δd 1n ×v |n ×v | ≡ δd 1d 1 , where δd 1 is completely determined by the constraint of orthogonality v · d 1 = 0, namely Then we obtain from Eqs. (B5). Additionally, we only keep data points of which the perturbations are not overly large: because the uncertainty of time measurement is assumed small. From Monte Carlo simulation, a large sample array of (δt i , δτ i , δv, δd 1 ) is generated. Fitting the data set into the inequality (16), C is found to be 1.24 (resp. 1.92) among 95% (resp. 99%) of 10 6 samples at ε cut = 0.5, see Fig. 8. The value of C is insensitive (≤ 0.1 absolute uncertainty) to cutoff ε cut and choice of separation distance δR.  (16) in Monte Carlo simulation (10 6 samples). The parameter ε cut is chosen to be 0.5, and the count is normalized. The value of C is found to be 1.24 (1.92) among 95% (99%) of samples.
One more important quantity is the root mean square (R·v) −1 as defined in Eq. (18). As shown in Fig. 9, the value of (R·v) −1 is found to be 12.3 (resp. 32.0) among 95% (resp. 99%) of 10 6 samples, which is independent of ∆R by its definition.

Cross section of certain network configurations
We generate 10 7 samples using Eqs. (B1) to Eqs. (B3), with the modification that and then check what percentage of the samples satisfy the condition (B4) for each station in different configurations, in the case of axions we also seek to satisfy the constraint (21). We also take the ratio σ/σ 0 of each configuration to the single station. This was repeated for N = 1, 2, 100 and for d = 1, 100 but it was found to have no effect on the percentage of trajectories satisfying the conditions. The configurations used are as described in subsection IV C where ∆R = d. Additionally for X-ray we also used a flat disk configuration, where for each trajectory we see what percentage of samples satisfy r z = 0 and r 2 x + r 2 y < R 2 . Results are summarized in Table II. The cross section ratio σ/σ 0 determines the event rate of a local flash as defined in Eq. (20). We conclude σ/σ 0 ∼ O(0.1) for all configurations listed as expected. Monte Carlo simulation (10 6 samples). The count is normalized. The value of (R ·v) −1 is found to be 12. 3 (resp. 32.0) among 95% (resp. 99%) of samples.

Uniqueness of solutions
We first perform a simulation check for the axion network. In this case, additional constraints (21) applies to every station. Both triangular and square configurations are examined. The time delay t i and bandwidth τ i are retrieved from Eqs. (B5). Feeding the simulated time parameters into Eqs. (12), the equation set is solved under constraints (21). For solutions that have no more than 10% difference in velocity from each other, they are treated as identical. With 10 6 samples simulated, we find 94% of trajectories are the unique solutions in Eqs. (12) in the triangular configuration, and the rate of uniqueness is 100% in the square configuration.
For the x-ray network, the tetrahedral configuration is examined under the same process without the constraints (21) for generality. With 10 6 samples simulated, we find the rate of uniqueness is 100%.

Time delay and bandwidth
In order to get the typical values of time delay and bandwidth we do a Monte Carlo simulation with 10 7 samples. During the simulation for each valid trajectory we record the values of t i and τ i of each station according to Eqs. (B5), and then we take the average value along with the standard deviation. The statistics for t 1 and τ 1 are not collected because t 1 is always 0 by setting (10). A sample simulated signal received in a triangular network is presented in Fig. 10. The average time delay |∆t| and average bandwidth share similar values: 0.2 s (resp. 0.1 ms) for radiation of axions (resp. X-rays). This is as expected from definitions (8) and (9). The results are summarized in Table II  Each signal corresponds to each station detecting the axion flux due to a nearby AQN. From these signals we can obtain the parameters t i and τ i in order to solve the system of Eqs. (12)