Measurement of the cosmic ray proton spectrum from 40 GeV to 100 TeV with the DAMPE satellite.

The precise measurement of the spectrum of protons, the most abundant component of the cosmic radiation, is necessary to understand the source and acceleration of cosmic rays in the Milky Way. This work reports the measurement of the cosmic ray proton fluxes with kinetic energies from 40 GeV to 100 TeV, with 2 1/2 years of data recorded by the DArk Matter Particle Explorer (DAMPE). This is the first time that an experiment directly measures the cosmic ray protons up to ~100 TeV with high statistics. The measured spectrum confirms the spectral hardening at ~300 GeV found by previous experiments and reveals a softening at ~13.6 TeV, with the spectral index changing from ~2.60 to ~2.85. Our result suggests the existence of a new spectral feature of cosmic rays at energies lower than the so-called knee and sheds new light on the origin of Galactic cosmic rays.

INTRODUCTION

It is widely believed that the remnants of explosive stars in the Milky Way may play a substantial role in producing energetic cosmic ray (CR) particles (, ). The energy spectra of CRs are expected to be single power laws (PLs) until the energies exceed the maximum acceleration limits of the sources, based on the conventional Fermi acceleration models (). The diffusive transportation of CRs in the interstellar turbulent magnetic field results in a softening of the accelerated spectrum, by a PL form (for rigidities above a few tens of gigavolts) according to the boron-to-carbon ratio (). This general picture of CR production and propagation has been supported by measurements of CR energy spectra and composition ratios, as well as diffuse gamma rays ().

However, such a simple picture has been challenged by some recent high-precision measurements. Remarkable spectral hardenings of the energy spectra of CRs were revealed in the Advanced Thin Ionization Calorimeter (ATIC) (, ), Cosmic Ray Energetics And Mass (CREAM) (, ), Payload for Antimatter-Matter Exploration and Light-nuclei Astrophysics (PAMELA) (), and Alpha Magnetic Spectrometer (AMS)–02 () observations. The spectral hardenings suggest extensions of the traditional CR source injection, acceleration, and/or propagation processes, e.g., (–). The spectral behaviors of CRs at higher energies (more than tera–electron volts) are essential to understand the nature of the spectral hardenings, as well as the origin and propagation of CRs. Moreover, the extension of the spectra to peta–electron volt energies according to the PL indices measured at sub-tera–electron volt energies seems to be in conflict with the all-particle spectrum of CRs (). Such a puzzle may be solved if a significant spectral softening presents well below the so-called knee at several peta–electron volts. The precise measurements of the energy spectra of CRs above tera–electron volts are thus motivated by the test of potential new spectral features. The recent CREAM and NUCLEON data show hints that the energy spectra of CR nuclei may become softer above rigidities of 10 to 20 TV (, ). However, the result of the proton plus helium spectrum from the air shower experiment ARGO-YBJ shows a single PL form for energies between 3 and 300 TeV (). The CREAM result is mainly limited by its low statistics. The NUCLEON data, with again relatively low statistics at energies above tens of tera–electron volts, also suffer from sizeable systematic uncertainties that need to be properly included. The indirect measurements, on the other hand, suffer from poor composition resolution. Although the magnetic spectrometers can measure CRs very accurately, they are unable to reach energies well beyond tera–electron volts in the foreseeable future. Therefore, the calorimeter-based direct measurement experiments, with high statistics up to ∼100 TeV and well-controlled systematic uncertainties, are most suitable to solve the above problems.

RESULTS

In this work, we present the measurement of the proton spectrum with the DArk Matter Particle Explorer (DAMPE; also known as “Wukong” in China). DAMPE is a calorimetric-type, satellite-borne detector for observations of high-energy electrons, gamma rays, and CRs (, ). From top to bottom, the instrument consists of a plastic scintillator strip detector (PSD) (), a silicon-tungsten tracker converter (STK) (), a bismuth germanate (BGO) imaging calorimeter (), and a neutron detector (NUD) (). The PSD measures the charge of incident particles and serves as an anticoincidence detector for gamma rays. The STK reconstructs the trajectory and also measures the charge of the particles. The BGO calorimeter measures the energy and trajectory of incident particles and provides effective electron/hadron discrimination based on the shower images. The NUD provides additional electron/hadron discrimination. These four subdetectors enable good measurements of the charge (∣Z∣) with a resolution (Gaussian SD) of about 0.06e and 0.04e for the PSD () and the STK, respectively, the arrival direction with an angular resolution of better than 0.5° above 5 GeV, the energy with a resolution of higher than 1.5% for >10-GeV electrons/photons () and about 25 to 35% for protons up to 10 TeV (), and the identification of incoming particles with a proton rejection capability of about 3 × 104 when keeping 90% of electrons (). The DAMPE detector was launched into a 500-km Sun-synchronous orbit on 17 December 2015. The on-orbit calibration results demonstrate that DAMPE operates stably in space ().

The data used in this work cover the first 30 months of operation of DAMPE, from 1 January 2016 to 30 June 2018. The fraction of live time is about 75.73% after excluding the time when the satellite passes the South Atlantic Anomaly region, the instrument dead time, the time for on-orbit calibration, and the period between 9 September 2017 and 13 September 2017 during which a big solar flare occurred and may have affected the baseline of the detector. We select protons using the charge measured by the PSD (see Materials and Methods for details about the event selection). Figure 1 illustrates the reconstructed PSD charge spectra for low-Z nuclei for deposited energies of 447 to 562 GeV (left), 4.47 to 5.62 TeV (middle), and 20 to 63 TeV (right), together with the Monte Carlo (MC) simulations of protons and helium nuclei with GEANT v4.10.03 (). Note that small corrections from the reconstructed charge to the true particle charge () have not been applied in this work. The proton and helium peaks are separated in this plot. The contamination of the proton sample due to helium nuclei is found to be less than 1% for deposited energies below 10 TeV and about 5% around 50 TeV, as an effect of the energy-dependent charge selection (see Materials and Methods). Given the excellent electron-proton discrimination capability of DAMPE (), the contamination due to residual electrons is estimated to be about 0.05% in the whole energy range analyzed in this work.

Fig. 1

The combined signal spectra of PSD for protons and helium nuclei.

The left panel is for BGO deposited energies between 447 and 562 GeV, the middle panel is for BGO deposited energies of 4.47 to 5.62 TeV, and the right panel is for BGO deposited energies between 20 and 63 TeV. The on-orbit data (black) are shown, together with the best-fit templates of simulations of protons (blue), helium nuclei (green), and their sum (red). The vertical dashed lines show the cuts to select proton candidates in this deposited energy range.

The proton spectrum in the energy range from 40 GeV to 100 TeV is shown in Fig. 2 and tabulated in Table 1. Error bars represent the 1σ statistical uncertainties of the DAMPE measurements, and the shaded bands show the systematic uncertainties associated with the analysis procedure (inner band) and the total systematic uncertainties including those from the hadronic models (outer band). Previous measurements by space detectors PAMELA () and AMS-02 (), and balloon-borne detectors ATIC-2 (), CREAM (), and NUCLEON () are overlaid for comparison. The DAMPE spectrum is consistent with those of PAMELA and AMS-02. At higher energies, our results are also consistent with that of CREAM, ATIC-2, and NUCLEON when the systematic uncertainties are taken into account.

Fig. 2

Proton spectrum from 40 GeV to 100 TeV measured with DAMPE (red filled circles).

The red error bars show the statistical uncertainties, the inner shaded band shows the estimated systematic uncertainties due to the analysis procedure, and the outer band shows the total systematic uncertainties including also those from the hadronic models. The other direct measurements by PAMELA () (green stars), AMS-02 () (blue squares), ATIC-2 () (cyan diamonds), CREAM I + III () (magenta circles), and NUCLEON-KLEM () are shown for comparison. For the PAMELA data, a −3.2% correction of the absolute fluxes has been included (, ). The error bars of PAMELA and AMS-02 data include both statistical and systematic uncertainties added in quadrature. For ATIC-2, CREAM, and NUCLEON data, only statistical uncertainties are shown.

Table 1

Fluxes of CR protons measured with DAMPE, together with 1σ statistical and systematic uncertainties.

The systematic uncertainties include those associated with the analysis procedure σana (e.g., the event selection, the background subtraction, and the spectral deconvolution) and the energy responses due to different hadronic models σhad.

〈E〉 (GeV)	Emin (GeV)	Emax (GeV)	F ± σstat ± σana ± σhad (GeV−1 m−2 s−1 sr−1)
49.8	39.8	63.1	(2.97 ± 0.00 ± 0.14 ± 0.20) × 10−1
78.9	63.1	100.0	(8.43 ± 0.00 ± 0.40 ± 0.56) × 10−2
125.1	100.0	158.5	(2.38 ± 0.00 ± 0.11 ± 0.16) × 10−2
198.3	158.5	251.2	(6.64 ± 0.00 ± 0.31 ± 0.44) × 10−3
314.3	251.2	398.1	(1.89 ± 0.00 ± 0.09 ± 0.12) × 10−3
498.1	398.1	631.0	(5.39 ± 0.01 ± 0.25 ± 0.36) × 10−4
789.5	631.0	1000	(1.60 ± 0.00 ± 0.07 ± 0.11) × 10−4
1251	1000	1585	(4.81 ± 0.01 ± 0.23 ± 0.33) × 10−5
1983	1585	2512	(1.45 ± 0.01 ± 0.07 ± 0.13) × 10−5
3143	2512	3981	(4.45 ± 0.02 ± 0.21 ± 0.44) × 10−6
4981	3981	6310	(1.36 ± 0.01 ± 0.06 ± 0.13) × 10−6
7895	6310	10,000	(4.06 ± 0.04 ± 0.19 ± 0.40) × 10−7
12,512	10,000	15,849	(1.20 ± 0.02 ± 0.06 ± 0.12) × 10−7
19,830	15,849	25,119	(3.35 ± 0.07 ± 0.17 ± 0.33) × 10−8
31,429	25,119	39,811	(9.03 ± 0.26 ± 0.48 ± 0.89) × 10−9
49,812	39,811	63,096	(2.47 ± 0.11 ± 0.15 ± 0.24) × 10−9
78,946	63,096	100,000	(6.50 ± 0.40 ± 0.50 ± 0.64) × 10−10

DISCUSSION

The features of the proton spectrum measured by DAMPE in the energy range from 40 GeV to 100 TeV give fundamental information about the origin and propagation of Galactic CRs. A spectral hardening at a few hundred giga–electron volt energies is shown in our data, in agreement with that of PAMELA () and AMS-02 (). As discussed in several papers [() and references therein], the hardening can be due to either details of the acceleration mechanism, effects in the propagation in the Milky Way, or the contribution of a new population of CRs (e.g., a nearby source). Furthermore, the DAMPE measurement gives strong evidence of a softening at about 10 TeV. It is worth reminding that a maximum of the large-scale anisotropy has been observed just at that energy [see, e.g., (, )]. We fit the spectrum with energies between 1 and 100 TeV with a single PL model and a smoothly broken PL (SBPL) model, respectively, and find that the SBPL model is favored at the 4.7σ confidence level compared with the single PL one (see Materials and Methods for details of the fit). For the SBPL model fit, the break energy is TeV, the spectral index below the break energy is 2.60 ± 0.01, and the change of the spectral index above the break energy is −0.25 ± 0.07. Results recently published by CREAM () and NUCLEON () experiments also indicate a spectral softening at rigidities of ∼10 TV. However, these results are limited by low statistics or the lack of careful studies of systematic uncertainties.

The spectral hardening and softening are not compatible with the paradigm of a unique PL spectrum up to the all-particle knee at peta–electron volt energies, thus implying a deep revision of CR modeling in the Galaxy. For instance, the 10-TeV softening might be due to the exhaustion of the contribution of a given CR population. Either a local source on top of a PL background () or various types of sources () can be compatible with this scenario. Note that the spectral softening should not correspond to the knee of protons; otherwise, the expected all-particle spectrum would undershoot the observational data, either for mass-dependent or for charge-dependent knees of various species (). Therefore, the current DAMPE measurement of the proton spectrum, together with other measurements from the space and ground-based experiments, puts a severe constraint on the models of Galactic CRs.

MATERIALS AND METHODS

MC simulations

Extensive MC simulations were carried out to estimate the selection efficiencies, background contaminations, and the energy response matrix of hadronic cascades in the detector (particularly the BGO calorimeter). The GEANT v4.10.03 () was adopted for these simulations. There are two typical hadronic interaction models in the GEANT simulation tool, the QGSP_FTFP_BERT and FTFP_BERT models. Comparisons of the shower development (longitude and transverse distributions) between simulations and the beam test and on-orbit data show that the FTFP_BERT model matches better with the data. Therefore, we adopted the FTFP_BERT model as the benchmark of the MC simulations for incident energies less than 100 TeV. For higher energies, we used the FLUKA tool, which links the DPMJET model for the simulation (). The FLUKA-based results were also used as a cross-check and an estimate of the systematic uncertainties of the hadronic interaction models through comparison with the GEANT results.

An isotropic flux with the E−1.0 spectrum was generated for the detector simulation. We simulated protons, helium nuclei, and electrons in this analysis. After the charge selection, heavier nuclei are negligible for the proton analysis. In the analysis, the spectra were reweighted to E−2.7 for protons and helium nuclei and to E−3.15 for electrons. The final proton spectrum we measured is not exactly the same as that of E−2.7. However, changing the reweight spectrum with indices from 2.5 to 3.1 has little influence on the results.

Proton event selection

The events with energy deposit in the BGO calorimeter larger than 20 GeV were selected in this analysis to suppress the effect of the geomagnetic rigidity cutoff. The detailed event selection method is described as follows.

Preselection

DAMPE has four different triggers implemented on orbit: the unbiased trigger, the minimum ionizing particle (MIP) trigger, the low-energy (LE) trigger, and the high-energy (HE) trigger (). The events were required to satisfy the HE trigger condition to guarantee that the shower development starts before or at the top of the calorimeter. We further required that there are one or more hits in each sublayer of the PSD and at least one good track (defined below) in the STK.

To check the MC trigger efficiency with the flight data, the events with coincidence of signals from the first two BGO layers, tagged as unbiased triggers, were used. The unbiased trigger events were prescaled by 1/512 at low latitudes (≤20°) and 1/2048 at high latitudes. The HE trigger efficiency is estimated as , where NUnb is the number of events that pass the unbiased trigger condition and the proton selection (without the requirement of the HE trigger) and NHE∣Unb is the number of events that further pass the HE trigger.

Track selection

The track reconstruction algorithm may give more than one track due to the back-scattering particles. Here, we defined the “good” track as the track that has at least four hits in both the xz and yz layers, the reduced χ2 value of the fit is smaller than 25, and the angular deviation from the BGO shower axis is within 5°. In the presence of multiple good tracks in the STK, we selected the “best” one through a combined assessment of the length of the track and the match between the candidate track and the shower axis in the calorimeter. Specifically, the selected tracks were required to be the longest among all the good tracks, among which the one closest to the shower axis was lastly selected. We then applied the geometry cut on the selected track, which is required to pass through all the sublayers of the PSD and the calorimeter from top to bottom.

To validate the STK track efficiency, we selected a proton sample based on the BGO-reconstructed tracks. The track efficiency was estimated as the ratio of the number of events passing the above STK track selection to the total number of events.

Charge selection

The ionization energy loss in the PSD was used to measure the charge of the incoming particle. The charge was measured independently by the two PSD layers, which were averaged to get the final value. The measured charge was corrected for the path length of the particle in the PSD bars using the track information. The PSD and STK alignments were performed to maximally profit the charge reconstruction capabilities of the detector (, ). The charge measurements of the MC simulations showed an energy-dependent difference from that of the flight data, due primarily to the back-scattering particles. To reduce the effect on the charge selection efficiency and the helium background estimate, we applied energy-dependent corrections of the charge measurements for the MC simulations. We first parameterized the charge distributions of protons and helium nuclei in different deposited energy bins with a Landau-Gaussian convolution function and fit the function to the flight data and MC data separately. Then, the MC charges were shifted and shrank according to the best-fitting parameters to match with the flight data. Proton candidates were selected through a cut of the PSD charge. This cut depends on the BGO deposited energy (Edep) as

Note that this charge selection generally follows the logarithmic dependence of the ionization energy loss in the PSD with particle energy. It enables us to have a very small contamination from helium, i.e., less than 1% for deposited energies below 10 TeV and about 5% around 50 TeV. However, this selection results in an energy-dependent selection efficiency for protons, ranging from 94% for incident energies of 0.5 TeV to about 75% above 50 TeV. The charge selection efficiency for protons as a function of incident energy, derived from the GEANT4 FTFP_BERT MC simulations, is shown in Fig. 3A. This efficiency was used in deriving the effective acceptance shown in Fig. 3D.

Fig. 3

Some key information for the proton spectrum measurement.

(A) The charge selection efficiency of protons versus incident energies for the GEANT FTFP_BERT model. (B) The fraction of helium (red open circles) and electron (blue filled dots) backgrounds in the proton candidate events as a function of deposited energy. (C) Probability distribution of deposited energies in the BGO calorimeter for different incident energies, for the GEANT FTFP_BERT model. The color represents the fraction of events in each energy bin. (D) Effective acceptance of protons versus incident energies for the GEANT FTFP_BERT model.

The MC charge selection efficiency of each layer of the PSD was also validated by the flight data. For instance, to estimate the efficiency of the first PSD layer, we used the second PSD layer and the first STK layer measurements of the charge to select the sample and calculated how many of the events have correct charge measurement in the first PSD layer. The differences between MC simulations and the flight data were adopted as systematic uncertainties of the charge selection efficiencies.

Background estimate

The background for protons includes misidentified helium nuclei and a tiny fraction of electrons. The electrons were rejected thanks to different developments between hadronic showers and electromagnetic ones in the BGO calorimeter. The fraction of residual electrons in the proton sample was estimated to be about 0.05% for deposited energies larger than 20 GeV, using the template fit of the shower morphology parameter [ζ as defined in ()]. Helium nuclei are the main source of background for protons. We used the template fit of the PSD charge spectra to estimate the helium backgrounds. The templates were built on the basis of MC simulations (see Fig. 1). The fraction of helium contamination as a function of deposited energy in the BGO calorimeter is shown in Fig. 3B. It is ≲1% for deposited energies below 10 TeV and increases up to ∼5% around 50 TeV.

Energy measurement and spectral deconvolution

The energy of an incident particle was measured by the BGO calorimeter. Because of the limited thickness of the BGO calorimeter (∼1.6 nuclear interaction length) and the missing energy due to muon and neutrino components in hadronic showers, the energy measurements of CR nuclei are biased with some uncertainties. Therefore, MC simulations are required to estimate the energy response of the calorimeter. The energy resolution for protons was found to be about 25 to 35% for incident energies from 100 GeV to 10 TeV (). The linear region of the energy measurement can extend to incident energies of ∼100 TeV, thanks to the maximum reachable energy of 4 TeV for the dynode-2 readout device of each BGO bar (). For very few highest-energy events, the saturation of the maximum energy deposited BGO bar has been corrected on the basis of the simulated transverse shower profile. The Engineering Qualification Model of DAMPE was extensively tested using test beams at the European Organization for Nuclear Research (CERN) in 2014–2015 (, ). The test beam momenta are 5, 10, 150, and 400 GeV/c, respectively. The comparison between the test beam data and MC simulations shows a good agreement with each other (). The on-orbit calibration of the energy measurement was performed by means of the MIP signals in each BGO crystal ().

A deconvolution of the measured energy distribution into the incident energy distribution was applied. The number of events in the ith deposited energy bin, Ndep,, can be obtained via the sum of the number of events in all the incident energy bins, Ninc,, weighted by the energy response matrixwhere M is the probability that an event in the jth incident energy bin is detected in the ith deposited energy bin. We used MC simulations to derive the energy response matrix, applying the same selections as described above. Figure 3C shows the energy response matrix for different incident energies, for the FTFP_BERT model. The color represents the relative probability that a proton with Einc deposits Edep energy in the calorimeter. Equation 2 is solved with a Bayesian method to derive the incident event distribution ().

Acceptance and absolute flux

The effective acceptance is defined as the product of the geometric factor and selection efficiencies (including energy, trigger, track, and charge selections). The effective acceptance for the ith incident energy bin is calculated aswhere Agen is the geometrical factor of the MC generation sphere and Ngen, and Npass, are the numbers of generated events and those passing the selections, respectively. All efficiencies and the effective acceptance were obtained via the MC simulations. For the selection efficiencies, we compared the MC simulations and the flight data with selected control samples (see the “Proton event selection” section), and the differences were adopted as an estimate of the systematic uncertainties of the effective acceptance. Figure 3D shows the effective acceptance as a function of incident energy.

The absolute proton flux F in the incident energy bin [E, E + ΔE] is thenwhere ΔE is the width of the energy bin and Texp is the exposure time.

Systematic uncertainties

Several types of systematic uncertainties were investigated in this analysis, including the event selection, the background subtraction, the spectral deconvolution procedure, and the energy response. The systematic uncertainties related with the event selection were estimated through comparisons between MC simulations and the flight data. The total uncertainty of the selection efficiencies iswhere σtrigger ≈ 2.5%, σtrack ≈ 3.5%, and σcharge ≈ 1.8% are the corresponding systematic uncertainties of the trigger, track selection, and charge selection efficiencies, respectively.

The uncertainties due to the spectral deconvolution were estimated to be ≲1%, through regeneration of the response matrix and varying the spectral index from 2.5 to 3.1 when reweighting the simulation data. The systematic uncertainties due to the helium background subtraction were estimated by varying the charge selection condition (Eq. 1) by ±10% and repeating the analysis. The background subtraction gives ∼0.1% systematic uncertainties below 40 TeV and increases to ∼5% at higher energies. The above systematic uncertainties were added in quadrature to give the total systematic uncertainties associated with the analysis procedure (σana as given in Table 1).

The uncertainties of the fluxes due to different hadronic interaction models were estimated to be about 7% for energies less than 400 GeV via comparisons of the HE trigger efficiency and the energy deposit fraction between the 400-GeV test beam data and the GEANT FTFP_BERT simulation. Note that the 5- and 10-GeV test beam data are out of the interested energy range of the current study, and the 150-GeV test beam is limited by statistics and thus not used. For higher energies, we used the difference between the GEANT FTFP_BERT model and the FLUKA model to estimate such systematic uncertainties, which vary from 7 to 10%. A further check of the DPMJET model with the CRMC () interface gives negligible difference compared with the FLUKA model. Last, the uncertainties associated with the absolute energy scale are about 2% (), which are not corrected in this work.

Figure 4 summarizes the energy-dependent relative uncertainties of the proton fluxes. The statistical uncertainties refer to the Poisson fluctuations of the detected numbers of events in each deposited energy bin. To get the statistical uncertainties of the unfolded proton fluxes, an error propagation from the detected events to the unfolded fluxes is necessary to properly take into account the bin-by-bin migration due to the unfolding procedure (). To obtain a proper estimate of the full error propagation, we ran toy-MC simulations to generate fake observations in each deposited energy bin following the Poisson distribution and obtained the proton fluxes through the unfolding procedure. The root mean squares of the resulting proton fluxes were taken as the 1σ statistical uncertainties. We found that the systematic uncertainties due to the hadronic interaction models dominate in the whole energy range. The STK track efficiency uncertainties are the subdominant ones. In Table 1, the statistical uncertainties, the systematic one associated with the analysis procedure (added quadratically), and the systematic ones due to hadronic models are presented separately.

Fig. 4

Statistical and systematic uncertainties of the proton flux measurements.

Comparison of different spectral models

To quantify the spectral behaviors of the proton spectrum, we used a PL functionor an SBPL oneto fit the data.

We first focused on the high-energy part of the spectrum from 1 to 100 TeV. To properly account for the systematic uncertainties, we adopted a set of independent nuisance parameters w, which are multiplied on the input model (). The χ2 function is defined aswhere E, F, and σstat, are the median energy, flux, and statistical uncertainty of the measurement in the ith energy bin, F(E) is the model predicted flux in corresponding energy bin, S(E; ) is a piecewise function defined by its value w in corresponding energy range covered by the jth nuisance parameter, and is the relative systematic uncertainty of the data in such an energy range. The last term in the right-hand side is a Gaussian prior of the nuisance parameters. The energy range of [1,100] TeV is logarithmically divided into m pieces. According to the energy dependence of the systematic uncertainties in our work (Fig. 5), we adopted m = 4 in the energy band of the fit, which corresponds to two nuisance parameters per decade of energies.

Fig. 5

Comparison between the best fitting of the proton spectrum from 1 to 100 TeV with the SBPL function (solid line) and the DAMPE data.

The fit to the PL model gives χ2/df = 28.6/4, where df represents the number of degrees of freedom. For the SBPL model, we fixed the smoothness parameter s to be 5.0 due to a lack of good constraint on it and got χ2/df = 2.5/2. The results slightly differ for different values of s. The reduction of the χ2 value is about 26.1 for two more free parameters, suggesting a significance of ∼4.7σ in favor of the SBPL model compared with the PL model. The fit of the SBPL model gives GeV−1 m−2 s−1 sr−1, γ = 2.60 ± 0.01, TeV, and Δγ = −0.25 ± 0.07. The comparison of the best-fit result (solid line) with the data is shown in Fig. 5. We also tested the fitting with m = 3 (or 5) and found that the fitting parameters change very little, while the significance in favor of the SBPL model becomes 6.9σ (or 4.2σ).

We also fit the low-energy part of the spectrum from 100 GeV to 6.3 TeV using the SBPL model to address the spectral hardening feature. Again, four nuisance parameters were assumed, and s was fixed to be 5.0, which is close to that obtained by fitting the AMS-02 spectrum (). The fitting model parameters are GeV−1 m−2 s−1 sr−1, γ = 2.772 ± 0.002, Eb = 0.48 ± 0.01 TeV, and Δγ = 0.173 ± 0.007. As comparisons, the fit to the PAMELA data () gives TeV, γ = 2.850 ± 0.016, and Δγ = 0.18 ± 0.06, and the fit to the AMS-02 spectrum () gives TeV, , and , respectively. Our low-energy spectrum is slightly harder than that measured by PAMELA and AMS-02. The break energy inferred from the DAMPE data is roughly consistent with that of AMS-02 but is slightly higher than that of PAMELA. The Δγ values of these three results are consistent with each other. Nevertheless, we should note that the fitting functions and energy ranges adopted in (, ) are not exactly the same as those in our work and the comparison of the detailed numbers should be done with caution.

Note added: During the final stage for the publication of this paper, the CALET collaboration reported new measurements of the proton spectrum from 50 GeV to 10 TeV and confirmed the spectral hardening feature found by other measurements (). The lack of measurements above 10 TeV of CALET, however, hampers a cross-check of our main finding.