Earth Shielding and Daily Modulation from Electrophilic Boosted Dark Particles

Boosted dark particles of astrophysical origin can lead to nonstandard nuclear or electron recoil signals in direct detection experiments. %It has been shown that this interpretation successfully explains the excess of keV electron recoil events recently observed by the XENON1T experiment, and that a daily modulation of the signal in the detector is expected. We conduct an investigation of the daily modulation feature of a potential future signal of this type. In particular, we perform simulations of the dark particle interactions with electrons in atoms building up the Earth on its path to the detector, and provide in-depth predictions for the expected daily changes in the signal for various direct detection experiments, including XENONnT, PandaX, and LUX-ZEPLIN.


I. INTRODUCTION
Dark matter (DM) is certainly one of the greatest outstanding puzzles in modern particle physics.An enormous scientific effort has been undertaken, both on the theoretical and experimental sides, to shed more light on its nature, with great progress achieved in probing the parameter space of various particle physics models.Direct detection experiments offer a particularly promising way to search for DM, since in many models the DM particle is expected to undergo measurable recoils of nuclei and/or electrons in the detector.
A signal of this type has been hinted by the XENON1T experiment, which observed an excess of low-energy electron recoil events [1].One of the beyond Standard Model (SM) interpretations involved a boosted dark matter particle scattering on electrons [2,3].Although this effect was recently ruled out by the XENONnT experiment [4], and was most likely the result of beta decays of tritium, it still remains interesting to explore the possibility of detecting such signals in the future in XENONnT itself and other experiments.
To keep our analysis general, we will not require for the boosted particle to be the DM, thus we will refer to it simply as a boosted dark particle (BDP).If the BDP χ is much heavier than the electron, the observed electron energy deposition signal implies χ velocities of O(10 −1 ) c.Such fast-moving BDPs cannot come from the Milky Way halo and, instead, must be of astrophysical origin, produced, e.g., via semiannihilation χ + χ → χ + X (where X is a SM particle or a new particle eventually decaying to SM particles) [5], or via annihilation of a heavier dark sector particle, ψ, ψ + ψ → χ + χ [6].Either the Galactic Center (GC) or the Sun can be the dominant source of the BDP flux.In [7] specific particle physics models for χ and ψ are discussed, and ψ is shown to satisfy the requirements for a DM candidate, including an annihilation cross section leading to the observed DM relic density.
Searches for such BDPs have been proposed for large volume neutrino experiments, e.g., Super-Kamiokande [7][8][9][10], ProtoDUNE [11,12], IceCube [13,14], and DUNE [15][16][17][18] (see also [19][20][21][22] for related work).Here we focus on electrophilic BDPs, and our results are applicable to direct detection experiments like XENONnT, PandaX [23] and LUX-ZEPLIN [24].For electrophilic BDPs, when the BDPelectron scattering cross section is sizable, the electron ionization signal in direct detection experiments is expected to exhibit daily modulation due to the Earth shielding effect [2].This can be used to distinguish the BDP signal from various backgrounds.The information on the phase of the modulation can reveal the direction of the BDP flux, which would be of high importance for experimental analyses.
In this paper we extend the analysis of the daily modulation of the BDP signal of astrophysical origin, and explicitly account for the distribution of various elements inside the Earth, calculating their contributions to the BDP-electron ionization cross section.Our results apply to any direct detection experiment measuring electron recoil energies.The software used in this research is publicly available.The ionization form factor is calculated using AtomIonCalc 1 , which is refined from the software DarkARC [25,26].The software realEarthScatterDM 2 is used to simulate the BDP propagation inside the Earth, and was independently developed for this research.

II. BOOSTED DARK PARTICLE MODEL AND SOURCES
The particle physics model we consider here is a simple extension of the SM including four new particles: the DM ψ, the BDP χ, and the dark mediators V and Z .The dark mediator Z enables the annihilation ψ + ψ → χ + χ, the cross section for which is given by [7] where g ψ and g χ are the couplings of ψ and χ to Z , respectively, and Γ Z is the width of the Z .There is a wide range of parameter values for which this yields the correct DM relic density, i.e., when σ χ χv ≈ 3 × 10 −26 cm3 /s.However, there can also exist other DM annihilation channels involving SM particles, in which case σ χ χv can be much smaller.Some representative benchmark points are provided in [7].
The BDP interacts with electrons through the dark mediator V, as described by the Lagrangian terms If the mediator mass m V is much larger than the ∼ keV momentum transfer, in the parameter space of interest the cross section for BDPs scattering on free electrons simplifies to The benchmark points we consider here correspond to BDPelectron scattering cross sections of σ e = 10 −28 cm 2 , 10 −31 cm 2 , and 10 −33 cm 2 .In those scenarios, there is a wide range of values for the parameters g χ , g e and m V consistent with existing experimental bounds (see [2] for details).On the other hand, an electrophilic BDP can couple to a proton at the loop level through a mixing induced by charged leptons [27].However, we will not consider the resulting BDP-nucleus scattering channel for two reasons.First, the ionization cross section for the BDP, which is similar to an elastic scattering process on an electron at rest, typically dominates over the scattering cross section on a nucleus for m χ below 100 MeV.If we take oxygen as a benchmark, which dominates the inside of the Earth, the ratio is approximately 28 for m χ = 100 MeV.Furthermore, the coupling between a leptophilic BDP and a nucleon can be more suppressed if the massive mediator between the BDP and the lepton is a scalar or an axial vector.As discussed in [27], while maintaining the same parametric dependence in the BDP-electron scattering cross section, the BDP-proton interaction is only introduced at the 2-loop level for the scalar mediator, and it is absent up to 2-loops for the axial-vector mediator.Since including the BDP-nucleus scattering would not qualitatively change our results, we will neglect it.
The two main candidates for sources of a BDP flux are DM annihilation in the Galactic Center (GC) or halo, and DM capture with subsequent annihilation inside the Sun.Regarding the first possibility, the expected full-sky BDP flux from the GC can be estimated as [9] where σ χ χv denotes the thermally averaged annihilation cross section for ψ + ψ → χ + χ.For example, assuming thermally produced DM with mass m χ = 100 MeV, the expected BDP flux is Φ BDP GC 1.6 × 10 −2 cm −2 s −1 , saturated when ψ + ψ → χ + χ is the only annihilation channel, otherwise smaller.If the DM is produced nonthermally, then the annihilation cross section can be larger.
As for the second possibility, if the DM particles scatter off nuclei, they can be captured by the Sun and accumulate in its core.As discussed in [7], the Sun reaches capture-annihilation equilibrium for typical values of the DM scattering and annihilation cross sections, and the BDP flux becomes fully determined by the DM capture rate, taking the form [9] 3 where σ nucl ∼ 10 −42 cm 2 is on the order of the upper limit for the DM-nucleon scattering cross section set by spindependent DM direct detection experiments [29], and it was assumed that, at leading order, the scattering cross section is velocity-independent 4 .

III. MODEL OF THE EARTH
The Earth consists essentially of two parts: the core and the mantle.The eight most abundant atomic elements in the core and mantle [30,31] are shown in Table I.The remaining elements contribute a mass fraction below 1%.Due to the lack of precise information regarding the density of each element in terms of the distance from the Earth's center, we assume that a given element's mass fraction is constant in the core and mantle, and we take the value in each region to be the average value in Table I.The total density profile as a function of radius is taken from [32] and is shown in Fig. 1.The Earth is assumed to be isotropic, despite the complexity of its composition.
In the next section, we calculate the ionization form factor for all the elements in Table I [32], while the Earth's composition was adopted from [31].
the absolute abundance of elements at arbitrary radius r as demonstrated in Fig. 1, one can fully determine the scattering behavior of the BDP propagation in the Earth.

IV. DARK MATTER INDUCED IONIZATION
In this section we briefly summarize how BDPs ionize electrons bound inside atoms; for a more detailed discussion, see Appendices A, B, and C. The differential cross section for the ionization caused by an incoming BDP χ (with velocity v χ ) is given by where µ is the reduced mass of the BDP-electron system, a 0 = 1/(αm e ) is the Bohr radius, is the range for momentum transfer q, and F(q) is the BDP form factor, which for the model described by Eq. ( 2) is F(q) = 1.
The atomic form factor K(E R , q) for ionization describes the probability of obtaining a particular recoil energy of an ionized electron for a given momentum transfer q.We follow the calculation presented in [25,26].The wave functions of the electron initial states with quantum numbers (n, ) are taken to be the Roothan-Hartree-Fock (RHF) ground state wave functions whose radial part is described by a linear combination of Slater-type orbitals, The values of the parameters C j n , Z j , n j , as well as the binding energies for each element are provided in [33].The final state wave functions, which are asymptotically free spherical waves in a central potential, are given in [34].The atomic form factor K(E R , q) defined in [35,36] is related to the ionization response function where Θ is the Heaviside function.We have where E n B is the binding energy of the initial state electron, and k is the momentum of the final state ionized electron.We take into account contributions from all accessible states.A detailed calculation of K(E R , q) is presented in Appendices A and B.
For the energy regime of the BDP scenario considered here, the energy losses are dominated by the ionization process.Scattering with the valence and conducting electrons, due to their small binding energy, should recover the elastic scattering limit at E R ∼ O(1) keV.We thus treat the electron ionization in the tight binding approximation, where the electrons are assumed to have limited interactions with the neighboring atoms, so that the uncertainty of the molecular composition can be ignored.We also neglect the dissipation induced by the transitions of an electron among bound states (see [37,38]), since these are subdominant compared to the ionization when the typical recoil energy E R is much larger than the binding energy of valence electrons.Under those assumptions, we calculate the ionization form factor K(E R , q) for each of the elements listed in Table I and show the results in Fig. 2. For all the cases, as expected, the ionization form factors approach the kinetic region of elastic scattering, i.e., E R = q 2 /(2m e ), when E R is much larger than the binding energy.On the other hand, when E R is just enough to ionize an electron, q has a broader distribution.

V. PROPAGATION OF BOOSTED DARK PARTICLES INSIDE THE EARTH A. Overview of the Monte Carlo simulation
We assume that the BDPs are produced monochromatically and arrive at Earth from a fixed direction [2].Thus, the in-FIG.2: The atomic ionization form factor K(E R , q) for different atoms listed in Table I in the tight binding limit.The radial wave functions are determined using the RHF ground state wave functions in Eq. ( 7) with the coefficients C j n , Z j , n j and binding energies provided in [33].The q and E R distribution converges to E R = q 2 /(2m e ) in the large recoil energy limit, which is labeled by the white solid line.
coming BDP flux can be written as 5 5 We confirmed that the effects of the BDP interactions with the galactic medium and the Earth's atmosphere on its velocity distribution are negligible compared to the effect of the interactions inside the Earth as it travels to the detector.Given how the latter affects the BDP velocity distribution, where Φ 0 is the total initial flux directed towards the Earth.A schematic diagram of the model is shown in Fig. 3.
the former does not introduce a sizable modification to our assumption of a monochromatic energy spectrum for the BDP flux.The influence of the atmosphere to hadrophilic dark matter was discussed, for example, in [39].
FIG. 3: Schematic plot showing a BDP flux arriving from a particular direction with velocity v 0 χ .The polar angle θ is between the direction of the initial flux and the direction pointing from the Earth's center to the detector; the angle δ χ is the declination of the BDP flux direction in the equatorial coordinate system, ranging from −π/2 in the south to π/2 in the north.The blue and black solid lines denote the Earth's rotation axis and the incident direction of the flux, respectively.The light and dark orange regions correspond to the Earth's mantle and core, and the cyan cylinder denotes the detector.
In order to understand the propagation of the BDP inside the Earth, one first needs to consider the interaction between the BDP and the Earth's elements.According to Eq. ( 6), the probability distribution of the BDP final state after scattering is fully determined by the ionization form factor K(E R , q), where E R is the recoil energy and q is the momentum transfer.From Eq. ( 6), the mean free path of the BDP inside the Earth can be calculated as where the index a denotes the type of the Earth's element.n a (r) is the number density of element a at radius r, which can be calculated from Fig. 1. σ a ion (v χ , m χ ) is the ionization cross section between element a and a BDP with velocity v χ and mass m χ , obtained by integrating out the recoil energy E R and momentum transfer q in Eq. (6).
We have developed a Monte Carlo simulation to study the BDP propagation inside the Earth.The flow chart of the simulation is shown in Fig. 4, and a more detailed description can be found in Appendix D. We start with BDPs of mass m χ and velocity v 0 χ , evenly distributed on the plane perpendicular to v 0 χ .The main structure of the simulation is the iteration of scattering (i and f denote the initial and final state for each step, respectively).In each iteration, we first calculate the mean free path l ion fp (x i , v i , m χ ) using Eq.(10).Next, we use the exponential distribution e −l/l ion fp /l ion fp to sample l, which denotes the propagation distance for the BDP in this step of the iteration.The final position of the BDP, x f , can thus be easily calculated from x i , v i and l.Then, we sample the recoil energy E R and the momentum transfer q whose probability distribution is proportional to q × K(E R , q) according to the differential cross section in Eq. ( 6).Meanwhile the azimuthal angle β on the transverse plane with respect to the initial velocity is drawn from a flat distribution between 0 and 2π.The values of E R , q, and β fully determine the momentum transfer vector q, which is used to calculate the final velocity v f .Lastly, the pair ( x f , v f ) is used as the input for the next iteration as ( x i , v i ).
Additionally, in each iteration we check whether the trajectory crosses the mantle-core border.If it does, we recalculate the mean free path l ion fp and reset the starting point for this iteration at the spot where the crossing happens.Furthermore, the initial velocity v i remains unchanged.The location and velocity at each iteration are recorded.The simulation stops once the BDP exits the Earth or when its velocity is smaller than the threshold velocity, which is either the DM virial velocity or the minimum velocity to ionize an electron in xenon.For more details, please see Appendix D 2. Finally, we perform the simulation with different BDP's initial velocity directions to account for the effect of Earth's rotation, as demonstrated in Fig. 3.

B. Distortion of the velocity distribution
Due to propagation inside the Earth, the BDP velocity distribution is distorted when reaching the detector.The amount of distortion depends on the polar angle θ between the incoming BDP flux and the direction pointing from the Earth's center to the detector, as shown in Fig. 3. Before showing the results of the Monte Carlo simulation, we first present a qualitative estimate of the distortion of the BDP velocity distribution.
The distance traveled l inside the Earth depends on the depth of the detector d and the direction of the incoming BDP flux.In terms of θ, it can be written as where R D ≡ R E − d.In the limit d R E , l ranges from d to √ 2R E d on the near side ( π 2 < θ ≤ π) and from 2 ).The BDP kinetic energy, E kin ≡ m χ v 2 χ /2, is smeared due to dissipation from ionization.For each scattering, the typical energy loss in the elastic scattering limit is m e v 2 χ when the BDP is much heavier than electrons [2] (see Appendix C for a more detailed discussion).Thus the energy dissipation can be approximated in terms of the mean free path l ion fp as from which one can derive the dissipation of velocity as In the elastic scattering approximation, the mean free path can be written as l free fp (r) = n e (r)σ e −1 , where n e (r) = a n a (r)Z a is the electron density including the contributions of all elements inside the Earth and σ e is the scattering cross section between the BDP and a free electron.In Fig. 5, we compare the mean free path for ionization, l ion fp , with the one from elastic scattering, l free fp .At low BDP velocities, the finite binding energy suppresses the ionization.On the other hand, when v χ 10 −2 c, l free fp (r) serves as a good approximation for l ion fp (r).In this approximation, taking the electron number density as 1 × 10 24 /cm 3 near the Earth's surface, 1.3 × 10 24 /cm 3 at the mantle, and 3 × 10 24 /cm 3 at the core, the mean free path of the BDP in each region is l S fp ∼ 100 m × 10 −28 cm 2 /σ e , l M fp ∼ 75 m × 10 −28 cm 2 /σ e , and l C fp ∼ 33 m × 10 −28 cm 2 /σ e , respectively.
According to Eq. ( 13), one can define the effective distance at which the velocity distortion is significant, This can be used to classify the distortion of the velocity distribution into several cases: In Fig. 6, we show the results of our simulation for the BDP velocity distribution when it reaches the detector as a function of cos θ.Nine cases are presented, illustrating how the differential velocity distribution depends on various model parameters.The first row corresponds to different values of the BDP mass m χ .Equation (13) implies that the larger the mass, the less distorted the velocity distribution is after scattering inside the Earth.The second row corresponds to a variation of the initial velocity v 0 χ .The third row compares three cases with various BDP-electron scattering cross sections σ e , corresponding to scenarios with strong interaction, weak interaction, and extremely weak interaction, respectively.

VI. DAILY MODULATION OF IONIZATION SIGNALS
Due to Earth's rotation, the angle θ between the direction of the incoming BDP flux and the detector varies with a period of one day, ) where δ χ is the declination of the source of the BDP flux and δ D is the detector's declination projected onto the celestial sphere.The time at which the BDP flux reaches at upper culmination of the detector, denoted as t 0 , can be determined using 2πt 0 /(24h) ≡ (α χ − λ D ) in terms of Greenwich Mean Sidereal Time (GMST).Here, α χ represents the right ascension of the BDP flux, and λ D represents the longitude of the detector.As an example, assuming a BDP flux from the GC with δ χ,GC −29.00 • and α χ,GC 266.40 ), respectively.Since the GC is on the southern hemisphere and the three detectors we consider in this study are on the northern hemisphere, the detectors are on the far side of the Earth with respect to the BDP flux for the majority of the time.
Apart from BDP from the GC, one can also consider BDP from the Sun.In this case, the daily modulation is more conveniently described by the Coordinated Universal Time (UTC), shown in the bottom panel of Fig. 7.The value of δ χ,Sun varies from −23.5 • on December 21st to 23.5 • on June 20th.We take t 0 = 11.1 h according to the longitude of Gran Sasso.The signal rate for each experiment can be written as where the differential cross section is provided in Eq. ( 6).N d 4.2 × 10 27 ton −1 is the number of xenon atoms in the detector.The electron recoil energy spectrum varies with time.
The first column of Fig. 8 shows the velocity distribution for the three cases.As expected, the distortion of the flux varies when changing the ratio between Φ 0 and σ e .The three dashed lines correspond to the values of cos θ for the detector at Gran Sasso at t − t 0 = 0 h, 6 h, and 12 h.In the most dis-torted case, i.e., when σ e = 10 −28 cm 2 , the flux is completely shielded at t − t 0 = 12 h, while for σ e = 10 −33 cm 2 the distortion of the flux is negligible.The second column shows the recoil energy spectrum from Eq. ( 16) at the three different time t − t 0 .The distortion of the flux leads to a shift of the recoil energy events towards lower energy bins.In addition, a time-averaged spectrum is denoted by the orange solid line.Finally, the last column shows the time evolution of the event rate in the three bins: [0, 5 keV], [5 keV, 10 keV], [10 keV, 15 keV], and the sum of those three bins.
All those unique features in the experimental data can be used to extract the properties of the BDP in a systematic and comprehensive manner.
One way to examine the daily modulation signal is to perform a Fourier transform on the data [40].The signal can be parametrized as where T is the modulation period.For example, T is either one sidereal day or one solar day (when the BDP source is the GC or the Sun, respectively).A n is the amplitude in the Fourier series and t n is the relative phase.For the signal, t n converges to t 0 for a given recoil energy.A fit to t n provides information on the direction of the BDP flux.If we correlate the time series of the signals for three different detectors, one expects the differences in t 0 for each detector to be related to the differences in the detector locations.

VII. DISCUSSION
In this paper, we carried out a detailed analysis of the daily modulation of the signal expected from BDPs interacting with electrons.Such an effect can be searched for in terrestrial DM direct detection experiments such as XENONnT, PandaX and LUX-ZEPLIN.
In particular, we developed a Monte Carlo code to simulate the BDP's trajectory through the Earth to the detector.Considering a benchmark scenario with the BDP source located at the GC, we calculated the expected time variation of the signal in terms of the electron recoil event rate.Our predictions can be directly compared to current and future data.
It is worth noting that a different study regarding the daily modulation of the BDP signal has been carried out in [41], where the BDP is being considered to be produced by cosmic ray scattering.Instead of the BDP-electron scattering, the focus of that work is on the hadronic interaction of the BDP.The hadronic form factor suppresses the interaction probability at large momentum transfer, in which case the distortion of the flux becomes most pronounced in the intermediate energy regime.Combined with the difference in the initial velocity distribution of the BDP flux, this leads to substantially different predictions for the BDP energy spectrum in the detector after the signal propagated through Earth.
A possible future extension of our work includes calculating local geographic effects on the BDP signal.It would also Let us first consider the elastic scattering, assuming the electron is free and at rest.For a contact interaction with F(q) = 1, in the nonrelativistic approximation of the final state electron, one can derive the differential cross section as a function of the recoil energy, E R = k 2 /2m e , as which is a flat distribution for E R < 2µ 2 v 2 χ /m e .We now include the effect of the binding energy and consider the ionization process.This requires the velocity of the incoming BDP particle and its mass to be large so that the momentum transfer is enough to trigger the ionization.In this limit, the integration range (q − , q + ) in Eq. ( A12) is (E R /v χ , 2m χ v χ ) at leading order, which effectively becomes (0, ∞) for the integral.Thus below the energy cutoff in Eq. (C1), the ratio between the differential cross section for the ionization of an electron with (n, ), i.e.Eq. (A12), and that of a free electron scattering, i.e.Eq. (C1), can be written as This is defined to be the effective electron number for a given atomic level.The numerical results for xenon are shown in Fig. 9.The results converge to the number of the electrons for that shell.Notice that it requires the sum of the final state angular momentum number to a large number to properly mimic the final state wave function when E R is large.
This result is also consistent with the kinetic distribution of the ionization form factor in the large recoil energy limit.In Fig. 10, we show the result for xenon.In the limit of a large recoil energy, the form factor peaks at E R q 2 /(2m e ), which recovers the kinetic distribution of elastic scattering.The simulation starts with four input parameters: m χ , σ e , v 0 χ , and N, which correspond to the BDP mass, the BDPelectron elastic scattering cross section evaluated at q = αm e , the initial incident BDP velocity, and the number of simulation events, respectively.We consider a BDP flux from the GC or the Sun.We set the direction of the z-axis to always coincide with the direction of the incoming BDP flux, see Fig. 3 for an illustration.
To generate the initial BDP flux evenly distributed on the plane orthogonal to the z-direction, we first draw a random value from the uniform distribution [0, R 2 E ) where R E is the Earth's radius, then we define ρ xy as the square root of the previously generated number.Next, we draw a random azimuthal angle φ E from a uniform distribution [0, 2π).With this choice, we can calculate the position of each BDP particle entering the Earth in the Cartesian coordinate (x, y, z) as,

Propagation inside the earth
Next we consider the propagation of the BDP inside the Earth.The simulation procedure of the BDP propagation is shown on the flow chart in Fig. 4. For each iteration, we first calculate the mean free path, l ion fp , of the BDP particle traveling inside the Earth, using Eq.(10).Next, we determine the travel distance between two successive scatterings in the Earth's mantle or core according to an exponential probability distribution, Combining l with the final velocity calculated from the previous step, we obtain the endpoint of the trajectory in each iteration where the scattering occurs.It becomes subtle when the trajectory hits the mantle-core border before it ends.If this happens, we set the location where the trajectory hits the mantle-core border as the new starting point x i of this iteration while the velocity is left unchanged.

Reconstructing the BDP flux distribution
For each scattering, we determine the electron recoil energy E R and the momentum transfer q using the ionization form factor.According to Eq. ( 6), q × K(E R , q) describes the joint probability of E R and q in an ionization process.In Fig. 2,  K(E R , q) for all elements listed in Table I are demonstrated.It is worth noting that when the binding energy of an electron is much smaller than the BDP kinetic energy, E R and q become closely correlated, and the most probable values are those satisfying E R q 2 /(2m e ).We use the generalized acceptancerejection method [43] to acquire a (E R , q) pair corresponding to the probability distribution given by q × K(E R , q).
A further dynamical constraint on the (E R , q) pair is imposed, E R ≤ qv i − q 2 /(2m χ ), which is equivalent to the condition of q ∈ (q − , q + ) used in Eq. ( 6).
The magnitude of the BDP final velocity, v f , after the scattering is written as, The polar angle of the final velocity, α, respect to the direction of the initial velocity can be calculated as To fully determine the direction of the BDP final state after the scattering, we sample the azimuthal angle β following a uniform distribution [0, 2π).The final state of the BDP in each iteration is thus determined, including its location x f and velocity v f .These will be used as the inputs for the next iteration.There are two conditions for the simulation to stop.First, there is a minimal recoil energy required to ionize an electron in xenon.It can be used to set a lower bound for the BDP velocity as v ion min = 2E min R /m χ , with E min R ≡ 10 eV is set in this study.Thus the threshold velocity in our simulation is chosen to be the maximum value between the DM virial velocity, i.e. 10 −3 c, and v ion min .Second, the BDP may reach the surface r = R D before its velocity becomes smaller than the threshold velocity, in which case it is no longer relevant.Under both conditions, the simulation of the event will be stopped automatically.
After the simulation, we define the "hit events" as the ones which reach the surface r = R D .We collect the velocity and position of each event.If the BDP-electron interaction is strong, the BDP loses its energy rapidly in the Earth, and the BDP particles can only penetrate the r = R D sphere at most once, i.e. when they just enter the Earth.On the other hand, when the interaction is weak, most of the BDP particles travel across the r = R D sphere twice, this leads to a doubling of the number of events, to 2N.In this subsection, we explain how to convert the distribution of "hit events" to the BDP velocity distribution that can be used to calculate the event rate in a detector.
In Fig. 11, we show an example of the "hit event" distribution projected to the x − y plane in both the near side (z < 0) and the far side (z ≥ 0) of the Earth.On the near side, the events are almost equally distributed on the x − y plane, which is consistent with the initial setup in Sec.(D 1).However on the far side, the events are more densely distributed near x 2 + y 2 R E where R E is the radius of the Earth.
For the parameter space we are interested in, the transverse component of the BDP velocity is much smaller than the one along the z−axis after the propagation, thus the event rate in a detector can be approximately calculated using the modified BDP flux along the z−axis.With a proper normalization, the "hit event" rate per area on the r = R D sphere is simply related to the modified BDP flux by a factor of 1/cosθ.

a. Reconstruction of the general velocity distribution
In more general cases, the BDP can reach the detector from all directions.In this section, we study the conversion from the "hit event" rate per area to the general velocity distribution.
The number density of BDP particles with velocities within d 3 v χ is where n χ is the BDP number density and f ( v χ ) is the 3velocity distribution.For an infinitesimal area element d s, the rate of particles hitting this area with velocities within d 3 v χ is This can be used to calculate the the differential BDP event rate at a detector, where N d is the number of target atoms in the detector.
Appendix E: Constraints on non-distorted BDP The first results from the XENONnT experiment [4] for E R < 30 keV ruled out the previous excess of the electron recoil events seen by XENON1T [1], thus providing the most stringent constraints on the BDP scenario.
The BDP model discussed in the main text contains several parameters, including the mass m χ , the incoming velocity v 0 χ , the total flux Φ 0 , and the cross section for scattering on free electrons σ e .To simplify the analysis, we consider the BDP flux without the effect of Earth's shielding.Under this assumption, the product σ e ×Φ 0 is degenerate, and the incoming direction of v 0 χ is not relevant.The ratio of σ e and Φ 0 becomes significant once σ e is large enough to cause shielding associated with the propagation inside the Earth, as discussed in the main text.Therefore, one is left with three free parameters: σ e × Φ 0 , m χ and v 0 χ .For each pair of m χ and v 0 χ , we calculate the corresponding −2∆ ln L = −2 ln(L S +B /L B ) using the XENONnT data [4], where L is the likelihood function, and derive the 90% upper bounds on σ e × Φ 0 via requiring −2∆ ln L = χ is required for the electron recoil energy to be above the threshold.On the other hand, with both m χ and v 0 χ sufficiently large, the elastic scattering limit with dR/dE R scaling as (v 0 χ ) 2 is recovered, as discussed in the previous appendix.Based on the results in Fig. 12, we choose our benchmark parameter in the main text to be σ e × Φ 0 = 10 −36 s −1 .

DarkFIG. 4 :
FIG.4:The flow chart of the Monte Carlo simulation.The indices i and f are used to denote the initial and final states during each iteration.In each sampling of l, we check if the path passes through the mantle-core border; if true, we use the coordinates where the path hits the mantle-border as the new starting point x i , and we keep v i unchanged; this step is not shown in the flow chart.See Appendix D for a more detailed discussion.

FIG. 5 :
FIG. 5:The ratio of the mean free path of the elastic scattering l free fp and that of the ionization l ion fp , as a function of the BDP velocity.It converges to 1 for large v χ .We take m χ = 50 MeV as a benchmark.

FIG. 6 :
FIG. 6: The BDP velocity distribution upon reaching the detector, after its propagation through the Earth.The three rows correspond to different choices of the BDP mass, the initial BDP velocity and the cross section, respectively.The value of the initial BDP velocity is indicated by the cyan dotted line.We use the color bar to characterize the normalized differential flux distribution as a function of the BDP velocity.∆Φ is the flux within each bin of Log 10 [v χ /c].

21 FIG. 7 :
FIG. 7: Top panel: Value of cos θ as a function of the sidereal time for XENONnT, PandaX and LUX-ZEPLIN, respectively, assuming the BDP flux originates in the GC.Bottom panel: The value of cos θ as a function of UTC for the XENONnT experiment on four different days of the year, assuming the BDP flux arrives from the Sun.

FIG. 9 :
FIG. 9:The effective electron number for various xenon shells in the large velocity and heavy BDP limit.In particular, n n eff ≡ n n eff , as defined in Eq. (C2).

FIG. 10 :
FIG.10:The atomic ionization form factor K(E R , q) for xenon.The white solid line corresponds to E R = q 2 /(2m e ).

40 FIG. 11
FIG. 11: "Hit event" distribution projected on the x − y plane.We choose m χ = 50 MeV, σ e = 10 −31 cm 2 and v 0 χ = 0.02 c.The left panel shows the near side (z < 0) while the right panel shows the far side (z ≥ 0).The initial event number is N = 10 5 .Both x and y axes are equally divided into 100 bins.

FIG. 12 :
FIG.12:The 90% exclusion regions on the product σ e × Φ 0 for several values of m χ and a range of v 0 χ , plotted using the first results of the XENONnT experiment[4].