Search for topological defect dark matter with a global network of optical magnetometers

Ultralight bosons such as axion-like particles are viable candidates for dark matter. They can form stable, macroscopic field configurations in the form of topological defects that could concentrate the dark matter density into many distinct, compact spatial regions that are small compared with the Galaxy but much larger than the Earth. Here we report the results of the search for transient signals from the domain walls of axion-like particles by using the global network of optical magnetometers for exotic (GNOME) physics searches. We search the data, consisting of correlated measurements from optical atomic magnetometers located in laboratories all over the world, for patterns of signals propagating through the network consistent with domain walls. The analysis of these data from a continuous month-long operation of GNOME finds no statistically significant signals, thus placing experimental constraints on such dark matter scenarios.

The nature of dark matter, an invisible substance comprising over 80% of the mass of the universe [1,2], 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 [3].One of the leading hypotheses is that dark matter consists of ultralight bosons such as axions [4] or axion-like particles (ALPs) [5][6][7].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 manifest as stable, macroscopic field configurations in the form of topological defects [8][9][10] or composite objects bound together by self-interactions such as boson stars [11,12].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 [9,10,12].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 [8,9].We note that while some models suggest that axion domain walls cannot survive to the present epoch [13][14][15], there do exist a number of ALP models demonstrating the theoretical possibility that ALP domain walls or composite dark matter objects with similar characteristics [12,16,17] can survive to modern times [18][19][20] and have the characteristics of cold dark matter [9,10,21].
Since ALPs can interact with atomic spins [3], the passage of Earth through an ALP domain wall affects atomic spins similarly to a transient magnetic field pulse [9,12].Considering a linear coupling between the ALP field gradient ∇a(r, t) and atomic spin S, the interaction Hamiltonian 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 [22], 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 [3,5].
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, even though it couples to the particle spin rather than the magnetic moment, can be treated as a "pseudo-magnetic field" as it causes energy shifts of Zeeman sublevels.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 [3].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-toencounter.We assume that they are described by the Standard Halo Model for virialized dark matter [23].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 [24].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 Sec.I of the Supplementary Information and Refs.[9,12,22].
The Global Network of Optical Magnetometers for Exotic physics searches (GNOME) is a worldwide network searching for correlated signals heralding beyondthe-Standard-Model physics which currently consists of more than a dozen optical atomic magnetometers, with stations (each with a magnetometer and supporting devices) in Europe, North America, Asia, the Middle East, and Australia.A schematic of a domain-wall encounter with the GNOME is shown in Fig. 1 acquisition systems [25], 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 the Methods and Ref. [26] (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 [27] while retaining sensitivity to exotic spin-couplings between ALP dark matter and atomic nuclei.
If the ALP field only couples to electron spins, interaction between the ALP field and the magnetic shield will reduce ALP-induced signal amplitudes in each magnetometer by roughly the magnetic shielding factors of 10 6 -10 7 as discussed in Ref. [28].Therefore in the present work we only consider interactions between ALP fields and atomic nuclei.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 Sec.II of the Supplementary Information and Ref. [29]).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 Sec.II of the Supplementary Information.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 the Methods and Ref. [30].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 [12] and bursts of exotic low-mass fields (ELFs) from cataclysmic astrophysical events [31] are presently underway.
Here we report the 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 dark-matter objects based on similar ideas were carried out using atomic clock networks [10,23,32,33], and there are a number of experimental proposals utilizing other sensor networks [34][35][36][37].All of these networks are sensitive to bosonic dark matter with a scalar coupling to Standard Model particles [3].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 [3] 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.

SEARCH FOR ALP DOMAIN-WALL SIGNATURES
There have been four GNOME Science Runs between 2017 and 2020 as discussed in the Methods.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 [30].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 component of the relative velocity between the domain wall and the Earth that is perpendicular to the domain-wall plane, 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 veloc-ity 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. [30].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 [30].The domain-wall 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 the Methods and Ref. [30].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.The thresholds on these tests are chosen to guarantee an overall detection efficiency ≥ 95% for the search algorithm, considering both the incomplete velocity lattice coverage and the detection probability (see Fig. 6).
The search data are analyzed for domain-wall encounters using the algorithm presented in Ref. [30].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 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. Rare domain-wall-crossing events producing signals below a magnitude-to-uncertainty ratio of 12.6 are indistinguishable from the background.Therefore, we base constraints on ALP parameters on the absence of any detection above the "loudest event" in a manner similar to that described, for example, in Ref. [38].
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 [30]; we conservatively estimate the network sensitivity assuming the domain wall comes from the least-sensitive direction.Figure 3 shows the active time T (∆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 The significance 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. ( 5)].The right axis shows the normalized number of events over a period of a year.The event with greatest magnitude-to-uncertainty ratio is found at 12.6.
ratio of ζ = B p /B p .The active time, T (∆t, B p ), can be used to constrain ALP domain-wall parameter space as discussed in Sec.I of the Supplementary Information.
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 [38].Note that since T 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%.

CONSTRAINTS ON ALP DOMAIN WALLS
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 exis- tence 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 speed v, 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-to-uncertainty 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 separation L can be related through the ALP domain wall model described in Refs.[9,22].A full derivation of how observable parameters are related to ALP parameters is given in Sec.I of the Supplementary Information.
The colored region in Fig. 4a describes the symmetry breaking scales up to which GNOME was sensitive with 90% confidence.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-touncertainty 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 for different ratios ξ ≡ f SB /f int .Because B p ∝ m a [see Eq. (SI.10) in Sec.I of the Supplementary Information], 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.(11) and (12)].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 ALP-domainwall parameter space during Science Run 2 can be excluded.
Our experiment explores ALP parameter space up to f int ≈ 4×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 [39].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 [40,41].
The parameter space for f int and m a explored in this search is well outside typical predictions for the quantum chromodynamics (QCD) axion (see, for example, Refs.[42,43]).However, for ALPs a vast array of possibilities for the generation of ALP masses and couplings are opened by a variety of beyond-the-standard-model theories, meaning that the values of f int and m a explored in our search are theoretically possible (see, for example, the reviews [44,45] and references therein).
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.The intermittent operation of some magnetometers due to technical difficulties during Science Runs 1-3 made it difficult to search for signals persisting 1 hour.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 [46,47] 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 [12], relaxion halos [48], and bursts of ex-otic low-mass fields from black-hole mergers [31].
In terms of the data-analysis algorithm used to search for ALP domain walls, recent studies [49] have considered a possible back-action that the Earth may have on a domain wall when certain interactions are significant; namely up-to-quadratic coupling terms between a scalar field and fermions.In contrast to Ref. [49], the present work analyzes a completely different interaction, namely a linear coupling between a pseudoscalar field and fermion spins, which produces no significant backaction effect.Regardless, it would be worthwhile to consider interactions generating similar back-action effects of the Earth on domain walls and the ALP field in later analysis.Further, 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) −3 (where ∆t is given by Eq. ( 9)).This makes the algorithm computationally intensive.We are investigating a variety of analysis approaches, such as machine-learning-based algorithms, to address these issues.
Table I.Characteristics of the magnetometers active during Science Run 2. The station name, location in longitude and latitude, orientation of the sensitive axis, type of magnetometer (NMOR [50,51], rf-driven [26], or SERF [52]), 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 Sec.II of the Supplementary Information.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 [29].sensors, including accelerometers, gyroscopes, and unshielded magnetometers, to measure local perturbations that could mimic a dark matter signal.Suspicious data are flagged [26] and discarded during the analysis.

Location
The number of active GNOME magnetometers during the four Science Runs and the combined network noise as defined in Ref. [30] is shown as a function of time in Fig. 5.Although Science Run 4 was carried out over a longer period of time than Science Run 2, it featured poorer noise characteristics and consistency of operation compared 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.The number of active magnetometers during Science Runs 1 and 3 was often less than four, which in insufficient to characterize a domainwall crossing.We thus present analysis efforts on Science Run 2 data.
Here we provide more details on the analysis procedure.It is composed of three stages to identify events likely to be produced by ALP domain-wall crossings: preprocessing, velocity scanning, and post-selection [30].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 firstorder 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 uncertainty resulting from the calibration error is later used to determine agreement with the domain-wall model, but not in the magnitude-to-uncertainty ratio estimate resulting from the model, since the calibration error affects signal and noise in the same way.
Second, in the velocity-scanning stage, the data from the individual magnetometers are time-shifted according to different relative velocities between Earth and the ALP domain walls.In order to sample 97.5% of the velocity probability distribution, a scan of the speeds from 53.7 to 770 km/s with directions covering the full 4π solid angle is chosen, therefore the domain walls can take any orientation with respect to the Earth movement.Note that this distribution considers just the observable perpendicular component of the relative domain-wall velocity and neglects the orbital motion of the Earth around the Sun.For low relative velocities, both the time between signals at different magnetometers as well as the signal duration diverge.So the velocity range is determined by the chosen 97.5% coverage and the maximum relative speed of domain walls travelling at the galactic escape speed.
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 [30] 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 pvalue 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.
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 [53].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 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%.This results in an overall detection efficiency ≥ 95% for the search algorithm, considering both the incomplete velocity lattice coverage and the detection probability.
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.Because 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 [54], 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, S = − 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. ( 5), 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 direction m, needed to produce ζ = 1 is calculated as in Ref. [30], 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. [30].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. (6).
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 [38], The domain-wall thickness is determined by the ALP mass, and is on the order of the ALP reduced Compton wavelength λ a [22], 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 [24], 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 [9], The average domain-wall separation is then approximated by which results in the average domain-wall encounter rate, We assume the typical relative domain-wall speed to be equal to the galactic rotation speed of Earth.
The ALP parameter space is constrained by imposing that r ≥ R 90% .The experimental constraint on the coupling constant is written as (see Eq. (SI.13) in Sec.I of the Supplementary Information), 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 [23].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 that T varies depending on the target m a and ξ.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 [55].
The following Lagrangian terms are considered in natural units ( = c = 1 here and throughout the supplementary information) for a new complex scalar field φ, where λ is a unitless constant and S 0 is a constant with units of energy [9].Other references may use a minus sign in front of the S N 0 term, which results in a similar potential, up to a phase.The end result of Eq. (SI.1) is that the axion potential will have a maximum at zero, while the minus-convention will have a vacuum at zero.
The axion field is obtained by reparameterizing the complex field φ in Eq. (SI.1) in terms of the real field S (Higgs) and a (axion), The second term in Eq. (SI.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. (SI.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 (SI. 3) The gradient of the field is then This has the full width at half maximum, Using the domain-wall solution [Eq.(SI. 3)] and integrating the energy density of the domain wall over x yields the surface tension (energy per unit area) [9], 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 .The sensitive axis of a magnetometer is either defined by a leading magnetic field applied to the atoms, or, in the case of zero-field magnetometers, by the axis orthogonal to the plane defined by the propagation directions of the pump and probe light.Variations of the magnetic field along the sensitive axis are measured.The spin coupling from Eq. (SI.8) can induce a similar shift in energy levels to a magnetic field.The maximum energy shift is determined by plugging the largest gradient from Eq. (SI.4) into Eq.(SI.8), ∆E = i∈e,p,n,. . .
where i labels the species of fermion, is the projected spin coupling, η = cos θ for the angle between the axion gradient and sensitive axis θ, and f 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, g F,j B j σ j η j = 4 µ B m a ξ ≡ B p , (SI.10) 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. [30]).
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 where v 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 speed v 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 (SI.12) 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. [38].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.(SI.5) and Eq.(SI.10).One finds a sensitivity bound for f int in terms of m a and ξ,

II. 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. [28].For this reason, GNOME magnetometers are most sensitive to interactions of the ALP field with nuclear spins.

A. 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 [52,56,57].For SERF magnetometers a weighted average of the ground-state Zeeman sublevels over both ground-state hyperfine levels is optically pumped and probed.
Table SI.1 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. [29].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: F . (SI.15) 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: . (SI.17) 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. [29], a reasonable estimate for K, Rb, and Cs nuclei can be obtained from the nuclear shell model 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 [58].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 [29].
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.[59,60] and to the results of detailed nuclear shell-model calculations where available [61][62][63].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 SI. 1.Further details are discussed in Ref. [29].

B. 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 [64] described by a density matrix in the Zeeman basis given by [65,66] ρ = Ce β•F , (SI.20 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 . (SI.21) In the low-spin-polarization limit, abundance, so spin-exchange collisions average over both ground-state hyperfine levels of all three species.Taking into account the relative abundances of the different atomic species at the cell temperature of ≈ 150 • C (≈ 9% 39 K, ≈ 65.5% 85 Rb, ≈ 25.5% 87 Rb), we find that g hf ≈ 0.193, σ p hf = −0.073+0.010 −0.000 , and thus σ p hf / g hf = −0.38 +0.05 −0.00 for the Hefei magnetometer.

Lewisburg magnetometer
The Lewisburg GNOME station employs a SERF magnetometer in a closed-loop, two-beam configuration.The vapor cell contains only 87 Rb atoms.The Lewisburg SERF magnetometer operates with a spin-polarization P ≈ 0.5, outside the low-spin-polarization regime.Discussions of the high-polarization regime are given in Refs.[66,68].For a nucleus with I = 3/2, the effective Landé g-factor is given by g hf = g s 1 + P 2 6 + 2P

Figure 1 .
Figure 1.Visualization of an ALP domain-wall crossing.(a) The image shows the Earth together with the position and sensitive axes of the GNOME magnetometers during Science Run 2. Position and sensitive axes are show as red arrows.The crossing direction of the domain wall is represented as a black arrow (see TableI).(b) Simulation of the signals expected to be observed from a domain-wall crossing at the different magnetometers conforming the network.

Figure 2 .
Figure2.Significance of the Search events.The blue dashed line represents the cumulative number events expected from the background in the twenty-three days of data from Science Run 2. 10.7 years of time-shuffled data are used to evaluate the background.Such duration is an arbitrary choice but is sufficiently long to characterize the background.The number of candidate events measured in the background data is re-scaled to the duration of Science Run 2. The solid green line represents the cumulative number of events measured in Science Run 2. The red crosses indicate at which magnitude-to-uncertainty ratio new events are found in the search data.The upper axis indicates the statistical significance in units of Gaussian standard deviations of finding one event in the search data.The significance 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. (5)].The right axis shows the normalized number of events over a period of a year.The event with greatest magnitude-to-uncertainty ratio is found at 12.6.

Figure 3 .
Figure3.Sensitivity of the GNOME network to domain walls.Amount of time T , indicated in color, that the GNOME had a normalized pseudo-magnetic field magnitude sensitivity above B p (i.e., the domain wall would induce a magnitude-touncertainty ratio of at least one) for domain walls with a given duration ∆t (defined as the FWHM of a Lorentzian signal) throughout Science Run 2. The upper axis shows the range of ALP masses to which GNOME is sensitive [see Eq. (9)].The characteristic shape of the sensitive region is a result of the filtering and averaging of the raw data as described in the Methods.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 active GNOME magnetometers recording data and their background noise change.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.

Figure 4 .
Figure 4. Bounds on the ALP parameter space.The bounds are drawn from the presented analysis of Science Run 2 with 90% confidence level.The relationship between ALP theory parameters and measured quantities are discussed in Sec.I of the Supplementary Information.(a) In color, upper bound on the the interaction scale for axion-nucleon coupling, fint, to which the GNOME was sensitive as a function of ma and the ratio between the symmetry-breaking and interaction scales ξ ≡ fSB/fint.Dashed horizontal lines highlight the cross section used in (b) with the respective color.(b) Cross-sections of the excluded parameter volume in (a) for different ratios ξ.We note that domain walls may not be the only form of dark matter so ρDW < 0.4 GeV/cm 3 .If the domain-wall energy density is significantly smaller, this would affect the bounds shown here.

Figure 5 .
Figure 5. Summary of the GNOME performance during the four Science Runs from 2017 to 2020.The raw magnetometer data are averaged for 20 s and their standard deviation is calculated over a minimum of one and a maximum of two hours segments depending on the availability of continuous data segments.For each binned point, the combined network noise considering the worse case domain-wall crossing direction is evaluated as defined in Ref. [30].(a) One-day rolling average of the number of active sensors.(b) Multi-colored solid line represents the one-day rolling average of the combined network noise and the multi-colored dashes show the noise of the individual sampled segments.The data are preprocessed with the same filters used for the analysis.The number of magnetometers active is indicated by the color of the line and dashes.

Figure 6 .
Figure 6.Summary of the true-positive analysis results.(a)shows the detection probability for domain-wall-crossing event with randomized parameters (as discussed in the text) as a function of p-value and directional-consistency thresholds.The inserted events have a magnitude-to-uncertainty ratio between 5 and 10.The black line indicates the combination of parameters corresponding to a 97.5% detection probability.The white dot indicates the particular thresholds chosen for the analysis.(b) Shows the mean detection probability reached for different magnitude-to-uncertainty ratios for the chosen thresholds.

2 [F
(F + 1) + I(I + 1) − S e (S e + 1)] .(SI.16) σ p,F = S p • I I(I + 1) σ N,F , = S p (S p + 1) + I(I + 1) − ( + 1) 2I(I + 1) σ N,F , = σ p σ N,F , (SI.18) 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 is T .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. (SI.13) defines an exclusion region.