Earth-Catalyzed Detection of Magnetic Inelastic Dark Matter with Photons in Large Underground Detectors

Inelastic dark matter with moderate splittings, $\mathcal{O}({\rm few} \; {\rm to} \; 150)$ keV, can upscatter to an excited state in the Earth, with the excited state subsequently decaying, leaving a distinctive monoenergetic photon signal in large underground detectors. The photon signal can exhibit sidereal-daily modulation, providing excellent separation from backgrounds. Using a detailed numerical simulation, we examine this process as a search strategy for magnetic inelastic dark matter with the dark matter mass near the weak scale, where the upscatter to the excited state and decay proceed through the same magnetic dipole transition operator. At lower inelastic splittings, the scattering is dominated by moderate mass elements in the Earth with high spin, especially $^{27}$Al, while at larger splittings, $^{56}$Fe becomes the dominant target. We show that the proposed large volume gaseous detector CYGNUS will have excellent sensitivity to this signal. Xenon detectors also provide excellent sensitivity through the inelastic nuclear recoil signal, and if a future signal is seen, we show that the synergy among both types of detection can provide strong evidence for magnetic inelastic dark matter. In the course we have calculated nuclear response functions for elements relevant for scattering in the Earth, which are publicly available on GitHub.

One of the most intriguing ideas for detecting inelastic dark matter is through the upscatter of the excited state and its subsequent decay into a photon.This was first attempted in 'luminous dark matter' [47] to explain the DAMA annual modulation, and subsequently as dark matter detection through 'two easy steps' [12].The parameter space of these proposals, however, has been ruled out by direct-detection bounds on the nonobservation of the nuclear recoil from upscattering of dark matter into its excited state [15].In 2019, we (with Harnik) [48] explored using the upscatter and subsequent decay to set bounds on large inelastic splittings, of order hundreds of keV, where the nuclear recoil sensitivity is severely limited by kinematics: only dark matter with the very highest velocities were capable of scattering off the moderately heavy elements (xenon and iodine) used in the large detect-detection experiments.We found upscatter off heavy elements in the Earth (dominantly lead) provided a sufficient flux of excited states such that Borexino was capable of setting bounds by using the sidereal-daily modulation intrinsic to this signal [48].The focus of the previous paper was dark matter inelastic scattering off nuclei through a contact interaction, such as Z-exchange, and subsequent decay to a photon, which arises for example in a narrowly-split model of Higgsinos [16,49].Later work also explored the possibility that neutrinos, similarly, could upscatter and subsequently decay into photons [50,51].
As discussed in [48], large gaseous detectors, proposed for directional detection of WIMP dark matter, could provide sensitivity also in a lower range of inelasticity where direct detection experiments have increasingly better sensitivity to the rate of nuclear recoil from upscattering.The inelasticity range we found in [48] was from approximately 40-150 keV, where it was possible that the photon signal in the gaseous detectors could surpass the-then bounds from xenon experiments.However, the model employed to analyze the sensitivity, namely upscatter through a contact interaction separate from decay through a magnetic dipole transition operator, becomes somewhat baroque in this inelasticity range.In particular, the magnetic dipole transition can easily dominate over the contact interaction in the upscatter rate.If this occurs, the sensitivity will be modified, since the decay length is no longer independent of the upscatter rate in the Earth.Thus, in this paper we focus on dark matter that can upscatter off nuclei into an excited state exclusively through a magnetic dipole transition operator, that then decays via the same operator back into dark matter and a photon.In the theory where the magnetic dipole transition occurs as a result of a Dirac fermion is split into two Majorana states, this is called Magnetic Inelastic Dark Matter (MIDM) [52,53].However, our results also apply to theories where the dark matter and the excited state are a Dirac fermions with "Dirac-to-Dirac" magnetic dipole transition interaction.These classes of models provide a concrete scenario in which as the strength of the magnetic dipole transition is varied, and therefore the resulting inelastic scattering rate off nuclei and the decay rate of the excited state are properly related.The physics that we explore contains similarities to our previous work [48]: dark matter χ upscatters to χ 2 off an element N in the Earth χ + N → χ 2 + N .The excited state continues in the same general direction as the dark matter itself, for tens to thousands of kilometers, before decaying inside a large underground detector.This type of signal can exhibit a strong sidereal-daily modulation because when the detector is on the Cygnus-facing side of the Earth, there is far less "target material" available to upscatter.
The signal we will consider therefore consists of single photons with an energy equal to the mass splitting δ of inelastic dark matter, and with a rate that can exhibit a strong modulation with a period of a sidereal day.The shape and phase of the modulation is predicted and depends on the geographical location of the detector, allowing for excellent signal-to-background discrimination.Our focus in this paper will be on using a large proposed gaseous detector, CYGNUS [54], which has been proposed to reach 1000 m 3 in volume.The ostensible objective of a directional detector is to observe scattering off the nuclei that make up the gas of the detector.Our signal, however, is the spontaneous decay of an excited state into dark matter and a photon, and so depends only on the volume of the detector, not the mass of the gas in the detector.All we require is that there is a decent efficiency to absorb the photon with few to hundreds of keV energy in order to generate a scintillation signal.Indeed, the ideal mode of operation of a large gaseous detector would be to inject the lowest pressure gas that has decent efficiency to scintillate, and choosing a gas that has exceptionally-low radioactive backgrounds (and/or can be made exceptionally pure), such as helium.In our analysis, we will consider the signal both in the presence of radioactive backgrounds at a level that the CYGNUS collaboration reports can be straightforwardly achieved [54], as well as a hypothetical mode of operation where the detector has been fully optimized for the inelastic photon signal, i.e., background-free.These will provide the realistic range of sensitivity depending on the efforts to reduce backgrounds and the time in a low-pressure, radiopure gas mode of operation.
We summarize the notation used in the paper in Table 1.

Dark Matter with a Magnetic Dipole Transition
Dark matter may be electrically neutral and yet possess higher-dimensional interactions with the photon.For fermionic dark matter, which is the focus of this paper, the leading interaction at dimension-5 is a magnetic dipole moment operator, that we will show, as a byproduct of our analysis, is strongly constrained from bounds on direct detection.Here Ψ is necessarily a Dirac fermion, Σ µν is a commutator of gamma matrices defined in Eq. (A.9), and the coefficient of the magnetic dipole moment μΨ has the dimensions of [mass] −1 .
An alternative possibility is that a magnetic dipole moment for dark matter is either absent or suppressed, and there is instead a magnetic dipole transition to an excited state in the dark sector.
A well-known model that has a vanishing magnetic dipole moment but allows for a magnetic dipole transition is Magnetic Inelastic Dark Matter that consists of a Dirac fermion split into two Majorana fermions [52,53], χ 2 and χ, with a small mass difference δ ≡ m χ 2 − m χ .We take δ > 0, and so χ is the dark matter while χ 2 is an excited state in the dark sector.As we detail in Appendix A, the magnetic dipole transition between two Majorana fermions arises from the operator where the magnetic dipole transition coefficient μχ has the units of [mass] −1 .
An alternative possibility is that there are two Dirac fermions with magnetic dipole transition operators between them.This arises in the Standard Model, for example, where the neutral, isospin-0 baryons Λ 0 s and Σ 0 have both magnetic dipole moments and magnetic dipole transitions among them.It is possible that dark sectors with dark baryons exist where the magnetic dipole moments vanish, and yet magnetic dipole transitions are present [55].The operators for the Dirac-to-Dirac magnetic dipole transitions have a nearly identical form to the analogous Majorana-to-Majorana transition, as shown in Appendix A. The two transitions, Ψ D1 → Ψ D2 and Ψ D1 → Ψ D2 in Eq. ( 3), are independent.If the magnetic dipole moment of dark matter vanishes (μ 11 = 0), we can map the theory space of the Majorana-to-Majorana magnetic dipole transition precisely onto the Diracto-Dirac magnetic dipole transition.A simple way to understand the correspondence is to first imagine that dark matter were entirely composed of Ψ D1 (and not ΨD1 ), which means the only magnetic dipole transition is Ψ D1 → ΨD2 .The decay of ΨD2 → Ψ D1 proceeds through the same interaction.As one shifts from a universe in which dark matter is made entirely of Ψ D1 into, for example, a symmetric abundance of Ψ D1 and ΨD1 , these two "species" of dark matter have two independent magnetic dipole transition operators that permit the two species to scatter into the excited states of ΨD2 and Ψ D2 , respectively.If the dark sector preserves CP, μ12 = μ21 , and so each "half" of the dark matter (particles or antiparticles) has the same magnetic dipole transition moment to an excited state as if there were only particles, or equivalently, if the particle transition were Majorana-to-Majorana. Hence, our use of Eq. ( 2) for our calculations in the paper applies to both a Majorana-to-Majorana model as well as a Dirac-to-Dirac model, so long as the magnetic dipole moment of dark matter vanishes.
The magnetic dipole transition operator could arise from a perturbative or nonperturbative ultraviolet (UV) completion.A perturbative generation of the operator through an electroweak loop, integrating out some charged states of mass m ± , is expected to lead to μχ ∼ eg 2 /(16π 2 m ± ), where g is the electroweak coupling.Instead, if the dark matter is a composite of a new stronglycoupled sector, the scale suppressing the operator is anticipated to be of order the mass of the dark matter itself, and so there is no loop suppression; for example, this is the case for the magnetic dipole transition of Σ 0 to Λ 0 s .
Given that this interaction involves integrating out heavy electrically charged particles, there is an upper bound on size of the interaction by simply requiring that no new charged particles are in conflict with the present LEP and LHC bounds.The robust bound on charged particles comes from LEP II [56,57], where m ± ≳ 100 GeV.This bound is essentially independent of lifetime.There are stronger bounds from the LHC in specific model-dependent scenarios.We will generally consider m χ = 1 TeV for most of our analyses, and so given that the dark matter is the lightest particle in the dark sector, the charged states are necessarily heavier, m ± > m χ .We normalize the operator in Eq. ( 2) as and so in practice the ratio m χ /m ± is absorbed into our definition of g M .Given that we take m ± > m χ = 1 TeV, this is fully sufficient to allow us to consider g M being either perturbative, g M ∼ 1/(16π2 ), or nonperturbative, g M ∼ 1.Our sensitivity estimates, and bounds from existing experiments, will be presented as bounds on g M .
The magnetic dipole transition operator in Eq. ( 2) allows for several interesting processes including upscattering off nuclei and excited state decay.(Scattering off electrons will be negligible for the dark matter masses we are interested in.)Our focus will be on upscatter of dark matter off nuclei in the Earth, followed by excited state decay into a monoenergetic photon, discussed in the following sections.There is also the possibility that the excited state could downscatter off nuclei, although this process is suppressed compared with spontaneous decay. 1

Excited state decay
The heavier dark matter state can decay through emission of a photon, χ 2 → χγ.The width of the excited state is where this decay rate applies to both the Majorana-to-Majorana as well as the Dirac-to-Dirac magnetic dipole transition-mediated decay.The lifetime of the excited state is τ 2 , and so when it moves with speed v in the Earth frame, the typical distance ℓ χ 2 it will travel is The decay length is illustrated in Figure 1 for both perturbative (left) and nonperturbative (right) values of g M .Importantly, splittings δ below O(few) keV will be below the threshold of foreseeable experimental searches, and above O(500) keV is also problematic due to the initial upscatter becoming kinematically disallowed (see [48]).Distances of order km and the Earth's radius, R E , are illustrated for reference.When ℓ χ 2 ≳ R E , the number of viable targets will drop significantly  6), as a function of δ, for m χ = 100 GeV, 1 TeV, and 10 TeV, for weakly-coupled (left) and strongly-coupled (right) theories; the extent of the shaded regions denote variation of v in the range 100 − 550 km/sec.as most decays take place outside the volume of the Earth.On the other hand, when ℓ χ 2 ≲ 1 km, the overburden becomes nearly isotropic in the vicinity of the detector, suppressing the modulating nature of our signal.Therefore the region between the gray dashed lines, where km ≲ ℓ χ 2 ≲ R E , is optimal.

Up-scattering
The magnetic dipole transition operator leads to dark matter scattering off Standard Model (SM) particles.A detailed and thorough analysis of this process is given in [58].For direct detection, the scattering of dark matter off a nucleon N = (p, n) occurs through the insertion of two photon currents, here N ′ denotes the same nucleus as N but with a different momentum.The first term is just what one would expect for the electromagnetic coupling of a nucleon, where p µ ≡ p µ N + p µ N ′ , q µ ≡ p µ N ′ − p µ N , and F 1 (q 2 ) and F 2 (q 2 ) are the conventional electromagnetic form factors of the nucleon with F 1 (0) = Q N and F 2 (0) = κ N .The full magnetic dipole moment of the nucleon is g N μN s N , with μN ≡ e/(2m N ) (the nucleon magneton) and s N = 1/2, and so we have The second term arises from the magnetic dipole transition (MDT) Since the speed distribution for dark matter in our halo has a cutoff at an escape velocity v esc /c ≪ 1, we can take the non-relativistic expansion of the spinors u i and ūi in Eqs. ( 8), (11) and then match to the set of inelastic scattering operators presented in Ref. [14].The only operators generated are O 1,4,5,6 which are defined as Here, ⃗ v ⊥ ≡ ⃗ v + ⃗ q/(2µ N ), ⃗ S χ is the dark matter spin vector, and ⃗ S N is the nucleon spin vector.As emphasized in [59], these constitute a complete set of Galilean and Hermitian invariants. 2 The operators O 1 and O 4 are the typical spin-independent and spin-dependent dark-matter-nucleon couplings, respectively.The operator O 5 is the coupling of dark matter spin with ⃗ q × ⃗ v ⊥ , which can be thought of as an analogue of the spin-orbit coupling ⃗ S χ • ⃗ L N , treating the orbital angular momentum of the nucleon as ⃗ L N = ⃗ r × ⃗ p with ⃗ p = m N ⃗ v and ⃗ r = ⃗ q/m 2 N .The operator O 6 , analogous to a tensor force, will turn out not to contribute to the scattering in our case, as we will see below.The amplitude for dark matter scattering off nucleons through the magnetic dipole transition operator is with the c N i given by The operator coefficients are given in terms of the nucleon charge Q N = (1, 0) and magnetic dipole moments g N ≃ (5.59, −3.83) for the (proton, neutron) respectively.It will also be useful to write the c i in terms of isospin quantities: Going from the DM-nucleon coupling Eq. ( 16) to the DM-nucleus coupling requires a model for how nucleons behave inside nuclei.We follow [60] and assume that the DM-nuclear interaction can be described by the individual DM-nucleon interactions summed over all nucleons in the nucleus.The insertion of the nucleon operators into the nucleus means that the operators O i must be evaluated for bound nucleons and the DM as a plane wave.These operators can be multipoledecomposed and the final spin-averaged, non-relativistic matrix element squared can be factorized as where j N is the spin of the nucleus, being the minimum dark matter speed necessary to scatter off the nucleus with recoil energy E R = q 2 /(2m N ), and µ is the dark-matter-nucleus reduced mass.In general the sum over k depends upon eight DM response functions, R k , and nuclear response functions W k .The DM response functions depend upon the details of the particle physics through the c i and the nuclear response functions are determined by expectation values of the allowed six nuclear operators in the nucleus. 3he DM response functions for inelastic dark matter with general couplings were presented in [14].For the case of the MIDM with the c i as given in Eqs. ( 17)-( 20) only four DM response functions are non-zero. 4These correspond to k = M, Σ ′ , ∆, ∆Σ ′ and are Here the τ, τ ′ are isospin labels and c τ i are given in Eqs. ( 21), (22).
The determination of the corresponding nuclear response functions requires knowledge of the nuclear wave functions, since the W k are expectation values of operators in nuclear states, The technical details of this procedure can be found in Appendix B, which follows [58][59][60].The outcome is that the nuclear response functions are simple polynomials (times an exponential) of the exchanged momentum q.For instance, the only non-zero response function for 24 Mg is W 00 M , which is

DM response function, nuclear response function
Table 1: Notation used in the paper.
Once the DM and nuclear response functions have been determined, the differential scattering cross section is given by where cos θ cm is the cosine of the center-of-mass scattering angle.One can compute the total scattering rate of MIDM with N T scattering targets per unit mass as where v min and v max bound the kinematically-allowed speeds (see e.g.[15] for details).For simplicity we will assume the velocity distribution, in galactic coordinates, is given by a Maxwell-Boltzmann distribution f gal (v MB ) with escape speed of v esc = 550 km/s and velocity dispersion of 220 km/s.To transform from the galactic frame velocity ⃗ v MB to the Earth frame velocity ⃗ v, we follow the series of Galilean boosts and rotations outlined in [48].We take the local DM density to be ρ χ = 0.3 GeV/cm 3 .The cross section dσ Z,A MDT /dE R is proportional to Eq. ( 30), and is easily derived from the relation (see e.g.Appendix D of [58] for a full derivation), For a given element the rate is calculable using the product of response functions in Eq. (23).
As an example, we show in Figure 2 the differential scattering rates for two isotopes that will be of significant relevance for us: 27 Al and 131 Xe.We show the total differential scattering rate (that includes interference among contributions), as well as the individual contributions from the different coefficients (which are individually squared to get the sub-component rates) c 1 , c 4 , and c 5 , for two specific values of the inelasticity δ = 5, 50 keV.Because both elements have non-zero spin (5/2 and 3/2 respectively), c 4 is non-zero.One of the critical observations in Figure 2 is that the scattering rates of both of these elements receive a large contribution from c 5 .For 27 Al, c 5 provides the dominant contribution for most recoil energies, while c 4 (and the interference between c 4 and c 5 ) also make a significant contribution to the scattering rate.For 131 Xe, c 5 dominates for E R ≲ 50 keV, whereas the response is dominated by the spin-dependent contribution c 4 for E R ≳ 50 keV.Notice also that in both cases the spin-independent contribution, c 1 , is negligible.The recoil energy-dependent (and δ-dependent) nuclear responses are one of the characteristics of theories that scatter through a magnetic dipole transition.

Signal Rate Calculation
The calculation of the expected event rate from dark matter upscattering off a nucleus in the Earth and subsequently decaying into a photon in a large volume underground detector was presented in [48].Here we review the inputs to the calculation as they pertain to dark matter scattering through the magnetic dipole transition operator.
The event rate involves integrating over all possible scattering sites in the Earth and requires knowing the elemental composition of the Earth.Unlike our previous calculation [48], the nontrivial dependence of the magnetic dipole transition cross section on several nuclear response functions means that we need to know the specific isotopic abundances of materials in the Earth, rather than the more straightforward atomic abundances.We denote the number density for a particular element with atomic number Z and mass number A as n Z,A (r), as a function of the radius r.For the radial dependence of the abundances, we employ a simplified three-layer model of the Earth Table 2: For each of the elements used to calculate the rate, we show the spin/parity, isotopic fraction (by number density), magnetic dipole moment (in units of the nuclear magneton, e 2mp ), and the number densities per cubic kilometer in the core, mantle, and crust.Additional information for more elements is provided in the ancillary materials.This data was compiled from Refs [61,[63][64][65], which report uncertainties not greater than 10% for the isotopes relevant to this work.consisting of a core, mantle, and crust.Thus, our number densities are piecewise-constant functions of radius.In particular, we take the Earth to be a sphere of radius 6371 km with the core and mantle having outer radii of 3483 km and 6341 km, respectively.
We extract the isotopic abundances for each element from an online database hosted by Brookhaven National Laboratory [61] and remove any isotope whose lifetime is shorter than the age of the Earth.Where possible, we cross-checked the results with [62].We then convolve these isotopic fractions, denoted f , with the known terrestrial abundances for the elements in the core [63], mantle [64], and crust [65] to arrive at the isotopic number density in each of the three regions of the Earth.We have collected these and other important quantities for stable isotopes on Earth in an ancillary file available on GitHub [66].
As we have seen, the magnetic dipole transition scattering cross section is a recoil energydependent combination of several nuclear response functions.In our numerical calculations the full dependence is taken into account.However, as a qualitative guide for which elements in the Earth will provide the dominant contributions to the rate, we present in Figure 3 the number density n T for each element (aggregated over all isotopes) (left panel), and the spin-weighted quantity, S × f × n Z,A , for each isotope (right panel).These are approximate measures of the density or number of spin-independent and spin-dependent targets, respectively.From these plots we expect that the rate will be dominated by contributions from 24 Mg, 25 Mg, 27 Al, 29 Si, and 56 Fe.Note that 57 Fe is another element that gives a non-zero contribution to the rate, especially for larger inelastic splittings.However, we found a high level of sensitivity of the nuclear response functions to choices about the details of the nuclear physics model for this particular nucleus, and so we did not include this in our calculations, making our bounds conservative.We present, in Table 2, the parameters of the five elements that dominate the signal rate and are used in our calculations.
In the centre-of-mass frame the scattering of the DM is isotropic, see Eq. (30), and the kinematics is straightforward.The outgoing speed of the excited DM state after a nuclear scatter  is with µ the reduced mass of the DM and the nuclear target, and similarly for µ 2 with m χ replaced by m χ 2 .In the Earth's frame (equivalently the lab frame), the scattering angle is related to that in the center-of-mass frame by The direction that the DM must scatter to reach the detector is determined by the relative position of the scatter site and the detector.This one direction in the lab frame could arise from two possible scattering angles in the center-of-mass frame, both of which result in χ 2 arriving at the detector.However, these correspond to two different speeds for χ 2 in the Earth frame, following the upscatter process χ + N → χ 2 + N .These two different speeds are For heavy dark matter, m χ > m N , the scattering is preferentially forward and the lab scattering angle θ lab is kinematically limited, with the maximum5 possible scattering angle denoted θ lab max .The cosine of this maximum scattering angle is Focusing on the case of heavy dark matter, we use Eq. ( 36) to express the outgoing speeds in the lab frame in a convenient form In order for χ 2 to arrive at the detector, the direction from the scattering site (at position ⃗ r s ) to the detector (at position ⃗ r D ) must lie within a cone subtended by opening angle θ lab max , and the fraction of the cone that the detector covers is given by [R D /(|⃗ r s −⃗ r D |θ lab max )] 2 , where we assume the detector is spherical with radius R D .Note that the amount of available scattering material grows as |⃗ r s − ⃗ r D | 2 while the likelihood of scattering towards the detector scales as |⃗ r s − ⃗ r D | −2 .Thus, all scatter sites within a decay length of the detector are of (approximately) equal importance.However, at distances larger than the typical decay length ℓ χ 2 , more distant scatter sites become less important.
The two different lab-frame speeds v lab out,± lead to two different decay lengths and consequently have different probabilities to decay in the detector.The probability that an excited state with lifetime τ 2 moving at speed v lab out will travel a distance L = |⃗ r s − ⃗ r D | and decay inside the detector is We illustrate the lab-frame kinematics for one choice of the outgoing speed in Figure 4.
Furthermore, for the two possible center-of-mass scattering angles there are two different values of the exchanged momentum, q ± , resulting in two different scattering cross sections.When calculating the total rate we sum over both of these possibilities.Integrating over all scattering angles in the lab frame requires introducing a Jacobian to relate the center-of-mass frame, where the scattering is isotropic, to the lab frame, Combining all of the effects discussed above we arrive at the final expression for the signal rate for scattering off one particular isotope in the Earth, The total rate will be a sum over all of the relevant isotopes in the Earth.In what follows we limit ourselves to a subset of all elements and focus on those which will give the dominant contribution to the rate, which are those summarized in Table 2.We are always interested in the limit where the detector size is smaller than the decay length R D ≪ ℓ χ 2 .Hence the rate scales as the volume of the detector and is independent of its geometry.

Rates at a typical detector
The rate of scattering is sensitive to the path length (L) in Eq. ( 39) that the dark matter must travel to reach the detector.This means that the signal rate at underground detectors is (sidereal) time-dependent [48], with a maximum rate at the time that the Cygnus constellation is furthest below the horizon.This is because the DM flux appears in Earth frame to be coming from the direction of the constellation Cygnus, and the rate is maximized when there is the greatest number of scattering targets within a typical decay length in front of the underground detector.The significant time-dependence of the rate -sidereal-daily modulation -depends upon the latitude of the laboratory, the distribution of elements within the Earth (Table 2), and the dark matter parameters.In our analysis we consider two locations for detectors: Gran Sasso, Italy at 42.6 • North and the Stawell Underground Physics Laboratory (SUPL), Australia at 37.1 • South.The Gran Sasso underground lab is at a depth of approximately 1.4 km; it contained the (now decommissioned) Borexino experiment and is a possible future location for a detector belonging to the worldwide CYGNUS directional dark matter experiment [67].SUPL, whose depth is around 1 km, is the current home of the SABRE-South experiment [68,69] as well as another possible home for a CYGNUS detector.The Cygnus constellation is below the horizon for about 2.5 hours for Gran Sasso and 18.5 hours at SUPL.
Using as inputs the information about the distribution of isotopes in the Earth, Table 2, along with the scattering cross sections and lifetime of the excited state discussed in Section 2, we can  calculate the signal rate using Eq.(39).For both Gran Sasso and SUPL we assume the underground detector has a volume of 1000 m 3 , as proposed for a future CYGNUS detector [54].As an illustration of the expected event rate for a de-excitation photon to be produced in a detector at each location, we show the rate in Figure 5 as a function of time throughout a single day, which we choose to be January 1.The phase of the daily modulation in the figure is arbitrary, with the peak chosen to line up with dark matter "noon" (which is not the same as local terrestrial noon).The solar time corresponding to the peak rate varies throughout the solar year, due to the difference between a sidereal day (approximately 23 hours and 56 minutes long) and a solar day.
In Figure 5(a) we show the rate at Gran Sasso for DM of mass m χ = 1 TeV, inelastic splitting δ = 15 keV, and magnetic dipole transition coefficient g M = 2 × 10 −3 .For these DM parameters, the typical decay length of the excited state is considerably larger than the size of the Earth; see Figure 1.This means the event rate grows linearly with the path length of the dark matter through the Earth.We show the rate for various elements but it is dominated at this latitude by scattering off elements in the crust, in particular aluminum.In Figure 5(b) we show the signal expected at SUPL for the same DM model.Over much of day the rate is again dominated by scattering off aluminum but the rate is an order of magnitude larger than at Gran Sasso due to the increased path length that DM takes through the Earth, since the Cygnus constellation is in the northern hemisphere while SUPL is in the southern hemisphere.At the time of highest rate, when Cygnus is lowest in the sky, some of the DM that de-excites in the detector has passed through the core of the Earth and the rate from iron briefly dominates (since there is very little aluminum in the core of the Earth, see Table 2).Note that in comparing the rate shape at Gran Sasso to that at SUPL, the peak is broader at SUPL because Cygnus is below the horizon for a longer fraction of the day at SUPL in the southern hemisphere.In Figure 5(c) we show the rate for DM with the same mass and magnetic dipole transition coupling but with an inelastic splitting of δ = 100 keV.With these parameters the decay length is shorter but is still comparable to the radius of the Earth; see Figure 1.For such a large splitting, the only possible scattering (among the elements we consider) is off iron, and again the rate at SUPL is larger than at Gran Sasso due to the increased overburden.
In Figure 6 we present a series of plots that show how the total daily rate, integrated over a sidereal day, varies as both mass splitting δ and size of magnetic dipole transition coefficient g M are changed.There is a intricate interplay between the scattering cross sections, which scales as g 2 M [Eq.(30)], the width of the excited state, which scales g 2 M δ 3 [Eq.( 5)], and the location in the Earth of the elements where the upscatter transition χ + N → χ 2 + N occurs.At large mass splitting and large coupling, where the distance traveled by χ 2 is < ∼ O(1) km, there is little difference between the rates expected at SUPL and Gran Sasso.At large mass splitting (when the scattering is mostly forward, off iron, and DM needs a high speed to upscatter) and small coupling the rate at SUPL is greater than that at Gran Sasso due to the increased overburden.Furthermore, the rate at SUPL scales as g 2 M whereas at Gran Sasso it scales as ∼ g 2.4 M , since there is a smaller fraction of overburden that is optimal for χ 2 decaying within Gran Sasso.At lower δ where scattering can occur off more elements the signal scales as a higher power of g M .The rapid decrease in signal at low mass splitting is due to the χ 2 lifetime becoming longer than the diameter of Earth.

Sensitivity to MIDM
We have seen so far that the total signal rate of photons from MIDM decay can be quite large, with many events occurring per day and a substantial sidereal-daily modulation effect, over the range O(few) keV ≲ δ ≲ 150 keV and for g M ≪ 1.Now we move to estimate the sensitivity to this signal in future experiments, as well as existing constraints on these parameters.In Section 4.1, we estimate the sensitivity of a future CYGNUS (or similar) gaseous detector to this monoenergetic photon signal.In Section 4.2, we recast present direct detection constraints on the same parameter space.

Sensitivity of CYGNUS
We calculate the sensitivity of a search for the photon signal arising from MIDM at the position of Gran Sasso as well as SUPL, which are two proposed locations under consideration for the eventual CYGNUS detector [67].We assume a fiducial volume of 1000 m 3 and a total experimental integration time of t int = 6 years, as projected in [54].The results are given in Figure 7 and explained below.The total signal rate (i.e.averaging modulation out) is the integral of Eq. (39) over the total integration time and is shown by the red contours in Figure 7.In practice, we simply sum over the rate in Eq. ( 39) in 1-hour bins to obtain the total rate/day; see e.g. Figure 6.Multiplying by the expected number of days of planned CYGNUS observation (N days = 6 • 365.25 = 2191.5),we also obtain the total signal rate.To be conservative we computed the daily rates using the position of Earth on January 1, even though the DM flux is expected to be highest in June and lowest in December or January (see e.g.[71]).
The dominant background in CYGNUS originates from radioactivity of the detector and surrounding materials, with a rate estimated to be of order 10 4 electron events/keV/year and approximately flat with energy [54].The energy resolution is of order 10% at few keV and is anticipated to be CYGNUS detector located at Gran Sasso (left) and SUPL (right), using t int = 6 year with the background rate projection of [54] (dark-blue shaded region), and t int = year with zero background (light-blue shaded region).The light gray (dark gray) region is the excluded by LZ [70] by requiring the number of signal events N sig ≲ 16.6 (356).The red contours give the total rate in the detector volume, and the dashed gray lines correspond to decay lengths roughly of order the overburden at each location; see text for details.
following [54] (see also the recent review of directional DM detection [72]). 6A significant rate of single photons originating from the decay of χ 2 might be interpreted as an excess in electronlike events, and therefore one might expect sensitivity to this signal to require the total signal rate exceed the expected background rate.However, the identifiable pattern of sidereal-daily modulation should allow for significantly-improved discrimination from background, as described below.
Following our previous work [48], we perform a simplified modulation analysis by dividing a sidereal day into two equal-size bins, labeled "signal on" and "signal off", corresponding to the halfday with highest and lowest rates, respectively.The reach of CYGNUS can then be estimated by requiring that the rate during "signal-on", Γ on , does not exceed the expected statistical fluctuation of events during "signal-off", whose rate is Γ off ; for a confidence interval of nσ we require where N off ≡ Γ off t off and t off = (N days /2) days is the total time during "signal-off".Note that Γ off contains a contribution from both the DM signal, during the "signal-off" period, and any additional intrinsic background in the experiment.The resulting 1σ sensitivity estimation is given in Figure 7 by the dark-blue region bounded by a solid blue line.Note that a dedicated search using a more robust modulation analysis (taking into account the predicted peak shape, for example) may significantly improve sensitivity relative to our estimates.
Note that the CYGNUS collaboration have proposed several ways to reduce their intrinsic backgrounds even further, with the eventual target of a zero-background experiment [54].We estimate the resulting sensitivity in this case using a total integration time of t int = 1 year; the result is given by the light-blue shaded regions in Figure 7.We observe that after only 1 year of zero-background CYGNUS running, one would already obtain substantial improvement in sensitivity relative to either 6 years of running with expected backgrounds, or existing direct detection experiments (see Section 4.2).This emphasizes the importance of further background rejection techniques in the development of CYGNUS-like experiments.
In Figure 7, we observe that the sensitivity diminishes at both small and large g M .At small g M , the total signal rate simply drops and the modulating signal becomes more difficult to observe.On the other hand, at large g M , the sensitivity is cut abruptly once the decay length of the excited state becomes shorter than (roughly) the typical overburden of rock that the DM sees along its path to the detector.Gran Sasso is located at latitude 43 • N and 1.4 km depth, and its typical overburden is of order ℓ ob ∼ 2 km, whereas SUPL is located at 37 • S at a depth of 1 km which implies a larger average overburden of ℓ ob ∼ 10 km; see [48] for details about the dependence of overburden on latitude.Therefore we expect the modulation to be significantly diminished once the most energetic DM particles have decay lengths ℓ χ 2 shorter than ℓ ob ; we illustrate this using the gray dashed contours, which match roughly the shape of the blue contours as expected.
Note that, as explained above, the sensitivity is also reduced at small δ once the decay length becomes comparable to the size of the Earth; solving Eq. ( 6) for ℓ χ 2 = R E gives g M ∼ 10 −3 (250 keV/δ) 3/2 v/450 km/sec for m χ = TeV.However, due to several mitigating factors below δ ∼ 100 keV, including the increased range of kinematically-allowed velocities and possible upscattering targets, the actual reduction in sensitivity in Figure 7 occurs at δ ≲ 30 − 50 keV rather than O(100)s keV.
Finally, in Figure 7 we also observe a distinct difference in the sensitivity contours between the two locations we considered, Gran Sasso and SUPL.There is a decrease in sensitivity for SUPL when compared to Gran Sasso at large δ, most pronounced for the zero-background sensitivity contour (light blue).This is due the shape of the modulation signal in SUPL, which is broader than at Gran Sasso; see Figure 5.This leads to sizable signal rate during the "signal-off" time in our simplified on-off analysis.This decrease is despite the relatively large total signal rate as illustrated by the red contours, and motivates an optimized modulation analysis.In addition, below the light gray region (discussed in the next section), there is a significant increased sensitivity to modulation at SUPL compared with Gran Sasso, broadly for δ < ∼ 100 keV and g M < ∼ 0.01.Computing the characteristic decay length of the excited state in this region, we find from Eq. ( 6) that ℓ χ 2 ≳ R E (see also Figure 1).Since SUPL's location in the southern hemisphere implies the dark matter wind is coming from directions that are below the horizon for most of the day, there is a significant increase in the upscatter rate arising from the much larger density of iron in the core of the Earth.Once the (δ, g M ) values give shorter decay lengths, the excited state upscattering in the core does not typically live long enough to reach the location of SUPL, causing a return to sensitivities that more closely resemble those at Gran Sasso, which is not able to benefit from the iron core upscatter enhancement in the rates.

Constraints from direct detection
Direct detection experiments searching for the nuclear recoil signal arising from the inelastic scattering of χ off nuclei are also highly sensitive; we focus here on LZ [70], which has the strongest limit to date on the parameter range considered in this work.
The rate for the nuclear recoil signal can be obtained by using Eq. ( 31) and integrating over the various isotopes of xenon with the kinematically-allowed range of velocities permitted by an inelastic collision [15].The cross section itself, dσ/dE R , is proportional to Eq. ( 30), and for a given xenon isotope is calculable using the product of response functions in Eq. ( 23).We have used the results of [60], that provided the nuclear structure functions for all relevant isotopes of xenon, with mass numbers A = 128, 129, 130, 131, 132, 134, and 136.Using the natural abundances of each isotope [61] as a proxy for xenon in the detector, we sum the contributions of each to the rate in Eq. ( 31), while fixing the total LZ fiducial target mass to 5.5 t [70].The total rate R is then computed as the integral of dR/dE R over E R in the sensitive range of the detector; for LZ, the detection efficiency is ≳ 50% for 5 keV ≲ E R ≲ 55 keV, which we take as the relevant range for our estimation. 7Finally, we estimate the expected number of events from MIDM in the detector during the LZ runtime of t int = 60 live days [70] as N sig ≃ R • t int .If N sig exceeds the statistical fluctuation in the number of observed events in this energy range, the corresponding parameters are therefore constrained.
The most recent dark matter search results obtained by the LZ experiment yielded ∼ 11 events in the nuclear recoil band [70].This number corresponds to the expected background under the assumption that the magnetic dipole transition inelastic nuclear scattering, χ + N → χ 2 + N , yields an S1 and S2 response consistent with the LZ experiment's calibrated nuclear recoil band.Using a 90% confidence interval around the observed 11 events gives N sig ≲ 16.6 events; the resulting constraint in the g M − δ plane is shown by the light-gray shaded region in Figure 7.We observe that, even in this case, the sensitivity of CYGNUS using current background estimations is still competitive with direct-detection constraints.Significant gains are possible for a southernhemisphere detector (right panel) and/or if the background estimates improve in the future.
It is interesting to consider that the LZ nuclear recoil band was established through neutron scattering off xenon.The magnetic dipole transition off nuclei proceeds, however, through photon exchange whose cross section is dominated by the O 5 nuclear response for the recoil energies of relevance to the bound.This interaction is significantly enhanced at low q 2 , and is quite unlike the contact interaction, O 1 , that neutron scattering captures in the calibration of LZ.The central value and spread of the nuclear recoil band has not been calibrated to a dark matter candidate that may be affected by the presence of the electrons.This suggests that it is at least possible the S1 and S2 response obtained from LZ's neutron calibration may not accurately represent the the nuclear recoil band for MIDM scattering.The most conservative interpretation of the LZ data for placing a limit on MIDM scattering would be to ignore the nuclear recoil band, and use all observed background events to set a limit; in the recent run, this total was 333 events (see Table 1 of [70]).The nuclear recoil calibration-independent bound can be set by using a 90% confidence interval around 333 events (i.e.356 events, assuming Gaussian distribution) as shown by the dark gray shaded region in Figure 7.This is definitively the upper bound for any possible scattering off xenon that leaves an electronic signal in the detector.Depending on the details of the S1 and S2 response to the specific nuclear responses involved for MIDM scattering, the true limit will lie somewhere between the dark and light shaded gray regions.We suspect that, given the typical momentum exchange even at the lowest recoil energies is |q| ≃ √ 2m Xe E R ≃ 30 MeV (for E R ≃ 3 keV), which is large compared to the inner electron binding energies, the S1 and S2 response is likely closer to the neutron scattering band.However, we believe it is important to not absolutely dismiss a potential difference, and so we show both the dark gray and light gray regions in Figure 7 to remind the reader of this subtlety.

Discussion
We have demonstrated that magnetic inelastic dark matter can be well-tested by a large underground gaseous detector using sidereal-daily modulation as the operative method to find the excited state decay to a photon line signal and distinguish it from background.There is clearly an advantage if the detector is located in the southern hemisphere, e.g., at SUPL, since this gives the largest rates and the maximum potential for discovery as we show in Figure 7.While all elements in the Earth contribute the total rate, we find the specific elements 27 Al and 56 Fe dominate the contribution to the upscatter process for the inelasticities that can be probed: few keV ≲ δ ≲ 150 keV.Smaller inelasticities δ ≲ few keV suffer an increasingly larger fraction of excited states decaying on distances much larger than the Earth, and thus suppressing the total rate.Larger inelasticities δ ≳ 150 keV suffer from a much smaller fraction of elements in the Earth that permit upscattering, as well as excited states decaying on distances that can be small compared with the detector overburden (∼ 1 km), diminishing the modulation of the signal.The "sweet spot" for the photon line signal that arises from magnetic inelastic dark matter fortuitously occurs in an energy regime that large gaseous detectors are capable of detecting.
That said, the existing xenon experimental results already probe a substantial portion of the parameter space.Exactly how much parameter space depends in part on whether the S1 and S2 response that the xenon experiments have calibrated using neutron scattering overlaps well with the nuclear recoil response arising from the magnetic inelastic transition.The reason for a possible difference has to do with the dominance of the O 5 nuclear response that is strongly enhanced at low q 2 due to photon exchange, which is quite unlike the calibration performed with neutron scattering that is dominated by very short range (pion exchange and contact) interactions.Given that the typical momentum exchange is larger than 30 MeV even at the lowest recoil energies that the experiments measure, we suspect that the xenon detector response is likely similar to neutron scattering, but we emphasize that this is an important issue to revisit.
The large gaseous detector that provided our probe of the photon signal for magnetic inelastic dark matter, CYGNUS, remains at the proposal stage.Xenon experiments, on the other hand, are well underway probing ever-smaller cross sections that include magnetic inelastic dark matter through upscattering off nuclei.Our view of the interplay among these experiments is that if xenon experiments were to find evidence for a nuclear recoil signal above background, the entire community would be immediately interested in understanding the nature of this dark matter signal and its interactions with the SM, as well as the astrophysical properties of the local dark matter halo, such as the local density and velocity distribution.Our study demonstrates that large gaseous detectors are complementary in their ability to probe the specific theory that was the focus of this paper: magnetic inelastic dark matter.Given that the inelasticity could be small, and thus may be difficult to discern in the nuclear recoil energy spectrum itself, our technique of using upscattering off the elements in the Earth and subsequent decay to a photon would provide outstanding information to distinguish among the various possible dark matter theory explanations.
While our main results in Figure 7 used the dark matter mass m χ = 1 TeV, here let us comment on what happens as one deviates from this value.In the case where m ± ∼ m χ , we can lower m χ to a few hundred GeV before running into LHC bounds (and at m ± = m χ = 100 GeV, the robust LEP II bounds) on charged particles.In this regime, lowering the dark matter mass will result in the standard increase in the rates from the increase in the dark matter number density holding the local dark matter energy density fixed, i.e., n χ ∝ ρ χ /m χ .The second effect is related to the strength of the upscattering and subsequent decay.These rates, however, depend on only on the magnetic dipole transition coefficient, μχ (when m χ ≫ m N , so that the reduced mass that enters the upscatter cross section is well-approximated by just the nucleus mass), that has the mass dependence μχ ∝ 1/m χ .Combining these effects gives the approximate scaling of the rates, σ χ n χ v ∝ g 2 M /m 3 χ .One can therefore read off the bounds and sensitivities to g M for other dark matter masses m χ ≳ O(100s) GeV by reinterpreting the y-axis to be i.e., probing smaller (larger) values of g M for smaller (larger) masses m χ .If instead we imagine m ± ≫ m χ , then the y-axis becomes where this expression continues to assume m χ ≫ m N .If instead we consider changing m χ to be near or below m N , there are several nonlinearities that prevent a simple scaling: there is a reduction in the rate due to the dependence of the reduced mass on m χ , and the sidereal daily modulation is diminished significantly since the scattering is no longer dominantly in the forward direction.
One of the auxiliary results from our analysis is that we also bound the magnetic dipole moment of dark matter (in the limit δ → 0); the limit can be readily obtained from Figure 7.For the purposes of this paragraph, let's write the magnetic dipole moment as in Eq. ( 1), i.e. μΨ = g M e/(4m χ ), using the same dimensionless coefficient g M .Then, using the bound we obtained from the LZ nuclear recoil band (i.e., assuming this applies for a magnetic dipole moment interaction), we find The direct detection bound scales proportional to σ χ n χ ∝ g 2 M /m 3 χ when m χ ≫ m N .For an electroweak dark matter candidate that acquires a magnetic dipole moment through a loop of charged particles (e.g., a one-loop diagram with W and a set of electrically charged excited states in the dark sector), we estimate g M ∼ g 2 mχ 16π 2 m ± , where g ≃ 0.65 is the electroweak coupling.In the case where the charged states in the dark sector are only slightly heavier than the dark matter itself, m ± ≃ m χ , the LZ bound corresponds to In the case of strongly-coupled dark matter, for two cases of g M we obtain It is striking that the xenon experiments are now probing dark matter masses over 100 TeV when dark matter has a magnetic dipole moment comparable to the neutron!
We also have calculated the nuclear responses for a much larger set of elements than originally presented in [59] by using the BIGSTICK nuclear physics code.Some of our results are new, while others we were able to cross-check with Refs.[59,73]; we find agreement in all of the overlapping results (see Figure 8).We have provided all of our numerical results for the response functions as auxiliary files in [66].
Finally, we remark that there is an interesting parameter space that may provide an opportunity for xenon experiments themselves to be sensitive to the photon line signal.Once the dark matter mass is small, m χ ≲ 10 GeV, the nuclear recoil constraints become rapidly weaker because nuclear scattering leaves an energy deposition in the detector that is nearly always below the energy threshold for detection.However, for 1 GeV ≲ m χ ≲ 10 GeV, it is possible for magnetic inelastic dark matter to still have a substantial rate to upscatter in the Earth, with the excited state decaying into a photon line signal.So long as few keV ≲ δ ≲ 10 keV, the photon line would show up as part of the electron background.While there is unfortunately only a very small modulating event fraction, the total event rate is, nevertheless, potentially large enough that this could be seen above backgrounds.Hence, our signal could provide a unique opportunity to access magnetic inelastic dark matter using existing xenon experiment data, that we plan to return to in future work [74].
The many-body matrix element of a one-body operator (as above) can be written as a product of the one-body density matrix times the matrix elements of the one-body operator: (in spin and isospin), and Ψ J;τ |α|,|β| are the One-Body Density Matrix (OBDM) elements for the ground-to-ground state transition.The composite nuclear states are represented by the collective quantum numbers α and β.We assume the underlying single-particle basis is given by states of the harmonic oscillator, in which case, Therefore, the problem of calculating the matrix elements ⟨α|O|β⟩ reduces to the determination of a finite set of polynomials p(y).The quantum numbers on which p(y) depends are: • N α : Initial and final principal quantum number, N α = 2(n α − 1) + ℓ α = (0, 1, 2, ...).
• J: The rank of the operator in the multipole expansion.
In [76], it is shown that these polynomials have the form where the J i , K i , I i are integers.We need only determine the coefficients in these polynomials once for each operator we use, as they do not depend on the element except in the selection of quantum numbers.We do so by using the Mathematica package described in Ref. [60].
The One Body Density Matrices, the Ψ in Eq. (B.3), can be determined for each element and each quantum number by using a model for the nuclear physics interactions.To do this we employ a modern automated code, BIGSTICK [77], which in turn uses input files from other codes, Oxbash [78] and NuShell [79].The calculation of the OBDM elements in BIGSTICK requires two files: a model file (typically with extension .sps)and an interaction file (typically with extension .int).The model files represent a model for a nuclear core with filled energy levels, plus valence protons and neutrons.In the model sp.sps, for example, the core is an 16 O (8 protons and 8 neutrons), corresponding the the magic number Z = N = 8 of filled states for both the proton and neutron.This corresponds to filling of the lowest-lying orbitals: 1s 1/2 (2 states), with orbital angular momentum ℓ = 0 and total j = 1/2; plus the lowest-lying ℓ = 1 orbitals (6 total states) labeled 0p 1/2 (j z = 1/2) and 0p 3/2 (j z = 3/2).The next available states are 0d 5/2 , 0d 3/2 , and 1s 1/2 (total of 10+2 = 12 states) for the valence protons and neutrons.Once these are all filled, we hit the next magic number Z = N = 20, fill a new core (now Ca) and begin again with the next set of orbitals.
The interaction files quantify the interactions among the valence nucleons in the background of the potential in the filled core (e.g. 16O), and so a given interaction model is only relevant in the context of a particular model of the core.We tabulate several model/interaction file combinations which are available for use in BIGSTICK, along with their corresponding valence states and range of applicability in N and Z, in Table 3.Note that BIGSTICK does not support different interaction models for protons and neutrons and above Z = 20 the number of protons and neutrons in stable isotopes become increasingly unequal.Thus, we expect the errors to grow as we move to heavier nuclei.We combine the polynomials for the matrix elements and the one body density matrices using the Mathematica package of Ref. [60] to determine the final nuclear response functions, W k .Note that doing so requires a little translation to convert the output of BIGSTICK to the input of the Mathematica package.
For some elements and response functions we can compare the output of our approach to results presented in Refs.[59,73]; see Figure 8.

B.1 Nuclear Response Functions
Below we present the form of the nuclear response functions for the four operators relevant for MIDM, for each of the elements studied in this paper.Note that higher-precision coefficients, including all the nuclear response functions (W M , W Σ ′′ , W Σ ′ , W Φ ′′ , W Φ′ , W ∆ , W Φ ′′ M , W ∆Σ ′ ), for these and additional elements are collected in an ancillary data file on GitHub [66].The file contains information for the elements 12 C, 14 N, 24 Mg, 25 Mg, 27 Al, 28 Si, 29 Si, 30 Si, 54 Fe, 56 Fe, 57 Fe, 58 Fe, 58   3: A summary of BIGSTICK models/interactions which can be used, along with their range of applicability in Z and N .We list all the available options that (a) treat p and n symmetrically, and (b) do not skip over any state (valence states are ordered by energy). 60Ni, 70 Ge, 72 Ge, 73 Ge, 74 Ge, 87 Sr, and 88 Sr.We expect the nuclear physics uncertainties to be sub-percent for light elements (Z < ∼ 20) and grow to O(10% − 50%) for elements as heavy as iron, and beyond.O from our calculation (solid lines) to previous works (dashed lines).For 27 Al (upper-left panel), we compare to Ref. [73], while for 29 Si (upper-right) and 56 Fe (lower-left), we compare to Ref. [59]; line colors correspond to the response functions for the five relevant MIDM operators for all isospin combinations (see legend).

Figure 1 :
Figure 1: Decay length of excited state, Eq. (6), as a function of δ, for m χ = 100 GeV, 1 TeV, and 10 TeV, for weakly-coupled (left) and strongly-coupled (right) theories; the extent of the shaded regions denote variation of v in the range 100 − 550 km/sec.

Figure 3 :
Figure 3: Number density (left) and spin-weighted isotopic fraction times number density for elements with non-zero densities on Earth (right) for each atomic number Z [63-65].

Figure 4 :
Figure 4: Illusration of the dark matter upscatter and lab frame kinematical quantities.The green region represents the Earth, while the gray region represents a 2-d projection of the cone subtended by opening angle θ lab max .In this figure, only one (of two) choices of v lab out,± was taken.

Figure 5 :
Figure 5: Daily modulation rates, using Eq.(39), at an underground detector placed in Gran Sasso (a) or SUPL (b) for dark matter of mass 1 TeV and splitting δ = 15 keV.In (c) the splitting is δ = 100 keV and the only scattering takes place off iron.Solar hour was normalized so that the peak of the scattering rates occurs in the middle of the plot; the actual solar hour corresponding to the peak scattering rate will vary from solar hour 0 to 24 over the course of a year.

Figure 6 :
Figure6: Total daily rate, using Eq.(39), for scattering off different elements at Gran Sasso (left) and SUPL (right) as a function of mass splitting δ and size of the magnetic dipole transition coefficient g M .

1 Figure 7 :
Figure7: Sensitivity estimate for the modulating MIDM signal in a 1000 m 3 CYGNUS detector located at Gran Sasso (left) and SUPL (right), using t int = 6 year with the background rate projection of[54] (dark-blue shaded region), and t int = year with zero background (light-blue shaded region).The light gray (dark gray) region is the excluded by LZ[70] by requiring the number of signal events N sig ≲ 16.6 (356).The red contours give the total rate in the detector volume, and the dashed gray lines correspond to decay lengths roughly of order the overburden at each location; see text for details.

Figure 8 :
Figure 8: Direct comparison of the response functions W τ τ ′O from our calculation (solid lines) to previous works (dashed lines).For 27 Al (upper-left panel), we compare to Ref.[73], while for29 Si (upper-right) and56 Fe (lower-left), we compare to Ref.[59]; line colors correspond to the response functions for the five relevant MIDM operators for all isospin combinations (see legend).