Direct Measurement of the Cosmic-Ray Carbon and Oxygen Spectra from 10 GeV$/n$ to 2.2 TeV$/n$ with the Calorimetric Electron Telescope on the International Space Station

In this paper, we present the measurement of the energy spectra of carbon and oxygen in cosmic rays based on observations with the Calorimetric Electron Telescope (CALET) on the International Space Station from October 2015 to October 2019. Analysis, including the detailed assessment of systematic uncertainties, and results are reported. The energy spectra are measured in kinetic energy per nucleon from 10 GeV$/n$ to 2.2 TeV$/n$ with an all-calorimetric instrument with a total thickness corresponding to 1.3 nuclear interaction length. The observed carbon and oxygen fluxes show a spectral index change of $\sim$0.15 around 200 GeV$/n$ established with a significance $>3\sigma$. They have the same energy dependence with a constant C/O flux ratio $0.911\pm 0.006$ above 25 GeV$/n$. The spectral hardening is consistent with that measured by AMS-02, but the absolute normalization of the flux is about 27% lower, though in agreement with observations from previous experiments including the PAMELA spectrometer and the calorimetric balloon-borne experiment CREAM.

In this paper, we present the measurement of the energy spectra of carbon and oxygen in cosmic rays based on observations with the Calorimetric Electron Telescope (CALET) on the International Space Station from October 2015 to October 2019. Analysis, including the detailed assessment of systematic uncertainties, and results are reported. The energy spectra are measured in kinetic energy per nucleon from 10 GeV/n to 2.2 TeV/n with an all-calorimetric instrument with a total thickness corresponding to 1.3 nuclear interaction length. The observed carbon and oxygen fluxes show a spectral index change of ∼0.15 around 200 GeV/n established with a significance > 3σ. They have the same energy dependence with a constant C/O flux ratio 0.911 ± 0.006 above 25 GeV/n. The spectral hardening is consistent with that measured by AMS-02, but the absolute normalization of the flux is about 27% lower, though in agreement with observations from previous experiments including the PAMELA spectrometer and the calorimetric balloon-borne experiment CREAM.

INTRODUCTION
Direct measurements of charged cosmic rays (CR) provide information on their origin, acceleration and propagation in the Galaxy. Search for possible chargedependent cutoffs in the nuclei spectra, hypothesized to explain the "knee" in the all-particle spectrum [1][2][3][4], can be pursued at very large energy by magnetic spectrometers with sufficient rigidity coverage (MDR) or by calorimetric instruments equipped with charge detectors capable of single element resolution. Furthermore, recent observations indicating a spectral hardening in proton and He spectra [5][6][7][8][9][10] as well as in heavy nuclei spectra [11][12][13][14] around a few hundred GeV/n, compelled a revision of the standard paradigm of galactic CR based on diffusive shock acceleration in supernova rem-nants followed by propagation in galactic magnetic fields, and prompted an intense theoretical activity to interpret these unexpected spectral features [15][16][17][18][19][20][21][22][23][24][25][26]. The CALorimetric Electron Telescope (CALET) [27][28][29] is a spacebased instrument optimized for the measurement of the all-electron spectrum [30,31], which can also measure individual chemical elements in CR from proton to iron and above in the energy range up to ∼1 PeV. CALET recently confirmed the spectral hardening in the proton spectrum by accurately measuring its power-law spectral index over the wide energy range from 50 GeV to 10 TeV [32]. In this Letter, we present a new direct measurement of the CR carbon and oxygen spectra from 10 GeV/n to 2.2 TeV/n, based on the data collected by CALET from October 13, 2015 to October 31, 2019 aboard the Inter-national Space Station (ISS).

CALET INSTRUMENT
CALET consists of a CHarge Detector (CHD), a finely segmented pre-shower IMaging Calorimeter (IMC), and a Total AbSorption Calorimeter (TASC). CHD is comprised of two hodoscopes made of 14 plastic scintillator paddles each, arranged in orthogonal layers (CHDX, CHDY). The CHD can resolve individual chemical elements from Z = 1 to Z = 40, with excellent charge resolution. The IMC consists of 7 tungsten plates inserted between eight double layers of 1 mm 2 cross-section scintillating fibers, arranged in belts along orthogonal directions and individually read out by multianode photomultiplier tubes. Its fine granularity and imaging capability allow an accurate particle tracking and an independent charge measurement via multiple samples of the particle ionization energy loss (dE/dx) in each fiber. The TASC is a homogeneous calorimeter made of leadtungstate (PbWO 4 ) bars arranged in 12 layers. The crystal bars in the top layers are read out by photomultiplier tubes, while a dual photodiode/avalanchephotodiode (PD/APD) system is used for each channel in the remaining layers. A dynamic range of more than six orders of magnitude is covered using a front-end electronics with dual gain range for each photosensor. The total thickness of the instrument is equivalent to 30 radiation lengths and 1.3 nuclear interaction lengths. A more complete description of the instrument can be found in the Supplemental Material (SM) of Ref. [30]. CALET was launched on August 19, 2015 and installed on the Japanese Experiment Module Exposure Facility of the ISS. The on-orbit commissioning phase aboard the ISS was successfully completed in the first days of October 2015, and since then the instrument has been taking science data continuously [33].

DATA ANALYSIS
We have analyzed flight data (FD) collected in 1480 days of CALET operation. The total observation live time for the high-energy (HE) shower trigger is T = 3.00 × 10 4 hours, corresponding to 84.5% of total observation time. Raw data are corrected for non-uniformity in light output, time and temperature dependence, gain differences among the channels. The latter are individually calibrated on orbit by using penetrating proton and He particles, selected by a dedicated trigger mode [34,35]. After calibrations, each CR particle track is reconstructed and a charge and an energy are assigned for each event.
Monte Carlo (MC) simulations, reproducing the detailed detector configuration, physics processes, as well as detector signals, are based on the EPICS simulation package [36,37] and employ the hadronic interaction model DPMJET-III [38]. An independent analysis based on FLUKA [39,40] is also performed to assess the systematic uncertainties.
The CR particle direction and its entrance point in the instrument are reconstructed by a track finding and fitting algorithm based on a combinatorial Kalman filter [41], which is able to identify the incident track in the presence of a background of secondary tracks backscattered from TASC. The angular resolution is ∼ 0.1 • for C and O nuclei and the spatial resolution on the determination of the impact point on CHD is ∼220 µm.
The identification of the particle charge Z is based on the measurements of the ionization deposits in CHD and IMC. The particle trajectory is used to identify the CHD paddles and IMC fibers traversed by the primary particle and to determine the path length correction to be applied to the signals to extract the dE/dx samples. Three independent dE/dx measurements are obtained, one for each CHD layer and the third by averaging the samples (at most eight) along the track in the top half of IMC. Calibration curves of dE/dx are built by fitting FD subsets for each nuclear species to a function of Z 2 by using a "halo" model [42]. These curves are then used to reconstruct three charge values (Z CHDX , Z CHDY , Z IMC ) from the measured dE/dx on an event-by-event basis [43]. For high-energy showers, the charge peaks are corrected for the systematical shift to higher values (up to 0.15 e) with respect to the nominal charge positions, due to the large amount of shower particle tracks backscattered from TASC whose signals add up to the primary particle ionization signal. A charge distribution obtained by averaging Z CHDX and Z CHDY is shown in Fig. S1 of the SM [44]. The charge resolution σ Z is ∼ 0.15 e (charge unit) for CHD and ∼ 0.24 e for IMC, respectively, in the elemental range from B to O.
The shower energy E TASC of each event is calculated as the sum of the energy deposits of all the TASC channels, after stitching the adjacent gain ranges of each PD/APD. The energy response of TASC was studied in a beam test carried out at CERN-SPS in 2015 with accelerated ion fragments of 13, 19 and 150 GeV/c momentum per nucleon [45]. The MC simulations were tuned using the beam test results as described in the Energy measurement section of the SM [44].
Carbon and oxygen candidates are identified among events selected by the onboard HE shower trigger, based on the coincidence of the summed signals of the last two IMC layers in each view and the top TASC layer (TASCX1). Consistency between MC and FD for triggered events is obtained by an offline trigger with higher thresholds (50 and 100 times a minimum ionizing particle (MIP) signal for IMC and TASC, respectively) than the onboard trigger removing possible effects due to residual non-uniformity of the detector gain.
In order to reject possible events triggered by particles entering the TASC from lateral sides or with significant lateral leakage, the energy deposits in the first TASC layer (TASCX1) and in all the lateral bars are required to be less than 40% of E TASC . Late-interacting events in the bottom half of TASC are rejected by requiring that the energy deposit in the last layer is < 0.4 × E TASC , and the layer, where the longitudinal shower development reaches 20% of E TASC , occurs in the upper half of TASC. Events with one well-fitted track crossing the whole detector from CHD top to the TASC bottom layer and at least 2 cm away from the edges in TASCX1 are then selected. The fiducial geometrical factor for this category of events is SΩ ∼510 cm 2 sr, corresponding to about 50% of the total CALET acceptance. Carbon and oxygen candidates are selected by applying window cuts, centered on the nominal charge values (Z = 6, 8), of half-width 0.4 e for Z CHDX and Z CHDY , and 2σ Z for Z IM C , respectively. Particles undergoing a charge-changing nuclear interaction in the upper part of the instrument (Fig. S2 of the SM [44]) are removed by the three combined charge selections and by requiring the consistency, within 30%, between the mean values of dE/dx measurements in the first four layers in each IMC view.
Distributions of E TASC for C and O selected candidates are shown in Fig. S3 of SM [44], corresponding to 6.154×10 5 C and 1.047×10 6 O events, respectively. In order to take into account the relatively limited energy resolution (Fig. S4 of the SM [44]) energy unfolding is necessary to correct for bin-to-bin migration effects. In this analysis, we used the Bayesian approach [46] implemented in the RooUnfold package [47] in ROOT [48]. Each element of the response matrix represents the probability that primary nuclei in a certain energy interval of the CR spectrum produce an energy deposit in a given bin of E TASC . The response matrix (Fig. S5 of the SM [44]) is derived using MC simulation after applying the same selection procedure as for FD. The energy spectrum is obtained from the unfolded energy distribution as follows: where ∆E denotes energy bin width, E the particle kinetic energy, calculated as the geometric mean of the lower and upper bounds of the bin, N (E) is the bin content in the unfolded distribution, ε(E) the total selection efficiency (Fig. S6 of the SM [44]), U () the unfolding procedure, N obs (E TASC ) the bin content of observed energy distribution (including background), N bg (E TASC ) the bin content of background events in the observed energy dis-tribution. Background contamination from different nuclear species misidentified as C or O is shown in Fig. S3 of the SM [44]. A contamination fraction N bg /N obs < 0.1% is found in all energy bins with E TASC < 10 3 GeV, and between 0.1% and 1% for E T ASC > 10 3 GeV.

SYSTEMATIC UNCERTAINTIES
In this analysis, dominant sources of systematic uncertainties include trigger efficiency, energy response, event selection, unfolding procedure, MC model. HE trigger efficiency as a function of E TASC was inferred from the data taken with a minimum bias trigger. HE efficiency curves for C and O are consistent with predictions from MC simulations, as shown in Fig. S7 of the SM [44]. In order to study the flux stability against offline trigger efficiency, the threshold applied to TASCX1 signal was scanned between 100 and 150 MIP signal. The corresponding systematic errors range between -4.2% (-3.1%) and 3.7% (7.3%) for C (O) depending on the energy bin. The systematic error related to charge identification was studied by varying the width of the window cuts between 0.35 e and 0.45 e for CHD and between 1.75 σ Z and 2.2 σ Z for IMC. That results in a flux variation depending on the energy bin, which is less than 1% below 250 GeV/n and few percent above. The ratio of events selected by IMC charge cut to the ones selected with CHD in different E TASC intervals turned out to be consistent in FD and MC. Possible inaccuracy of track reconstruction could affect the determination of the geometrical acceptance. The contamination due to off-acceptance events which are mis-reconstructed in the fiducial acceptance was estimated with MC to be ∼ 1% at 10 GeV/n and decrease to less than 0.1% above 60 GeV/n. To investigate the uncertainty in the definition of the acceptance, restricted acceptance (up to 20% of nominal one) regions were also studied. The corresponding fluxes are consistent within statistical fluctuations. A different tracking procedure, described in Ref. [49], was also used to study possible systematic uncertainties in tracking efficiency. Results are consistent with those obtained with the Kalman filter algorithm, hence we consider negligible this source of systematic error. The uncertainty in the energy scale is ±2% and depends on the accuracy of the beam test calibration. It causes a rigid shift of the measured energies, affecting the absolute normalization of the C and O spectra by +2.6% −2.8% , but not their shape. As the beam test model was not identical to the instrument now in orbit, the difference in the spectrum obtained with either configuration was modeled and included in the systematic error. Other energy-independent systematic uncertainties affecting the normalization include live time (3.4%, as explained in the SM of Ref. [30]) and long-term stability of the charge measurements (< 0.4%). The uncertainties due to the unfolding procedure were evaluated by using different response matrices, computed by varying the spectral index (between -2.9 and -2.5) of the generation spectrum of MC simulations, and the Singular Value Deconvolution method, instead of the Bayesian approach, in RooUnfold software [47]. Since it is not possible to validate MC simulations with beam test data in the high-energy region, a comparison between different MC models, i.e. EPICS and FLUKA, was performed. We found that the total selection efficiencies for C and O determined with the two models are in agreement within < 1.5% over the whole energy range, but the energy response matrices differ significantly in the low and high energy regions. The resulting fluxes show maximum discrepancies of 9% (7.8%) and 9.2% (12.2%), respectively, in the first and last energy bin for C (O), while they are consistent within 6.6% (6.2%) elsewhere. This is the dominant source of systematic uncertainties. Materials traversed by nuclei in IMC are mainly composed of carbon, aluminum and tungsten. Possible uncertainties in the inelastic cross sections in simulations or discrepancies in the material description might affect the flux normalization. We have checked that hadronic interactions are well simulated in the detector, by measuring the survival probabilities of C and O nuclei at different depths in IMC, as described in the SM [44]. The survival probabilities are in agreement with MC prediction within < 1% (Fig. S8 of the SM [44]). Background contamination from different nuclear species estimated with FLUKA and EPICS simulations differ by less than 1%. The energy dependence of all the systematic uncertainties for C and O is shown in Fig. S9 of the SM [44]. The total systematic error is computed as the sum in quadrature of all the sources of systematics in each energy bin.

RESULTS
The energy spectra of carbon and oxygen and their flux ratio measured with CALET in an energy range from 10 GeV/n to 2.2 TeV/n are shown in Fig. S10, where current uncertainties that include statistical and systematic errors are bounded within a gray band. CALET spectra are compared with results from space-based [14,51,52,56,57] and balloon-borne [50,[53][54][55] experiments. The measured C and O fluxes and flux ratio with statistical and systematic errors are tabulated in Tables I,  II and III of the SM [44]. Our spectra are consistent with PAMELA [56] and most previous experiments [51][52][53][54][55], but the absolute normalization is in tension with AMS-02 [14]. However we notice that C/O ratio (Fig. S10 (c)) is consistent with the one measured by AMS-02. In Fig. S11 of the SM [44], it is shown that CALET and and (c) ratio of carbon to oxygen fluxes, as a function of kinetic energy E. Error bars of CALET data (red) represent the statistical uncertainty only, while the gray band indicates the quadratic sum of statistical and systematic errors. Also plotted are other direct measurements [14,[50][51][52][53][54][55][56][57]]. An enlarged version of the figure is available as Fig. S10 in the SM [44].
. AMS-02 C and O spectra have very similar shapes but they differ in the absolute normalization, which is lower for CALET by about 27% for both C and O. Figure 2 shows the fits to CALET carbon and oxygen data with a double power-law function (DPL, Eq. S1 in SM [44]) above 25 GeV/n. A single power-law function (SPL, Eq. S2 in SM [44])) fitted to data in the energy range [25,200] GeV/n and extrapolated above 200 GeV/n is also shown for comparison. The effect of systematic uncertainties in the measurement of the energy spectrum is modeled in the χ 2 minimization function with a set of 6 nuisance parameters as explained in detail in the SM [44]. The DPL fit to the C spectrum yields a spectral index γ = −2.663 ± 0.014 at energies below the transition region E 0 = (215 ± 54) GeV/n and a spectral index increase ∆γ = 0.166 ± 0.042 above, with χ 2 /d.o.f. = 9.0/8. For oxygen, the fit yields γ = −2.637 ± 0.009, E 0 = (264 ± 53) GeV/n, ∆γ = 0.158 ± 0.053, with χ 2 /d.o.f. = 3.0/8. SPL fits to CALET carbon and oxygen data above 25 GeV/n are shown in Figs. S12 and S13 of the SM [44]. They give γ = −2.626 ± 0.010 with χ 2 /d.o.f. = 27.5/10 for C, and γ = −2.622 ± 0.008 with χ 2 /d.o.f. = 15.9/10 for O, respectively. A frequentist test statistic ∆χ 2 is computed from the difference in χ 2 between the fits with SPL and DPL functions. For carbon (oxygen), ∆χ 2 = 18.5 (12.9) with 2 d.o.f. (i.e. the number of additional free parameters in DPL fit with respect to SPL fit) implies that the significance of the hardening of the C (O) spectrum exceeds the 3σ level. We also checked that the spectral hardening is not an artifact of the energy binning and unfolding, by increasing the bin width by a factor 2.5 as shown in Fig. S14 of the SM [44]. The resulting flux difference is negligible when compared with our estimated systematic uncertainties. In order to study the energy dependence of the spectral index in a model independent way, the spectral index γ is calculated by a fit of d[log(Φ)]/d[log(E)] in energy windows centered in each bin and including the neighbor ±3 bins. The results in Fig. 3 show that carbon and oxygen fluxes harden in a similar way above a few hundred GeV/n. The carbon to oxygen flux ratio is well fitted to a constant value of 0.911±0.006 above 25 GeV/n (Fig. S15 of the SM [44]), indicating that the two fluxes have the same energy dependence.

CONCLUSION
With a calorimetric apparatus in low Earth orbit, CALET has measured the energy spectra of carbon and oxygen nuclei in CR and their flux ratio from 10 GeV/n to 2.2 TeV/n. Our observations allow to exclude a single power law spectrum for C and O by more than 3σ; they show a spectral index increase ∆γ = 0.166 ± 0.042 for C and ∆γ = 0.158 ± 0.053 for O above 200 GeV/n, and the same energy dependence for C and O fluxes with a constant C/O flux ratio 0.911 ± 0.006 above 25 GeV/n. These results are consistent with the ones reported by AMS-02. However the absolute normalization of our data is significantly lower than AMS-02, but in agreement with previous experiments. Improved statistics and refinement of the analysis with additional data collected during the lifetime of the mission will allow to extend the measurements at higher energies and improve the spectral analysis, contributing to a better understanding of the origin of the spectral hardening.   Contamination of each nuclear species with Z > 4 is estimated by rescaling its ET ASC distribution measured in FD by the ratio, estimated with MC simulations, of its reconstruction efficiency to the probability of being misidentified as C or O. Background due to proton and helium is computed by normalizing their ET ASC distributions from MC to the number of events expected from previous flux measurements.

ENERGY MEASUREMENT
Differently from the case of electrons, the energy released in TASC by interacting CR nuclei is only a fraction of the primary particle energy, the electromagnetic component of the hadronic cascades, originating from the decays of π 0 secondaries produced in the showers. Though a significant part of the hadronic cascade energy leaks out of the calorimeter because of its limited thickness (1.3 λ I ), the energy deposited in the TASC by the electromagnetic shower core scales linearly with the incident particle energy, albeit with large event-to-event fluctuations. As a result, the energy resolution is poor by the standards of total containment hadron calorimetry in experiments at accelerators. Nevertheless, it is sufficient to reconstruct the steep energy spectra of CR nuclei with a nearly energy independent resolution. The TASC response to nuclei was studied at CERN SPS in 2015 using a beam of accelerated ion fragments with A/Z = 2 and kinetic energy of 13, 19 and 150 GeV/n [S1]. In Fig. S4a, the E TASC distributions for C nuclei at 150 GeV/n is shown as an example. C nuclei in the beam are selected with CHD and the HE trigger is applied. The resulting distribution looks nearly gaussian, the energy released in the TASC is ∼20% of the particle energy and the resolution σ E is close to 30%. The mean energy deposited in TASC by different nuclear species in the beam (selected with CHD) is plotted as a function of the kinetic energy per particle in Fig. S4b. The energy response of TASC is linear up to the maximum available particle energy of 6 TeV (obtained with a primary beam of 40 Ar nuclei). The energy response derived from MC simulations was tuned using the beam test results. Correction factors are 6.7% for E TASC < 45 GeV and 3.5% for E TASC > 350 GeV, respectively, while a simple linear interpolation is used to determine the correction factor for intermediate energies.
For flux measurement, energy unfolding is applied to correct E TASC distributions of selected C and O candidates for significant bin-to-bin migration effects (due to the limited energy resolution) and infer the primary particle energy. In this analysis, we apply the iterative unfolding method based on the Bayes' theorem [S2] implemented in the RooUnfold package [S3, S4]. The response matrix is derived using MC simulations of the CALET flight model after applying the same selection as for FD and used in the unfolding procedure (Fig. S5). Each element of the matrix represents the probability that primary nuclei in a certain energy interval of the CR spectrum produce an energy deposit in a given E TASC bin.   . The data are fitted to a gaussian (blue line); the mean value is 20% of the kinetic beam particle energy, and the energy resolution (defined as the standard deviation to mean ratio) is ∼30%. (b) Energy linearity of TASC as measured at CERN SPS with beams of 19 (blue dots) and 150 (black squares) GeV/n ion fragments with A/Z = 2. The red line represents a linear fit to the data. The fitted slope is 0.194 ± 0.005, indicating that on average ∼20% of the particle energy is deposited in TASC.  . The FD to MC efficiency ratio is shown in the lower canvas. The HE trigger efficiency was measured directly from the data by using dedicated runs where in addition to HE, a low-energy (LE) trigger was active. The trigger logic is the same for both trigger (i.e. coincidence of the pulse heights of the last two pairs of IMC layers and the top TASC layer) but lower discriminator thresholds are set for the input signals in case of LE trigger, allowing to trigger also penetrating nuclei with Z > 2. The ratio of the number of events counted by both triggers to those registered by the LE trigger only is an estimate of the HE trigger efficiency in each bin of deposited energy. To compare with simulations, we apply exactly the same method to MC samples, where both trigger modes are modeled.

SURVIVAL PROBABILITIES
In order to check that hadronic interactions in the detector are well simulated, we have measured the survival probabilities of C and O nuclei at different depths in IMC and compared with the ones expected from MC simulations. For each HE-triggered event, selected as C or O by means of CHD only, several measurements of the charge are obtained by dE/dx samples along the particle track in pairs of adjacent scintillating-fiber (SciFi) layers. The coupled SciFi layers are spaced by about 2.3 cm and interleaved with W plates (except for the first pair), aluminum honeycomb panels and CFRP (carbon fiber reinforced polymer) supporting structures. From the charge distribution of each pair of Scifi layers, the number of events that did not interact in the material upstream the fibers can be estimated by selecting the nuclei with a charge value consistent with the one measured in CHD. The reduction in the number of not-interacted events, normalized to the ones selected in CHD, allows to measure the survival probability as a function of the material thickness traversed by the particle in the IMC, as shown in Fig. S8. In MC data, the survival probabilities can be calculated by using the true information on the point where the first hadronic interaction occurs in the detector. The survival probabilities measured in FD are in good agreement (within < 1%) with MC predictions, as shown in the bottom panels of Fig. S8.     [S5]. Error bars of CALET data represent the quadratic sum of statistical and systematic uncertainties. The factor (1.27±0.03) is obtained by a simultaneous minimization of the sum of squared residuals between CALET and AMS data points for C and O.

SPECTRAL FIT METHOD
We have tested two models for the C and O energy spectra: a double power-law (DPL) in energy and a single power-law (SPL) where C is a normalization factor, γ the spectral index, and ∆γ the spectral index change above the transition energy E 0 . We have fitted the data using a chi-square minimization. The effect of systematic uncertainties in the measurement of the energy spectrum is modeled introducing in the χ 2 function a set of nuisance parameters. The χ 2 function is defined as i , E i are the measured flux, the statistical and systematic errors, and the energy of the i-th bin (geometric mean of bin edges), respectively; N = 22 is the number of data points and the fit starts from the fourth point. The elements of the vector p are the free parameters of the model function Φ used in the fit: p = {C, γ, ∆γ, E 0 } for Eq. S1, and p = {C, γ} for Eq. S2. β is a set of independent nuisance parameters β j with j = 1, ..., N np introduced to properly account for systematic uncertainties in the fit. The energy range to be fitted (from 25 to 2200 GeV/n) is logarithmically divided into N np intervals. According to the energy dependence of the systematic fractional errors (Fig. S9), we chose N np = 6; therefore the width of each interval covers three consecutive energy bins of the spectrum. A nuisance parameter β j is assigned to each interval. In the fit, a gaussian constraint (the last term in Eq. S3) is applied to each parameter β j , with zero mean and standard deviation equal to 1. In the χ 2 function, the systematic errors σ sys i are multiplied by a piece-wise function S(E i , β), which assumes the value of β j of the corresponding energy interval within which E i falls. The results of the fits with the two models are shown in Figs. S12 and S13 for C and O spectrum, respectively. We also performed the fits by varying the number of nuisance parameters, with similar results and significance. FIG. S12. Fit to the carbon energy spectrum with (a) a DPL function (Eq. S1) and (b) a SPL function (Eq. S2), in the energy range [25,2000] GeV/n. Error bars of CALET data points represent the sum in quadrature of statistical and systematic uncertainties. The DPL fit yields γ = −2.663 ± 0.014, E0 = (215 ± 54) GeV/n, ∆γ = 0.166 ± 0.042, with χ 2 /d.o.f. = 9.0/8. The SPL fit yields γ = −2.626 ± 0.010 with χ 2 /d.o.f. = 27.5/10. The difference ∆χ 2 = 18.5 between the fits with the two models, with two additional free parameters in DPL fit with respect to SPL fit, allows to exclude the single power law hypothesis at the 3.9σ level. FIG. S13. Fit to the oxygen energy spectrum with (a) a DPL function (Eq. S1) and (b) a SPL function (Eq. S2), in the energy range [25,2000] GeV/n. Error bars of CALET data points represent the sum in quadrature of statistical and systematic uncertainties. The DPL fit yields γ = −2.637 ± 0.009, E0 = (264 ± 53) GeV/n, ∆γ = 0.158 ± 0.053, with χ 2 /d.o.f. = 3.0/8. The SPL fit yields γ = −2.622 ± 0.008 with χ 2 /d.o.f. = 15.9/10. The difference ∆χ 2 = 12.9 between the fits with the two models, with two additional free parameters in DPL fit with respect to SPL fit, allows to exclude the single power law hypothesis at the 3.2σ level. FIG. S14. CALET carbon (a) and oxygen (b) spectra derived with two different energy binnings. In the reference spectra (red squares) the energy range from 10 GeV/n to 2.2 TeV/n is divided into 22 bins (Tables I and II); all bins are chosen to have relative width commensurate with the TASC energy resolution σE, with the exception of the last bin whose width is 2 × σE. To study possible binning related effects in the spectra, C and O fluxes are also derived by dividing the energy range into 8 bins (blue dots), with relative bin width ∼ 2.5 × σE, but the last bin whose width is 3.5 × σE. The wide-bin spectra have the same shape as the original ones. The flux differences are well within the error bars, representing the quadratic sum of statistical and systematic uncertainties.  10.0-12.6 (2.471 ± 0.008 +0.106