Search for topological defect dark matter using the global network of optical magnetometers for exotic physics searches (GNOME)

Samer Afach, 2 Ben C. Buchler, Dmitry Budker, 2, 4 Conner Dailey, ∗ Andrei Derevianko, Vincent Dumont, Nataniel L. Figueroa, 2 Ilja Gerhardt, Zoran D. Grujić, 9 Hong Guo, Chuanpeng Hao, Paul S. Hamilton, Morgan Hedges, Derek F. Jackson Kimball, Dongok Kim, 15 Sami Khamis, Thomas Kornack, Victor Lebedev, Zheng-Tian Lu, Hector Masia-Roig, 2, † Madeline Monroy, 13 Mikhail Padniuk, Christopher A. Palm, Sun Yool Park, ‡ Karun V. Paul, Alexander Penaflor, Xiang Peng, Maxim Pospelov, 21 Rayshaun Preston, Szymon Pustelny, Theo Scholtes, 22 Perrin C. Segura, § Yannis K. Semertzidis, 15 Dong Sheng, Yun Chang Shin, Joseph A. Smiga, 2, ¶ Jason E. Stalnaker, Ibrahim Sulai, Dhruv Tandon, Tao Wang, Antoine Weis, Arne Wickenbrock, 2 Tatum Wilson, Teng Wu, David Wurm, Wei Xiao, Yucheng Yang, Dongrui Yu, and Jianwei Zhang Johannes Gutenberg-Universität Mainz, 55128 Mainz, Germany Helmholtz-Institut Mainz, GSI Helmholtzzentrum für Schwerionenforschung, 55128 Mainz, Germany Centre for Quantum Computation and Communication Technology,

Results are reported from the first full-scale search for transient signals from exotic fields of astrophysical origin using data from a newly constructed Earth-scale detector: the Global Network of Optical Magnetometers for Exotic physics searches (GNOME). Data collected by the GNOME consist of correlated measurements from optical atomic magnetometers located in laboratories all over the world. GNOME data are searched for patterns of signals propagating through the network consistent with exotic fields composed of ultralight bosons such as axion-like particles (ALPs). Analysis of data from a continuous month-long operation of the GNOME finds no statistically significant signals consistent with those expected due to encounters with topological defects (axion domain walls), placing new experimental constraints on such dark matter scenarios.

I. INTRODUCTION
The nature of dark matter, an invisible substance comprising over 80% of the mass of the universe [1][2][3], is one of the most profound mysteries of modern physics. Although evidence for the existence of dark matter comes from its gravitational interactions, unraveling its nature likely requires observing non-gravitational interactions between dark matter and ordinary matter [4]. One of the leading hypotheses is that dark matter consists of ultralight bosons such as axions [5][6][7] or axion-like particles (ALPs) [8][9][10]. Axions and ALPs arise from spontaneous symmetry breaking at an unknown energy scale f SB , which, along with their mass m a , determines many of their physical properties.
ALPs can form stable, macroscopic field configurations in the form of topological defects [11][12][13][14] or composite objects bound together by self-interactions such as boson stars [15][16][17]. Such ALP field configurations could concentrate the dark matter density into many distinct, compact spatial regions that are small compared to the galaxy but much larger than the Earth. In such scenarios, Earth-bound detectors would only be able to measure signals associated with dark matter interactions on occasions when the Earth passes through such a dark-matter object. It turns out that there is a wide range of parameter space, consistent with observations, for which such dark-matter objects can have the required size and abundance such that the characteristic time between encounters could be on the order of one year or less [13,17,18]. This opens up the possibility of searches with terrestrial detectors. Here we present the results of such a search for ALP domain walls, a class of topological defects which can form between regions of space with different vacua of an ALP field [11][12][13][14]. 1 Since ALPs can interact with atomic spins [4], the passage of Earth through an ALP domain wall affects atomic spins similarly to a transient magnetic field pulse [13,17]. Considering a linear coupling between the ALP field gradient ∇a(r, t) and atomic spin S, the interaction Hamil- * Currently: University of Waterloo, Department of Physics and Astronomy, N2L 3G1,Ontario, CA † Electronic address: hemasiar@uni-mainz.de ‡ Currently: JILA, NIST and University of Colorado, and Department of Physics, University of Colorado, Boulder Colorado 80309-0440, USA § Currently: Department of Physics, Harvard University, Cambridge, MA 02138 ¶ Electronic address: jsmiga@uni-mainz.de 1 A different type of bosonic field configuration is the Q-ball [19,20] which is spherically symmetric having nearly constant field magnitude in the bulk with a rapid fall off at the boundaries. It is another possible dark matter candidate. We note that the signal pattern produced in GNOME for sufficiently large Q-balls with couplings to atomic spins can be similar to that of ALP domain walls [17].
tonian can be written as where is the reduced Planck's constant, c is the speed of light, r is the position of the spin, t is the time, and f SB /ξ ≡ f int is the coupling constant in units of energy described with respect to the symmetry-breaking scale f SB [21], where ξ is unitless. In most theories, the coupling constants f int describing the interaction between Standard Model fermions and the ALP field are proportional to f SB ; however, f int can differ between electrons, neutrons, and protons by model-dependent factors that can be significant in some cases [4,8].
In analogy with Eq. (1), the Zeeman Hamiltonian describing the interaction of a magnetic field B with an atomic spin S can be written as where γ is the gyromagnetic ratio. Since Eqs. (1) and (2) have the same structure, the gradient of the ALP field can be treated as a pseudo-magnetic field. An important distinction between the ALP-spin interaction [Eq. (1)] and the Zeeman interaction [Eq. (2)] is that while γ tends to scale inversely with the fermion mass, no such scaling of the ALP-spin interaction is expected [4]. The amplitude, direction, and duration of the pseudomagnetic field pulse associated with the transit of the Earth through an ALP domain wall depends on many unknown parameters such as the energy density stored in the ALP field, the coupling constant f int , the thickness of the domain wall, and the relative velocity v between Earth and the domain wall. The dynamical parameters, such as the velocities of the dark matter objects, are expected to randomly vary from encounter-to-encounter. We assume that they are described by the Standard Halo Model for virialized dark matter [22]. Furthermore, the abundance of domain walls in the galaxy is limited by physical constants, m a and f SB , as these determine the energy contained in the wall and the total energy of all domain walls is constrained by estimates of the local dark matter density [23]. The expected temporal form of the pseudo-magnetic field pulse can depend both on the theoretical model describing the ALP domain wall as well as particular details of the terrestrial encounter such as the orientation of the Earth. The relationships between these parameters and characteristics of the pseudo-magnetic field pulses searched for in our analysis are discussed in Appendix A and Refs. [13,17,21].
The Global Network of Optical Magnetometers for Exotic physics searches (GNOME) is a worldwide network searching for correlated signals heralding beyond-the-Standard-Model physics which currently consists of more than a dozen optical atomic magnetometers, with stations in Europe, North America, Asia, the Middle East, and Australia. A schematic of a domain-wall encounter with the GNOME with stations used in this study is shown in Fig. 1. The measurements from the magnetometers are recorded with data-acquisition systems, synchronized to the Global Positioning System (GPS) time, and uploaded to servers located in Mainz, Germany, and Daejeon, South Korea. Descriptions of the operational principles and characteristics of GNOME magnetometers are presented in Sec. IV and Ref. [24], details on the GNOME data acquisition system are discussed in Ref. [25] (see also Table I).
The active field sensor at the heart of every GNOME magnetometer is an optically pumped and probed gas of alkali atoms. Magnetic fields are measured through variations in the Larmor spin precession of the optically polarized atoms. The vapor cells containing the alkali atoms are placed inside multi-layer magnetic shielding systems which reduce background magnetic noise by orders of magnitude [26] while retaining sensitivity to exotic spin-couplings between ALP dark matter and atomic nuclei. 2 Since all GNOME magnetometers presently use atoms whose nuclei have a valence proton, the signal amplitudes measured by the GNOME due to an ALP-spin interaction are proportional to the relative contribution of the proton spin to the nuclear spin (as discussed in Appendix B and Ref. [28]). This pattern of signal amplitudes [Eq. (1)] can be characterized by a pseudo-magnetic field B j measured with sensor j: where is the normalized pseudo-magnetic field describing the effect of the ALP domain wall on proton spins and µ B is the Bohr magneton. The ratio between the Landé gfactor and the effective proton spin (g F,j /σ j ) accounts for the specific proton-spin coupling in the respective sensor. This ratio depends on the atomic and nuclear structure as well as details of the magnetometry scheme, see Appendix B. Since each GNOME magnetometer measures the projection of the field along a particular sensitive axis, the factor η j is introduced to account for the directional sensitivity. This factor, given by the cosine of the angle between B p and the sensitive axes, takes on values between ±1.
In spite of the unknown properties of a particular terrestrial encounter with an ALP domain wall, the GNOME would measure a recognizable global pattern of the associated pseudo-magnetic field pulse amplitudes described by Eq. (3), as illustrated in Fig. 1b. The associated pseudo-magnetic field pulses would point along a common axis, have the same duration, and exhibit a characteristic timing pattern. The data-analysis algorithm used in the present work to search for ALP domain walls is described in Sec. IV and Ref. [29]. The algorithm searches for a characteristic signal pattern across the GNOME having properties consistent with passage of the Earth through an ALP domain wall. Separate analyses to search for transient oscillatory signals associated with boson stars [17] and bursts of exotic low-mass fields (ELFs) from cataclysmic astrophysical events [30] are presently underway.
Here we report the first results of a dark matter search with the GNOME: a search for transient couplings of atomic spins to macroscopic dark-matter objects, and therefore demonstrate the ability of the GNOME to explore parameter space previously unconstrained by direct laboratory experiments. Searches for macroscopic darkmatter objects based on similar ideas were carried out using atomic clock networks [18,[31][32][33], and there are a number of experimental proposals utilizing other sensor networks [34][35][36][37][38]. All of these networks are sensitive to bosonic dark matter with a scalar coupling to Standard Model particles [4]. The GNOME is sensitive to a different class of dark matter: bosons with pseudoscalar couplings to Standard Model particles. Pseudoscalar bosonic dark matter generally produces no observable effects in clock networks [4] but does couple to atomic spins via the interaction described by Eq. (1). Thus the GNOME is sensitive to a distinct, so far mostly unconstrained, class of interactions as compared to other sensor networks.

II. RESULTS
There have been four GNOME Science Runs between 2017 and 2020 as discussed in Sec. IV. Here we analyze the data from Science Run 2, which had comparatively good overall noise characteristics and consistent network operation (as seen in Fig. 5). Nine magnetometers took part in Science Run 2 that spanned from 29 November 2017 to 22 December 2017. The characteristics of the magnetometers are summarized in Table I.
Before the data are searched for evidence of domainwall signatures, they are preprocessed by applying a rolling average, high-pass filters, and notch filters to the raw data. The averaging enhances the signal-to-noise ratio for certain pulse durations, avoids complications arising from different magnetometers having different bandwidths, and reduces the amount of data to be analyzed. The high-pass and notch filters reduce the effects of longterm drifts and noisy frequency bands. We refer to the filtered and rolling-averaged data set as the "search data." The search data are examined for evidence of collective signal patterns corresponding to planes with uniform, non-zero thickness, crossing Earth at constant velocities. The imprinted pattern of amplitudes depends on the domain-wall crossing velocity [29]. As noted in Sec. I, we assume that the domain-wall-velocity probability density function follows the Standard Halo Model for virialized dark matter. The signature of a domain-wall crossing the magnetometer network depends on the perpendicular component to the domain-wall plane of the relative velocity between the Earth and the domain wall, v ⊥ . A lattice of points in velocity space is constructed such that the search algorithm covers 97.5% of the velocity probability density function. The algorithm scans over the velocity lattice and, for every velocity, the data from each magnetometer are appropriately time-shifted so that the signals in different magnetometers from a hypothetical domain-wall crossing with the given velocity occur at the same time. For each velocity and at each measurement time, the amplitudes measured by each magnetometer are fit to the ALP domain-wall crossing model described in Ref. [29]. As a result, estimations for signal magnitude and domain-wall direction, along with associated uncertainties, are obtained for each measurement time and all lattice velocities. The magnitude-to-uncertainty ratio of an event is given by the ratio between the signal magnitude and its associated uncertainty.
The search algorithm uses two different tests to evaluate if a given event is likely to have been produced by an ALP domain-wall-crossing: a domain-wall model test and a directional-consistency test [29]. The domainwall model test evaluates whether the event amplitudes measured by the GNOME magnetometers match the signal amplitudes predicted by the ALP domain-wall crossing model, and is quantified by the p-value as discussed in Sec. IV and Ref. [29]. The directional-consistency test checks the agreement between the direction of the scanned velocity and the estimated domain wall direction, and is quantified by the angle between the two directions normalized by the angle between adjacent lattice velocities.
In order to evaluate the detection probability of the search algorithm, a well-characterized data set that includes domain-wall-crossing signals with known properties is required. For this purpose, we generate a background data set by randomly time-shuffling the search data so that the relative timing of measurements from different GNOME stations is shifted by amounts so large that no true-positive events could occur. By repeating the process of time shuffling, the length of the background data can be made to far exceed the search data. This method is used to generate background data with noise characteristics closely reproducing those of the search data [39]. A set of pseudo-magnetic field pulses matching the expected amplitude and timing pattern produced by the passages of Earth through ALP domain walls are inserted into the background data to create the test data. The test data are used to evaluate the truepositive probability for the search algorithm for different thresholds on the p-value and the normalized angular difference. The detection probability as a function of the thresholds is studied in order to find a combination of thresholds that ensures the detection of 97.5% of the inserted domain walls (see Sec. IV and Fig. 6). This results in an overall detection efficiency ≥ 95% for the search algorithm, considering both the incomplete velocity lattice coverage and the detection probability.
The search data are analyzed for domain-wall encounters using the algorithm presented in Ref. [29]. The cumulative distribution of candidate events as a function of their magnitude-to-uncertainty ratio is shown as a solid green line in Fig. 2. The background data are analyzed in the same way as the search data to estimate the expected number of background events during Science Run 2. Time-shuffling is used to create 10.7 years of background data. The number of candidate events measured in the background data re-scaled to the duration of Science Run 2 is shown as a dashed blue line in Fig. 2. The significance of a candidate event in the search data is given by the probability of detecting one or more background events at a magnitude-to-uncertainty ratio above that of the candidate event [see Eq. (10)]. The probability of finding a number of candidate events in the background data follows Poisson statistics. The candidate event in the search data with the largest magnitude-to-uncertainty ratio (= 12.6) had a significance of less than one sigma. Therefore, we find no evidence of an ALP domain-wall crossing during Science Run 2.
In order to evaluate the domain-wall characteristics excluded by this result, the observable domain-wall crossing parameters above 12.6 magnitude-to-uncertainty ratio during Science Run 2 are determined. The GNOME has nonuniform directional sensitivity [29]; we conservatively estimate the network sensitivity assuming the domain wall comes from the least-sensitive direction. Figure 3 shows the active timeT (∆t, B p ), i.e., how long the network was sensitive to domain-walls as a function of pseudo-magnetic field-pulse-magnitude sensitivity, B p , and pulse duration, ∆t. A signal with pseudo-magnetic field magnitude B p produces a magnitude-to-uncertainty ratio of ζ = B p /B p . The characteristic shape of the sensitive region is a result of the filtering and averaging of the raw data as described in Sec. IV. Averaging reduces the sensitivity of the search data to short pulse durations and high-pass filtering suppresses sensitivity for long ∆t. The GNOME sensitivity varies in time as the number of : Amount of timeT that the GNOME is sensitive to domain walls with a given duration ∆t and normalized pseudomagnetic field magnitude sensitivity B p throughout Science Run 2, the magnitude of an event needed to induce a signal with a magnitude-to-uncertainty ratio of one [see Eq. (11)].
Only the worst-case direction is considered. The plot assumes the parameters of the analysis: 20 s averaging time, 1.67 mHz first-order zero-phase Butterworth filter, and 50 Hz and 60 Hz zero-phase notch filters with a Q-factor of 60.
active GNOME magnetometers recording data and their background noise change. The active time,T (∆t, B p ), can be used to constrain ALP domain-wall parameter space as discussed in Appendix A. If one assumes a probability distribution for the number of domain-wall encounters, an upper bound on the rate R C of such encounters can be calculated with a confidence level C. We assume a Poisson probability distribution for the domain-wall crossings. Since the excess number of events in the search data as compared to the background data was not statistically significant, the upper bound on the observable rate is given by the probability of measuring no events during the effective time [40]. Note that sinceT depends on the parameters of the domain-wall crossing, our constraint on the observed rate depends on the ALP properties. We choose the confidence level to be C = 90%.

III. DISCUSSION
The analysis of the GNOME data did not find any statistically significant excess of events above background during Science Run 2 that could point to the existence of ALP domain walls, as seen in Fig. 2. The expected rate of domain-wall encounters, r, depends on the ALP mass, m a , the domain-wall energy density in the Milky Way, ρ DW , the typical relative domain-wall speedv given by the Standard Halo Model, and the symmetry breaking scale, f SB . The region of parameter space to which GNOME is sensitive is defined by the ALP parameters expected to produce signals above 12.6 magnitude-touncertainty ratio with rates r ≥ R 90% during Science Run 2 (see Fig. 3). Based on the null result of our search, the sensitive region is interpreted as excluded ALP parameter space.
The ALP parameters and the phenomenological parameters describing the ALP domain walls in our galaxy, namely the thickness ∆x, the surface tension or energy per unit area σ DW , and the average separationL can be related through the ALP domain wall model described in Refs. [13,21]. A full derivation of how observable parameters are related to ALP parameters is given in Appendix A and summarized below.
The domain-wall thickness is determined by the ALP mass, and is on the order of the ALP reduced Compton wavelength λ a [21], The constant prefactor of 2 √ 2 is obtained by approximating the spatial profile of the field-gradient magnitude as a Lorentzian and defining the thickness as full width at half maximum. For a given relative velocity component perpendicular to the domain wall v ⊥ , the signal duration is We assume that domain walls comprise the dominant component of dark matter. Thus with the energy density ρ DW ≈ 0.4 GeV/cm 3 in the Milky Way [23], the energy per unit area (surface tension) in a domain wall, σ DW , determines the average separation between domain walls, L. The surface tension σ DW is related to the symmetry breaking scale [13], The average domain-wall separation is then approximated byL which results in the average domain-wall encounter rate, The colored region in Fig. 4a describes the symmetry breaking scales up to which GNOME was sensitive with 90% confidence (as detailed in Appendix A). The parameter space is spanned by ALP mass, maximum symmetry breaking scale, and ratio between symmetry breaking scale and coupling constant. The shape of the sensitive area shown in Fig. 4a is determined by the event with the largest magnitude-to-uncertainty and the characteristics of the preprocessing applied to the raw data. Figure 4b shows various cross-sections of Fig. 4a for different ratios between the symmetry breaking scale and the coupling constant indicated by the dashed lines. The upper bound of f SB that can be observed by the network is shown in Fig. 4b , there is a sharp cut-off for low ALP mass where the corresponding field magnitude falls below the network sensitivity. Even though B p increases for large m a , the mean rate of domain-wall encounters decreases with increasing mass [see Eqs. (8) and (9)]. Correspondingly, the upper limit for the symmetry-breaking scale f SB is ∝ 1/ √ m a . Given that no events were found, the sensitive region of ALPdomain-wall parameter space during Science Run 2 can be excluded.
Our experiment explores ALP parameter space up to f int ∼ 3×10 5 GeV (see Fig. 4). This goes beyond that excluded by previous direct laboratory experiments searching for ALP-mediated exotic pseudoscalar interactions between protons which have shown that f int 300 GeV over the ALP mass range probed by the GNOME [41]. Although astrophysical observations suggest that f int 2 × 10 8 GeV, there are a variety of scenarios in which such astrophysical constraints can be evaded [42,43].
Future work of the GNOME collaboration will focus both on upgrades to our experimental apparatus and new data-analysis strategies. One of our key goals is to improve overall reliability and the duration of continuous operation of GNOME magnetometers 3 . Additionally, magnetometers varied in their bandwidths and reliability, and stability of their calibration. These challenges were addressed in Science Run 4 through a variety of magnetometer upgrades and instituting daily worldwide test and calibration pulse sequences. However GNOME suffered disruptions due to the worldwide COVID-19 pandemic. We plan to carry out Science Run 5 in 2021 to take full advantage of the improvements. Furthermore, by upgrading to noble-gas-based comagnetometers [44,45] for future science runs (Advanced GNOME), we expect to significantly improve the sensitivity to ALP domain walls. Additionally, GNOME data can be searched for other signatures of physics beyond the Standard Model such as boson stars [17], relaxion halos [46], and bursts of exotic low-mass fields from black-hole mergers [30]. In terms of the data-analysis algorithm used to search for ALP domain walls, in future work we aim to improve the efficiency of the scan over the velocity lattice. The number of points in the velocity lattice to reliably cover a fixed fraction (e.g., 97.5%) of the ALP velocity probability distribution grows as (∆t) [where ∆t is given by Eq. (6)]. This makes the algorithm computationally intensive. We are investigating a variety of new analysis approaches, such as machinelearning-based algorithms, to address these issues.

IV. METHODS
The GNOME consists of over a dozen optical atomic magnetometers, each enclosed within a multi-layer magnetic shield, distributed around the world [26]. GNOME magnetometers are based on a variety of different atomic species, optical transitions, and measurement techniques: some are frequency-or amplitude-modulated nonlinear magneto-optical rotation magnetometers (NMOR) [47,48], some are rf-driven optical magnetometers [24], while others are spin-exchange-relaxation-free magnetometers (SERF) [49]. A detailed description and characterization of six GNOME magnetometers is given in Ref. [24]. A summary of the properties of the GNOME magnetometers active during Science Run 2 is presented in Table I.
Each GNOME station is equipped with auxiliary sensors, including accelerometers, gyroscopes, and unshielded magnetometers, to measure local perturbations that could mimic a dark matter signal. Suspicious data are flagged [24] and discarded during the analysis.
The number of active GNOME magnetometers during the four Science Runs and the combined network noise as defined in Ref. [29] is shown as a function of time in Fig. 5. The data in Science Run 4, although they extended for a longer period of time, featured poorer noise characteristics and consistency of operation to Science Run 2. Since many GNOME stations underwent upgrades in 2018 and 2019, further characterization of Science Run 4 data is needed, and results will be presented in future work.
Here we provide more details on the analysis procedure discussed in Sec. II. It is composed of three stages to identify events likely to be produced by ALP domain-wall crossings: preprocessing, velocity scanning, and postselection [29]. First, in the preprocessing stage, a rolling average and filters are applied to the GNOME magnetometer raw data which are originally recorded by the GPS-synchronized data-acquisition system at a rate of 512 samples/s [25]. The rolling average is characterized by a 20 s time constant. Noisy frequency bands are suppressed using a first-order Butterworth high-pass filter at 1.67 mHz together with the notch filters corresponding to 50 Hz and 60 Hz power line frequencies with a quality factor of 60. These filters are applied forward and backward to remove phase effects. This limits the observable pulse properties to a frequency region to which all magnetometers are sensitive. Additionally, it guarantees that the duration of the signal is the same in all sensors. We note that these filter settings may be changed in future analysis. The local standard deviation around each point in the magnetometer's data is determined using an iterative process. Outliers are discarded until the standard deviation of the data in the segment converges. The local standard deviation is calculated taking 100 downsampled points around each data point.
Additionally, auxiliary measurements have shown that the calibration factors used by each magnetometer to convert raw data into magnetic field units experience change over time due to, for example, changes in the environmental conditions. Upper limits on the calibration factor errors due to such drifts over the course of Science Run 2 have been evaluated and are listed in Table I. Calibration factor errors result in magnetic field measurement errors proportional to the magnetic field B j . The uncer-   Fig. 5). The station name, location in longitude and latitude, orientation of the sensitive axis, type of magnetometer (NMOR [47,48], rf-driven [24], or SERF [49]), and probed transition are listed. The bandwidth indicates the measured -3 dB point of the magnetometers' frequency response to oscillating magnetic fields. The calibration error takes into account potential temporal variation of the magnetometers' calibration over the course of Science Run 2, and is estimated based on auxiliary measurements. The rightmost column lists the estimated ratio between the effective proton spin polarization and the Landé g-factor for the magnetometer, σp/g, which depends on the atomic species and the magnetometry scheme as described in Appendix B. The σp/g value is used to relate the measured magnetic field to the signal expected from the interaction of an ALP field with proton spins. The indicated uncertainty describes the range of values from different theoretical calculations [28]. The corresponding time-shifted data along with their local standard deviation estimate are fetched from each magnetometer's rolling averaged full-rate data at a rate of 0.1 samples/s. This reduces the amount of data to process, while keeping the full timing resolution.
The step size used in the speed scan is chosen so that a single step in speed corresponds to time-shift differences of less than the down-sampled sampling period. For each speed, a lattice of directions covering the full 4π solid angle is constructed. The angular difference between adjacent directions is informed by sampling rate and speed [29] such that, as for the speed scan, a single step in direction results in time-shift differences of less than the down-sampled sampling period. With the settings used, the velocity-scanning lattice consists of 1661 points. This number scales with the down-sampled sampling rate cubed.
After time-shifting, the pulses produced by a domainwall crossing appear simultaneously as if all the magnetometers were placed at the Earth's center. This process results in a time-shifted data set for each lattice velocity on which for each time point a χ 2 -minimization is performed to estimate the domain-wall parameters. An ALP domain-wall-crossing direction and magnitude, B p , with the corresponding p-value quantifying the agreement is obtained. The p-value is evaluated as the probability of obtaining the given χ 2 value or higher from the χ 2 -minimization. The p-value is calculated using the quadrature sum of the standard deviation of the data and the uncertainty due to drifts of the calibration factors. All data points in every time-shifted data set are considered potential events, characterized by time, p-value as well as direction and magnitude B p with their associated uncertainties. The magnitude-to-uncertainty ratio of an event ζ is the ratio between this B p and its associated uncertainty.
Third, in the post-selection stage, two tests are carried out to check if a potential event is consistent with an ALP domain-wall crossing. The domain-wall model test evaluates if the observed signal amplitudes are consistent with the expected pattern of a domain-wall crossing from any possible direction. It is quantified by the aforementioned p-value. The directional-consistency test is based on the angular difference between the estimated domainwall crossing direction and the direction of the velocity corresponding to the particular time-shifted data set being analyzed. In a real domain-wall crossing event, these two directions should be aligned.
To evaluate the consistency of a potential event with a domain-wall crossing, we impose thresholds on the p-value and the angular difference normalized with respect to the angular spacing of the lattice of velocity points for that speed. The thresholds are chosen to guarantee a detection probability of 97.5% with the minimum possible false-positive probability. The false-positive analysis is performed on the background data. The true-positive analysis is performed on test data consisting of background data with randomly inserted domain wall signals as we describe below.
A single signal pattern may appear as multiple potential events in the analysis, whereas we are seeking to characterize a single underlying domain-wall crossing event. For example, a signal consistent with a domainwall crossing lasting for multiple sampling periods would appear as multiple potential events in a single timeshifted data set. Furthermore, even if such a signal lasts for only a single sampling period, corresponding potential events appear in different time-shifted data sets. Since it is assumed that domain-wall crossings occur rarely, such clusters of potential events are classified as a single "event." In order to reduce double-counting of these events, conditions are imposed. If potential events passing the thresholds occur at the same time in different time-shifted data sets or are contiguous in time, the potential event with the greatest magnitude-to-uncertainty ratio is classified as the corresponding single event.
The true-positive analysis studies the detection probability as a function of the thresholds. Multiple test data sets are created featuring domain-wall-signal patterns with random parameters by inserting Lorentzianshaped pulses into the background data of the different GNOME magnetometers. The domain-crossing events have magnitudes of B p randomly selected between 0.1 pT and 10 4 pT and durations randomly selected between 0.01 s and 10 3 s. The distributions of the these randomized parameters are chosen to be flat on a logarithmic scale. Additionally, the signals are inserted at random times with random directions. In order to simulate calibration error effects, the pulse amplitudes inserted in each magnetometer are weighted by a random factor whose range is given in Table I. The crossing velocity is also randomized within the range covered by the velocity lattice. For each inserted domain-wall-crossing event, the p-value, the normalized angular difference, and the magnitude-to-uncertainty ratio are computed. Figure 6a shows the detection probability as a function of the threshold on the lower-limit of the p-value and the threshold on the upper-limit of the normalized angular difference. We restrict the analysis in Fig. 6a to events inserted with a magnitude-to-uncertainty ratio between 5 and 10. This enables reliable determination of the true-positive detection probability without significant contamination by false positive events since the background event probability above ζ = 5 is below 0.01% in a 10 s sampling interval. Since the detection probability increases with the signal magnitude, we focus on the events below ζ = 10. The detection probability is then the number of detected events divided by the number of inserted events. The black line marks the numerically evaluated boundary of the area guaranteeing at least 97.5% detection. All points along this black line will yield the desired detection probability, so the particular choice is made to minimize the number of candidate events when applying the search algorithm to background data. These values are 0.001 for the p-value threshold and 3.5 for the directional-consistency threshold (represented as a white dot in Fig. 6a). Figure 6b shows that the detection probability is greater than 97.5% for events featuring a magnitude-to-uncertainty ratio above 5 and guarantees ≥ 95%.
Since the noise has a nonzero probability of mimicking the signal pattern expected from an ALP domain-wall crossing well enough to pass the p-value and directional consistency tests, we perform a false-positive study on background data. The analysis algorithm is applied to T b =10.7 years of time-shuffled data in order to establish the rate of events expected solely from background. Be-cause of the larger amount of background data analyzed, lower rates and larger magnitude-to-uncertainty ratios are accessible as compared to the search data. Based on the false positive study, the probability of finding one or more events in the search data above ζ, is [50], where T = 23 days is the duration of Science Run 2 and n b (ζ) is the number of candidate events found in the background data above ζ. The significance is then defined as, where erf −1 is the inverse error function. The significance is given in units of the Gaussian standard deviation which corresponds to a one-sided probability of P .
After characterizing the background for Science Run 2, the search data are analyzed. The results are represented as a solid green line in Fig. 2. For ζ > 6 only a few events were found. The event with largest magnitudeto-uncertainty ratio, ζ max was measured at 12.6 followed by additional events at 6.2 and 5.6. From Eq. (10), the significance associated with finding one or more events produced by the background featuring at least ζ max is lower than one sigma. This null result defines the sensitivity of the search and is used to set constraints on the parameter space describing ALP domain walls.
The observable rate of domain-wall crossings depends on how long GNOME was sensitive to different signal durations and magnitudes. For the evaluation of this effective time, the raw data of each magnetometer are divided into continuous segments between one to two hours duration depending on the availability of the data. The preprocessing steps are applied to each segment. Then the data are binned by taking the average in 20 s intervals. To estimate the noise in each magnetometer, the standard deviation in each binned segment is calculated to define the covariance matrix Σ s . The domain-wall magnitude, crossing with the worse case directionm, needed to produce ζ = 1 is calculated as in Ref. [29], for each bin. The matrix D ∆t contains the sensitivity axes of the magnetometers, the factor σ p /g, and the effects of the preprocessing as a function of the signal duration as described in Ref. [29]. Such prepocessing effects rely on a Lorentzian-shaped signal and give rise to the characteristic shape of Fig. 3. The effective time,T , is defined as the amount of time the network can measure a domain wall with duration ∆t and magnitude B p producing ζ ≥ 1. Monte Carlo simulations analyzing segments with inserted domain-wall encounters on raw data show a good agreement with the sensitivity estimation in Eq. (11).
Assuming that the domain-wall encounters follow Poisson statistics, a bound on the observable rate of events above ζ max with 90% confidence is set as [40], The ALP parameter space is constrained by imposing that r ≥ R 90% . The experimental constraint on the coupling constant is written as (see Eq. (A.13) in Appendix A), The signal duration can be written in terms of the mass of the hypothetical ALP particle and the specific domainwall crossing speed, ∆t = 2 √ 2 vmac . When calculating the constraints on f int , we fix the domain-wall crossing speed to the typical relative speed from the Standard Halo Model,v = 300 km s −1 [31]. In contrast to the signal duration, the pseudo-magnetic field signal depends on all parameters of the ALPs, the mass and the ratio between the coupling and symmetry breaking constants, B p = 4mac 2 ξ µ B ζ . Figure 4 is given by Eq. (13) taking ζ = 12.6. The shape of the constrained space is given by the fact thatT varies depending on the target m a and ξ. A detailed derivation of Eq. (13) is given in Appendix A.
Domain walls form when a field can monotonically vary across vacuum states; two degenerate vacua or possibly the same state if, for example, the field takes values on the 1-sphere. This is the case for ALPs which arise from the angular part of a complex scalar field [51].
The following Lagrangian terms are considered in natural units ( = c = 1 here and throughout this appendix) for a new complex scalar field φ, 4 where λ is a unitless constant and S 0 is a constant with units of energy [13]. The axion field is obtained by reparameterizing the complex field φ in Eq. (A.1) in terms of the real field S (Higgs) and a (axion), The second term in Eq. (A.1) will break the U (1) symmetry of the complex field into a discrete Z N symmetry, φ → exp(2πik/N )φ for integer k. This corresponds to the axion shift-symmetry, a → a + 2πS0 N k. The Higgs mode obtains a vacuum expectation value S → S 0 , and the axion field has degenerate vacua or ground energy states at a = πS0 N (2k + 1) for integer k. One can define the symmetry-breaking scale as f SB = S 0 /N . Reparameterizing the complex scalar field in Eq. (A.1) and setting the Higgs mode to the vacuum expectation value, the axion Lagrangian is where the axion mass is m a = N S 0 √ 2λ. This can be seen by matching the second-derivative of the cosine term at the minima to a scalar mass term.
For simplicity, a static domain wall in the yz-plane separating domains of −πf SB and +πf SB is considered. Solving the classical field equations, one finds The gradient of the field is then .
This has the full width at half maximum, Using the domain-wall solution [Eq. (A. 3)] and integrating the energy density of the domain wall over x yields the surface tension (energy per unit area) [13], Interactions observable by magnetometers involve coupling between the axion field gradient and the axialvector current of a fermionic field. For a fermion field ψ, the interaction is The axial-vector current is related to the spin S, so that the interaction Hamiltonian becomes i.e., for spin-1 / 2 particles, 1/ S = 2. Optical magnetometers operate by measuring the change in atomic energy levels between two energy states with magnetic quantum numbers differing by ∆m F . A leading magnetic field is applied to the atoms, and variations in the magnetic field along the leading field are measured. The spin coupling from Eq.
where i labels the species of fermion, σ (i) = is the projected spin coupling, η = cos θ for the angle between the axion gradient and sensitive axis θ, and f (i) int is the interaction coupling to particle i. In general, we will combine the terms into an effective ratio 2σj fint , where j now labels the magnetometer. Comparing this to the energy shift due to a magnetic field, ∆E B = g F,j µ B ∆m F B j , one obtains a relationship for a normalized pseudo-magnetic field, for ξ ≡ fSB fint , g F,j being the g-factor for the magnetometer j, and S (i) = 1/2 since we only consider coupling to spin-1 / 2 particles. Here, the normalization is such that B p is the same for all magnetometers, though each individual sensor may observe a different pseudo-magnetic field, B j .
There are two factors that must be considered when determining if axion domain walls are observable by the network: the magnitude of the signal B p and the rate of signals. Domain walls are assumed to exist in a static (or virialized) network across the galaxy through which Earth traverses. For a domain-wall velocity v, the duration of a signal is ∆t = ∆x/v. Filters and bandwidth limitations generally reduce the magnitude by a factor dependent on the signal duration, which affects the sensitivity of the network (see Appendix in Ref. [29]).
Meanwhile, if the domain walls induce a strong enough signal to be observed, but are so infrequent that one is unlikely to be found over the course of a measurement, then the network is effectively insensitive. If the energy density of domain walls across the galaxy is ρ DW , then the average rate of domain walls passing through Earth is given by wherev is the typical relative speed. The physical parameters describing the ALP domain walls (m a , f SB , and f int ) must be related to the parameters observable by the network (B p and ∆t). The energy density ρ DW and the typical relative speedv are assumed according to the observed dark matter energy density and the galactic rotation velocity, respectively.
In order to determine if a set of physical parameters is observable, the likelihood that no events are found must be constrained. This constraint defines the confidence level of the detection. The probability of observing k events given that one expects to observe µ events is given by the Poisson probability mass function, However, the network also has some detection efficiency < 1, so there could be multiple domain-wall-crossing events, but no detection. In particular, the chance of missing no events given that there were k events is (1 − ) k . For an event rate r and measurement time T , the probability that no events are detected is then A bound on the event rate R C at confidence level C is then given by demanding the probability of observing at least one event 1 − e − R C T ≥ C, likewise, one would expect to observe event rates The physical parameter space of the ALPs is constrained by demanding that r ≥ R C . Similar arguments for defining constraints can be found, e.g., Ref. [40]. The total time that the network is sensitive to the measurable parameters,T (∆t, B p ), may be less than the total measurement time. These parameters are related to the physical parameters via Eq. (A.5) and Eq. (A.10). One finds a sensitivity bound for f int in terms of m a and ξ, where B p is the sensitivity of the network and ζ is the magnitude-to-noise ratio induced by a signal with magnitude B p . The main calculation from the network data needed for this sensitivity bound isT . This is calculated by measuring the sensitivity of the network over time for different signal durations ∆t and integrating over the time during which the network is sensitive to B p = ζB p . Finally, if, after analyzing the data, no domain walls are found, Eq. (A.13) defines an exclusion region.

Appendix B: Conversion between magnetic field units and proton spin coupling
The amplitude of a signal appearing in the magnetic field data from a GNOME magnetometer due to interaction of atomic spins with an ALP field a via the linear coupling described by Eq. (1) varies based on the atomic species. In every GNOME magnetometer, the atomic vapor cell is located within a multi-layer magnetic shield made of mu-metal and, in some cases a ferrite innermost layer. Interactions of the ALP field with electron spins in the magnetic shielding material can generate a compensating magnetic field that could partially cancel the energy shift due to the interaction of the ALP field with atomic electrons in the vapor cell, as discussed in Ref. [27]. For this reason, GNOME magnetometers are most sensitive to interactions of the ALP field with nuclear spins.

Deriving spin-projection
All GNOME magnetometers active during Science Run 2 measure spin-dependent interactions of alkali atoms whose nuclei have valence protons. Thus the GNOME is primarily sensitive to spin-dependent interactions of ALP fields with proton spins. Consequently, the expected signal amplitude measured by a GNOME magnetometer due to the pseudo-magnetic field pulse from passage of Earth through an ALP domain wall must be rescaled by the ratio of the proton spin content of the probed ground-state hyperfine level(s) to their gyromagnetic ratio. Some GNOME magnetometers optically pump and probe a single ground-state hyperfine level, while others rely on the technique of spin-exchange relaxation free (SERF) magnetometry in which the spinexchange collision rate is much faster than the Larmor precession frequency [49,52,53]. For SERF magnetometers a weighted average of the ground-state Zeeman sublevels over both ground-state hyperfine levels is optically pumped and probed. Table II shows the relevant factors needed to convert the magnetic field signal recorded by GNOME magnetometers into the expected pseudo-magnetic field due to interaction of an ALP field with the proton spin. Detailed calculations are carried out in Ref. [28]. The relationship of the expectation value for total atomic angular momentum F to the nuclear spin I can be estimated based on the Russell-Saunders LS-coupling scheme: where S e is the electronic spin and L is the orbital angular momentum. GNOME magnetometers pump and probe atomic states with L = 0, which simplifies the above equation to: For L = 0 the projection of I on F is given by The above relations define the fractional spin polarization of the nucleus relative to the spin polarization of the atom: The next step is to relate σ N,F to the spin polarization of the valence proton σ p,F for a particular F . As discussed in Ref. [28], a reasonable estimate for K, Rb, and Cs nuclei can be obtained from the nuclear shell model II: Fractional proton spin polarization σp,F , Landé g-factors gF , and their ratios for the ground state hyperfine levels used in GNOME, and the weighted average of these values across both hyperfine levels ( σp hf / g hf ) applicable to SERF magnetometers in the low-spin-polarization limit. The estimates are based on the single-particle Schmidt model for nuclear spin [54] and the Russell-Saunders scheme for the atomic states. Uncertainties in the values for σp,F describe the range of different results from calculations based on the Schmidt model, semi-empirical models [55,56], and large-scale nuclear shell model calculations where available [57][58][59]. The uncertainties in σp,F and σp hf are one-sided because alternative methods to the Schmidt model generally predict smaller absolute values of the proton spin polarization. See Ref. [28] for further details.

Atom (state)
σp,F gF σp,F /gF σp hf / g hf by assuming that the nuclear spin I is due to the orbital motion and intrinsic spin of only the valence nucleon and that the spin and orbital angular momenta of all other nucleons sum to zero. This is the assumption of the Schmidt or single-particle model [54]. In the Schmidt model, the nuclear spin I is generated by a combination of the valence nucleon spin (S p ) and the valence nucleon orbital angular momentum , so that we have where it is assumed that the valence nucleon is in a welldefined state of and S p , and σ p is defined to be the fractional proton spin polarization for a given nucleus [28].
For comparison between GNOME magnetometers using different atomic species, it is essential to evaluate the uncertainty in the estimate of σ p,F based on the Schmidt model. To estimate this uncertainty, we compare calculations of σ p,F based on the Schmidt model to the results of the semi-empirical calculations described in Refs. [55,56] and to the results of detailed nuclear shell-model calculations where available [57][58][59]. Conservatively, we assign the uncertainty in σ p,F to be given by the full range (maximum to minimum) of the values of σ p,F calculated by these various methods. It turns out that in each considered case, the estimate based on the Schmidt model gives the largest value of |σ p,F |, causing the theoretical uncertainties in estimates of σ p,F to be one-sided as shown in Table II. Further details are discussed in Ref. [28].

SERF magnetometers
SERF magnetometers operate in a regime where the Larmor frequency is small compared to the spin-exchange rate, so that the rapid spin-exchange locks together the expectation values of the angular momentum projection M F in both ground-state hyperfine levels of the alkali atom. Because the Landé g-factors g F in the two groundstate hyperfine levels have nearly equal magnitudes but opposite signs, the magnitude of the effective Landé gfactor in a SERF magnetometer, g hf , is reduced compared to that in optical atomic magnetometers where a single ground-state hyperfine level is probed.
To calculate the effective Landé g-factor averaged over hyperfine levels, g hf , for a SERF magnetometer, it is instructive to consider the equation describing the magnetic torque on an alkali atom, where g s ≈ 2 is the electron g-factor and we have ignored the contribution of the nuclear magnetic moment.
In the SERF regime, where the alkali vapor is in spinexchange equilibrium, the populations of the Zeeman sublevels correspond to the spin-temperature distribution [60] described by a density matrix in the Zeeman basis given by [61,62] ρ = Ce β·F , (B.7) where C is a normalization constant and β is the spintemperature vector defined to point in the direction of the spin polarization P with magnitude β = ln 1 + P 1 − P . (B.8)