Simulation study on cosmic ray background at large zenith angle based on GRANDProto35 coincidence array experiment

Xiang-Li Qian?Xu Wang ?Hui-Ying Sun ?Zhen Wang ?Olivier Martineau-Huynh

Abstract Neutrino detection in the 100 PeV energy region is the ultimate means of studying the origin of ultra-highenergy cosmic rays, in which the large radio detection array giant radio array for neutrino detection (GRAND)project aims to use to decipher this century-old problem.The GRANDProto35 compact array is a microform of 35 radio prototype detectors for the GRAND experiment,which verifies the reliability of GRAND performance through operation, and data analysis of the prototype detectors. As radio detectors are a novel development in recent years, and their indexes need to be verified by traditional detectors, the GRAND Cooperation Group designed and constructed the GRANDProto35 coincidence array composed of radio detectors and scintillation detectors. This study simulated the changes in detection efficiency, effective area, and event rate of cosmic rays with zenith angle based on this coincidence array. The study found that the 1017 eV energy region is sensitive to GRANDProto35 detection. When the energy exceeded 1017 eV, the array detection efficiency could reach more than 95% and the effective area was up to ～2 × 106 m2.A simulation study on cosmic ray events with large zenith angles showed that the event rate detected by the array decreased significantly with increasing zenith angle, and the event rate of cosmic rays was approximately 0.1 per day for a zenith angle of 75°.This serves as the background pollution rate for neutrino observation caused by largeangle cosmic-ray events, providing an important reference for further experiments. The study results will be verified after the joint operation of the coincidence array.

Keywords GRANDProto35 · GEANT4 · Scintillation detector · Cosmic ray

1 Introduction

Cosmic rays were discovered in 1912 by the Austrian physicist Victor Franz Hess using an ionization chamber in a high-altitude balloon experiment,for which he shared the 1936 Nobel Prize in Physics. Although more than 100 years have passed since, the origin, acceleration, and propagation of cosmic rays are yet to be solved.Therefore,they were listed as one of the 11 great scientific questions for the new century by the American National Research Council in 2002 [1],and identified as a twenty-first-century research frontier by the National Science and Technology Council in 2004 [2]. In fact, since the cosmic ray was discovered,it has always been one of the frontiers in highenergy astrophysics. The range of energies starts at approximately 107eV and reaches 3 × 1020eV for the most energetic cosmic ray ever detected. As a first approximation, the energy spectrum of cosmic rays dN/dE can be fitted over long ranges with a power law E-αwith a spectral index of α=-2.7. The cosmic ray flux drops off with more than 30 orders of magnitude as the energy increases by 11 orders of magnitude. The overall cosmicray energy spectrum and composition are shown in Fig.1,showing four obvious spectral breaks, which are usually referred to as the knee at ～4 PeV [3],the second knee at 100 PeV [4], the ankle at EeV [5], and the Greisen–Zatsepin–Kuzmin cutoff at dozens of EeV [6, 7]. It is now accepted that the transition from galactic to extragalactic origin of cosmic rays occurs at energies of the knee. Significantly, a primary cosmic ray with energy of approximately 3×1020eV has been detected [8], which is far higher than the highest-energy (1012～1013eV) protons that have been produced in terrestrial particle accelerators.What are the ultra-high-energy cosmic-ray (UHECR)sources? How do cosmic rays accelerate cosmic rays?What happens during the journey of cosmic rays to the Earth? These questions are of great interest in the field of higher-energy cosmic rays.

Neutrino detection provides the only direct method of studying the origin and acceleration of cosmic rays.Direct detection of charged cosmic ray particles is not sufficient to find the sources, as their directions change as they travel through interstellar magnetic fields. Ultra-high-energy(UHE) gamma rays and neutrinos are produced by UHECRs interacting inside their sources.However,gamma rays interact with the cosmic microwave background,causing an energy cascade down to GeV～TeV,where it is more difficult to disentangle them from gamma rays produced in unrelated phenomena. During propagation, neutrinos travel unimpeded, their energies unaffected by interactions; therefore, a technique using only UHE neutrinos is the most promising way to probe sources of UHECRs [10]. In 2013, 2 events around 1 PeV deposited energies and an additional 26 events in the 30–400 TeV energy range were observed in the IceCube experiment [11, 12], which ushered in the era of high-energy neutrino astronomy. To date, IceCube has observed approximately 80 neutrinos,the arrival directions of which are consistent with being isotropic,indicating extragalactic origin [13]. Due to the limitations of the effective area of IceCube, neutrinos above 10 PeV energy have not been found,whereas UHECRs produce neutrinos in the 100 PeV energy region. Therefore, experiments with higher sensitivity are needed to decipher their origin. The Giant Radio Array for Neutrino Detection (GRAND) experiment was designed for this purpose.It is a large neutrino observation array planned in China, aiming to study and discover the sources of UHECRs [14, 15]. For the experiment, radio detectors are used.Because the radio detector is novel,and its performance parameters need to be tested by mature traditional detectors, the GRAND Cooperative Group established the coincidence array GRANDProto35, consisting of 21 scintillation unit detectors and 35 antenna detectors.The present study explores its performance using a simulation method.

Fig.1 (Color online) All-particle energy spectrum of cosmic rays and the energy spectrum of each component [9]

The structure of this paper is as follows: Section 2 describes the experimental instruments, Sect. 3 introduces the simulation process of the atmospheric shower and detector, Sect. 4 provides the simulation results, and finally, there is a summary of this study.

2 Experiment

The GRAND experiment is to be built at the foot of a high mountain,which will as a target for the observation of neutrinos. It is designed to detect radio signals generated by extensive air showers of high-energy cosmic rays and neutrinos. GRAND is expected to be completed around 2030 and will consist of 200,000 radio antennas operating in the 50–200 MHz range of the low-frequency band, with a total area of approximately 200,000 km2, significantly larger than the total detection area of IceCube. The GRANDProto35 prototype coincidence array, which is the first-stage program of the GRAND experiment, consists of a tri-polarization antenna array and a plastic scintillation detector array based on the detection project of ‘‘the first universe light’’ 21 CentiMeter Array (21CMA) [16] and Tianshan Radio Experiment for Neutrino Detection(TREND) [17], which is located in the Tianshan Mountains in the Xinjiang Province of China. The tri-polarization radio antenna in the coincidence array is mainly used to detect radio radiation signals, whereas the scintillation detector array, with its high sensitivity and high detection efficiency of cosmic-ray observation, conducts joint observation with a radio antenna array and cross-validation of the observation results and radio observations, so as to improve the radio detection method as well as the efficiency of radio self-triggered detection of cosmic rays.

Figure 2 shows the coincidence array distribution of the GRANDProto35 experiment. The antenna array and scintillator array are arranged according to a rectangle based on the 21CMA array.The whole array is 800 m long from east to west and 4200 m long from north to south. Due to the geomagnetic field effect, the radio radiation signals generated by cosmic rays perpendicular to the geomagnetic field direction will be more concentrated. Therefore, this array arrangement method is more conducive to detecting air shower events from the north. The radio antenna has three pairs of mutually perpendicular butterfly antenna arms to sample the electric field components in the three polarization directions. The scintillation detector is mainly composed of a plastic scintillator, photoconductive box,and photomultiplier tube (PMT). The internal structure is shown in Fig.3.The plastic scintillator is an EJ-200 model with dimensions 707 mm × 707 mm × 20 mm. A Hamamatsu PMT R7725 with a diameter of 2’’ (5.04 cm) is installed on the bottom of the photoconductive box and an interface is reserved for testing next to the PMT.To further verify the detection efficiency of the radio antenna, the scintillation unit detector is placed with a tilt to the north,with zenith angle 50°.The DAQ chain and trigger logic of the scintillator array are fully independent of the radio array;analog signals collected by the PMTs are digitized in real time. If a pulse is observed simultaneously in the signals from three or more scintillators, the signals are written to disk.Offline comparison between scintillator and radio data will quantify the radio-detection efficiency and contamination by background events. Therefore, the performance of the scintillation detector and the ability of the array to detect cosmic rays are very important to the whole coincidence array. It is necessary to conduct a detailed simulation study.

3 Simulation

3.1 Air shower simulation

High-energy elementary particles and atomic nuclei enter the atmosphere and interact with nuclei of the air,usually several tens of kilometers high. In such collisions,many new particles are usually created and the colliding nuclei evaporate to a large extent. The new particles may collide with further nuclei of the air and eventually produce a large number of secondary particles, usually including hadrons, electrons, photons, and muons, and spread over several square kilometers of area according to a certain distribution. The number of secondary particles increases rapidly as the shower or cascade of particles moves downward in the atmosphere. The development of this ensemble of charged and neutral particles in the atmosphere is called an extensive air shower(EAS).CORSIKA(COsmic Ray SImulations for KAscade) is a detailed Monte-Carlo (MC) program to study EASs in the atmosphere initiated by photons, protons, nuclei, and various other particles. It has been widely used by international experimental groups to simulate the air shower process.For our simulation, we adopted CORSIKA version 7.6400 [18], within which the EGS4 model was used for electromagnetic interaction,QGSJET II-03 for high energy(＞80 GeV)hadronic interaction,and the FLUKA software package for low-energy hadronic interaction [19]. For primary cosmic rays, the proton component was mainly considered,with a simulated energy spectrum from 1016to 1018eV, zenith angle from 45° to 80°, and azimuth angle from- 10°to 10°,such that the primary cosmic-ray events occurred in the same direction as the coincidence array.Considering the CPU consumption and statistical error of samples in high-energy event simulation, 1×105primary protons were generated and the first interaction occurred at an altitude of several tens of kilometers. An energy cut of 0.1 GeV was set for hadrons and muons and 0.001 GeV for electrons and photons.These secondary particles interacted with the scintillation detector.For the simulation of the unit detector and the whole array, we developed a G4GRANDProto35 simulation program based on GEANT4, a simulation software package widely used in high-energy physics experiments [20, 21].

Fig.2 (Color online)Coincidence array distribution of the GRANDProto35 experiment composed of 35 radio antennas and 21 scintillation detectors

Fig.3 (Color online) Structure of scintillation detector

3.2 Unit detector simulation

The G4GRANDProto35 simulation program simulated the entire scintillator array based on the detector distribution,as shown in Fig.2.In fact,as long as one unit detector simulation is optimized correctly, the other detectors can be completed by replication for the entire array.In the unit detector simulation, a 1:1 real simulation was carried out according to the structure shown in Fig.3.

The luminescence properties of the EJ-200 plastic scintillator (ELJEN Technology [22]) used by the scintillation detector are shown in Table 1. In the G4GRANDProto35 program, a real simulation was conducted based on the table. When charged particles pass through a scintillation crystal,scintillation light is released isotropically with random polarization. According to the relation between the optical refraction and the refractive index of the scintillator material, when emitted photons strike the smooth interface of the crystal, light with an incident angle greater than the critical angle (38°) generates total reflection in the scintillator, and light with anincident angle less than the critical angle can exit the scintillator and enter the photoconductive box. To reduce the impact of total reflection and increase the light output,the surface of the scintillator in the direction of the PMT is sanded (rough surface), so that a photon with a large incident angle can also emit the scintillator after several diffuse reflections. Thus, the scintillator has one polished surface, one rough surface, and four machine-cut edges.Tyvek 1082D1http://www.dupont.com, a high-reflectivity material, was affixed to the outside surface of the scintillator and the inner wall of the photoconductive box to significantly improve the efficiency of the PMT for receiving photons.

Table 1 Parameters of EJ-200 plastic scintillator

Both the simulation of the scintillator interface and the simulation of photon propagation in the photoconductive box mainly involve optical properties. Therefore, the UNIFIED model for optical simulation in GEANT4 is adopted in this simulation [23]. The UNIFIED model assumes that the surface of the medium is composed of a number of micro-facets and defines the whole plane as the average surface. Each micro-facet has its own normal direction, and all the micro-facet normals are averaged to obtain the average surface normal of the whole plane. The angle between a micro-facet normal and the average surface normal, α, is assumed to follow a Gaussian distribution of standard deviation of σα. The UNIFIED model allows σαto indicate the surface roughness values over a wide range. The model also allows a reflection coefficient,RC,to be specified to simulate an external diffuse reflector.When a beam of photons is incident on the medium surface, four possible reflection types may occur: specular reflection (of a micro-facet and of the average surface),backward reflection, and Lambertian reflection. Correspondingly, the parameters for the description of reflectivity include Csl(specular lobe constant of a micro-facet),Css(specular spike constant of the average surface), Cbs(backscatter spike constant), and Cdl(diffuse lobe constant). These four parameters satisfy the following equation:

Figure 4 is a schematic diagram of photon propagation in the UNIFIED model. The values of various optical model parameters in this study are shown in Table 2.

Similar to the optical simulation of the scintillator interface, the UNIFIED optical simulation model is also used for the Tyvek reflection material. The reflectivity of the Tyvek reflecting material is measured and found to reach 95% in the sensitive region of the PMT [24], which significantly increases the receiving efficiency of the PMT.After diffuse reflection through the photoconductive box,some of the photons emitted from the scintillator reach the PMT. In the specific simulation process, the quantum efficiency and receiving efficiency of PMT R7725 are parameterized, and then each photon is digitized. Finally,the charge and time information are output and provided for later array simulation analysis.

Fig.4 (Color online)Schematic diagram of photon propagation in the UNIFIED model

Table 2 Parameters of UNIFIED model for scintillator at different interfaces

4 Simulation results

The reliability of the developed simulation software required testing before the physical results of the CORSIKA sample simulation based on the software were obtained. Therefore, the unit detector was simulated and the experimental data were compared first to verify the software reliability.Based on this,the scintillator array was simulated and the array performance simulation results were given.

4.1 Simulation results of unit detector

The unit scintillation detector mainly measures two physical quantities: charge and time-to-digital converter(TDC) time. Charge is mainly used to measure the energy and core position of the air shower, while TDC time determines the direction of the air shower.Therefore,these two quantities were mainly tested for laboratory testing of the unit detector reliability.

A coincidence detection experiment was carried out with a scintillation detector in the laboratory. A small scintillation detector with a side length of 5 cm was placed directly under the standard unit detector as coincidence.The count was triggered whenever the two detectors were triggered almost simultaneously. The coincidence ensured that the standard unit detector detected the secondary particles of cosmic rays and eliminated the noise interference of electronics and other sources. Meanwhile, the number of photoelectrons and TDC time resolution information were recorded.The time resolution is the time interval distribution of the coincidence events recorded when the two detectors were triggered. Figure 5 is a comparison between the number of photoelectrons and time resolution data detected by the unit detector and the MC simulation data, which are consistent with the experimental data. The optimal value of the photoelectrons generated in the unit detector is ～190, and the time resolution is ～ 1 ns. The simulation results are consistent with the experimental results,indicating that the simulation results of the detector are reliable and that the array simulation, as well as the study on physical expectations, can be furthered.

4.2 Simulation results of array detection efficiency

The array detection efficiency is an important feature for measuring the performance of an array. It is defined as the probability that a cosmic-ray event is detected when the air shower core is inside the detector array and the trigger condition of the array is met.Owing to the limited number of array detectors and the 300 m east–west spacing of the scintillation detector,the array trigger rate is relatively low.In the simulation,if the number of photoelectrons collected by the unit detector exceeds 100, this is defined as the triggered hits of the detector. We focus on the detection efficiency of the scintillator array under the conditions of different triggered hits of the detector. Figure 6 shows the variation of array detection efficiency with cosmic ray energy for triggered hits of 1, 2, 3 and 4, respectively.When the triggered hits are 1 or 2, the shower direction of the primary cosmic rays cannot be determined;therefore,it cannot be used as the trigger condition of the experimental array.The detection efficiency under this condition is given here only for a comparison with the detection efficiency under the normal trigger condition.The triggered hits must be greater than or equal to 3 as the real triggering condition.According to the figure,above 1017eV of energy,the detection efficiency under several trigger conditions reached more than 95%, and for energy of 1017.2eV, the efficiency reached 99%, indicating that this energy region is a sensitive region for detection of GRANDProto35.This result also shows that the layout of the scintillator array is appropriate for detecting high-energy cosmic rays.

4.3 Simulation results of array effective area

The effective area is one of the most important variables for characterizing the detection ability of an experimental array. The detection efficiency may characterize the threshold energy of the array and the sensitive region for energy segment detection, while the effective area can show the detection ability of the experimental detector visually, and can quantitatively estimate the expected physical results.The effective area can be calculated using the following formula:

where n is the number of triggered events,N is the number of dropped events, Asis the sampling area, and θ is the zenith angle of the shower events [25]. According to the formula,when Asis the array area and the zenith angle θ is 0°(vertical incidence),the effective area is just the product of detection efficiency and array area. Equation 2 can be used to calculate the effective area of the detector array at different energies of primary cosmic rays. Because we are interested in large-zenith-angle events, in the simulation,the zenith angle for the primal cosmic ray is set from 45°to 80° and the azimuth angle from -10°to 10°, with 0°pointing due north. This direction was chosen because of the layout of the array (see Fig.2), for the purpose of increasing the effective detection area of the unit detector for an air shower.In the simulation,the placement mode of the unit detector was adjusted accordingly. The horizontal and zenith angles of the unit detector were simulated.Figure 7 shows the variation of effective area with energy under two different placement modes for this coincidence array. The effective area when the detector was placed at 50°was clearly better than when it was placed horizontally,as the direction of the simulated primary cosmic ray was greater than 45°.When the normal direction of the detector was closer to this direction, the effective section for the detection of secondary particles was increased, consistent with expectations. When the energy was 1017eV, the effective area reached ～2×106m2, equivalent to the area of the scintillator array. This indicates that all cosmic rays with energy around 1017eV can be detected by the detector array and that the 1017eV energy region is the sensitive region for the array. As shown in the figure, the effective area increased with the cosmic ray energy and was larger than the area of the array,because some high-energy events that dropped outside the array were also detected by the array.

4.4 Estimation of large-zenith-angle event rate

Fig.5 (Color online)Comparison between simulation and experimental results of the unit detector: a for number of photoelectrons; b for time resolution information

The dominant neutrino detection principle in GRAND is described by:a ντmakes a τ by interacting underground or inside a mountain, and the τ exits the surface of the Earth and propagates into the atmosphere. The τ decays into a hadronic mode with a high probability, thus producing an EAS [26]. The ντ-→τ conversion probability strongly depends on the thickness of the surface layer or mountain.The detection efficiency can be significantly enhanced by placing the detector in a mountainous terrain. The GRANDProto35 coincidence scintillator array is used to place the unit detector toward the mountain with a zenith angle of 50° to increase the detection efficiency. We simulate the cosmic-ray event rate detected by the scintillator array with different zenith angles, and also under two different placement modes.The results,presented in Fig.8,show that the event rate decreases gradually with increasing zenith angle.When the zenith angle is 75°,the detector is placed horizontally and with a zenith angle of 50 degrees, achieving detection rates of cosmic rays of 0.01 and 0.1, respectively. On the other hand, the background pollution of neutrinos is mainly caused by the cosmic-ray air shower. The larger the cosmic-ray zenith angle, the thicker the atmospheric layer, and the greater the absorption of secondary particles. At this point, the detector may detect higher-energy primary cosmic rays, which reduces the pollution of the cosmic-ray background to the number of events for neutrino observation. In fact, owing to the mountain block to the original cosmic ray shower, the actual observation background will be less than this value.We look forward to further testing this result through experiments.

Fig.6 (Color online) Variation of array detection efficiency with cosmic-ray energy under different triggered hits

Fig.8 (Color online) Variation of event rate with zenith angle

5 Conclusion

The primary physical objective of the large radio detection array GRAND is to detect neutrinos in the 1018eV energy region and to determine the origins of UHECRs. The GRANDProto35 compact array is the microform of the 35 prototype antenna detectors for the GRAND experiment. To test this novel radio detector, a traditional EAS array consisting of 21 scintillation detectors was formed by the collaborative group, and to further understand the performance of the array, G4GRANDProto35 simulation software was developed to simulate the array. First, the reliability of the simulation software was verified by the charge and time resolution simulation of the unit detector. Based on this, simulated calculations of the array detection efficiency and effective area were performed. The simulation results show that the 1017eV energy region is the sensitive region for GRANDProto35 detection, and that when the energy is greater than 1017.2eV,the detection efficiency of cosmic rays is higher than 99%,with an effective area of approximately 2 × 106m2. Because GRANDProto35 detects neutrinos with the mountain as the target, and large-zenith-angle events are detected, the scintillation detector is placed with a zenith angle of 50° in the experiment and simulation to increase the detection efficiency for large-angle cosmic-ray events.This study simulated the variation in the cosmic-ray event rate with zenith angle, and concluded that when the zenith angle of cosmic rays was 75°and the detector was placed at 50°, the cosmic-ray detection rate was approximately 0.1.It is expected that after the operation of the GRANDProto35 coincidence array, the simulation results of this study, and thus the reliability of the G4GRANDProto35 software developed herein, can be further tested.

Author contributionsAll authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Xiang-Li Qian, Xu Wang,Hui-Ying Sun,Zhen Wang and Olivier Martineau-Huynh.The first draft of the manuscript was written by Xiang-Li Qian and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.