ENERGY SPECTRUM AND CHEMICAL COMPOSITION OF ULTRAHIGH ENERGY COSMIC RAYS FROM SEMI-RELATIVISTIC HYPERNOVAE

Wang et al. (2007, 2008) have suggested that hypernova remnants can accelerate heavy or intermediate-mass nuclei to ultrahigh energies above 10¹⁹ eV. Here for simplicity we only focus on the semi-relativistic part of the ejecta, since only this part is relevant to the acceleration of UHECRs. Particles are accelerated in the shock region where the semi-relativistic ejecta is freely expanding before being decelerated by the swept-up circum-stellar medium. The size of this free-expansion phase region for ejecta of a particular velocity β_shc and kinetic energy E_K is $R_{\rm HN}=E_k v_w/\Gamma ^2c^2\dot{M} \simeq 2\times 10^{17} E_{k,51} \Gamma ^{-2} \dot{M}_{-5}^{-1} v_{w,3} \, {\rm cm}$, where Γ = (1 − β²_sh)^−1/2 is the bulk Lorentz factor of the ejecta, ${\dot{M}}=10^{-5}\dot{M}_{-5}\ {M_\odot \ \rm yr^{-1}}$ is the wind mass loss rate (whose average value is 3 × 10⁻⁵ M_☉ yr⁻¹ for W-R stars), and v_w = 10³v_{w, 3} km s⁻¹ is the wind velocity (Willis 1991; Chevalier & Li 1999). We assume that a fraction of _B of shock internal energy goes into the magnetic field and hence the magnetic field energy density is B²/8π ≈ 2Γ²β²_Bρ_wc², where $\rho _w=\dot{M}/4\pi R_{\rm HE}^2v_w$ is the mass density of the stellar wind at radius R_HN. Thus, by equating the adiabatic cooling time with the acceleration time, the maximum energy of accelerated particles is

$$\begin{eqnarray} E_{\rm max}&\simeq& Z e \Gamma BR_{\rm HE}\beta _{\rm sh}\nonumber\\ &=& 3.7\times 10^{20} \left(\frac{Z}{26}\right) \epsilon _{B,-1}^{1/2}\Gamma ^2\beta _{\rm sh}^2\dot{M}_{-5}^{1/2} v_{w,3}^{-1/2} \,{\rm eV}, \end{eqnarray} \tag{ 1 }$$

where Z is the nuclear charge number of the particle and e is the charge of electrons. Note that, in the semi-relativistic hypernova scenario, the adiabatic cooling of UHE nuclei is much more efficient than other cooling processes, such as synchrotron cooling, photopion production, and photo-disintegration processes (Wang et al. 2008), so we do not consider them here.

Hypernova can also provide sufficient energy in UHECRs (Wang et al. 2007; Liu et al. 2011). According to Katz et al. (2009), the energy production rate³ in UHECRs per logarithmic energy interval inferred from the measured flux by PAO is ∼10^43.5 erg Mpc⁻³ yr⁻¹. Hypernova ejecta with Γβ ∼ 0.5 can produce cosmic rays with energies ∼10²⁰ eV and the kinetic energy in such velocity ejecta is typically ∼10⁵¹ erg (e.g., Wang et al. 2007; Chakraborti et al. 2011). Assuming that half of the kinetic energy goes into cosmic rays and that 10% of this energy is distributed in every logarithmic energy interval, the required local event rate of hypernovae is then ∼600 Gpc⁻³ yr⁻¹, which is consistent with the observed rate (e.g., Soderberg et al. 2006; Liang et al. 2007).

It is usually assumed that shock acceleration leads to a power-law energy spectrum for particles. The differential number of particles accelerated by one hypernova can be described by

$$\begin{equation} \frac{dN_i}{dE_i}=\frac{U_iE_i^{-p}}{\int _{E_{i,{\rm min}}}^{E_{i,{\rm max}}}dE_i E_i^{1-p}} {\rm exp}\left(-\frac{E_i}{E_{i,\rm max}}\right){\rm exp}\left(-\frac{E_{i,\rm min}}{E_i}\right), \end{equation} \tag{ 2 }$$

where the subscription i represents the species of the nuclei (e.g., i = H, He, and etc.), E_{i, max} and E_{i, min} are respectively the maximum and minimum energy of cosmic rays injected into the intergalactic space by the hypernova, and U_i is the total energy of nuclei of species i released by one hypernova into the intergalactic space. According to Equation (1), hypernova remnants can accelerate iron nuclei to energies about ∼5 × 10²⁰ eV, so we fix E_{max, Fe} = 10^20.5(Z/26) eV in the following calculation. E_{i, min} is the low energy cutoff in the spectrum, since nuclei with energy below E_{i, min} ≈ Z_ieB_hostR_host cannot escape from the confinement by the magnetic fields in the host galaxies, where B_host is the typical magnetic field of the host galaxy and R_host is the scale height of the host galaxy. Taking our own galaxy as an approximation, we have E_{min, Fe} ≈ 10¹⁹(Z/26)B₋₆(h_z/0.3 kpc) eV, where h_z is the scale height of the Galaxy. For a single velocity ejecta,⁴ the power-law index is usually assumed to be p ≃ 2 for both non-relativistic shocks and semi-relativistic shocks (Kirk & Schneider 1987). In the following calculation, we will assume the spectrum to be power law with exponential cutoffs above E_max and below E_min.

Given that both E_{i, max} and E_{i, min} are rigidity dependent (i.e., E_{i, max(min)}∝Z_iE_{H, max(min)}), the total mass of nuclei of species i released into the intergalactic space by one hypernova is

$$\begin{equation} M_{i,\rm CR}=A_i\int _{E_{i,\rm min}}^{E_{i,\rm max}}\frac{dN_i}{dE_i}dE_i\propto \frac{A_i}{Z_i}U_i. \end{equation} \tag{ 3 }$$

Approximating Z_i = A_i/2 (except for H) and assuming that the same proportions of particles of each species get accelerated (i.e., the value of M_{i, CR}/M_i is independent of i), we have

$$\begin{equation} U_{\rm H}:U_{i}:U_{j}= 2M_{\rm H}:M_{i}:M_{j}.\quad (i,j\ne H, i\ne j). \end{equation} \tag{ 4 }$$

With the above equation, one can transform the abundances of nuclei of each species to the abundances of these nuclei at certain given energy, as described by Equation (2).

In the hypernova scenario, cosmic-ray particles originate from the circum-stellar wind material or hypernova ejecta material. The W-R stellar wind contains abundant intermediate-mass elements such as C and O, which dominate the chemical enrichment of interstellar medium (ISM, e.g., Abbott 1982; Bieging 1990). The mass ratio of different elements for typical W-R stellar wind is M_He:M_C:M_O:M_X = 0.32:0.39:0.25:0.04 (Bieging 1990), derived based on observation of emission line of these elements in W-R stars and the stellar evolution model (see, e.g., Willis 1982; Prantzos et al. 1986; van der Hucht et al. 1986), where X denotes the elements other than He, C, and O. In the following calculation, we will use this composition as a representative case for the cosmic rays from the stellar wind.

On the other hand, the hypernova ejecta contains both intermediate-mass elements and heavy elements. It consists of products of nuclear reaction in the interior of the W-R star, such as C, O, Mg, Fe, as well as some heavier elements produced by the explosive nucleosynthesis. Chemical composition of the ejecta can be derived by modeling the spectra and light curves of hypernovae (e.g., Nakamura et al. 2001; Fryer et al. 2006) or through numerical simulation of SN explosion based on some specific stellar model (e.g., Georgy et al. 2009). Best fit for the early spectra and light curve of SN 1998bw is obtained in the model named "CO138E50" by Nakamura et al. (2001), where the input explosion kinetic energy and ejecta mass are, respectively, E_K = 5 × 10⁵² erg and M_ej = 10 M_☉. This model yields a chemical composition⁵ of M_C:M_O:M_Ne:M_Mg:M_Si:M_S:M_Ca:M_Fe= 0.006:0.71:0.037:0.034:0.083:0.041:0.007:0.09. In the following calculation, we will use this composition as a representative case for cosmic rays from the hypernova ejecta. In both scenarios, the amount of hydrogen is negligible, because the progenitor of Type Ib/c SN is stripped of H envelope, as indicated by lack of H line in the SN spectrum (Iwamoto et al. 1998; Nakamura et al. 2001). So during the explosion stage, there is little amount of H in both the wind and the SN ejecta.

While UHECRs are propagating in the intergalactic space, there are three main attenuation processes due to interaction with CMB or CIB photons: the photodisintegration process, the Bethe–Heitler process, and the photopion production process. The photodisintegration process causes nuclei to lose their nucleons and thus change their species, but does not reduce their Lorentz factor, while the latter two processes can reduce the Lorentz factors of nuclei. Thus, we need to consider the evolution of the mass number A(t) and Lorentz factor γ_N(t) of a nucleus with time jointly. In the hypernovae scenario, given that E_max/A ≲ 10¹⁹ eV, the corresponding maximum energy of each nucleon is below the threshold energy for photopion production with CMB photons even at high redshifts. Although these nuclei can interact with more energetic CIB photons, the photopion energy loss rate with such photons is much lower than the energy loss via the Bethe–Heitler process by interacting with much denser CMB photons, so we neglect the photopion energy loss in the following calculation. In this work, we adopt the CIB model suggested by Razzaque et al. (2009) and Finke et al. (2010), and use the online data of the CIB intensity at different redshifts.⁶

When an UHE nucleus is propagating in the intergalactic space, its Lorentz factor evolves with time as

$$\begin{equation} {-}\frac{d\gamma _N(t)}{dt}=\gamma _N(t)H(z)+\dot{\gamma }_{N,\rm BH}(t), \end{equation} \tag{ 5 }$$

where γ_N(t)H(z) represents the adiabatic energy losses due to cosmological expansion and $\dot{\gamma }_{N,\rm BH}(t)$ is energy loss rate due to the Bethe–Heitler process. Here $H(z)=H_0\sqrt{\Omega _m(1+z)^3+\Omega _{\Lambda }}$ is the Hubble constant at time t = t(z), where H₀ = 71 km s⁻¹ Mpc⁻¹, Ω_m = 0.27, and Ω_Λ = 0.73. In an isotropic photon background, the Bethe–Heitler energy loss rate of a nucleus of the Lorentz factor γ_i is given by

$$\begin{eqnarray} \dot{\gamma }_{N,\rm BH}(t)&=&\frac{Z^2}{A}\dot{\gamma }_{p,\rm BH}(t)=\frac{Z^2}{A}\frac{c}{2\gamma _N}\int _{\epsilon _{\rm th}}^{\infty }d\epsilon _{\gamma }\sigma _{\rm BH}(\epsilon _{\gamma })f_{\rm BH}(\epsilon _{\gamma })\epsilon _{\gamma }\nonumber\\ &&\times\,\int _{\epsilon _{\gamma }/2\gamma _i}^{\infty }d\epsilon \frac{n_{\gamma }(\epsilon,z)}{\epsilon ^2}, \end{eqnarray} \tag{ 6 }$$

where $\dot{\gamma }_{p,\rm BH}(t)$ is the Bethe–Heitler energy loss rate for a proton of the same Lorentz factor γ_N, Z is the nuclear charge of the nucleus, σ_BH and f_BH are the cross section and the fraction of energy loss in one interaction for the Beth–Heitler process (Chodorowski et al. 1992), and _th is the threshold energy for the interaction. Here _γ is the photon energy in the rest frame of the nucleus while is the photon energy in the lab frame. n_γ(, z) is the number density of background (CMB+CIB) photons of energy at redshift z.

The nucleus suffers from loss of nucleons by the photo-disintegration process, so the mass number evolves as

$$\begin{equation} {-}\frac{dA}{dt}=R_A(t,\gamma _i), \end{equation} \tag{ 7 }$$

where R_A is the reaction rate, which is

$$\begin{equation} R_A(t,\gamma _N)=\frac{c}{2\gamma _N^2}\int _{\epsilon _{\rm th}}^{\infty }d\epsilon \sigma _{{\rm dis},A} (\epsilon)\epsilon \int _{\epsilon /2\gamma _i}^{\infty }d\epsilon _{\gamma }\frac{n_{\gamma }(\epsilon _{\gamma },z)}{\epsilon _{\gamma }^2}, \end{equation} \tag{ 8 }$$

where σ_{dis, A} is the total photo-disintegration cross section for a nucleus of mass number A. The photo-disintegration cross section is dominated by the giant dipole resonance (GDR) up to 30 MeV, with the threshold energy between about 10 MeV and 20 MeV for all nuclei (in the nucleus rest frame). The cross section in this energy range can be modeled by a Gaussian form (e.g., Puget et al. 1976; Stecker & Salamon 1999; Hooper et al. 2008) or a Lorentzian form (e.g., Khan et al. 2005; Anchordoqui et al. 2007; Hooper et al. 2007). From 30 MeV to 150 MeV, the quasi-deuteron process becomes dominant and the cross section can be approximated as a plateau (e.g., Puget et al. 1976; Stecker & Salamon 1999; Ahlers & Taylor 2010). In this work, we use the tabulated cross section data generated by the code TALYS⁷(Goriely et al. 2008) which considered all the individual nucleon emission channels for nuclei with A ⩾ 5. For nuclei with A < 5, we adopt the Gaussian form near the threshold along with a plateau at higher energies to describe the total cross section.

We use the Runge–Kutta method to solve Equations (5) and (7) jointly so that we can trace the evolution history of the Lorentz factor γ_N(z; z_s, E_s, A_s) and mass number A(z; z_s, E_s, A_s) of a specific nucleus injected from the source at redshift z = z_s, where E_s and A_s are the initial energy and mass number of the nucleus. The evolution of the mass number is assumed to develop along the Puget–Stecker–Bredekamp chain (see e.g., Puget et al. 1976; Stecker & Salamon 1999; Ahlers & Taylor 2010). We also record the redshifts, Lorentz factors, and number of secondary protons⁸ that are disintegrated from the parent nuclei so that the contribution of secondary protons to cosmic-ray spectrum can be properly included. The nuclei are injected from the maximum redshift z_max = 6 to the minimum redshift z_min = 0.001 (corresponding to a distance D = 4 Mpc). The number density of hypernovae at redshift z is assumed to follow the star formation history (SFH), given by (Hopkins & Beacom 2006; Yüksel et al. 2008)

$$\begin{equation} n(z)\propto \left\lbrace \begin{array}{@{}lll}(1+z)^{3.4}, \quad z<1 \\\\ (1+z)^{-0.3}, \quad 1<z<4 \\\\ (1+z)^{-3.5}, \quad z>4. \end{array} \right. \end{equation} \tag{ 9 }$$

We collect all the nuclei and secondary products that arrive at the Earth (z = 0) and put them into the corresponding energy bin and species bin according to their Lorentz factors γ(z = 0) and mass numbers A(z = 0). Then we can get the energy spectrum of cosmic rays in each species bin or get the all-particle energy spectrum by adding up all species particles at a certain energy.

The cosmic-ray composition at the sources strongly affect the final composition and energy spectrum of cosmic rays arriving at the Earth. First, we consider the chemical composition at the source equal to that of the W-R stellar wind, where the mass ratio among dominant elements is M_He: M_C: M_O = 0.32: 0.39: 0.25 with negligible H. The initial spectrum at the source is assumed to be a power law with p = 2. The result of the final spectrum and composition are presented in Figure 1. One can see that with the typical W-R stellar wind composition, the cosmic-ray flux drops too fast at high energies to account for the observations by PAO. The apparent reason for the fast drop is that the low E_max for intermediate-mass nuclei such as C and O results in an exponentially cutoff at energy ≲ 10²⁰ eV. However, the situation will not be ameliorated even if we raise E_{max, O} to higher energies, e.g., to 3 × 10²⁰ eV. This is because that lower energy background photons, which are more abundant, are involved in the photodisintegration interaction for intermediate-mass nuclei due to their higher Lorentz factors, compared with heavy nuclei of the same energy. Therefore the attenuation lengths for the intermediate nuclei with energies larger than 10²⁰ eV are so small (e.g., Puget et al. 1976; Allard et al. 2006) that even those who emitted by the nearest sources at z_min = 0.001 (i.e., 4 Mpc) will be effectively attenuated by the background photons and hence the flux of cosmic rays at such energies is severely suppressed after propagation. Only if E_{max, O} is extremely high (e.g., ∼3 × 10²¹ eV) so that secondary protons of energies up to (2–3) × 10²⁰ eV can be produced in photodisintegration process, and then the spectrum becomes flat at the highest energy end (Allard et al. 2007). However, such a high E_{max, O} can hardly be achieved in the hypernova scenario.

Zoom In Zoom Out Reset image size

Presence of heavier nuclei in the source composition would be beneficial to fit the observation data, since heavy nuclei have higher steepening energies due to interactions with background photons. Hypernova ejecta produced after the explosion can provide heavy elements besides intermediate-mass elements. Figure 2 shows the propagated spectrum for cosmic-ray composition at the sources equal to the hypernova ejecta composition in the "CO138E50" model for SN 1998bw (Nakamura et al. 2001). With the contribution by heavy nuclei at the sources, the propagated spectrum is in better agreement with the PAO data than the wind composition case.

Zoom In Zoom Out Reset image size

The fit can be further improved if some non-standard effects are taken into account. In panel (a) of Figure 3, we increase the initial iron abundance by a factor of three, i.e., iron nuclei constitute a fraction of ∼25% of the total mass at the sources. This could be achieved since the enrichment of heavy elements is not uniformly distributed spatially, as already seen in some Galactic supernova remnants (e.g., Hwang et al. 2000). A larger explosion energy would also lead to more heavy elements synthesized during the explosion (e.g., Metzger et al. 2011). As iron nuclei have larger steepening energy in the spectrum due to interactions with background photons, a harder spectrum is expected at the highest energy if more iron nuclei are present. This is consistent with earlier results in Allard et al. (2008). Another possible effect is a hard injection CR spectrum. Although in conventional shock acceleration theory, the power-law index is p ≃ 2, some observations (e.g., Horns & Aharonian 2004; Green 2009; Reynolds et al. 2011; Abdo et al. 2011) and theoretical calculations (e.g., Ellison et al. 1996; Malkov 1997; Vainio & Schlickeiser 1999; Vainio et al. 2003; Tammi & Duffy 2009) have suggested harder cosmic-ray spectra with indices p < 2. In panel (b), we show the propagated spectrum for an injection spectrum with power-law index p = 1.6. One can see that the fit to the data gets better than the case of p = 2. Besides, the source density distribution in the nearby universe could be nonuniform. Since our Galaxy locates inside the Local Group and Local Supercluster, we examine the effect of local overdensity of sources on the propagated spectrum. Panel (c) shows the spectrum for a local overdensity in the source number by a factor of two relative to the average within the size of 30 Mpc (Blanton et al. 2001). Due to the attenuation by background photons, higher energy cosmic rays are mainly contributed by closer sources. Thus a local excess of the source number density increases the flux of higher energy cosmic rays and thus hardens the propagated spectrum at the highest energy end.

Zoom In Zoom Out Reset image size

We should note that, although the deviation between the propagated spectrum and observational data is relatively large in the wind composition scenario, the situation can also be improved if some of the above effects are taken into account. Panel (d) of Figure 3 presents the result after considering a harder injection cosmic-ray spectrum of p = −1.6 as well as a local overdensity in the source number by a factor of two. The theoretical energy spectrum agrees reasonably well with the observational data.

In Figure 4, we present how the average mass number of cosmic rays evolves with energy for different scenarios discussed above. Note that the starting energy of this figure is set to 3 × 10¹⁸ eV since our model does not account for the CR spectrum data below this energy, for which a Galactic component contribution is needed. One can see that, for the hypernova ejecta composition scenarios, the average mass number increases gradually with energy, which is consistent with the finding by PAO that the composition of UHECRs becomes increasingly heavy with energy. In the wind composition scenarios, the average mass number increases more slowly with energy. However, since the hadronic interaction models at such high energies are not well understood, such a composition is still consistent with the measurements within the uncertainties of theoretical expectations.

Zoom In Zoom Out Reset image size