Dark Matter from Monogem

As a supernova shock expands into space, it may collide with dark matter particles, scattering them up to velocities more than an order of magnitude larger than typical dark matter velocities in the Milky Way. If a supernova remnant is close enough to Earth, and the appropriate age, this flux of high-velocity dark matter could be detectable in direct detection experiments, particularly if the dark matter interacts via a velocity-dependent operator. This could make it easier to detect light dark matter that would otherwise have too little energy to be detected. We show that the Monogem Ring supernova remnant is both close enough and the correct age to produce such a flux, and thus we produce novel direct detection constraints and sensitivities for future experiments.


I. INTRODUCTION
Identifying dark matter (DM) is one of the most important problems in particle physics and cosmology today, because although it makes up the vast majority of the matter in the universe, its particle properties remain a mystery [1][2][3].Decades of searching with direct detection experiments have produced strong constraints, but no unambiguous detection [4][5][6].This lack of a discovery has motivated searches for DM beyond the traditional Weakly Interacting Massive Particle (WIMP) paradigm [7][8][9][10][11].
We present a new mechanism for producing highvelocity DM: the upscattering of DM by nuclei in supernova shocks.Supernova shocks can expand at velocities in excess of 0.01c (e.g.[60]), an order of magnitude larger than typical DM velocities at Earth's position in the Milky Way.If DM is light enough, collisions can boost it up to roughly double this expansion velocity, producing a small flux of light but high-momentum particles.We focus here on the specific case of the Monogem Ring supernova remnant.As shown below, the ratio of the Monogem Ring's distance from Earth to its age is comparable to the upscattered DM velocity, meaning such a flux would be visible today.Note that this mechanism is entirely distinct from Fermi acceleration, which, as mentioned above, could accelerate millicharged DM in supernova shocks [55][56][57].We do not require DM to be charged, only to have nonzero elastic scattering cross section with nuclei.
Given that we consider DM interacting with large, but still nonrelativistic, velocities, this mechanism is ideal for probing DM with velocity-dependent interactions with nuclei.Specifically, nonrelativistic effective field theory (NREFT) [61,62] provides numerous effective operators beyond the traditional spin-independent/spin-dependent dichotomy, many of which produce cross sections that scale strongly with velocity.In this work, we produce new constraints on several of these operators, constraining new parameter space at low mass and large coupling.Furthermore, we show that detectors being deployed in the coming year will be sensitive to supernovaupscattered DM over a wide range of parameter space for these operators.
This paper is organized as follows.In Sec.II, we introduce nonrelativistic effective field theory, and discuss the velocity scaling of the operators presented.In Sec.III, we describe supernova expansion into the surrounding interstellar medium (ISM).In Sec.IV, we compute the flux of supernova-upscattered DM reaching Earth.In Sec.V, we compare DM-induced event rates to experimental data.In Sec.VI, we present our results.In Sec.VII, we give our conclusions.

II. NONRELATIVISTIC EFFECTIVE FIELD THEORY
Dark matter direct detection searches are typically performed in the context of spin-independent (SI) or spin-dependent (SD) cross sections.These are velocityand momentum-independent operators that arise naturally from nonrelativistic reductions of many generic UVcomplete models (see, e.g., Ref. [63]).However, there exist models in which momentum-dependent operators can dominate over the constant operators, motivating a more general treatment of possible interactions [64][65][66].The now-standard framework for the nonrelativistic effective field theory of DM interactions (NREFT, also called nonrelativistic effective operators) was laid out in Ref. [61], in which the various possible responses are combinations of the following four Galilean-invariant quantities: where ⃗ S χ and ⃗ S N are respectively the DM and nucleon spins, µ N is their reduced mass, ⃗ q is the exchanged momentum in a collision, and ⃗ v N is the DM-nucleon relative velocity.Ref. [61] tabulates 11 such operators, including the standard spin-independent response as O 1 and the spin-dependent response as O 4 .A few additional operators that do not appear in the above reference can also be constructed, as detailed in Ref. [62], resulting in a final list of 15 linearly independent operators appropriate for DM of spin 0 or 1 2 .In Table I, we list the form of each of these operators, as well as the velocity scaling of the total elastic scattering cross section each operator produces.In the third column, v is the DM-nucleus relative velocity.Additional operators can be constructed for higher-spin DM [67,68], but we restrict our work to the operators listed here.
The DM-nucleon velocity can be separated into a sum of the DM-nucleus velocity, and the relative velocity of each nucleon with respect to the nucleus.The internal dynamics can then be projected onto nuclear basis states.
As such, each of the operators listed here couples to a different combination of nuclear basis states.To turn the DM-nucleon interactions into DM-nucleus matrix element, this change of basis must be performed, and evaluated for each individual isotope using a nuclear shell model.This allows the cross section to be factored into a DM response function, dependent on the DM-nucleus relative velocity v, and a nuclear response function which suppresses interactions for finite momentum exchanges.Nuclear response functions were computed and tabulated in the form e −2y times a polynomial in y, with y ∝ q 2 , by [61] and [69] for a majority of elements that we will consider here.
From Table I, we see that almost all of the effective operators listed produce strongly velocity-dependent cross sections, with O 15 scaling as strongly as v 6 .It is valuable to probe velocity-dependent operators at velocities significantly higher than the galactic escape velocity, as the interaction strength may be orders of magnitude larger than it would be at halo velocities.In our analysis, we focus on O 15 , and on two operators, O 3 and O 6 , that scale as v 4 .

III. SUPERNOVA EJECTA A. Sedov-Taylor Solution
The early stages of the expansion of a supernova shock into a uniform medium are described by the Sedov-Taylor blastwave solution [70][71][72][73].In the free expansion or ejecta-dominated phase, the mass of the swept up ISM is small compared to the ejecta mass, and has little effect on the expansion of the shock.The Sedov-Taylor phase begins when enough ejecta has been swept up to significantly impact the expansion of the shock, usually around 1000 years after the explosion.This causes the expansion to slow down more rapidly.
The Sedov-Taylor solution allows us to compute the radius and expansion velocity of the shock as a function of time, for a given explosion energy, ejecta mass, and ambient ISM density.We adopt the parametrization given in Refs.[74,75], convenient for the discussion of supernova shocks: the radius is given by and the shock velocity by Here λ F E and λ ST are parameters describing the free expansion and Sedov-Taylor phases, respectively.For a Type Ia supernova expanding into a uniform medium, λ ST = 2/5, and λ F E = 4/7.
is a scale radius, and , with M ej denoting the ejecta mass, E SN the explosion energy, and n 0 the ambient ISM density.The factor of 1.27m p is the average ISM mass per nucleus.For a Type II supernova, expanding into the progenitor star's presupernova wind, λ ST = 2/3, and λ F E = 6/7.The scale radius is R 0 = Mej Vw Ṁ , where V w is the wind velocity and Ṁ is the mass loss rate.Following [74], we assume V w = 10 km/s and Ṁ = 10 −5 M ⊙ /year.In this case . The density of the presupernova wind around the supernova remnant, which we will use in the following Section, is

B. Ejecta Composition
We are also interested in the elemental composition of the ejecta, as different nuclei can in principle have very different scattering cross sections with DM.Ref. [76] performed simulations to model nucleosynthesis in stars ranging from 15 to 25 M ⊙ , and reported the abundance of different nuclear isotopes in the ejecta relative to solar abundances.In Table II, we report the mass fractions of the most abundant nuclei in supernova ejecta, averaged from the results of five simulations reported by Ref. [76] at different progenitor masses.The mass fractions are dominated by hydrogen and helium, as in the Sun, but the abundances of many heavier elements are enhanced by around a factor of 10 compared to solar abundances.

C. The Monogem Ring
The Monogem Ring is one of the closest supernova remnants to Earth, with an apparent diameter of 25°on the sky.Its distance from Earth is believed to be approximately D = 300 parsecs, based on supernova energetics, as well as the parallax distance to both the (arguably) associated pulsar PSR B0656+14 and a seemingly interacting star cluster in the Gemini H α ring [77][78][79].The most recent Sedov-Taylor calculation gives the Monogem Ring an age of 68,000 years, an explosion energy of 8.38 × 10 50 erg, and a surrounding ISM density of 3.73 × 10 −3 cm −3 [79].These fits do not, however, provide a value for the ejecta mass.Inserting these fit values into Eqns. 2 and 3 gives a present day radius and velocity that agree with the reported values to within a factor of 2, the values being nearly independent of the assumed ejecta mass.

IV. DARK MATTER UPSCATTERING IN SUPERNOVA SHOCKS A. Dark Matter Velocity Distribution
We assume that all of the ejecta produced by the supernova is located in a thin shell at R s (t) traveling at the shock velocity V s (t).While the density distribution within a supernova remnant is model dependent, it peaks sharply near the shock in the widely used Sedov model [73], and even more sharply in the more recent Chevalier model [80], with the result that the majority of the mass is indeed concentrated within ∼10% of the shock radius.The number of DM particles that this shell encounters per unit time is where we adopt the standard value ρ χ = 0.3 GeV cm −3 (the distance to Monogem is small compared to the scale radius of the Milky Way halo).In the limit of an optically thin shell, meaning that the probability of a dark matter particle colliding with a nucleus is ≪ 1, the probability that a dark matter particle is struck by a nucleus as the shell passes by it is given by where m i is the mass of a nucleus of species i, and the sum is over the nuclei in the ejecta.Strictly speaking, the above formula only accounts for dark matter being struck by the ejecta.If we add to the shell the mass of the presupernova wind that it sweeps up as it expands, treating this ambient material as composed purely of hydrogen, the probability becomes (7) where n(r) is the density of the presupernova wind obtained from Eq. 4, and i = 1 denotes hydrogen, so that the delta function picks out only the hydrogen contribution.
To turn this into a differential rate, we replace σ χi with dσχi dv , where we use dE = m χ vdv to write This allows us to write where the factor of m χ has canceled out, though the differential cross section still implicitly depends on the mass.To avoid confusion, we note that the factor of V s here represents the rate at which the ejecta is passing through the DM, while the factor of v arises from the transformation of the differential cross section above.The flux at Earth can then be obtained by integrating this flux over all time, ensuring that at a given time t, only particles upscattered to a velocity v = D/(Age − t) contribute: particles at any other velocity would have either already passed the Earth, or not reached us yet.Accounting for the geometric factor of 4πD 2 , this yields,

B. Numerical Implementation
We compute the DM flux at Earth (Eq.10) in a general NREFT framework using a modified version of the publicly available Capt'n General code [81,82].This code was originally developed to model DM capture in the Sun in effective theories.In this work, we have modified it to model DM being upscattered by stellar material to a velocity v, based on the equations in the preceding Sections, rather than being downscattered to below stellar escape velocity.Interactions could in principle be described by an arbitrary combination of the operators listed in Table I, for both isoscalar and isovector interactions, though we consider only individual operators and only the isoscalar case for this work.
For the case of spin-independent scattering (O 1 ), we validate our results by independently computing the flux at both Monogem and Earth, and find no significant difference between the two implementations.
Fig. 1 shows an example DM flux at Earth for m χ = 1 GeV.The different steps are due to scattering with different nuclei, which leads to a maximum DM velocity that depends on the nucleus mass.The large peak at low velocity comes from scattering with hydrogen, while the subsequent small peaks are due to helium, carbon, nitrogen, etc.The general trend of rising flux with higher velocity is due to the differential cross section's strong dependence on momentum transfer, while the downward slope after each step is due to the range of supernova ages (and thus ejecta velocities) that contribute to the upscattered flux.Note that the sharpness of the peaks is due to our assumption that all the ejecta travels at the same velocity, and is not an artifact of the plot's velocity resolution.When setting limits, the sharpness of these peaks is blunted by the detector energy resolution and the range of allowed recoil energies, so such a sharply peaked velocity distribution should not lead to unphysically narrow recoil spectra.

V. DATA ANALYSIS A. Direct Detection Rates
Given a DM flux, mass, and differential cross section, it is possible to compute the recoil spectrum in a detector.The DM flux Φ Earth (v), with units of (cm −2 s −1 )(cm/s) −1 , plays the role that the quantity ρχ mχ vf (v) would play in traditional direct detection formalism.In other words, the detection rate can be given by where N is the nucleus making up the detection material.
To compute the resulting recoil spectrum in the NREFT framework, we use a slightly modified version of the publicly available WIMpy NREFT code [83].This code computes recoil spectra for numerous nuclei and for an arbitrary combination of effective operators, and has been used by experimental collaborations such as DEAP [84].

B. Detectors Used
Because Φ Earth is proportional to the differential cross section, the recoil spectrum in Eq. 11 scales with two factors of the cross section.This means that increasing the cross section rapidly increases the detectability of the supernova-upscattered DM.For this reason, the best existing detectors to search for supernova upscattered DM are surface detectors, which have little shielding other than the atmosphere.In this section, we use data from the athermal phonon detector operated by the SuperCDMS Collaboration [17], because of its low energy threshold, location on the Earth's surface, and its relatively large exposure compared to other surface detectors such as the CRESST surface run [13].(This SuperCDMS detector does have some copper shielding, which will be discussed below.)However, as the surface detector target material is silicon, it has no sensitivity to purely spin-dependent operators such as O 6 .For this reason, we also derive limits from PICO-60 [85], which is particularly sensitive to spin-dependent interactions.Although PICO-60 is deep underground, the nuclei that make up most of Earth's crust do not interact via O 6 , so the shielding is primarily due to nitrogen in the atmosphere, and to subdominant components of Earth's crust like aluminum and sodium.
We furthermore find that a future underground detector with low threshold, low background, and large exposure would be able to probe a large region of DM parameter space.Specifically, we focus on the SuperCDMS germanium HV detectors, planned for deployment at SNO-LAB.The projected exposure for these detectors is 36 kg-yr [18], compared to the gram-day scale exposure of the existing surface run [17].For the Ge HV detectors, we assume a detection threshold of 40 eVnr as reported in Ref. [86], and take the projected background from this same reference, with a rate of 1-10 events/kg/yr/keVnr for most of the ROI.
Finally, we can also compute projected limits for a large, future xenon experiment like DARWIN.If no candidate events are seen in a 200 ton-year exposure [87] and assuming the same energy range and efficiency as XENON1T [88], such a detector could also set limits on supernova-boosted DM.

C. Dark Matter Propagation to the Detector
The couplings that we aim to probe are large enough that DM particles would scatter in the Earth's crust or atmosphere before reaching a detector, deflecting them away from their original trajectory and causing them to lose energy.A commonly used formalism to model this attenuation is to assume that DM travels in a straight line to the detector, modeling its energy loss as a continuous process, with the energy-loss rate where the summation is over nuclei in the crust or atmosphere.This formalism is not strictly correct when modeling light DM, as particles may scatter through large angles, making the straight-line assumption inaccurate.However, Ref. [89] showed that this straight-line approximation produces reasonably accurate or even conservative results when compared to more complete numerical simulations.As a result, the publicly available VERNE code [90,91] has been used for multiple low-mass direct detection analyses, where the scattering is essentially isotropic [16,17,24].In fact, because we consider cross sections that scale strongly with velocity, assuming that all particles lose the same amount of energy is conservative, as particles that lose less than the average amount of energy will interact with a larger cross section once reaching the detector.Ref. [92], and limits from CDMS (orange) are from Ref. [93].
We assume that any flux arriving from below the horizon is completely stopped by the Earth, and that Monogem is below the horizon for exactly half the day (it is actually below the horizon for very slightly less than half the day, so this is mildly conservative).We further assume that all particles arriving from above the horizon arrive at an angle 45°from vertical, as this is roughly Monogem's median zenith angle at the latitudes of the detectors we consider.This is again conservative, as the particles that traverse less shielding would arrive with much larger velocities and thus cross sections at the detector.
For deep-underground detectors like PICO-60, the Su-perCDMS SNOLAB experiment, and DARWIN, we neglect atmospheric shielding, considering only elements in the Earth's crust.We take the crust composition from Ref. [89].For the SuperCDMS surface run, we consider ∼14 meters water equivalent (mwe) of atmospheric shielding (i.e. the 10 mwe thickness of the atmosphere divided by cos 45°).The surface run also has 5 cm of copper around it, but unfortunately we have not been able to find differential cross sections for copper tabulated for the effective operators in question.Instead, we have computed cross sections for several different isotopes of iron, nickel, germanium, silicon, and aluminum, and opted to conservatively model copper as 56 Fe, as this nucleus had the largest cross section of any of the nuclei we considered.
For a more complete discussion of attenuation in the analogous case of cosmic ray-upscattered dark matter, see Refs.[28,31,96].

D. Limit Setting
Following Ref. [17], we set limits from the SuperCDMS surface run using the optimal interval method [97].This method assumes that all observed events could in principle be dark matter, so no background subtraction or fitting is performed.Roughly speaking, a dark matter mass-cross section pair is ruled out when the corresponding recoil spectrum is significantly larger than the observed event rate over a properly optimized energy interval.
Fig. 2 shows the surface run data compared to the predicted DM spectrum for a m χ = 1 GeV and c 2 15 m 4 v = 7 × 10 28 , a coupling that is just barely ruled out by the data.As with the velocity distribution, the recoil spectrum rises with energy because of the dependence on momentum transfer.Increasing the coupling increases the height of the recoil spectrum, but also pushes it to lower energy due to attenuation.If the coupling is made too large, the recoil spectrum will be either pushed entirely below the detector's energy threshold, or hidden by the high background at low energy.
Ref. [85] reports a total of 3 DM candidate events in the full PICO-60 data set.We determine the 90% CL by computing the range of mass and cross section that would produce an expectation value of 6.7 events in the full PICO-60 exposure.
The projection for DARWIN is determined in a similar way: assuming a 200 ton-year exposure [87], and the same efficiency and energy range as XENON1T [88], we compute the region of parameter space with an expectation value of at least 2.3 events, the 90% CL if zero candidate events are observed.
Finally, for the SuperCDMS SNOLAB projections, we perform a profile likelihood analysis of artificial data generated from projected backgrounds.For each DM mass and coupling, we generate 20 artificial data sets based on the projected background reported in Ref. [86].For each data set we perform a profile likelihood analysis, in the form of an unconstrained fit over the six projected background spectra plus the signal spectrum.We rule out a mass-coupling pair if at least 90% of the data sets rule it out at the 90% level.

VI. RESULTS
Fig. 3 shows our limits and projections for three effective operators: O 3 , O 6 , and O 15 .The total cross sections for O 3 and O 6 scale with v 4 , with O 6 being representative of spin-dependent interactions, while the cross section for O 15 scales with v 6 .
Our limits from the SuperCDMS surface run reach down to masses of about 70 MeV, about 20 MeV lower than limits based on the virialized Standard Halo Model (SHM) dark matter using the same detector.The limit from PICO-60 extends to just below 0.5 GeV, compared to the minimum mass of ∼3 GeV excluded by the SHM PICO-60 analysis.The large gap in our surface run exclusion region for O 3 , from about 0.6-1.2GeV is due to the shape of the observed data.The slope of the observed recoil spectrum is such that the signal spectrum for this mass range is always hidden by the data.Lighter particles lose energy very inefficiently during attenuation, while heavier particles carry more energy to begin with, meaning both lighter and heavier particles can be excluded for some range of couplings.
The projections show similar behavior.The projections for SuperCDMS extend down to 80-90 MeV, compared to a minimum mass of 400-500 MeV for the SHM analysis.The projections for DARWIN reach a minimum mass of 1.35-1.4GeV, down from ∼5 GeV for a SHM analysis.
Finally, cosmological and astrophysical limits have been set on effective operators or, more generally, velocity-dependent scattering.Ref. [95] set limits on effective operators based on DM scattering with hydrogen and helium leading to damping in the CMB at small scales.Ref. [94] set limits on general velocity-dependent scattering based on the Milky Way satellite galaxy population, which would be suppressed in the presence of such DM-nucleus interactions.Though these limits are not explicitly limits on effective field theory parameters, we recast them as limits on the corresponding operators.We show constraints from both these references in Fig. 3.
In principle, more detailed modeling of the ejecta's velocity structure could impact our limits.However, we expect this to be a small effect, especially at the low masses we are most interested in.As a cross check, we have computed how our limits would change if Monogem's expansion were treated as that of a Type Ia supernova rather than Type II, modifying the equations of Sec.III A according to Ref. [74,75].Qualitatively, the largest change is that the shock is expanding into the ambient ISM rather than a dense presupernova wind, resulting in higher velocity but less swept up hydrogen (in both cases we assume 5 M ⊙ of ejecta, admittedly too large for a Type Ia supernova, but our interest here was just on the impact of the velocity evolution).The limits we obtained in that case were slightly stronger for m χ ≳ 1 GeV, but virtually unchanged for lower masses.Because our results are not strongly affected by such a large change, we expect that the effect of a wider distribution of ejecta velocities would likewise be small.Further, a more realistic distribution of ejecta velocities would not necessarily lead to a less sharply peaked DM velocity distribution, because the shock, by its very nature, produces a sharp jump in both density and velocity.

VII. CONCLUSIONS
Supernovae are among the most energetic events in the Universe.In this work, we considered the Monogem Ring supernova remnant as a source of upscattered dark matter, as has been previously done for cosmic rays, the Sun, blazars, and supernova neutrinos.Exploiting the highbut still nonrelativistic-velocity of such dark matter, we derive limits and projected sensitivities to light dark matter that scatters via velocity-dependent effective operators.Our results extend to lower mass than existing nonrelativistic effective field theory searches, and cover a wide range of couplings.Our work provides a new avenue for future experiments to search for light dark matter that would otherwise be doomed to lurk below threshold.
helpful discussions, and especially to Bradley Kavanagh for updates made to, and help with, the WIMpy NREFT code.We also thank the anonymous referee for their helpful suggestions, which improved the paper.

FIG. 3 .
FIG.3.Limits on three non-relativistic operators due to the high-velocity dark matter component from the Monogem Supernova remnant.Results of this paper are shown in red and purple, where the lower edges are limited by the experiment sensitivity, and the upper edges by the overburden attenuation of DM flux on its way to the detector.Standard Halo Model (SHM) limits from CRESST-II (blue) are taken from Ref.[92], and limits from CDMS (orange) are from Ref.[93].The black cosmological limits are based on Milky Way satellites[94] (dashed) and CMB damping[95] (solid).
The authors are supported by the Arthur B. McDonald Canadian Astroparticle Physics Research Institute, NSERC, and the Province of Ontario through an Early Research Award.Computations were performed on equipment funded by the Canada Foundation for Innovation and the Ontario Government and operated by the Queen's Centre for Advanced Computing.Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science, and Economic Development, and by the Province of Ontario.

TABLE I .
List of effective dark matter-nucleon interaction operators, and the zero-momentum scaling of the corresponding cross section with the DM-nucleus relative velocity v.The operators we consider in this work are denoted in purple text.