High-energy emissions from neutron star mergers

In 2017, LIGO-Virgo collaborations reported detection of the first neutron star merger event, GW170817, which is accompanied by electromagnetic counterparts from radio to gamma rays. Although high-energy neutrinos were not detected from this event, mergers of neutron stars are expected to produce such high-energy particles. Relativistic jets are launched when neutron stars merge. If the jets contain protons, they can emit high-energy neutrinos through photomeson production. In addition, neutron star mergers produce massive and fast ejecta, which can be a source of Galactic high-energy cosmic rays above the knee. We briefly review what we learned from the multi-messenger event, GW170817, and discuss prospects for multi-messenger detections and hadronic cosmic-ray production related to the neutron star mergers.


Introduction
Binary neutron star (BNS) mergers have been actively discussed as sources of multi-messenger astrophysics for a long time. A close BNS can merge within the Hubble time because the orbit of the BNS decreases due to emission of gravitational waves (GWs). This is confirmed by the observation of the binary pulsar [1].
When the BNS merges, a compact object, either a black hole or massive neutron star, is left. The remnant object accrets surrounding material and releases a large amount of gravitational energy, which is expected to launch relativistic jets. If the jets accelerate electrons by dissipating their kinetic energy, they emit gamma-rays, which can be observed as a short gamma-ray burst (SGRB) [2][3][4]. The r-process elements are also expected to be produced by the neutron star mergers because the neutron stars consist of neutron-rich material [5]. According to the numerical relativity simulation, the mergers create massive ejecta [6], and decay of the radioactive nuclei powers the optical/infrared transients (kilonova/macronova; hereafter macronova) [7][8][9][10]. Optical and infrared observations of SGRB afterglows give some hints of macronovae (e.g., [11]). The jets and ejecta of neutron star mergers form forward shocks through interaction with ambient matter. Electrons are accelerated at the shocks, which produce a broad band afterglow emission through the synchrotron emission [12][13][14]. The BNS mergers can be sources of hadronic high-energy particles. If the jets of SGRBs contain protons, they are accompanied with highenergy neutrinos [15]. The kinetic energy of the jets dissipates in the dissipation region through some plasma processes, such as shocks [16] or magnetic reconnection [17]. * e-mail: szk323@psu.edu The protons are accelerated to very high energy there and interact with the target photons observed as SGRBs, producing pions that decay to the neutrinos. Also, the ejecta of macronovae can accelerate cosmic-ray protons beyond the knee, because they are faster than supernova ejecta which accelerate the Galactic cosmic rays below the knee.
The multi-messenger observation of GW170817 confirmed most of the pictures above [18]. This event was detected by the GWs, radio waves, optical, ultraviolet (UV), infrared (IR), X-rays, and MeV gamma-rays. However, GeV and TeV gamma-rays and neutrinos are not detected, despite that the BNS mergers are expected to emit these high-energy particles. The GW170817 observations give us the physical quantities of the macronova ejecta, which enables us to discuss the hadronic high-energy processes related to BNS mergers in more quantitative manner. In this paper, we briefly review what we learned from GW170817, and discuss future prospects for high-energy neutrino detections and hadronic cosmic-ray production from neutron star mergers. The main topics of this paper (Section 3 and 4) are published in [19][20][21].

GW170817
In 2017, LIGO-Virgo collaborations reported multimessenger detection of a BNS merger event (GW170817; [18]). First, LIGO-Virgo detected GW signals from a inspiraling binary of total mass of 2.8 M and mass ratio of ≥ 0.7. This tells us that the binary consists of two neutron stars [22]. From the GW data analysis, the luminosity distance is estimated to be 40 Mpc. The localization of the event is around 30 degree 2 , which is much smaller than the previous binary black hole merger events owing to the observation by three detectors. 1.7 sec after the merger, the Fermi Gamma-ray Burst Monitor (GBM) and the SPectrometer on board INTE-GRAL Anti-Coincidence Shield (SPI-ACS) detected a short gamma-ray burst, GRB 170817A [23]. This prompt gamma-ray emission confirmed the BNS merger paradigm of SGRBs. We learned that at least some fraction of SGRBs occur through BNS mergers. However, the luminosity of GRB 170817A is L iso 1.6 × 10 46 erg s −1 , which is much lower than the typical SGRBs occurred at cosmological distance, L iso ∼ 10 51 erg s −1 . The emission mechanism of such low-luminosity prompt gamma-rays is still controversial. One possibility is that GW 170817 is an off-axis event of the classical SGRB (e.g., [24,25]). Another is the shock breakout from the macronova ejecta (e.g., [26,27]).
Eleven hours later, the optical counterpart is identified at NGC 4993, a galaxy located around 40 Mpc away from the Earth [28][29][30][31]. The distance is consistent with that obtained from the GW analysis. Initially, the transient evolves rapidly in both luminosity and color, compared to supernovae. The color becomes redder in later time (e.g. [32]). This rapid and redward evolutions are consistent with the theoretical modeling of a macronova powered by radioactivity of r-process elements. The modeling tells us that two-component models fit the data well: the fast-light component with a low opacity (M ∼ 0.01M − 0.03M , V ∼ 0.3c, κ ∼ 1 cm 2 g −1 ) and the slow-heavy component with a high opacity (M ∼ 0.03M − 0.05M , V ∼ 0.1c − 0.2c, κ ∼ 10 cm 2 g −1 ) (e.g., [33][34][35][36][37]). To produce the massive, fast, and low-opacity ejecta, the remnant central object should be a temporal hypermassive neutron star. If the hypermassive neutron star exists for a long time, the strong X-ray and gamma-rays should be detected weeks to months after the merger [38]. Since we did not detect such signals, the hypermassive neutron star should collapse to a black hole after the macronova ejecta is produced.
The X-ray and radio counterparts are detected 9 days and 16 days after the merger, respectively [39][40][41][42][43][44]. The spectrum of the afterglow is consistent with the synchrotron emission from electrons with a single power-law distribution of the index 2.2. The light curve shows the gradual increase in both radio and X-ray bands, L ∝ t 0.7 for more than 100 days after the merger (e.g. [45,46]). This ruled out the classical top-hat jet model seen from off-axis, and a radial or azimuthal structure is demanded. To explain this feature, two possibilities are intensely discussed: the quasi-spherical cocoon (or ejecta) with a radial structure [45,47] and the azimuthally structured jet [46,48]. Resolving the emission region is a smoking gun for distinguishing these two models. 230 days after the merger, the VLBI radio observation detected a superluminal motion of the emission region, which undoubtedly indicates that the emission region has a relativistic velocity [49]. Also, they found that the emission region is so compact that even the VLBI observations cannot resolve it. These results favor the structured jet model rather than the quasi-spherical cocoon model. This supports that GW170817 is a canonical SGRB seen from off-axis. The light curve of the afterglow started to fade around 150 days with a decreasing rate of ∼ t −2.2 [50]. Such a rapid decreasing also favors the structured jet model, because the quasi-spherical cocoon models predict slower decreasing rates.
The higher-energy gamma rays and neutrinos from GW170817 are not detected, despite the intense search by several collaborations [51][52][53][54]. We discuss the neutrino emissions and implications of this result in the next section.

High-energy neutrinos from neutron star mergers
From GW170817, we learned that BNS mergers indeed produce relativistic jets. These jets are expected to produce high-energy neutrinos if they contain protons. Since the emission is strongly beamed toward the jet direction, detection from off-axis events are very challenging. Here, we discuss the detectability of the neutrinos from on-axis events.
After BNS mergers, the jets launched from the central engine interact with the ejecta of macronovae, and the outcomes can be classified into two cases: the successful SGRBs with late-time activities (Section 3.1) and the failed SGRBs with choked jets (Section 3.2). For the successful jet case, the jets successfully penetrate the ejecta of macronovae, resulting in the canonical SGRBs. For the failed SGRBs case, the jets fail to penetrate the ejecta, and bright gamma-rays are not detected from this system. Only the neutrinos are able to come out from the inside of the ejecta. We call these neutrinos "trans-ejecta neutrinos". See [19] and [20] for details of the neutrinos from the successful jet and the choked jet cases, respectively.
We calculate the neutrino fluence using phenomenological formula. We consider a single power-law proton spectrum with a spectral index s = 2: We appropriately take into account the energy-dependent cross-section of photomeson production [55], proton cooling processes (synchrotron, adiabatic cooling, and photomeson production for the successful and choked jets, Bethe-Heitler process for the choked jets), pion cooling processes (synchrotron and adiabatic cooling for the successful and choked jets, proton-pion inelastic collisions for the choked jets), and muon cooling processes (synchrotron and adiabatic cooling for the successful and choked jets). During the propagation from the source to the Earth, the neutrino oscillation changes the flavor ratio, which is calculated using the tri-bimaximal mixing matrix (e.g., [56]) is the neutrino fluence without the oscillation and d L is the luminosity distance.
We discuss the prospects for the neutrino detection coincident with GWs by IceCube and IceCube-Gen2. The expected number of neutrino detection is estimated to be where A eff is the effective area and δ is the declination angle.The effective area for ν µ -induced events with IceCube is given by [57]. For lower energy of 1 PeV, the effective area for up-going + horizontal events is larger than that for down-going events, because the atmospheric muons are shielded by the Earth. For the higher energies, the neutrinos are also blocked by the Earth, so the effective area for the down-going events is larger. For IceCube-Gen2, we use 10 2/3 larger effective area than that for IceCube. Although such a simple scaling might not be a good approximation for the down-going events, this treatment suffices for a demonstration purpose. The detection probability of k neutrino events is given by the Poisson distribution.

Neutrinos from SGRBs
The observed SGRBs are followed by afterglows. The classical afterglow theory based on the forward shock model predicts a decreasing light curve with a single power-law function [58][59][60]. SGRBs are often followed by the afterglow with flat light curves (extended emission for t ∼ 10 2 −10 3 sec, and plateau emission for t ∼ 10 3 −10 4 sec [61][62][63]). X-ray flares are also observed during the afterglow of SGRBs [64]. Since the classical forward shock models have difficulty to explain the time variability of these emissions, they are likely to originate from prolonged central engine activities [65]. The energy fluence of the late-time emissions are comparable to that of the prompt emission [63], so they are important components of this system with regard to energetics. The late time emissions are also expected to produce high-energy neutrinos. Since these components can have a lower Lorentz factor and a lower break energy than those for the prompt emission, the photon density can be higher, leading to the efficient neutrino production. We consider broken power-law spectra for the target photons with indices α = 0.5 and β = 2.0 below and above the break energy, respectively. We calculate the neutrino spectra from the prompt emission, the extended emissions (two cases), the plateau emission, and the X-ray flare with the model parameters tabulated in Table 1. These parameters are obtained from the observations, although some of them are not constrained very well. The resulting physical quantities are tabulated in Table 2. The late-time emissions can accelerate protons up to several EeVs to 100 EeV, depending on the component. Figure 1 shows the neutrino spectra with the baryon loading factor ξ acc = 10. We can see two breaks in the spectra. The first break appears due to the change of the photon spectral index. The second break comes from the pion synchrotron cooling. Owing to the higher luminosity and the lower Lorentz factor, the extended emissions have denser photon fields. Therefore, the extended emission is the most efficient neutrino emitter of the four. We discuss the neutrino detection probability coincident with GWs. Hereafter, we focus on the neutrinos from the extended emissions. Since the neutrino fluences strongly depend on the Lorentz factor, we consider the distribution of the Lorentz factor, assuming that the lognormal function describes the distribution function: where F 0 is the normalization factor, Γ 0 is the mean Lorentz factor, and σ Γ is the dispersion of the distribution. We fix the other parameters to those of EE-mod and EE-opt shown in Table 1. We consider the merger events at 300 Mpc, which is the detection horizon of the design sensitivity for advanced LIGO. Using the local event rate obtained by the SGRB observations, 4 − 10 Gpc −3 yr −1 [66,67], we expect 2-5 extended emissions within 10 years if half of SGRBs are followed by extended emissions Table 1. Used parameters for each component of SGRBs. Γ is the Lorentz factor of the jet, L * γ,iso and E * γ,iso are the isotropic equivalent luminosity and energy fluence in the observed energy band, r diss is the dissipation radius, E γ,pk is the observed break energy of the photon spectrum, and the last column shows the energy band of the SGRB observations. This  Table 2. Resulting physical quantities for each component of SGRBs. B is the magnetic field in the dissipation region, L γ,iso and E γ,iso are the total isotropic equivalent luminosity and energy fluence, E p,M is the maximum energy for protons, E ν,µ is the critical neutrino energy above which the muon cooling is effective, E ν,π is the break energy of neutrino spectrum due to pion cooling. This  (cf. [68]). Assuming that all the merger events at 300 Mpc are detected by GWs, we tabulate the resulting probabilities of neutrino detection coincident with GWs within 10 years of operation in Table 3. For optimistic cases, we can highly expect the coincident neutrino detection with GWs even with the current facilities. For the moderate case, the neutrino detection is probable with IceCube-Gen2, while it is challenging with IceCube.

Trans-ejecta neutrinos
The optical/UV/IR counterparts of GW170817 confirmed that neutron star mergers produce massive ejecta. The ejecta is produced immediately after the merger, while the launch of the relativistic jets can be delayed with a time lag of t lag 1 sec. Thus, the jets interact with the ejecta, forming a cocoon surrounding the jets [69][70][71][72]. If the luminosity of the jets is low or the duration of the jet launch is short, the jets are choked inside the ejecta. This choked jet system is expected for a wide parameter range [73]. In this system, the photons are completely absorbed by the ejecta, while the neutrinos can penetrate the ejecta. Using such trans-ejecta neutrinos, we can discuss the physical conditions of the choked jet system without the electromagnetic signals.
In the choked jet systems, there are two possible dissipation site: the internal shocks and collimation shocks [74]. During the jet interacting with the ejecta, the cocoon collimates the jets by pushing the jets inward, which forms the collimation shock. If the central engine creates the strong velocity fluctuation of the jets, the internal shocks can be formed below the collimation shock [16]. We draw a schematic picture of this system in Figure 2. The jet head, the interaction region of the jet and ejecta, also has strong forward and reverse shocks. However, we cannot expect particle acceleration there, because the density of the shock upstream is too high (see the next paragraph).
The typical size of the choked jet system is ∼ 10 10 cm due to the upper limit of the time lag, t lag 1 sec. This is much smaller than the typical emission region of SGRBs (see Table 1). Hence, the dissipation region is very dense. If the shock upstream is too dense, the shock is mediated by radiation, causing a gradual velocity change [75]. This prevents the particles from being accelerated. The necessary condition for particle acceleration at the shock is given by [76] τ u = n u σ T l u 1, where τ u is the optical depth for the upstream, n u is the density at the upstream, σ T is the Thomson cross sectin, and l u is the length of the upstream fluid. To satisfy this condition, the Lorentz factor of the jets should be Γ j 200 for the internal shocks and Γ j 500 for the collimation shocks with a typical parameter set of the system. We set the ejecta mass M ej = 0.01 M , the ejecta velocity V ej = 0.33c, the lag time t lag = 1 sec, the jet opening angle θ j = 0.3 rad, the duration of the jet launch t dur = 2 sec, and the kinetic luminosity L k,iso = 10 51 erg s −1 .
The neutrinos are produced at the shock downstream. The downstream of the collimation shock has a Lorentz factor of a few, leading to the very high baryon density and strong magnetic field there. This causes the strong cool-ing of pions, both by synchrotron and pion-proton inelastic collisions. The cutoff energy in the neutrino spectrum is typically less than 0.3 TeV, which is a too low energy to be detected by IceCube. On the other hand, the downstream of the internal shock has a relatively high Lorentz factor, Γ ∼ 300, so they can emit high-energy neutrinos of 100 TeV. Thus, we focus on the neutrinos from the internal shocks when we discuss the neutrino detectability. Figure 3 shows the neutrino fluences from the internal shocks of the choked jet systems for the optimistic (model A) and moderate (model B) cases, whose parameters are tabulated in Table 4. The target photons are provided from the downstream of the collimation shock, where the photon distribution is the Planck function owing to the high optical depth. The photon density at the downstream of the internal shocks is so high that this system can be calorimetric. This leads to a flat neutrino spectrum below the cutoff energy caused by the pion cooling around 100 TeV, although the muon cooling causes a slightly softer neutrino spectrum than that for protons. The optimistic and moderate cases differ in the Lorentz factor and the jet luminosity, which mainly change the cutoff energy and normalization of the fluence, respectively. We also plot the neutrino spectrum for a case with the successful jet case (model C). The trans-ejecta neutrinos can be produced when the jet head is inside the ejecta even for the successful jet case. Such trans-ejecta neutrinos can be detected as precursor neutrinos of the SGRBs.
Using the fluences shown in Figure 3, we estimate the detection probability of neutrinos coincident with the GWs. The upper two parts of Table 5 give the detection probability for a single merger event at a given distance. If the merger happens at 40 Mpc, the neutrinos are detectable with IceCube for the optimistic case, and IceCube-Gen2 is likely to detect the neutrinos even for the moderate case. On the other hand, if the merger happens at 300 Mpc, the coincident detection is challenging with IceCube-Gen2 even for the optimistic case. The lower part provides the neutrino detection rate per year. Here, we use a neutron star merger rate obtained by the LIGO/Virgo collaborations, ∼ 1.5 × 10 3 Gpc −3 yr −1 , and consider the uniformly distributed population in the local universe. For the optimistic case, we can highly expect the coincident detection of GWs and neutrinos by IceCube with a few years of operation. Even with the moderate case, the coincident detection is likely for several years of operation with IceCube-Gen2.

Implications from GW170817
Although our models predict detectable neutrino fluences for some optimistic cases, such high-energy neutrinos are not detected from GW170817. However, our models are not constrained by this result. First, GW170817 turned out to be the off-axis SGRBs. For the off-axis events, the neutrino fluence decreases with (θ v /θ j ) 2 for θ v < 2θ j and (θ v /θ j ) 3 for θ v > 2θ j [24,77], where θ j and θ v are the jet opening angle and the viewing angle. From the VLBI observation, these angles are estimated to be θ v ∼ 20 degree and θ j 5 degree, leading to L v /L on 1/32, where L v and   L on are the isotropic equivalent luminosities from off-and on-axis observers, respectively. Thus, the neutrino fluence is too low to be detected by the current facilities. Second, GW170817 occurred at the southern hemisphere, where the sensitivity of IceCube is lower at E ν 1 PeV [57]. This makes it difficult to detect neutrinos of 100 TeV emitted from the choked jets and extended emissions for the optimistic model. KM3NeT will be useful to detect such neutrinos in the southern hemisphere. Finally, GW170817 may not have the efficient neutrino production sites. The extended emission is not detected from the event, and the relativistic jet is observed from this event. Hence, it is possible to have neither the extended emissions nor the choked jets. Future on-axis events will provide detection of neutrinos or put strong constraints on the physical parameters of the choked jets and late-time activities. Table 4. The used parameters for the internal shock models. L k,iso is the isotropic equivalent kinetic luminosity of the jets, Γ j is the Lorentz factor of the downstream of the internal shocks, t dur or t bo is the duration of the jet launch (for models A and B) or the breakout time from the ejecta (for model C), ξ acc is the baryon loading factor, and Γ rel−is is the Lorentz factor of the internal shock. This  Table 5. Detection probability of neutrinos coincident with GWs by IceCube and IceCube-Gen2 (Gen2). "up+hor" means the up-going and horizontal events. "down" means the down-going events. Since the effective area of the down-going events with IceCube-Gen2 is very uncertain, we avoid discussing it. This

Super-knee cosmic rays from neutron star merger remnants
Neutron star mergers produce fast and massive ejecta. This ejecta interacts with the ambient medium, forming a forward shock that accelerates the cosmic rays. The UV/optical/IR counterparts of GW170817 provide the mass (0.01 M − 0.05 M ) and velocity (0.1c − 0.3c) of the ejecta of BNS mergers. Also, the GW observation gives rough estimate of the event rate, ∼ 1.5 × 10 3 Gpc −3 yr −1 . This enables us to estimate the cosmic-ray production at the neutron star merger remnants (NSMRs). Since the ejecta of NSMRs are faster than that of supernova remnants (SNRs), NSMRs can produce higher energy cosmic rays than SNRs. In this section, we discuss whether the NSMRs can account for the cosmic-rays above the knee. See [21] for details of this section. At a NSMR, the balance between acceleration and age gives the maximum energy of cosmic-rays at a given time, which is expressed as where Z i is the charge of the particle species i, e is the elementary charge, B is the magnetic field, R ej is the ejecta radius, and V ej is the ejecta velocity. The time evolution of the velocity and radius of the ejecta is ballistic before the deceleration time, and given by the Sedov-Taylor solution after that. Then, the maximum energy of cosmic rays through an entire life of a NSMR is obtained at the deceleration time, which is estimated to be where we use a typical parameter set with the ejecta mass M ej 0.03 M , the initial ejecta velocity V ini 0.25c, the ambient density n amb 0.1 cm −3 , and the magnetic field B 0.4 mG at the deceleration time.
From the GW observation, the local BNS merger rate is estimated to be 1.5 × 10 −6 Mpc −3 yr −1 . Using the density of the Milky-way-size galaxies, ∼ 0.01 Mpc −3 yr −1 , the occurrence time of neutron star mergers in our Galaxy is estimated to be T mer 1.5 × 10 −4 yr −1 . On the other hand, the escape time of the Galactic cosmic rays are estimated to be T esc ∼ 20 − 400 Myr for particles of the rigidity R = 10 GV [78]. The rigidity dependence of the escape time is often assumed to be T esc ∝ R −δ . For δ 0.4, the escape time for the cosmic rays of R < 10 8 GV is always longer than the occurrence time. Hence, we can use the steady state assumption. Using the grammage obtained from the recent experiments [79] and one-zone approximation for the interstellar cosmic-ray density, the cosmic-ray intensity on Earth is estimated to be where X esc is the grammage, EQ E,inj is the injection term, M gas is the total gas mass in the cosmic-ray halo of the Milky-way Galaxy, δ 1/3 is the energy dependence of X esc , X esc ∝ E −δ , for the range of our interest. To calculate the injection term, we take into account the escape process from the NSMRs according to [80], where only the cosmic rays near the maximum energy can escape from the NSMR. Considering the time evolution of the ejecta radius and velocity discussed in the previous paragraph, the resulting escape spectrum has the same power-law index as that for the injection spectrum. The normalization of the injection term is given by E cr ≈ cr M ej V 2 ini /2, where we set the cosmic-ray production efficiency cr = 0.25. The composition ratio of the injection term is given by the model proposed by [81], in which the cosmic-ray injection efficiency is proportional to (A i /Z i ) 2 . Applying this model to the ambient medium of the solar abundance ratio, the abundance ratio of each element at the source is given by ( 17, 0.52, 0.024, 0.099, 0.027, 0.028, 0.14).
The resulting cosmic-ray intensity is shown in Figure 4.
The abundance ratio on Earth is written as ( f p , f He , f C , f O , f Ne , f Si , f Fe ) (0.10, 0.41, 0.028, 0.13, 0.037, 0.043, 0.26), which is different from that at the NSMRs due to propagation effect. The GeV-PeV and extragalactic components (see [82] for the extragalactic one) are also plotted in the figure. Our model can reproduce the observed flux well. The cosmic-rays from NSMRs are dominant for 10 7 GeV E p 10 9 GeV. Our model is also consistent with the hardening around 10 7 GeV recently reported by Icetop and Telescope Array Low-energy Extension (TALE) [83,84]. In addition, the light-component intensity for 10 7 GeV E p 3 × 10 8 GeV matches that reported by KASCADE-Grande [85]. The predicted composition ratio is also consistent with that obtained by the experiments, although the uncertainty is large.

Summary
We have briefly reviewed the multi-messenger event GW170817, future prospects for neutrino detections coincident with GWs, and super-knee cosmic rays from the remnants of neutron star mergers. The multi-messenger campaign of GW170817 confirmed that the BNS mergers are the progenitor of SGRBs caused by relativistic jets and macronovae powered by radioactivity of r-process elements.
The observed SGRBs are accompanied by the latetime emissions, which can emit neutrinos more efficiently than the prompt emissions. Also, if the jet duration is short or the jet luminosity is low, the ejecta of macronovae can choke the jets, causing the failed SGRBs. Such choked jets are also a strong neutrino source. Hence, the GWs can be accompanied by the high-energy neutrinos. We estimated detectability of neutrinos from these systems within the GW horizon of the design sensitivity of advanced LIGO. For both cases, the neutrino detection coincident with the GW is probable with IceCube-Gen2 with 10-year operation even with the moderate parameter set. IceCube is also possible to detect the neutrinos with 10 years of operation for the optimistic cases.
The neutron star mergers produce massive outflows, whose velocity is higher than that of the supernovae. Thus, the NSMRs can produce higher energy cosmic rays. We estimate the cosmic-ray spectrum from Galactic NSMRs, and find that NSMRs can be the dominant source of the cosmic rays from 10 PeV to 1 EeV. Our model can naturally explain the hardening feature around 10 PeV, and the   [82]). The gray band represents the experimental data [84,86]. The cyan region shows the intensity of the light composition [85]. Lower panel: the mean atomic number of the cosmic rays on Earth. The solid line shows our model prediction. This is consistent with the experimental results (yellow [87] and cyan regions [88]), although uncertainty is large. This figure is reproduced from [21]. spectrum of the light elements reported by KASCADE-Grande. S.S.K. thanks Kohta Murase and Peter Meszaros for fruitful discussion and continuous encouragement. This work is supported by JSPS oversea research fellowship and IGC fellowship.