Studies of the Earth shielding effect to direct dark matter searches at the China Jinping Underground Laboratory

Dark matter direct detection experiments mostly operate at deep underground laboratories. It is necessary to consider shielding effect of the Earth, especially for dark matter particles interacting with a large cross section. We analyzed and simulated the Earth shielding effect for dark matter at the China Jinping Underground Laboratory (CJPL) with a simulation package, CJPL Earth Shielding Simulation code (CJPL\_ESS), which is applicable to other underground locations. The further constraints on the $\chi$-N cross section exclusion regions are derived based on the studies with CDEX experiment data.


I. INTRODUCTION
Weakly interacting massive particles (WIMPs, denoted as χ) as a dark matter candidate, have not been observed in direct detection experiments, which provided exclusion limits [1][2][3][4][5]. Direct detection experiments are usually located in deep underground laboratories such as China Jinping Underground Laboratory (CJPL) [6], SNOLab [7], and LNGS [8], which means that the dark matter particle will go through a few kilometers of rock before reaching the laboratory. Generally, the transport of dark matter is not be affected by rocks, but the shielding effect may not be negligible for a large cross section. Therefore, the flux of dark matter at laboratory may be suppressed when assuming high cross sections. This phenomenon is called Earth shielding effect and has been studied in detail in Refs. [9][10][11][12].
At large scattering cross section with mean free path in rocks less than a few kilometers, the WIMPs will scatter with nucleus and lose kinetic energy before reaching the detectors at underground locations. Increasing scattering cross sections will reduce the WIMP fluxes above a given velocity. According to Refs. [9][10][11][12], the event rate of WIMP-nucleus (χ-N) scattering will decrease as the cross section increases at large cross section. The direct detection experiments at underground laboratory cannot detect dark matter with a too large χ-N cross section, due to the shielding effect of overburden. That means constraints of χ-N cross section is not an upper limit, but an exclusion region with both upper and lower bounds. References [9][10][11][12] have studied the shielding effect in detail, but the geometric model was the entire Earth with the laboratory placed on a spherical shell a few kilometers deep, where the undulating terrain above the laboratory has been ignored. In this way, these works may not be completely suitable for laboratories at traffic arXiv:2111.11243v2 [hep-ex] 9 Mar 2022 tunnel whose overburden are mountains, such as CJPL, LNGS. In this paper, we analyze and simulate the Earth shielding effect for dark matter at CJPL, and we also calculated the exclusion region of χ-N cross section with CDEX [1,[13][14][15][16][17][18][19][20] experiment data.

II. THEORY
In the direct search for dark matter, WIMPs are assumed to interact elastically with nucleons, losing a fraction of their kinetic energy and deviating from their original direction. While dark matter particles go through the rocks and scatter multiple times, their tracks are no longer straight lines but broken lines. Consequently, simulations that only have energy loss taken into account are not completely in line with the actual transport process, a complete Monte Carlo (MC) simulation based on χ-N scattering is necessary.
Mean free path (denoted as λ) is a key parameter for MC simulations. The overburden contains multiple nuclides, so the mean free path can be expressed as where ρ is the density of the transport medium, f i is the mass fraction of the i th nuclide in transport medium, m Ai is the mass of the i th nuclide, and σ χAi is cross section of χ-A i scattering, where A i refers to nucleus of the i th nuclide. For low mass dark matter whose mass m χ is less than 10 GeV/c 2 , the form factor F (q) ∼1, where q is equal to √ 2m A E R , and E R is nuclear recoil energy. In this way, σ SI χAi is equal to σ SI χN (0) A 2 i for χ-N spin-independent (SI) scattering, where σ SI χN refers to χ-nucleon SI-interactions cross section, and µ refers to reduced mass. The mean free path of χ at Jinping Mountain in SI-scattering process was shown in Fig. 1, where the blue line is that of the cross section at 10 −30 cm 2 , and the orange line is that of the cross section at 10 −31 cm 2 . At the same cross section, the mean free path decreases with the increase of χ's mass. Before every scattering, the free path of dark matter particle can be calculated by where ξ follows an uniform distribution between 0 to 1, denoted as ξ ∼ U(0,1). The χ-N scattering is treated as elastic scattering process. In the center of mass (c.m.) coordinate frame, the magnitude of χ's velocity remains the same in scattering process, which is expressed as v χc = v χc , where v χc refers to the velocity of χ in the c.m. frame before scattering, v χc refers to the velocity of χ in the c.m. frame after scattering. In a laboratory frame, the velocity of χ after The mean free path of χ at Jinping Mountain in SI-scattering process, the blue line is that of cross section at 10 −30 cm 2 , and the orange line is that of cross section at 10 −31 cm 2 . The density of the transport medium is set to be 2.7 g/cm 3 , and the composition of the transport medium is adopted as Table I. At the same cross section, the mean free path decreases with the increase of χ's mass.
scattering is where v c is the velocity of the c.m. frame in lab frame, n is the direction of χ after scattering in the c.m. frame. The direction of χ after scattering is isotropic in the c.m. frame, so the scattering angle cosα ∼ U(−1, 1). The scattering angle θ in lab frame can be calculated by (4)

III. SIMULATION
CJPL is located at the center of a traffic tunnel under Jinping Mountain in Sichuan province, southwest China, where the latitude and longitude is approximately (E101.7 • , N28.2 • ). The length of the tunnel is about 17.5 km, and the rock overburden at CJPL is about 2400 m vertically. The altitude of the laboratory is about 1600 m, and the highest point of Jinping Mountain is about 4000 m [6]. The topographic map of Jinping Mountain is shown in Fig. 2. The rocks of Jinping Mountain are mainly marble, whose density is about 2.7 g/cm 3 , and the main elements included are O, Ca, Mg, and C. On the basis of the test result of ore composition, the detailed composition of the rocks is shown in Table. I.
The Earth model was adopt following Refs. [21,22], and the radius (R ⊕ ) was set to be 6400 km. We spliced the Earth and Jinping Mountain as a geometric model and developed an Earth shielding effect simulation program for CJPL laboratory based on the geometric model, so-called CJPL laboratory Earth shielding effect simulation code (denoted as CJPL ESS).
In CJPL ESS simulation framework, coordinate system is established with the laboratory stationary and at the origin, as shown in Fig. 2. The initial dark matter particles are generated on a sphere concentric with the Earth with a radius of R 0 = R ⊕ +h, where h = 3 km, the center of which defined as point O e . To improve calculation time efficiency, nonuniform sampling is applied to the sampling of initial position p 0 of χ. For χ with given m χ and scattering cross section, the β, which is the angle between p 0 O e and y axis, was sampled following the probability density distribution function where N is the normalization factor, x = cosβ, λ 0 is the mean free path of χ in Jinping Mountain. While λ 0 R 0 , the initial particles are uniformly generated on the spherical surface. As λ 0 decreases, the initial particles are more likely to be generated in locations near the laboratory. To ensure that the simulation result is unbiased, the weight of each particle is set to be 1/f (x). The initial velocity is got by where v χ,gal is the velocity of χ in the Galaxy, the distribution of which follows Maxwell distribution [23], v lab is the velocity of the laboratory in the Galaxy, which is equal to v earth + v rot , sum of the velocity of Earth ( v earth ) and the velocity of the Earth's rotation at the laboratory ( v rot ). v lab is related to sidereal time and the latitude and longitude of the laboratory, which can be calculated as Ref. [24]. Particles are collected on a horizontal plane with a radius of r c centered on the laboratory, and the particle that went through the plane are counted to calculate the velocity distribution. r c is set to be 100 km while λ 0 > 100 km and set to be 0.5 km, while λ 0 < 50 km and decreases from 100 to 0.5 km as λ 0 from 100 to 50 km. The flow chart of CJPL ESS code is shown in Fig. 3, and the code can output the velocity distribution of χ at CJPL with a given m χ , σ χN and time (or a period of time). Parallel computing based on message passing interface and bias sampling method [25] have been used to speed up the program running time.
According to Eq. (4), the scattering angle for WIMP interaction on nucleus can be large, especially for light WIMP. Consequently, the particle will almost completely lose its initial direction after multiple scatterings. The track of two simulated events are shown in Fig. 4. The particle at m χ of 0.05 GeV/c 2 and SI cross section of 10 −30 cm 2 reached the lab after hundreds of scatterings, that is because the particle can only lose a small part of its total kinetic energy by each scattering, and the maximum rate of energy loss at each scattering is 4µ χN /m N . In this way, the direction of particle at lab has a very weak correlation with the initial direction, and the distance traveled is much greater than the displacement.
Using CJPL ESS code, we simulated the velocity distribution of χ at different masses and cross section at CJPL. The velocity distributions of χ normalized to 0.3 GeV/cm 3 at mass of 0.2 GeV/c 2 and different cross section are shown in Fig. 5. It is found that Earth shielding effect becomes noticeable at σ SI χN > 10 −36 cm 2 . The flux of χ was slightly enhanced at σ SI χN from 10 −36 cm 2 to 10 −31 cm 2 , because more particles that would not otherwise go through the laboratory after several times scatterings. As σ SI χN increases from 10 −31 cm 2 , the flux of χ with velocity greater than the cutoff velocity at lab is getting smaller and smaller.  The black dashed line is velocity distribution of χ without Earth shielding effect taken into account, and the solid lines are that at different SI cross section. In simulation, cutoff velocity of particle was set to be 100 km/s.

IV. RESULT
CDEX experiment is located at CJPL, and committed to searching dark matter with high purity germa-nium detector [1,[13][14][15][16][17][18][19][20]. CDEX-1B data and CDEX-10 data are used to Earth shielding effect analysis in this paper. CDEX-1B experiment has for its target a singleelement p-type point contact germanium (PPCGe) detector cooled by cold finger, with a fiducial mass of 939 g. A shielding system was built to shield the background radiation from outside the detector, constructed with 20 cm copper, 20 cm PE and 20 cm lead [17]. Different from CDEX-1B, CDEX-10 experiment used three tripleelement PPCGe strings (C10A, B, C) directly immersed in liquid nitrogen (LN 2 ), and the shielding system was constructed with LN 2 and 20 cm copper [1]. Compared to that of the overburden, the shielding effect of manually built shielding system to dark matter is negligible, and has been disregarded in simulation process.
In this work, both χ-N elastic scattering and χ-N inelastic scattering via Migdal effect [26][27][28] are taken into account. The event rate of χ-N elastic scattering is expressed as where E R is the nuclear recoil energy, N T is the number of target nuclei per unit detector mass, ρ χ is the density of dark matter, m χ is the mass of χ, v is the velocity of χ in lab frame, f v (v, σ) is velocity distribution of χ, which is the output of CJPL ESS code. Considering χ-N SI scattering and spin-dependent (SD) scattering, the differential cross section can be expressed as where µ χA refers to reduced mass between χ and nucleus A, σ SI χA and σ SD χA are SI and SD cross section between χ and A, F is the form factor. We will discuss the SIscattering and SD-scattering case separately later.
Following Ref. [28], the event rate of χ-N inelastic scattering via Migdal effect is calculated by where E det refers to the energy CDEX PPCGe detectors can detect, which is equal the summation of nuclear recoil energy and electron recoil energy, denoted as Q nr E R + E EM , where Q nr is the quenching factor and was treated following Ref. [28]. Considering SI scattering, the expected energy spectrum at χ of different mass and cross section have been calculated using Eqs. (7) and (9). The standard WIMP galactic halo assumption and conventional astrophysical models are used, with χ density ρ χ set to 0.3 GeV/cm 3 [23]. The energy spectrum at mass of 0. (eV), where E is expressed in keV.
6. As the SI cross section increasing from ∼ 10 −36 cm 2 to ∼ 10 −30 cm 2 , the count rate of χ-N scattering event increases first and then decreases, and decreases sharply at cross section from 1 × 10 −30 to 3 × 10 −30 cm 2 . Migdal effect has been considered in the exclusion line evaluation since it can significantly extend the lower sensitivity reach to light WIMPs at a given detector threshold. Energy loss due to Migdal effects due Earth attenuation is negligible, since the probability is extremely low compared to that for elastic scattering. For example, at m χ = 0.2 GeV/c 2 , the ratio for mean energy loss of each χ-N scattering via Migdal effect to that of χ-N elastic scattering is ∼ 20:1, the ratio of event rate is ∼ 1:20000, and the combined ratio of (energy loss × event rate) is ∼ 1:1000.
The exclusion region at 90% confidence level (CL) was calculated with CDEX-1B and CDEX-10 data, using binned Poisson method [34] and unified approach [35], respectively. The results are shown in Fig. 7 (a). Although the background and lower limit of the CDEX-1B and CDEX-10 are different, the upper limit of exclusion region of the two experiments are consistent, that's because the count rate of χ-N scattering decreases rapidly as σ SI χN increase to approach the upper limit. In this way, the upper limit of the exclusion region in Fig. 7 can basically represent the largest cross section for dark matter at the corresponding mass region that can be detected by all experiments at CJPL, within the χ-N SI scattering and context of standard spherical isothermal galactic halo model. In particular, it can be seen in Fig. 7(a) that the detailed topography of the Jinping Mountain overburden adopted in this work gives rise to substantial shifts to the upper exclusion bounds (red curve) at the (a) The exclusion region at 90% CL with Earth shielding effect considered, derived by CDEX data [1,17]. The blue region enveloped by blue solid line ("CDEX-1B ME") is the exclusion region on χ-N scattering via Migdal effect derived by CDEX-1B data, and the red region enveloped by red solid line ("CDEX-10 ME") is that derived by CDEX-10 data. Corresponding, the regions enveloped by dashed line ("CDEX-1B NR" and "CDEX-10 NR") are the exclusion regions on χ-N elastic scattering derived by CDEX data. The red region and blue region have the same upper limit. The parameter space has been further categorized into three regions by black solid line and black dashed line. As a comparison, the green line is the constraint calculated with flat topography. The dotted lines ("CDEX-1B ME w/o es" and "CDEX-10 ME w/o es") are the constraints without considering Earth shielding effect. The constraints labeled with "CDEX-1B NR" and "CDEX-1B ME w/o es" are from previous work [17,28]. (b) Summary of results from experiments with published exclusion regions, including CRESST [29], DarkSide [2], EDELWEISS [30], XENON-1T [29,31], CDEX, and constraints from the X-ray Quantum Calorimeter experiment (XQC) [32,33]. Among these experiments, XENON-1T and CRESST-III are located at LNGS, EDELWEISS-Surface and CRESST-Surface are located at the surface of the Earth, and the Earth shielding effect has been taken into account. The upper boundary of "XENON-1T ME" is lower than that of CDEX is because Ref. [31] adopted a more conservative calculation method. FIG. 8. The exclusion region at 90% CL for χ-N SD scattering of neutron-only case with Earth shielding effect considered, derived by CDEX data [1]. The momentum-transferdependent spin-structure function of 73 Ge was taken from Ref. [37]. The constraint of XENON-1T was taken from Ref. [31], where a more conservative calculation method has been adopted in evaluation of Earth shielding effect.
light mass region, relative to those (green curve) from models with flat topography [9]. The deviation is ∼45% for WIMP mass at 0.1 GeV/c 2 . Since it (red curve) is more precise applying detailed topography, we quote it as our official result. In addition, Fig. 7 (a) is further categorized into three regions. Region I is where the cross sections are so large that no WIMPs above the cutoff velocity (<1% in flux) would survive. Region III is where all (>99%) WIMPs will arrive laboratory without scattering with the ordinary matter nuclei as they travel through Earth's rock, and region II is between region I and III. Regions I and II are the regions where the Earth shielding effect works. In this way, at m χ < 1 GeV/c 2 , Earth shielding effect should be considered in calculation of exclusion curve.
The Earth shielding effect of χ-N SD scattering has also been simulated and evaluated. The odd-even nucleus with the largest mass fraction in Jinping Mountain and the Earth have been considered, and the spin expectation values, S n and S p , follow Ref. [36]. The constraints of neutron-only case was calculated and shown in Fig. 8. The proton-only case has not been considered, because the expectation value, S p , of 73 Ge approaches zero. The spin expectation values used in SD analysis are shown in Table II.

V. DISCUSSION
In this paper, we have introduced a Monte Carlo simulation package for Earth shielding effect. This package was based on χ-N elastic scattering and developed to be applied to simulate the Earth shielding effect of CJPL (or other underground laboratory that is under a mountain). Compared to previous works [9][10][11][12], the current work includes an effect due to topography of the underground locations. This improves on the accuracy of the exclusion bounds. To speed up the simulation, nonuniform sampling and bias sampling method have been adopted. The simulation package is also applicable to other underground locations, after revising the input parameters on topography and compositions. More optional physical processes will be added to this framework, including the inelastic scattering process, DM-electron scattering and so on.
In the framework of CJPL ESS, different form factor F (q) can be selected and adopted according to different physical process scenarios. The impact of form factor depends on the energy scale of the momentum transfer during the scattering. In this χ-N SI-scattering and SDscattering analysis with CDEX data, the approximation of omitting the form factor [F (q) ∼1] will bring an error of no more than 1%. This approximation cannot be applied to χ-N scattering process with large momentum transfer, such as cosmic ray boosted sub-GeV dark matter searches in our subsequent work [38], which has a great impact for different ways to handle the form factor.
The studies of Earth shielding effect allow us to derive the upper bounds of the exclusion region on the χnucleon scattering cross section, which will make the constraints from underground experiments more complete than the previous results which only presented one exclusion line. In this work, the data of CDEX experiment at CJPL are used in analysis and new constraints on χ-N SI scattering and SD scattering were achieved with the Earth shielding effect taken into account. Further more physics analysis with considering the Earth shielding effect will be carried out, such as modulation effect analysis, cosmic ray boosted sub-GeV dark matter searches [38] and so on.