{'fraction_non_alphanumeric': 0.03845862095122506, 'fraction_numerical': 0.027611317606007737, 'mean_word_length': 4.009118541033435, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Constraining UHECR source models by the TA SD energy spectrum

We have published the energy spectrum of ultra-high energy cosmic rays (UHECRs) for E>10^18.2 eV from our first 5-year observation with the surface detectors of Telescope Array. We found two breaking points in the power-law spectrum, a flattening known as the"ankle"or"dip"at E = 10^18.66 eV, and a steepening at E = 10^19.72 eV. The spectral features must be related to distribution of cosmic ray sources, injection spectrum, and energy loss processes during propagation in inter-galactic space. We constructed a model of source distribution of UHECRs and energy loss processes in the CMB photon field to calculate energy spectrum to be observed at the Earth. The model includes injection energy spectra at the sources, energy loss processes, and evolution of cosmic ray fluxes with red shift. We consider two cases that the UHECR sources are distributed uniformly, or same as galaxies in the local universe. We found that the spectrum can be well fitted with the model. We discuss to give a constraint to the model parameters, e.g. the injection power indices and cosmological evolution from the observed data.


Introduction
The Telescope Array (TA) [1] is a cosmic-ray observatory of the largest area in the northern hemisphere. TA experiment has 507 surface particle detectors (SDs) with 1.2 km spacing covering about 700 km 2 ground area. SDs are surrounded by 3 fluorescence detector (FD) stations which consist of 12 or 14 FDs. The duty cycle of the SD array is greater than 95 % throughout 5-year observation, whereas the FD duty cycle is about 10 % because observation is possible only in moon-less clear night. The TA SD energy spectrum given in Ref. [2] is obtained by the 4-year observation. This energy spectrum shows the ankle at 10 18.66 eV and the cutoff at 10 19.72 eV in 5.5 σ . Both of these features are firstly observed by High Resolution Fly's Eye (HiRes) [3] and also confirmed by Pierre Auger Observatory (PAO) [4]. We analyzed this energy spectrum with additional 1-year data.
Cosmic ray mass composition is quite important in interpreting cosmic ray spectral features. The result of the measurement of Xmax by TA FDs in stereo mode is consistent with proton models above 10 18.2 eV [5]. In other experiments, HiRes reported the consistency with proton models [6] but PAO reported the significant discrepancy between data and considered proton models [7].
Several models have been proposed to explain the spectral features like the flattening around 10 18.7 eV or the steepening at 10 19.72 . The dip model proposed by Berezinsky et al. [8] explains the ankle by the energy dependence of the e + e − pair production process when protons from extragalactic sources interact with cosmic microwave background (CMB) photons. The cutoff can be also explained by the interaction of protons with cosmic microwave background (CMB) photons. This process is known as GZK mechanism [9]. Therefore both of the ankle and the cutoff feature are given only by protons from extragalactic sources with these processes.
In this paper, we compared the observed TA SD energy spectrum with this source model in detail. In the next sec-tion, the model calculations and the fitting procedure are explained, and the fitted results are shown in section 3. The discussion and conclusions are described in section 4 and 5.

Methodology
We firstly assumed uniform source distribution as cosmic ray sources. In the latter part of this section, we also show the methodology when we assume local large scale structure (LSS) source distribution using galaxy distribution. This is a reasonable assumption if the astrophysical sources are UHECR sources. We also assumed that all the cosmic ray sources have the same power-law spectrum αE −γ , where α and γ g are the free parameters. The evolutions of source number densities are also applied uniformly as (1 + z) m , where m is the free parameter.
The energy loss processes of UHECR protons in the CM-B photon field is treated with the formulation in Ref. [10] for energy spectrum below 2 · 10 19 eV. In this formulation the energy loss time of protons is given by a polynomial of log E. The redshift energy loss is also taken into account as 1/E · dE/dz = 1/(1 + z). We employ CRPropa v2.0 [11] for energy spectrum with energies greater than 2 · 10 19 eV to simulate energy losses to consider the distribution of energy losses in each interaction.

Uniform source distribution
Using an uniform source distribution and the energy loss processes, the expected energy spectrum is calculated with the following procedure that is similar to the section 2.3 of Ref. [12] for energy spectrum below 2 · 10 19 eV. Firstly we start from the energy spectrum of a single source. We assume that each source isotropically emits cosmic rays.
where J 0 (E 0 ) is the source power above E 0 at the source and E 0 is the energy threshold at the source. The proper distance D(z) is given by where c is the speed of light and H(z) is the Hubble parameter at the redshift z.
The total integral flux F exp (E) within the redshift z max can be described by where n src (z = 0) is the UHECR source density at z = 0 and m is the assumed source number evolution. From the assumption, J 0 (E) can be given by where E max is the acceleration limit of cosmic rays at the source. E 0 can be calculated by the numerical integration of dE/dz from the energy E with the energy loss processes of protons which are explained above. F exp (E) is actually calculated using these calculation processes when z max = 5, , Ω Λ = 0.73 and E max = 10 21 eV. The differential energy flux f exp (E) = dF exp (E)/dE is also numerically calculated with given F exp (E).
Monte Carlo simulations of the trajectories of protons are carried out using CRPropa v2.0 to calculate energy spectrum above 2 · 10 19 eV. The uniform source distribution is set from 0 Mpc to 1400 Mpc, which is well beyond the GZK horizon. More than about 10 5 trajectories are obtained and we calculate differential energy flux from them. We connect this result to the flux from lower energies at 2 · 10 19 eV to obtain the final energy spectrum. We calculated expected number of events N exp in each energy bin from this calculated differential energy flux. In this calculation, we used the same acceptance and energy resolutions as in Ref. [2]. We compared N exp with the data.
The model fitting of the data is carried out for energies above 10 18.2 eV. The best-fit four model parameters α, γ g , m and the energy scale are determined by the fitting. The best-fit α and the energy scale are also determined for the fixed γ g and m in order to determine the confidence region of γ g and m.
The fitting procedures are the following. Here, we employed a binned likelihood analysis method considering Poisson probability distribution of number of events in each energy bin. We define the likelihood function where P(n; µ) = µ n exp(−µ)/n!, N i is the observed number of events in the i-th energy bin, E i is the i-th energy and N is the number of bins. In this case, N = 21. The likelihood ratio can be given by λ = ∏ i P(N i ; N exp (E i ))/ ∏ i P(N i ; N i ), and the quantity −2 ln λ follows χ 2 distribution. The degree of freedom of the χ 2 corresponds to N minus number of fitting parameters. The best-fit parameters are determined by minimizing this χ 2 . Here, we define p-values as the probability of obtaining larger χ 2 minimum than the observed one. We calculated p-values for the fixed parameters to draw the confidence region of the parameters.

LSS source distribution
As a model of UHECR sources along the large scale structure, we used about 110,000 galaxy samples of 2 Mass Extended Source Catalogue (XSCz) [13]. We selected the galaxies in the XSCz catalog with distances smaller than 250 Mpc and apparent magnitudes brighter than 12.5 in the Ks band (2.2 µm). We assumed an uniform matter distribution beyond 250 Mpc.
The galaxies within 5 Mpc were not included because they do not represent a proper statistical sample of LSS. We introduced a weighting factor for each selected galaxy to take into account faint galaxies below the limit 12.5 with the distribution of absolute magnitudes of galaxies as proposed in Ref. [14].
We calculated the distance dependence of ∑ i w i A i from the galaxies, where w i is the weight of each galaxy and A i is the relative TA SD exposure in the direction of the galaxy. Then we simulated the propagation from the source distribution ∑ i w i A i in the same way as Section.2.1. Fig.1 shows the calculated distance dependence of ∑ i w i A i from the galaxies.

Uniform source distribution
We show the expected energy spectrum with source parameters γ g = 2.36 and m = 4.5 determined by fitting the TA SD spectrum in Fig.2. The expected energy spectra are calculated when γ g = 2.00, 2.01, · · · 2.69 and m =   0.00, 0.01, · · · 9.99, and compared with the data in the same way.
The p-values calculated from the χ 2 are shown in Fig.3. The two parameters γ g and m are correlated, however we can give a constraint band to the parameters with the current statistics of TA as shown in Fig.3. γ g = 2.36 +0.08 −0.04 and m = 4.5 +0.6 −1.1 are obtained as a conclusion. Note that the number evolution (1 + z) m are known for several classes of astrophysical objects, e.g. m = 3 for Q-SOs and z < 1.3 [15], m = 4.8 for GRBs and z < 1 [16] and m = 5 for AGNs and z < 1.7 [17].

LSS source distribution
The blue line in Fig.2 is calculated in the same procedure as the pink line except for the source distribution. In LSS source distribution, there are some regions locally overdense because of local galaxy clusters like Virgo cluster as you can see in Fig.1. Therefore the flux at the high-energy end expected from the LSS source distribution is larger than that from the uniform sources as seen in Fig.2. The χ 2 of the LSS source distribution in Fig.2 is not much different from that of the uniform sources as a result because of large error bars of the data at the high energy end. The p-values calculated from the χ 2 are shown in Fig.4. γ g = 2.39 +0.08 −0.08 and m = 4.4 +0.9 −1.3 are obtained as a conclusion.

Discussion
The systematic uncertainty of TA SD energies is 22% [2]. This only changes the energy scale and does not affect the shape of the energy spectrum. We included the shift of the energy scale in the fitting as a free parameter as described above. Here, we define the shift of the energy scale as ∆ log E = log E − log E, where E is the shifted energy scale. The energy scale ∆ log E = −0.01 +0.02 −0.02 is obtained by the fitting if we assume uniform source distribution. The energy scale ∆ log E = 0.02 +0.04 −0.05 is determined by the fitting if we assume LSS source distribution. So the observed energy scale is consistent with these models. The systematic uncertainty of the energy scale is larger than the statistical fluctuation of that. So this can have a big influence on the conclusions in this paper.
In case of strong evolution models, the larger component of energy flux is produced from the distant sources. This situation can be strongly restricted from the observation of the diffuse gamma ray background. We calculated secondary photons in the propagation of cosmic ray protons with CRPropa v2.0. The sources with redshift below 0.7 are taken into account. This redshift corresponds to the event horizon of 10 18.2 eV cosmic ray protons. If m > 7, the energy spectrum of secondary photons exceeds the data points of the diffuse gamma ray background that is observed by Fermi Large Area Telescope [18].
If the source density has the dependence of the maximum acceleration energy E max or luminosity dependence, our conclusions can be different [19] [20]. If the source density n src (E max ) ∝ E −δ max , δ = 1.7 and γ g = 2.0, the energy spectrum becomes similar to the source with γ g = 2.7 and δ = 0 when we do not consider the evolution. If the clustering in the arrival directions is observed in the future, this dependence can be discussed experimentally.

Conclusions
We calculated the expectation of the cosmic ray energy spectrum of uniform source distribution and LSS source distribution with pure proton model. We fitted the TA S-D spectrum with this expectation. χ 2 /d.o. f . = 20.2/17 is obtained if we assume uniform source distribution. χ 2 /d.o. f . = 18.5/17 is determined if we assume LSS source distribution. Both source distributions can explain the data.
We investigated a compatibility of the TA SD spectrum and the expectation with source spectral index γ g and the evolution m as free parameters. We found that γ g = 2.36 +0.08 −0.04 and m = 4.5 +0.6 −1.1 are obtained in case of uniform source distribution. γ g = 2.39 +0.08 −0.08 and m = 4.4 +0.9 −1.3 are determined if we assume LSS source distribution. No evolution model is ruled out.
This result is determined mainly from the dip feature of the energy spectrum around 10 18.7 eV and TA SD needs more statistics to discuss the compatibility at the high energy end which is above 10 19.7 eV. We need more statistics in higher energies to rule out the uniform source distribution model.