COSMIC-RAY HELIUM HARDENING

Recently, the Cosmic Ray Energetics And Mass (CREAM) has directly observed the cosmic-ray (CR) compositions with high statistics in the wide energy range up to about 10¹⁴eV. Interestingly, CREAM shows N_p(E) ∝ E^{−2.66±0.02} for CR protons and N_He(E) ∝ E^{−2.58±0.02} for CR helium in the energy region 2.5 × 10¹²eV to 2.5 × 10¹⁴eV; that is, the spectrum of CR helium is harder than that of CR proton (Ahn et al. 2010). Although the difference in the spectral index Δs ≈ 0.08 appears to be small, the implications are of great importance as shown below. In addition, the spectral index becomes harder by ∼0.12 for CR protons and by ∼0.16 for CR helium at ≳2 × 10¹¹eV nucleon⁻¹ because the Alpha Magnet Spectrometer (AMS) shows N_p(E) ∝ E^{−2.78±0.009} for the CR proton (Alcaraz et al. 2000a) and N_He(E) ∝ E^{−2.74±0.01} for the CR helium (Alcaraz et al. 2000b) in the low-energy range 10¹⁰–10¹¹eV. These results have already been obtained by the Advanced Thin Ionization Calorimeter-2 (ATIC-2; Panov et al. 2009).

The Fermi Gamma-ray Space Telescope has recently observed the spectrum of CR electrons in the wide energy range from 7 × 10⁹eV to 10¹²eV (Ackermann et al. 2010). Fermi shows that the observed data can be fitted by a power law with the spectral index in the interval 3.03–3.13 and the spectral hardening at about 10¹¹eV, which may have the same origin as that of the CR nuclei. (For other models, see, e.g., Kashiyama et al. 2011; Kawanaka et al. 2010; Ioka 2010, and references therein.) Note that we do not discuss CR positrons in this Letter.

Supernova remnants (SNRs) are thought to be the origin of the Galactic CRs. The most popular acceleration mechanism at SNRs is the diffusive shock acceleration (DSA; Axford et al. 1977; Krymsky 1977; Bell 1978; Blandford & Ostriker 1978). In fact, Fermi and AGILE show that middle-age SNRs interacting with molecular clouds emit gamma rays (e.g., Abdo et al. 2009; Tavani et al. 2010) and the gamma-ray observations support that SNRs produce the bulk of Galactic CRs (e.g., Ohira et al. 2011; Li & Chen 2010).

According to DSA theory, the spectrum of accelerated particles at a shock does not depend on CR elements, but depends only on the velocity profile of the shock. Thus, naively, recent CR observations seem to show that the acceleration site of CR helium is different from that of CR protons (Biermann et al. 2010). However, the different site scenario cannot naturally explain why the ratio of CR helium and protons, N_He/N_p, at 10⁹eV is similar to the cosmic abundance (Y = 0.25). Furthermore, the difference in the spectral index Δs ≈ 0.08 means that N_He/N_p at 10¹⁴eV is about three times higher than that at 10⁹eV. This enhancement is amazing because the mean helium abundance in the universe is virtually maintained constant. The stellar nucleosynthesis never enhances the mean helium abundance, which is the essential reason that the Big Bang nucleosynthesis is indispensable for the cosmic helium abundance. To account for the enhancement, we should consider inhomogeneous abundance regions. We show that this leads to a different spectrum of CR protons and helium when escaping from SNRs.

In this Letter, considering the inhomogeneous abundance region, we provide a new explanation about the different spectra of CR protons and helium, even if they are accelerated simultaneously. Our idea uses the fact that CRs escaping from SNRs generally have a different spectrum than those of the acceleration site (Ptuskin & Zirakashvili 2005; Ohira et al. 2010; Caprioli et al. 2010). The runaway CR spectrum depends not only on the acceleration spectrum at shocks but also on the evolution of the maximum energy and the number of accelerated CRs (Ohira et al. 2010). We also suggest that the spectral hardening of CRs is caused by the decreasing Mach number in the high-temperature medium. Both the inhomogeneous abundance and the high temperature can be realized in the superbubbles with multiple supernovae. Our conclusions are summarized as follows.

• 1.  
  Runaway CR spectra depend not only on CR spectra inside the SNR but also on the evolution of the maximum energy and the number of accelerated CRs. Therefore, taking account of the inhomogeneous abundance region, the runaway CR spectra of different CR elements have different spectra (Sections 2 and 3.1).
• 2.  
  Our model is in excellent agreement with the observed spectra of CR protons and helium. The harder spectrum of CR helium is due to the enhancement of the helium abundance around the explosion center. On the other hand, the concave spectra of all CR elements are due to the decreasing Mach number in the hot gas with ∼10⁶K. The concave spectra may also be produced by the CR nonlinear effect, the energy-dependent effects on the accelerated CRs (on α or β), the propagation effect (γ), and/or multi-components with different spectral indices (Sections 3.2 and 4).
• 3.  
  Within the single-component scenario, the hard helium spectrum suggests that the origin of the Galactic CRs is SNRs in superbubbles, although we are not excluding the multi-component scenario (Section 5).
• 4.  
  Our model predicts that heavier (at least volatile) CR elements also have harder spectra than those of CR protons and have concave spectra (Section 5).

In this section, we briefly review the runaway CR spectrum (see the Appendix of Ohira et al. 2010). We here use a variable χ (for example, the shock radius or the SNR age) to describe the evolution of an SNR. Let F_SNR(χ, p) and p_max(χ) be the CR momentum spectrum [(eV/c)⁻¹] and the maximum four-momentum of CR inside the SNR at a certain epoch labeled by χ, respectively. CRs escape in order from the maximum energy CRs because the diffusion length of high-energy CRs is larger than that of low-energy CRs. Then, the number of runaway CRs between χ and χ + dχ is

$$\begin{equation} F_{\rm SNR}(\chi,p_{\rm max}) \frac{{d}p_{\rm max}}{{d}\chi } {d}\chi, \end{equation} \tag{ 1 }$$

which corresponds to the number of runaway CRs between p = p_max(χ) and p = p_max(χ) + dp, F_esc(p)dp. Hence, F_esc(p) is

$$\begin{equation} F_{\rm esc}(p) = F_{\rm SNR}\big(p_{\rm max}^{-1}(p),p\big), \end{equation} \tag{ 2 }$$

where p⁻¹_max(p) is the inverse function of p_max(χ). Assuming F_SNR(χ, p) ∝ χ^βp^−s and p_max(χ) ∝ χ^−α, we obtain the runaway CR spectrum as

$$\begin{equation} F_{\rm esc}(p) \propto p^{-(s+\frac{\beta }{\alpha })}, \end{equation} \tag{ 3 }$$

where α and β are parameters to describe the evolution of the maximum energy and the number of accelerated CRs, respectively. (We use α ∼ 6.5 and β ∼ 1.5 later.) Therefore, the runaway CR spectrum F_esc is different from that in the SNR, F_SNR ∝p^−s. Figure 1 shows the schematic picture of the runaway CR spectrum. In this Letter, we use the shock radius, R_sh, as χ.

Zoom In Zoom Out Reset image size

The evolution of the maximum energy of CRs at the SNR has not been understood. This strongly depends on the evolution of the magnetic field around the shock (e.g., Ptuskin & Zirakashvili 2003). Although some magnetic field amplifications have been proposed (e.g., Lucek & Bell 2000; Bell 2004; Giacalone & Jokipii 2007; Ohira et al. 2009b) and investigated by simulations (e.g., Niemiec et al. 2008; Riquelme & Spitkovsky 2009; Ohira et al. 2009a; Inoue et al. 2009; Gargaté et al. 2010), the evolution of the magnetic field has not been completely understood yet. Here we assume that CRs with the knee energy escape at R = R_Sedov, where R_Sedov is the shock radius at the beginning of the Sedov phase. Furthermore, we use the phenomenological approach with the power-law dependence (Gabici et al. 2009; Ohira et al. 2010),

$$\begin{equation} p_{\rm max}(R_{\rm sh}) = p_{\rm knee} Z \left(\frac{R_{\rm sh}}{R_{\rm Sedov}}\right)^{-\alpha }, \end{equation} \tag{ 4 }$$

where p_knee = 10^15.5eV/c is the four-momentum of the knee energy. Note that α does not depend on the CR composition because the evolution of the maximum energy depends only on the evolution of the magnetic field and the shock velocity.

The evolution of the number of CRs inside the SNR has also not been understood. This depends on the injection mechanism (Ohira et al. 2010) and the density profile around the SNR. We here adopt the thermal leakage model (Malkov & Völk 1995) as an injection model. For the total density profile, ρ_tot(R_sh) ≈ m_p(n_p(R_sh) + 4n_He(R_sh)), where n_p and n_He are the number densities of proton and helium and m_p is the proton mass. The shock velocity of the Sedov phase is

$$\begin{equation} u_{\rm sh}(R_{\rm sh}) \propto \rho _{\rm tot}(R_{\rm sh})^{-\frac{1}{2}} R_{\rm sh}^{-\frac{3}{2}}. \end{equation} \tag{ 5 }$$

In the thermal leakage model, the injection momentum of the element i is proportional to the shock velocity, p_inj,i ∝ u_sh, and the number density of CR with momentum p_inj,i is proportional to the density, p³_inj,if_i(p_inj,i) ∝ n_i(R_sh), where f_i is the distribution function of the CR element i. Hence, the number of the CR element i with a reference momentum p = m_pc, F_SNR,i(R_sh, m_pc) is

$$\begin{eqnarray} F_{{\rm SNR},i}(R_{\rm sh},m_{\rm p}c)&\propto & R_{\rm sh}^3f_i(m_{\rm p}c) \nonumber \\ \ms &\propto & R_{\rm sh}^3p_{{\rm inj},i}^{s_{\rm low}+2}f_i(p_{{\rm inj},i}) \nonumber \\ \ms &\propto & R_{\rm sh}^3 n_i(R_{\rm sh})p_{{\rm inj},i}^{s_{\rm low}-1} \nonumber \\ \ms &\propto & n_i(R_{\rm sh})\rho _{\rm tot}(R_{\rm sh})^{\frac{1-s_{\rm low}}{2}}R_{\rm sh}^{\frac{3(3-s_{\rm low})}{2}}, \end{eqnarray} \tag{ 6 }$$

where $f_i(p)p^2\propto p^{-s_{\rm low}}$ and s_low is the spectral index in the nonrelativistic energy region. For the nonlinear DSA, the spectral index in the nonrelativistic energy region is different from that in the relativistic energy region (Berezhko & Ellison 1999). To understand the essential feature of the runaway CR spectrum, we here consider only the test-particle DSA, that is, s_low = s. Because $n_i(R_{\rm sh})\rho _{\rm tot}(R_{\rm sh})^{\frac{1-s}{2}}$ is not always a single power-law form, the evolution of the number of accelerated CRs cannot be always described by a constant β.

3.1. Different Spectra of CR Protons and Helium

According to the test-particle DSA theory, the index s of the relativistic CR energy spectrum depends only on the velocity jump at the shock,

$$\begin{equation} s = \frac{u_1 +2u_2}{u_1-u_2} = 2\frac{M^{2}+1}{M^{2}-1}, \end{equation} \tag{ 7 }$$

where we use the Rankine–Hugoniot relation at the second equation and M is the Mach number. Then, the index of the runaway CR spectrum, s_esc, is

$$\begin{equation} s_{\rm esc} = s + \frac{\beta }{\alpha }, \end{equation} \tag{ 8 }$$

in Equation (3). Therefore, if β/α (in particular β, the index for the accelerated CR number evolution) is different, the runaway CR spectrum is different between the CR compositions. This is our main idea to explain the helium hardening observed by CREAM and ATIC-2. From Equation (6), β depends on the ambient number density n_i. Therefore, different density profiles make different runaway CR spectra (see Section 4 for more details). Figure 2 shows the schematic picture of our idea.

Zoom In Zoom Out Reset image size

3.2. Spectral Hardening of All CRs at the Same Energy per Nucleon

In this subsection, we discuss the spectral hardening of the observed CRs. The Galactic CR spectrum observed at the Earth, F_obs, is obtained by the simple leaky box model,

$$\begin{equation} F_{\rm obs} \propto F_{\rm esc}(p)/D(p) \propto F_{\rm esc}(p) p^{-\gamma }, \end{equation} \tag{ 9 }$$

where D(p) ∝ p^γ is the diffusion coefficient (e.g., Strong et al. 2007). Hence, the index of the observed spectrum is

$$\begin{equation} s_{\rm obs} = s + \frac{\beta }{\alpha }+\gamma. \end{equation} \tag{ 10 }$$

The deviation from a single power law means that at least one of s, α, β, and γ has an energy dependence or that the origin of low-energy CRs below 10¹¹eV is different from that of high-energy CRs above 10¹¹eV. Although the multi-component scenario may be the case because there are many types of SNRs, we discuss the single-component scenario in this Letter.

First, we discuss the energy dependence of s. From Equation (7), s depends on the shock radius because the Mach number M decreases with the shock radius. Then we can expect the spectral hardening of all CR compositions at the same rigidity cp/Ze, that is, at approximately the same energy per nucleon. From Equation (5), the Mach number is

$$\begin{equation} M \approx 10^3 \left(\frac{\rho _{\rm tot}(R_{\rm sh})}{\rho _{\rm tot}(R_{\rm Sedov})}\right)^{-\frac{1}{2}}\left(\frac{T}{10^4 {\rm K}}\right)^{-\frac{1}{2}} \left(\frac{R_{\rm sh}}{R_{\rm Sedov}}\right)^{-\frac{3}{2}}, \end{equation} \tag{ 11 }$$

where T is the surrounding temperature and we assume that the ejecta mass and the energy of supernova explosion are 1M_☉ and 10⁵¹erg, respectively. From Equations (4), (7), and (11), we can obtain s as a function of p (see Section 4).

Alternatively, the spectral hardening can also be interpreted as the CR nonlinear effect (e.g., Drury & Völk 1981; Malkov & Drury 2001). This issue will be addressed in the future work.

Next, we discuss the energy dependence of β, which is the parameter to describe the evolution of the number of accelerated CRs. In Section 3.1, we consider different power-law forms for n_p(R_sh) and n_He(R_sh) to make the different spectrum of the CR proton and helium. Therefore, ρ_tot(R_sh) ≈ m_p[n_p(R_sh) + 4n_He(R_sh)] is not a single power-law form, and β has an energy dependence (see Section 4).

The energy dependence of γ will be soon precisely determined by AMS-02 (Pato et al. 2010). We do not discuss the energy dependence of α because the complete physics of the CR escape and magnetic turbulence is beyond the scope of this Letter.

In this section, specifying model parameters, we calculate the Galactic CR spectrum. For simplicity, we here assume the number densities of proton and helium as follows:

$$\begin{eqnarray} n_{\rm p}(R_{\rm sh})&=&n_{\rm p,0}\nonumber \\ n_{\rm He}(R_{\rm sh})&=& \zeta n_{\rm p,0} \left(\frac{R_{\rm sh}}{R_{\rm Sedov}}\right)^{-\delta }, \end{eqnarray} \tag{ 12 }$$

where n_p,0 is the number density of proton at R_sh = R_Sedov and ζn_p,0 is the normalization factor of the helium density. We set ζ = 10^{6.5(δ/α)−1} so that the helium abundance is that of the solar abundance, n_He/n_p = 0.1 (i.e., Y ≈ 0.25), when cp_max = Z GeV with Equation (4). Note that the power-law dependence is a first-step approximation for the mean value. Then, from Equations (2), (4), (6), and (9), the observed spectra of CR proton and helium are

$$\begin{eqnarray} F_{\rm obs,p} &=& F_{\rm p,knee}\left\lbrace \frac{1+\zeta \!\left({p}/{p_{\rm knee}}\right)^{\frac{\delta }{\alpha }}}{1+\zeta }\right\rbrace ^{\frac{1-s(p)}{2}} \nonumber \\ && \times \left(\frac{p}{p_{\rm knee}}\right)^{-[s(p)+\frac{3\left\lbrace 3-s(p)\right\rbrace }{2\alpha }+\gamma]}, \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} F_{\rm obs,He} &=& \epsilon F_{\rm p,knee}\left\lbrace \frac{1+\zeta \left({p}/{Zp_{\rm knee}}\right)^{\frac{\delta }{\alpha }}}{1+\zeta }\right\rbrace ^{\frac{1-s(p)}{2}} \nonumber \\ && \times \left(\frac{p}{Zp_{\rm knee}}\right)^{-[s(p)+\frac{3\left\lbrace 3-s(p)\right\rbrace -2\delta }{2\alpha }+\gamma]}, \end{eqnarray} \tag{ 14 }$$

where F_p,knee and F_p,knee are the normalization factors of CR proton and helium, Z = 2 for helium, and s(p) is obtained from Equations (4), (7), and (11). In this model, all parameters are α, γ, δ, , T, F_p,knee.

Figure 3 shows the comparison of our model with observations. We take into account the solar modulation effects with the modulation potential Φ = 450 MV (Gleeson & Axford 1968). Our model is in excellent agreement with the observed spectra, with α = 6.5, γ = 0.43, δ = 0.715, = 0.31, and T = 10⁶ K. The different spectra of CR protons and helium originate from the different density profiles in Equation (12). Figure 4 shows the evolution of the maximum energy of CRs and the spectral index of CRs inside the SNR. In the early phase, the spectral index s is 2 and after then the spectral index decreases with the shock radius because the Mach number decreases with the shock radius. The change of the spectral index s is about 0.1 which is almost the same as the observed hardening. The observed hardening is not the result of the change of the injection history, β, but the result of a spectral change of CRs inside the SNR. The high temperature T ∼ 10⁶ K is necessary for the spectral hardening Δs ∼ 0.1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition, our model also results in a concave spectrum of CR electrons, as observed (Ackermann et al. 2010). However, the evolution of the injection efficiency of CR electrons has not been understood well. So, we need further studies to discuss the CR electron spectrum in detail.

To make the different spectrum, our model requires that the helium abundance around the explosion center is higher than that of the solar abundance. SNRs in superbubbles are one of such candidates. Higdon et al. (1998) show that supernova ejecta can dominate the superbubble mass within a core radius of one-third of the superbubble radius. In the stellar wind and the supernova explosion, the stellar hydrogen envelope has lower density and higher velocity than that of helium. Then we expect that the helium fraction in the center of superbubbles is higher than that in the outer region. Furthermore, to make the concave spectrum, our model requires an ambient medium with high temperature, T = 10⁶K. This is also consistent with superbubbles. According to the CR composition study, SNRs in superbubbles have been considered as the origin of Galactic CRs (e.g., Lingenfelter & Higdon 2007; Ogliore et al. 2009). Particle accelerations in superbubbles have also been investigated by intensive studies (e.g., Bykov & Fleishman 1992; Parizot et al. 2004; Dar & De Rújula 2008; Ferrand & Marcowith 2010). We here considered a spherically symmetric system. The off-center effects may be important for the initial phase and thereby for the high-energy spectrum, because the shock radius at the beginning of the Sedov phase, R_Sedov, is about 20 pc which is comparable to the typical size of OB association, 35 pc (Parizot et al. 2004), and the shock radius R_sh is about 200 pc at the end of the Sedov phase. This is an interesting problem for future work.

Note that the spectral hardening can also be made by the nonlinear model, the energy dependence of the CR diffusion coefficient, and/or multi-components with different spectral indices. So, the high temperature may not be absolutely necessary. The stellar wind of red giants is one of the candidates for the cold and helium-rich ambient medium. Still, the dominant core-collapse supernova is Type II (e.g., Smartt et al. 2009), which has no helium-rich wind; therefore, the superbubble scenario looks more likely as the origin of the Galactic CRs above 10¹¹eV. For the CRs below 10¹¹eV, the spectral difference between CR protons and helium may be caused by the solar modulation and the inelastic interactions (Putze et al. 2011).

The spatial variation of the helium ionization degree can also change the injection history. The injection efficiency of the large rigidity is thought to be higher than that of the low rigidity since particles with large rigidity can easily penetrate through the shock front from the downstream region. If the ionization degree increases with the SNR radius, the CR helium spectrum becomes harder than the CR proton one, β_He < β_p. However, the rigidity dependence of the injection efficiency has not been understood completely. Moreover, the injection from neutral particles should also be understood (Ohira et al. 2009b; Ohira & Takahara 2010).

According to our model, CR spectra of volatile elements heavier than helium are also harder than that of protons. Low-energy CRs of refractory elements are thought to result from suprathermal injection by sputtering off preaccelerated, high-velocity grains (Ellision et al. 1997). To be accelerated to the relativistic energy, the refractory elements should be sputtered because the grains cannot be accelerated to the relativistic energy. The SNR shock velocity is not fast enough to accelerate the refractory elements to the knee energy when the refractory elements are injected because the sputtering timescale is too long. Therefore, refractory CRs around the knee energy should be injected by the standard manner similar to volatile CRs. In this case, the refractory CRs also have harder spectra than protons, although we need further studies of the injection of the refractory CRs at the knee energy.

If CRs trapped inside the SNR and released at the end of the SNR's life outnumber runaway CRs (see Figure 3 in Caprioli et al. 2010), our scenario does not work for producing hard and concave spectra. In our model with α ∼ 6.5 in Equation (4), trapped CRs have energy below 1GeV when they are released, that is, p_max ≲ Zm_pc when R_sh ≳ 10R_Sedov, and are not relevant for our interest. Higher energy CRs escape from the SNR even after they are advected to the downstream since the CR diffusion is faster than the expansion of the SNR. Our case is similar to the right of Figure 7 in Caprioli et al. (2010) where trapped CRs are released below 100 GeV. The energy boundary between trapped CRs and runaway CRs depends on the evolution of the maximum energy (α).

We thank the referees T. Suzuki, T. Terasawa, and A. Bamba for comments. This work is supported in part by grant-in-aid from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, no. 21684014 (Y.O. and K.I.), nos. 19047004, 22244019, 22244030 (K.I.).