Searching for Magnetic Monopoles with the Earth's Magnetic Field

Magnetic monopoles have long been predicted in theory and could exist as a stable object in our universe. As they move around in galaxies, magnetic monopoles could be captured by astrophysical objects like stars and planets. Here, we provide a novel method to search for magnetic monopoles by detecting the monopole moment of the Earth's magnetic field. Using over six years of public geomagnetic field data obtained by the Swarm satellites, we apply Gauss's law to measure the total magnetic flux, which is proportional to the total magnetic charge inside the Earth. To account for the secular variation of satellite altitudes, we define an altitude-rescaled magnetic flux to reduce the dominant magnetic dipole contribution. The measured magnetic flux is consistent with the existing magnetic field model that does not contain a monopole moment term. We therefore set an upper limit on the magnetic field strength at Earth's surface from magnetic monopoles to be $|B_{\rm m}|<0.13$ nT at 95% confidence level, which is less than $2\times 10^{-6}$ of Earth's magnetic field strength. This constrains the abundance of magnetically-charged objects, including magnetic black holes with large magnetic charges.


Introduction
As observed by Dirac close to a century ago [1], the magnetic monopole is a fascinating physical object that could elegantly explain the quantization of electric charges in Nature. Ever since, physicists have been studying magnetic monopoles both from theoretical and experimental directions. On the theory side, the discovery of Polyakov-'t Hooft monopoles in 1974 [2,3] has been applied to Grand Unified Theories (GUT) [4], which predict the GUT monopole mass around 10 17 GeV. For heavier masses above the Planck mass scale, magnetically charged black holes have long been proposed, which can have masses proportional to their magnetic charges [5][6][7]. Various experimental methods have been adopted to search for magnetic monopoles, e.g., detecting the quantized jump in magnetic flux when monopoles pass a superconducting quantum interference device [8] and searching for the Cherenkov light generated when the accelerated monopoles pass the large IceCube detector [9]. Magnetic monopoles could also be captured by stars and planets including our Earth, and their annihilations can produce detectable neutrinos and/or heat [10][11][12] (see also [12][13][14][15] for further constraints). In this article, we present a new way to search for magnetic monopoles by measuring the magnetic monopole moment of the Earth.
Ever since the classical 1839 paper by Gauss [16] (see [17] for English translation), measurements indicate the Earth's magnetic field is dominantly dipole-like with subleading contributions from higher moments. Although the general spherical harmonics formula developed by Gauss contains the monopole moment [16], it has been assumed to be zero because the traditional Maxwell equations lack magnetically charged objects. As magnetic monopoles have a solid theoretical motivation and a chance to be captured inside the Earth, the Earth's magnetic monopole moment could be non-zero. Performing a measurement of Earth's monopole moment is therefore a straightforward and worthwhile method to search for magnetic monopoles.
To measure different magnetic moments including the monopole moment of the Earth, one could perform a global fit like the one used for the International Geomagnetic Reference Field [18]. In this article, we adopt a simpler approach and apply Gauss's law to Earth's magnetic field to obtain the total magnetic flux, proportional to the enclosed magnetic charge. Ideally, one needs to perform synchronous measurements of magnetic field vectors at all points on a sphere surrounding Earth to obtain a precise value for the magnetic flux. In practice, ground-based measurements only cover a small fraction of the total area, although the measurements can be made synchronously. For satellite-based measurements, the satellites visit different locations at different times but do provide good whole-sky coverage. Knowing that the secular variation of magnetic fields is in general much longer than the satellite orbiting time (order of hours for one orbit and order of days for whole-sky coverage), we will use the monopole moment. More specifically, we will analyze the publicly available data by the Swarm satellites [19,20] spanning more than 6 years of observation. For equations, we use natural units with = c = ε 0 = 1, although numerical values are expressed in the International System of Units (SI).
2 Measuring the monopole moment 2.1 Measuring magnetic charges via Gauss's law For a single object or a group of objects with a total magnetic charge of Q at the center of the Earth, the Earth's magnetic field has a monopole moment of B m (r) = Q h 4 π r 2r = Q 2 e r 2r , where e = √ 4πα with α ≈ 1/137 as the fine-structure constant and h = 2π/e ≈ 68.5 e ≈ 21 is the magnetic coupling. Q = 1 is the minimal magnetic charge, corresponding to the Dirac quantization condition with e h = 2π [1]. Numerically, B m ≈ 0.082 nT × R ⊕ r 2 Q 10 19 , where R ⊕ ≈ 6371.2 km is the average radius of the Earth. For comparison, the measured Earth surface intensity has a magnitude of up to ≈ 65000 nT.
To measure the magnetic charge, one could adopt Gauss's law¸B(r) · dA = Q h. This requires a full-sky measurement of the magnetic vector field. For convenience, one could choose the manifold to be a sphere of radius R centered on Earth. Then, dA = R 2n dΩ withn as a unit surface vector pointing outward and dΩ = sin θdθdφ in spherical coordinates. For magnetic monopole objects, Here, we have defined a solidangle averaged magnetic field B, which is simply the amplitude of the monopole magnetic field at radius R. Its sign matches the sign of the magnetic charge. All higher multiple moments beyond B m do not contribute to B.
In practice, the measurement of magnetic field is not performed at a uniform radius-the Swarm satellite orbits have a variation of O(1%) during one orbit and decay over time. Thus, it is not possible to integrate the magnetic flux along a perfectly spherical closed manifold, and the surface's normal vectorn will not match the radial coordinate unit vectorr. So, a numerical integration of´B(r, θ, φ) ·r dΩ will not be zero, even in the absence of a monopole term (for Swarm's orbital parameters, B −70 nT; see the Supplemental Material). To suppress this measurement-induced dipole contribution, we use the following modified Gauss law to measure the magnetic field from the monopole charge Here, r(θ, φ) is the radius of the magnetic measurement at different angular directions and R ref is a fixed reference radius. For the dipole component, this is formally equivalent to integrating on a perfectly spherical surface at r = R ref , son =r and the dipole component contributes zero to the above quantity. Note that the Earth's higher-moment magnetic fields have non-zero contributions to the quantity B because the the higher moments scale with higher powers of r. For instance, the quadrupole moment has a magnitude of O(10%) of the dipole moment, and contributes around 0.5 nT for B using Swarm's orbit. Therefore, the r 3 scaling in Eq. (1) is practically useful to improve the sensitivity of searching for the monopole moment because it reduces contributions to B from the dipole and higher moments while preserving the monopole signal.

Swarm satellites
The European Space Agency Swarm mission launched in late 2013. It consists of three satellites in low-Earth near-polar orbits each measuring the local magnetic field vector once per second. Two satellites Swarm A and Swarm C orbit in a similar plane with only about 1.5 degrees east-west separation and initial altitude 460 km, while Swarm B is separated from the other two and has initial altitude 530 km. The altitudes vary with time due to atmospheric drag (causing an overall decrease with time) and orbital eccentricity (changing the satellite radius by less than 1% within each orbit). Because each orbit takes around 90 minutes and the Earth completes a rotation every 24 hours, each satellite can be thought of as sweeping out measurements at all latitudes in about 32 different longitudes per day. The orbital inclinations of the satellites are 87.4 • for Swarm A and C and 88 • for Swarm B, so that no measurements are taken within 2 • of the geographic poles. We use the VirES architecture to access the data via the Python package viresclient [21]. Magnetic field measurements are used from the Swarm L1b 1 Hz data product. When selecting data, we take the following considerations. Due to a sensor failure in Swarm C, we only use the data from the other two satellites. Bad data is removed using the Flags F filter, and only non-zero vector data for which the magnetic activity level Kp 3 are included to reduce the impact of the Sun (following the procedure in Refs. [22,23]). Unlike other analyses, we do not impose a cut on the Sun's elevation angle above the horizon, as this would prevent us from obtaining full 4π sky coverage (particularly near the poles).

Data analysis and model comparison
To calculate the flux, the measurements are binned into a × a degree angular patches and dday-long time bins. If the satellites do not cover all a × a degree patches within a given d day time bin (for example, due to the selection cuts), that bin is not used. Because one satellite covers about 32 different longitudes per day, it takes at least d 360/(32a) days for one satellite to cover every patch, although generally it takes longer due to selection cuts and the precise orbital path taken. Because it takes about 45 minutes for a satellite to transit from the north to the south pole, it spends about a/2 minutes in each angular patch. With 1 Hz measurement cadence, it takes about 30a measurements per angular patch in a single orbit. Due to the orbital inclinations, the patches closest to the poles are imposed to have size 3 We index the individual measurements, angular patches, and time bins by i, j, and k respectively. To estimate the magnetic flux for a given time bin, we first rescale the measurements of the magnetic field in the radial direction B r to a common radius R ref = 6850 km, the average radius of both satellites for the mission duration thus far, as Then, we average every rescaled measurement of the magnetic field in the radial direction within each angular patch with N j the number of measurements in patch j. These patch averages are then summed together, weighted by the area of each patch, to obtain B = j B r j w j , with w j = (4π) −1´p atch j sin θdθdφ. To obtain an estimate of the statistical error, the squared standard error of the mean is calculated for each angular patch as We ignore the errors in each measurement, which are 0.5 nT, as trivial compared to the errors from having different measurements at different positions and times within each patch. These errors are weighted-summed in quadrature to obtain the squared error of the total flux σ 2 = j σ 2 j w 2 j . For each time bin k, this procedure yields an estimate for B k and its statistical uncertainty σ k . The error for each time bin scales as σ k ∝ a.
The approach outlined above contains several systematic errors, so B cannot be identified with the magnetic flux B. These include the imperfect spherical coverage which depends on the orbital path taken, the finite angular patch size used to numerically integrate (1), the finite time bin size when compared against the non-static nature of the fields, and the other magnetic field components besides the dominant dipole component which are not cancelled using the r 3 rescaling in (1). Many of these systematic errors can be reduced by subtracting the data measurement B dat k from the value B mod k predicted by a model of Earth's field (that does not contain a monopole moment term) evaluated at the same coordinates (r, θ, φ, t) of each satellite measurement.
Models for Earth's magnetic field contain contributions from the dominant dipole moment term, higher moments coming from both the core and the lithosphere, and external sources including the ionosphere and magnetosphere field. The internal components are expressed as a spherical harmonic expansion with time-dependent coefficients. We use the CHAOS-7 model [23], which is fitted through 20 July 2020.
To search for a monopole, we compare how well the data agrees with the model, the latter of which contains no monopole moment term. First, we calculate the difference in each time bin Where the χ 2 difference between the time-binned measurements B dif k and the mean B dif indicates additional variability, we enlarge the total error σ by the Birge ratio r B = χ 2 /(N k − 1) [24,25] or σ = max[1, r B ] σ dif , with N k the number of time bins. When r B > 1, this enlarged error is essentially equivalent to the standard error of the mean of B dif , i.e., (σ sem Note, σ sem is irreducible in the limit a → 0 provided a is small enough for each time bin's measured B dif k to stabilize numerically (but not so small that there are empty angular patches without measurements). It can be interpreted as including additional variability effects not included in the purely statistical error, including different orbital paths or acceptance cuts in each time bin. 1 In the absence of a signal, we set a bound on the monopole charge by injecting an artificial signal B m (r) into the data. A new B dif is calculated by comparing the data+injected signal to the model, and the bound on the injected signal size is set where B dif differs from zero at 95% confidence level.
1 If a significant signal were detected, the variability could also be due to the time evolution of the monopole moment as charged objects are captured by Earth, an effect we neglect in the absence of a signal. To provide a sense of how much each component of Earth's field contributes to our estimate of the flux, the various contributions are summarized in Table 1, showing the prediction for B mod (the time average of the model alone) from various components of the CHAOS-7 model [23] using a = 2, d = 180, and Kp 3 (see also the shaded bands in Fig. 1). The core component is separated into the dipole and higher moments. The second column shows the calculation under idealized conditions where all measurements are taken at the same time and on an equal-radius sphere, with (θ, φ) still matching the satellite measurements. These terms formally go to zero when a → 0, so they can be interpreted as the systematic error from the discretization of the angular integral. The crust and external components contribute negligibly (see also Fig. 4 in the Supplemental Material).
The third column uses the actual variation in the satellite altitude and rescales the field by (r/R ref ) 3 as in Eq. (1), but still removes time dependence. This shows the error introduced by the r 3 rescaling. Note, the dipole results match whether taken on a sphere or rescaled because the two operations are mathematically equivalent. The dipole dominates the full model's error because its overall contribution to the field is larger, and its value is consistent with zero. On the other hand, the higher core moments have much larger B mod when calculated at the satellite radius (compared with both the dipole and higher moment calculation on sphere) because higher moments are proportional to r −n with n > 3. Thus, higher moments provide the dominant contribution to B mod . Unlike the on-sphere calculation, there is relatively little suppression of B mod as a → 0 because of the saturated contribution from the higher moments when evaluated at the satellite radius. Meanwhile, the mean external component evaluated at 2 The lithosphere component is time independent by model assumption.
satellite radius turns out to have a large significance compared with its statistical error, but it contributes only a small amount to the full model prediction.
Finally, the last column shows the inclusion of time-dependence to the model, which is most relevant to our calculation. For most components, time dependence is a subdominant effect in that there is minimal change between the third and fourth columns. The exception is the external component, which is expected to have significantly more time variation than the terrestrial components on timescales relevant to the data.

Constraints on magnetic monopole objects
One can apply the constraints on the total net magnetic charge |Q net | < Q max into various fundamental theories that predict magnetic monopole objects with different charge-to-mass ratios, which we parametrize by q ≡ Q/Q ext with Q ext = e M/[cos(θ W ) and e the electric coupling constant. A magnetic black hole has q 1 (saturated to equality in the extremal limit) and M > M pl . Because magnetic black holes can efficiently Hawking radiate into electrons if their temperature is sufficiently large, they satisfy q ∼ 1 whenever M 10 17 g [7,12,15], but can take on any q 1 at larger masses. Conversely, a monopole particle has q > 1 according to the weak gravity conjecture [26]. While a GUT monopole with Q = 2 and mass M 10 17 GeV/c 2 has q 1300, a gravitating composite monopole object with a large magnetic charge could have much larger mass with both small and large q [27]. Therefore, we treat q and M as free model parameters to set limits.
The capture rate of magnetic monopoles is estimated to be C cap ≈ πR 2 ⊕ 4π F . Here, F ≈ local monopole energy density ρ as a fraction of the local dark matter density [12]. 3 Using the Earth's lifetime τ ⊕ ≈ 1.4 × 10 17 s, the number of captured monopoles is N = C cap τ ⊕ ≈ 3800 f (10 15 g/M ). When N ≥ 1, the net charge of captured monopoles is Q net √ N Q. Thus, the constraint Q net < Q max can be expressed as a limit on the local density of monopoles: This is depicted in Fig. 2. Also shown is the Parker bound [28,29] derived from M31 in [12], f 6 × 10 −3 q −2 , which disappears when q < 0.08. The limit presented here is complimentary to other limits-for example from gas heating and white dwarf destruction [15]-in that it is a direct measurement as opposed to an inference from difficult-to-model astrophysical systems.

Discussion
The Earth's magnetic field provides an interesting way to search for new physics. Using satellitebased measurements, we have shown it is possible to estimate a bound on the magnetic flux with a modified version of Gauss's law. Even without any knowledge of the higher multiple moments, the monopole moment can be constrained at the level of O(nT) (the value for B when r 3 rescaling as in (1) is employed), and comparing to a model of the higher moments allows for bounds an order of magnitude stronger. This suggests that future model fits should include the possibility of a monopole term. Combining a global fitting technique with data from more observatories would provide the best achievable bound with present technology. We encourage future work in this direction to build on the work presented here. Mars, which also has satellite-based measurements of its magnetic field [30], could be another interesting target for future work.  Fig. 1, but B is calculated without the r 3 rescaling in Eq. (1). All components other than the dipole have |B| < 3 nT, far off the range of the plot but a bit larger compared to |B| < 1 nT when r 3 rescaling is employed. The central value for B dif is close to zero as in Fig. 1, while the error bar size is larger.
ponents and therefore do not play a significant role in the analysis, although they are always included. In addition, in Table 2 we show B dif for several combinations of d, a, and Kp cut other than the benchmark combination used in the main text. When the amount of data per angular patch is large enough, the results are insensitive to the choice of a in the presented range when d and the Kp cut are fixed, indicating stability in the numerical integration. The errors (including Birge ratio factors) are also insensitive to a as explained in Section 2.3. When d varies (with other variables fixed), the result fluctuates slightly because different amounts of data are included (2362 days of data are included in total, which is not always divisible by d, and the remainder are not used). Also, when d becomes too small, not all time bins meet the 100% sky coverage criteria, resulting in further fluctuations and larger error estimates. For results quoted in Table 2, only the a = 1, d = 90, and Kp 2 cut is affected by this sky coverage criteria, and only one time bin is excluded on these grounds. We also find the results to be largely insensitive to variation of the Kp cut. On the other hand, we do not present results for a < 1 because B dif and its error becomes numerically unstable,  Additionally, the data from both satellites Swarm A and B are largely consistent with each other. Using a = 2 and d = 180, the calculated B dif in each time bin for each individual satellite are plotted in Fig. 5. There, rather than using the average satellite radius for both