Probing Earth-Bound Dark Matter with Nuclear Reactors

Strongly-interacting dark matter can be accumulated in large quantities inside the Earth, and for dark matter particles in a few GeV mass range, it can exist in large quantities near the Earth's surface. We investigate the constraints imposed on such dark matter properties by its upscattering by fast neutrons in nuclear reactors with subsequent scattering in nearby well-shielded dark matter detectors, schemes which are already used for searches of the coherent reactor neutrino scattering. We find that the existing experiments cover new parameter space on the spin-dependent interaction between dark matter and the nucleon. Similar experiments performed with research reactors, and lesser amount of shielding, may provide additional sensitivity to strongly-interacting dark matter.


Introduction
One of the first signs that the Standard Model of particles and fields (SM) needs to be extended comes from cosmology and astrophysics.Large classes of fairly diverse data point to the existence of the energy-density component in the Universe that can cluster and form cosmic structures.This component, a manifestation of some beyond-SM physics (BSM), is commonly referred to as dark matter (DM).Despite overwhelming evidence for the DM existence, its identity remains a mystery.A number of important classes of candidates were explored both experimentally and theoretically over the past decades.These include weakly interacting massive particles, or WIMPs, as well as super-weakly interacting bosonic fields such as axions.Despite tremendous efforts spent both theoretically and experimentally, no convincing signals of these particles have been found.In this paper, we will concentrate on one type of dark matter that interacts stronger than a WIMP, but nevertheless stays elusive to most of the terrestrial experiments.
If cross sections of DM scattering on atoms are very large, DM particles can quickly thermalize in Earth's atmosphere, rock and soil, making it energetically sub-threshold for most of the direct detection detectors, E χ,kin → kT Earth < E threshold .
(1.1) Indeed, the kinetic energy of kT Earth ∼ 0.03 eV is typically far below the sensitivity of modern devices.It is often the case that if strongly-interacting particles saturate the entire DM abundance, the limits from the balloon and surface direct detection experiments are too strong (see e.g.Ref. [1,2]).However, for a fractional DM abundance, characterized by f χ , f χ ≤ 1, large interaction cross sections are not excluded.(See, however, an opposing opinion [3] claiming that for some limited range of masses and cross sections f χ = 1 is still allowed.) There have been numerous ideas on how to probe f χ ≪ 1, l coll ≪ R ⊕ models of DM, where l coll stands for a typical DM-atom scattering length.If thermalization happens beyond 100 m -1 km length scale, then usual direct detection (DD) constraints apply up to O(1) factor.If however the thermalization is rapid, then one has to invent alternative probes.Several ideas have been exploited so far: 1. Effects of thermalized DM on collection of SM particles (e.g.particles in a beam, or very cold atoms such as in liquid helium) [4][5][6].
2. Existence of very fast component of DM that has more penetrating power (such as fast DM component arising from collisions with cosmic rays [7][8][9]).
4. De-excitation of nuclear isomers that is enhanced because of the relatively large momentum exchange, even with thermalized DM [12][13][14], and power measurement using novel quantum devices [15].
5. Possibility to use laboratory sources to accelerate DM and study its subsequent scattering.In particular, the use of underground accelerators was entertained in Refs.[16,17].
In this paper, we explore the idea very similar to the last one on this list.Namely, we analyze the existing constraints from coupling of nuclear reactors with well-shielded nearby DD devices.Experimentally, this is not a hypothetical scheme, but a mature effort directed to detect the coherent neutrino scattering on nuclei (CEνNS) at E ν < few MeV (see e.g.Refs.[18,19]).The essence of our idea is as follows: nuclear reactors, while operating, constantly generate the flux of MeV-energy neutrons.Prior to/during the collisions with a moderator (e.g.p in H 2 O), neutrons have a chance of interacting with DM, thereby accelerating it to velocities that are much larger than thermal.Subsequent scattering with an atom inside the detector can be visible.Thus, we suggest to use the nuclear reactors as the most powerful "O(MeV) neutron beam dumps".We find that the existing experiments at commercial nuclear reactors provide additional constraints on thermalized fraction of DM that are almost competitive with the existing DD constraints in the spin-independent (SI) interaction case, and probe a new parameter region not explored by the existing DD experiments in the spin-dependent (SD) interaction case, respectively.In the future, the sensitivity to thermalized DM can be further improved, providing an additional physics motivation for "coupling" DD devices to nuclear reactors.
We note that the research nuclear reactors, that are less powerful than the commercial ones, require less shielding around a neutrino/DM detector, and therefore could provide additional constraints in the large cross section/low scattering range regime.This paper is organized as follows.In the next section we review the accumulation and distribution of strongly-interacting DM in the Earth.We calculate the velocitized DM flux generated by the up-scattering on neutrons.In Section 3 we analyze possible signal, solving the problem of DM slow-down due to multiple collisions via numerical simulations.We reach our conclusions in Section 4.

Earth-bound dark matter flux accelerated by nuclear reactor
The Earth is exposed to a continuous flux of DM particles, and a component of the DM strongly interacting with the SM particles, if exits, can be efficiently captured by the Earth.As a result, the number density of the Earth-bound DM particles greatly exceeds the standard halo DM number density, and yet these are harder to detect by the standard DM direct detection experiments due to their small kinetic energy after thermalization.In this regard, nuclear reactors can play an interesting role as neutrons inside the reactors can upscatter the Earth-bound DM particles, enhancing their kinetic energy to be observed at nearby detectors.
To study the nuclear-reactor-generated flux of non-thermal DM, we need to combine the following ingredients: i. accumulation of dark matter inside the Earth, ii.upscattering of DM on fast neutrons inside the reactor, iii.subsequent propagation through shielding and eventual scattering inside a DM detector.In this section, we address these topics in turn.

Earth-bound dark matter
A strongly-interacting dark matter component (which we call χ), owing to their interactions with the Earth-nuclei, efficiently accumulates inside the Earth-volume.Their time evolution inside the Earth-volume (assuming no annihilation among them) is determined by the balance between their capture (Γ cap ) and thermal evaporation (Γ evap ) In the following, we briefly discuss both of these rates.
Starting with Γ cap , we first define the maximal capture rate as geometric capture rate (Γ geo ) which occurs when all of the transiting DM particles through the Earth get captured.Of course, the capture rate, Γ cap , is a certain fraction of the geometric capture rate which crucially depends on the DM mass as well as DM-nucleon scattering cross section.In the optically thick regime, i.e., for large DM-nucleon scattering cross section, this fraction, commonly known as capture fraction (f cap ), can be quite significant.For heavier DM (m χ ≫ m A with m A is the typical nuclear mass), it can even reach unity, whereas, for GeV-scale DM, f cap ∼ 0.1 [2].We use recent numerical simulations [2] to estimate the value of f cap which agrees reasonably well with the previous analytical estimate [4].To summarize, we use the following capture rate for the optically thick regime where f χ = ρ χ /ρ DM is the energy fraction of χ with ρ DM being the local DM energy density, m χ is the mass of χ, v gal is the average DM velocity inside the galaxy, and R ⊕ is the Earth radius.We take ρ DM = 0.4 GeVcm −3 , v gal = 220 km/s, and R ⊕ = 6371 km in our numerical computation.
Thermal evaporation, which is particularly important for light DM, occurs when the DM particles obtain too much kinetic energy via thermal kicks from the stellar nuclei and overcome the escape velocity of the stellar object.If the collision length of the DM particles is significantly smaller than the size of the Earth, which is of primary interest here, the evaporation effectively occurs from the "last scattering surface", and we estimate the evaporation rate by adopting the Jeans' formalism [4] where R LSS (v LSS ) denotes the radius (thermal velocity) at the last scattering surface.The thermal velocity at the last scattering surface (v LSS ) is related to the temperature at the last scattering surface as where σ χj denotes the scattering cross-section between DM and the j-th nuclei and n j (r) denotes the number density of the j-th nuclei.For the density profile of the Earth, we use Preliminary Reference Earth Model [20], and for the chemical compositional profile, we use Table I of Ref. [21]. 1 For the density, temperature, and compositional profile of the Earth's atmosphere, we use NRLMSISE-00 model [22].Note that, for our parameter range of interest, R LSS resides either at the Earth's crust or slightly above the crust, i.e., R LSS ≃ R ⊕ [4].To be more specific, the height of the atmosphere, relevant in our analysis, is at most ∼ 80 km, and is therefore negligible compared to the Earth's radius R ⊕ = 6371 km, justifying R LSS ≃ R ⊕ .
1 Most elements inside the Earth do not have a nuclear spin, and relatively rare isotopes such as 29 Si and 27 Al provide the dominant contribution to the last scattering surface inside the crust in the SD case.To evaluate Eq. (2.4), we have adopted the rescaling (2.20) given below (with mSi replaced by each element and the spin contents given in [2]).
Finally, G χ (r) = V ⊕ n χ (r)/N χ represents the spatial distribution of the DM particles within the Earth-volume.This is governed by the Boltzmann equation that combines the effects of gravity, concentration diffusion, and thermal diffusion [23,24] ∇n χ (r) In the above equation, we have used the hydrostatic equilibrium criterion as the diffusion timescales for the χ particles are short as compared to the other relevant timescales.
denotes the thermal diffusion co-efficient [24], g(r) is the acceleration due to gravity, and T (r) denotes the temperature profile of the Earth [25].
For a uniform DM distribution inside the Earth, G χ (r) = 1.We note that, in our parameter range of interest, the thermal diffusion co-efficient κ, which depends on the mass of the nuclei (m A ), remains almost constant for a given DM mass.Given such a weak dependence of κ with m A , we have taken for spin-independent (spin-dependent) interactions, while solving the Boltzmann equation to obtain G χ .Whereas, when we obtain the last scattering surface in Eq. (2.4), which is crucial to estimate the evaporation rate, we have included the more realistic chemical composition of the Earth.By solving Eq. (2.1) with the initial condition N χ = 0 at t = 0, we obtain where t ⊕ ≃ 1.4 × 10 17 sec is the age of the Earth.When evaporation is negligible, Eq. (2.6) simply reduces to N χ (t ⊕ ) ≈ t ⊕ Γ cap .DM number density at the Earth's surface, where nuclear reactors are, is therefore given by In the light mass region, m χ ≲ 1 GeV, this is suppressed by the efficient evaporation, while in the heavy mass region, m χ ≳ 3 GeV, the normalized surface density G χ starts to drop since most of the DM sinks to the center of the Earth. 2 As a result, the DM number density at the location of the reactors takes its maximal value at m χ ∼ 2 GeV.

Dark matter acceleration at nuclear reactor
Nuclear reactors utilize a fission chain reaction to extract the nuclear energy.Each fission reaction produces a few energetic neutrons, E n ≳ MeV, and the reactors slow them down, or moderate them, as the fission cross section is greatly enhanced at the lower neutron energy.To explore BSM physics, one can thus think of the nuclear reactor as a neutron beam dump experiment; energetic neutrons are injected into a moderator (typically, water or heavy water) and potentially produce and/or upscatter BSM particles before losing the kinetic energy.In particular, in the case of our interest, the reactor neutrons can upscatter DM particles accumulated inside the Earth, transferring them sufficient kinetic energy to be observable at detectors near the reactors.
235 U emits about 2.5 neutrons and 200 MeV per fission, released mostly as kinetic energy of the fission products.From this information, we can estimate the total number of neutrons inside the nuclear reactor per time as where P is the thermal power of the nuclear reactor, and we normalize it on a typical power of a large commercial reactor.One can obtain similar estimates for 238 U and 239 Pu.We have checked that the same estimate reproduces the antineutrino flux at the location of the CONUS detector [26,27].The fission neutron spectrum is well-fitted by the so-called Watt distribution, and the normalized spectrum is given (for 235 U) as [28] 1 We treat it as the neutron injection spectrum, and calculate the resulting flux of the DM, while taking into account the elastic scattering of neutrons by hydrogen in water.(We will not consider heavy water reactors, noting that one can easily extend the same treatment to D 2 O).The probability of a single neutron upscattering the DM while traveling the distance dx is where n χ is the DM number density at the reactor which we compute by following Sec.2.1, σ χn is the cross section between the neutron and DM, and E χ is the energy of the final state DM.In our study, we consider both SI and SD interactions, σ SI χn and σ SD χn , and the discussion here equally applies to both cases.The spectrum of the DM from a single neutron, before sufficiently moderated, is given by where E(x) is the neutron energy after traveling the distance of x with E(0) = E n being the initial neutron kinetic energy, −dE/dx is the energy loss rate of the neutron inside the reactor, and E th (x th ) is the threshold energy (travel distance) of the neutron we follow, which depends on the DM signal of our interest.Number density of dark matter n χ is to be taken as n χ (R ⊕ ) calculated in Sec.2.1.The energy loss function of neutrons, −dE/dx, is well studied in the context of, e.g., nuclear reactor engineering, but we may take a simplified treatment and estimate it as where n p is the number density of the target particle (or the proton), σ pn is the protonneutron cross section, and E loss = E/2 is the mean energy loss per scattering, assuming that the moderator is composed of water. 3The neutron-proton cross section is well-known, and we include it numerically in our computation.Combining these formulas, the DM number spectrum at the reactor is estimated as Notice that this formula closely resembles those appearing in the context of beam dump experiments.In the following, we assume that the DM-neutron scattering is isotropic, which results in where Θ is the Heaviside step function and we take the non-relativistic limit of the neutron and the DM.The accelerated DM spectrum at the source is now simplified as where we assume that the proton number density n p is constant inside the reactor.

Event rate at the detector
In the previous subsection, we calculated the accelerated DM spectrum at the source.
To calculate the event rate at the detector, we need to include DM slow-down in the shielding between the reactor and the detector.In a realistic situation, one would have to take into account reactor walls, concrete shielding around the reactor, as well as shielding surrounding the detector.Here we pursue a simplified approach, simply lumping all amount of shielding into a single parameter characterizing its overall thickness.
We start from the event rate of the DM reaching the detector without being attenuated by the shielding (zero collision case): where V det is the volume of the detector, n t is the number density of the target particles t inside the detector, t det is the exposure time, l det is the distance between the detector and the reactor, dσ χt /dE eve is the differential cross section between the DM and the target particles, and the subscript "0" stands for zero-scattering.We assume that the shielding has the thickness l s , the number density n s , and the scattering cross section σ χs .For the isotropic scattering between χ and t, we obtain where When the cross section becomes sizable, the probability of having zero-scattering is exponentially suppressed and multiple-scattering contributions become significant.Since the DM is significantly lighter than the particles inside the shielding in the case of our interest, m χ /m s ≲ 0.1, the DM particles retain sufficient kinetic energy even after experiencing multiple scatterings (see App.A for more details).To calculate the multiple-scattering contributions properly, we perform a Monte Carlo simulation.For simplicity, we assume that the shielding is spherical around the reactor core (whose finite size we ignore) with the radius l s , and that the shielding is composed solely of Si with 95 % of A = 28 and 5 % of A = 29, and its density 2.4 g/cm 3 (the typical value of concrete).For the SI case, we then have with m Si = 28m n and r Si = 4.03×10 −13 cm.For elastic cross sections, we use an ansatz that interpolates between a perturbative answer where nucleons of a given nucleus contribute to the scattering amplitude, and the strong coupling regime when the scattering occurs mostly at the nuclear surface.In the SD case, the shielding comes dominantly from 29 Si with the nuclear spin I = 1/2.We then take where we assume the isospin symmetric coupling and use the spin average value of the neutron and proton ⟨S n ⟩ + ⟨S p ⟩ = 0.172 quoted in [2,29].Throughout our computation, the DM is non-relativistic and the target particle is at rest in the laboratory frame.The scattering length is then given by l coll = 1/n s σ χs .In The DM energy spectrum emitted from the nuclear reactor before and after 5 m of the shielding in the SI case.The parameters are taken as P = 3.9 GW, m χ = 2 GeV, and σ SI χn = 10 −28 cm 2 , corresponding to l s /l coll = 8.07.
our Monte Carlo simulation, we take the step size as ∆l = min (0.1 × l coll , 0.01 × l s ).For each step, we generate a random variable uniformly distributing between 0 and 1.If it is smaller than 1 − e −∆l/l coll , we generate a scattering and change the direction of the DM particle isotropically.We stop each Monte Carlo step once the particle goes outside the sphere of radius l s , and in each process, we keep track of the number of scatterings that the DM experienced.Note that the numerical result depends only on the ratio l coll /l s .In the left panel of Fig. 1, we show the probability of a DM particle being scattered N times before exiting the shielding, P N , computed by our Monte Carlo simulation for different values of l s /l coll .It is clear from the plot that once l s /l coll exceeds unity, the probability of experiencing multiple times of scatterings increases drastically.Eventually we expect that the average number of scatterings before exiting the shielding scales as (l s /l coll ) 2 due to the random walk process.The DM spectrum after the shielding is given by where f N is the probability of the DM with the initial kinetic energy E ending up with the kinetic energy E χ after N -scatterings (see App. A).We may assume for simplicity that the flux is isotropic and all the particles are still radially outgoing after the shielding.The full event rate is then given by Here, note that, the subscript "eve" refers to quantities at the location of the detector.
To see the effects of multiple scatterings, we plot the DM spectrum emitted from the reactor before and after the shielding, with l s = 5 m, in the SI case in the right panel of Fig. 1.We take P = 3.9 GW, m χ = 2 GeV, and σ SI χn = 10 −28 cm 2 , which leads to l s /l coll = 8.07.Although most DM particles experience multiple scatterings and the DM spectrum is depleted, a small portion of DM particles still retains sufficient kinetic energy, primarily due to the scattering kinematics.Since the initial energy of DM particles is high enough compared to the threshold energy of the CEνNS experiments, the multiple scattering components greatly enhance the sensitivity of the neutrino experiments to our scenario, as we will see in the next section.

Constraints from existing nuclear reactors
Recently CEνNS has attracted a lot of attention, and extensive experimental programs world-wide aim at detecting the CEνNS signals both from neutrinos produced at fixed target experiments and from the reactor neutrinos.The latter type of experiment locates their detectors close to nuclear reactors, and therefore provides an ideal platform for probing the Earth-bound DM accelerated by the reactor neutrons.In the following, we focus on the CONUS experiment [26,27] as a representative of such an experiment.Other experiments that are/will be sensitive to the accelerated Earth-bound DM include Bullkid [30], CONNIE [31], Dresden-II [32,33], MINER [34], NEON [35], ν-cleus [36], νGeN [37], RED-100 [38], Ricochet [39], Texono [40], SBC [41], among others.
The CONUS experiment [26,27] in Brokdorf, Germany searches for the CEνNS events from reactor neutrinos.The commercial power plant in Brokdorf has the thermal power P = 3.9 GW.The CONUS detector is composed of Ge (for simplicity we assume 92 % of 74 Ge and 8 % of 73 Ge) and is located l det = 17.1 m away from the reactor core center, with shielding of l s = 5 m composed of concrete. 4We may use the analysis in [26] with the dataset of 248.7 kg×day, resulting in n t V det t det = 1.75 × 10 32 sec for the SI case and 1.4 × 10 31 sec for the SD case (originating from 73 Ge), respectively.For σ χt , the scattering cross section on target nuclei, we use with m Ge = 74m n and r Ge = 5.26 × 10 −13 cm for the SI case, and with I = 9/2 and ⟨S n ⟩ + ⟨S p ⟩ = 0.47 [2,29] for the SD case, respectively.Again, the formula in the SI case interpolates between perturbative and non-perturbative regimes.(A more accurate treatment can be implemented by explicitly introducing the nucleon-DM potential, and calculating the resulting scattering cross section.)In our case, the DM scatters off the nuclei, not electrons, inside the detector, and hence we need a conversion of the nuclear recoil energy to the electron-equivalent recoil energy (which is better calibrated).The quenching factor k, defined as the ratio of the ionization energy generated by nuclear recoils to the one produced by electron recoils with the same energy transfer, is rather uncertain in the energy range of our interest.We may take k = 0.18 [26] and assume that it is energy-independent for simplicity.The energy region of interest depends both on the detectors and different scientific runs in the analysis in [26], but we may take it as [0.3 keV ee , 1 keV ee ] to simplify our analysis, where the subscript "ee" stands for the electron-equivalent.Assuming that the collected charge output is linear in the electron recoil energy, we thus require [26] 5.6 keV 1.7 keV dE eve dN eve dE eve < 85, ( to derive the constraints on the Earth-bound DM from the CONUS experiment.Note that, even though the CONUS experiment focuses on O(1) keV recoil energy since their main target is the CEνNS events, our DM flux is extended up to O(1) MeV as shown in 1.
Therefore, an analysis dedicated to our DM flux with e.g. higher energy threshold may further improve the sensitivity.
Our results in the SI case are shown in Fig. 2 for f χ = 1 (the left panel) and f χ = 10 −3 (the right panel).The region enclosed by the red line is the full result including the multiple scatterings while the smaller region enclosed by the red line is the result with only the zero scattering component.The figures show that the multiple scattering contributions indeed enhance the sensitivity for the large values of the cross section.In comparison, we show the direct detection constraints from CRESST-III [42], CRESST (surface) [43], Darkside-50 [44], EDELWEISS (surface) [45], rocket based X-ray Quantum Calorimeter (XQC) [46], and balloon-based RRS experiment [47].We observe that the current parameter region of the accelerated Earth-bound DM explored by the CONUS experiment is already covered by the other direct detection experiments, both for f χ = 1, and for smaller values of f χ .
As it is evident from Fig. 2, the main obstacle for extending the constraints towards larger cross sections in the SI case is the amount of shielding employed in the searches at commercial reactors.The CONUS experiment loses the sensitivity for σ SI χn ≳ 10 −28 cm 2 due to the shielding.Research reactors, that typically have two-to-three orders of magnitude less power are not widely used so far for the detection of CEνNS, where the main limiting factor is the neutrino flux.For the thermalized DM searches, however, the research nuclear reactors may indeed provide additional sensitivity.Assuming that a DD experiment can operate at ∼ 3 m distance from the research reactor core (a factor of 10 closer than for Figure 2: Constraints on the spin-independent cross section σ SI χn from the CONUS experiment [26] for f χ = 1 (left) and f χ = 10 −3 (right).The bigger red shaded regions are the full result with the multiple-scattering components, whereas, the smaller red regions are the result with only the zero scattering component, respectively.Red dashed contours indicate possible reach when DD experiments are placed in proximity of the nuclear research reactors (see text for details).In comparison, the existing surface/underground direct detection constraints from CRESST-III [42], CRESST (surface) [43], Darkside-50 [44], EDELWEISS (surface) [45], rocket based X-ray Quantum Calorimeter (XQC) [46], and balloon-based RRS experiment [47] are shown in gray shaded regions.Cosmological constraints from CMB, Lyman-α measurements [48][49][50] as well as the constraints discussed in [4] do not cover any additional parameter space, and therefore are not shown for clarity.a typical placement of a detector at commercial reactor) with the thermal power of P = 10 MW, and employ only 1 m of shielding5 , we plot the projected sensitivity with the dashed line contour.One can see that the in f χ = 1 case, the sensitivity extends to the region constrained only by the rocket-based XQC experiment.At f χ = 10 −3 , new unconstrained parameter space can be covered.
The new sensitivity that we claim can be achieved is right in the range of masses and cross sections for a hypothetical exotic stable di-nucleon/sexaquark state [6].If stable, such particles may contribute to the DM abundance.The upscattering in the nuclear research reactors is perhaps one of the most promising ways to limit/study such model.
Next, our results in the SD case are shown in Fig. 3 for f  Figure 3: Constraints on the spin-dependent cross section σ SD χn from the CONUS experiment [26] for f χ = 1 (left) and f χ = 10 −3 (right), including the full multiple scattering contributions (the bigger red-shaded region) and the zero-scattering contribution only (the smaller red-shaded region).In comparison, the existing surface/underground direct detection constraints from CDMSlite [52], CRESST (surface) [43], rocket based X-ray Quantum Calorimeter (XQC) [46], and balloon-based RRS experiment [47] (adapted from [2]) are shown in gray shaded regions.Cosmological constraints from CMB, Lyman-α measurements [48][49][50] (not shown) assumes DM-proton scattering, and can be alleviated for neutron only scattering.which has a nuclear spin one, 6 while the dominant composition of the shielding in the CONUS experiment do not possess a finite nuclear spin.This indicates that the shielding to the DM particles in the SD case is effectively smaller, relative to the DD experiments at the Earth's surface, compared to the SI case.We note that our lower bound scales as f −1/2 χ since the event rate scales as f χ × σ 2 χn where one factor of σ χn comes from the up-scattering by fast neutrons inside the reactor, while the other factor of σ χn is from the scattering inside the detector.On the other hand, the bounds from the DD experiments scale as f −1 χ , therefore the parameter region newly explored by the CONUS experiment further enlarges for smaller values of f χ than shown in Fig. 3.
We note several possible avenues for future improvement.The big effort for detecting CEvNS at nuclear reactors is not optimized for searches of the Earth-bound DM.Therefore, given the extensive programs of the CEνNS experiments and with a further analysis of the events, we expect an improvement in the sensitivity in the near future.Depending on the thickness of shielding relative to scattering length, the recoil spectrum from DM could 6 Refs.[2,45,53] quoted ⟨Sn⟩ = ⟨Sp⟩ = 1/2 for 14 N.We note that a simple nuclear shell model instead predicts ⟨Sn⟩ = ⟨Sp⟩ = −1/6 which reproduces the measured magnetic moment well [54].This enhances the upper bound of the DD experiments, performed on the Earth's surface, shown in Fig. 3 by a factor of 9.In our computation of the DM surface density, we adopt the shell model value for the cross section with 14 N.  be considerably harder than the CEνNS recoil spectrum.This is because the momentum transfer to a nucleus in the detector, in case of CEνNS, is limited by neutrino energy q ∼ E ν , and for accelerated DM scattering, the same recoil momentum can reach O( m p E n ).Since reactors produce neutrinos and neutrons in the same energy range, the nuclear recoil energy in the detector for DM scattering compared to CEνNS can be enhanced by up to O(m n /E ν ) factor.
We also would like to point out that currently, there is some tension between the results of CEνNS searches with Dresden-II [32,33] and the rest of efforts, and especially CONUS.The results of [32,33] are interpreted by authors as CEνNS, albeit with the use of a somewhat large value for the quenching factor (i.e.assuming more efficient translation of nuclear recoil into eVee.)This experiment has some residual neutron background, that is absent in CONUS due to more extensive shielding in the latter case.We point out that while these claimed results are difficult to reconcile within the minimal CEνNS assumptions, the strongly interacting and accelerated DM scattering hypothesis may be consistent with both.Indeed, higher amount of shielding at CONUS may effectively extinguish the DM scattering, while Dresden-II may still be sensitive to it.This possibility may require a dedicated investigation by the experimental collaborations, with more realistic simulations of DM moderation than those performed in this paper.

Summary
We have investigated the use of DD experiments near nuclear reactors in the novel regime.We regard nuclear reactors as O(MeV) energy neutron beam dumps that have intensities of up to 10 20 neutrons per second.This is a much higher intensity that can be achieved with e.g.any proton beam of the same energy.Therefore, one can consider new physics scenarios, where nuclear reactors can provide additional sensitivity.
This paper is dedicated to an investigation of neutron-induced upscattering of dark matter particles with subsequent scattering in a nearby DM detector.Coupling of nuclear reactors and DD experiments is not a hypothetical scheme, but a series of concrete experiments aiming to detect the elastic scattering of antineutrinos.We show that the same experiments provide complementary sensitivity to strongly-interacting dark matter that may accumulate in large quantities inside the Earth in the spin-independent interaction case.Moreover, in the spin-dependent interaction case, we find that these experiments probe a new parameter region not covered by the other existing DD experiments.
We found that the thickness of shielding employed in a typical experiment at commercial nuclear reactors searching for CEνNS is a limiting factor that does not allow for an efficient transmission of neutron-upscattered DM from a reactor to a detector, especially in the spin-independent interaction case.We point out that the DD experiments at research reactors that require far less shielding may be used to extend sensitivity to so far unconstrained parts of the parameter space.
For general β, E N varies uniformly in between E i (1 − β) N to E i , and its probability distribution is given by [57] f where x = E N /E i and k = 1, 2, ..., N with k ≤ N .
In Fig. 4, we show the probability distribution for the resonant scattering case (β = 1) for N = 1, 2, and 6 collisions, and we compare it with a simplified approach where the random variable z is replaced by 1/2.As a concrete example, suppose, we take "x" number of DM particles with a monochromatic initial energy of 50 keV and count the DM particles whose energy falls above, say, 20 keV.In the z → 1/2 approach, it is evident that after first collision, all of the "x" number of DM particles are above threshold, whereas, after second collision, none of them are above threshold (right panel).So, either all of them (100%) or none of them (0%).However, in the realistic scenario, always a non-zero fraction of these DM particles are above threshold and this fraction falls off with larger number of collisions (left panel).For example, after first collision, only 60% of these particles will be above threshold, whereas, after second collision, only 23% of these particles will be above threshold, and so on.Therefore, it is evident that z → 1/2 approach is inaccurate as it either over-estimates or under-estimates the realistic fraction for any number of collisions.In the realistic approach, as shown in Fig. 1, a small but non-zero fraction of DM retains large kinetic energy even after O(10) scatterings, and this component extends the sensitivity to the larger cross section, as shown in Fig. 2.

− 1 ]Figure 1 :
Figure1: Left panel: The probability of DM particles being scattered N times before exiting the shielding, P N , for three different values of l s /l coll , computed with our Monte Carlo simulation.In this plot, we generated 10 7 events and the wiggles that show up around P N ∼ 10 −4 − 10 −5 are due to statistical fluctuations.Right panel: The DM energy spectrum emitted from the nuclear reactor before and after 5 m of the shielding in the SI case.The parameters are taken as P = 3.9 GW, m χ = 2 GeV, and σ SI χn = 10 −28 cm 2 , corresponding to l s /l coll = 8.07.
χ = 1 (the left panel) and f χ = 10 −3 (the right panel).The region enclosed by the red line is the constraint from the CONUS experiment, while the other existing DD experiments are shown in gray.In this case, we see that the CONUS experiment probes a new parameter region, not covered by the existing DD experiments.Note that the atmosphere is composed mainly of 14 N

- 1 ]Figure 4 :
Figure 4: Probability distribution function of final kinetic energy E N after N number of collisions for the resonant scattering (DM mass = target mass) scenario.(Left) probability distribution of z is used, (Right) z is replaced by its average value 1/2.