<figure><img src="image_1.png" /></figure>

HAL
open science

# Non-linear acceleration at supernova remnant shocks and the hardening in the cosmic ray spectrum

S. Recchia, S. Gabici

► To cite this version:

S. Recchia, S. Gabici. Non-linear acceleration at supernova remnant shocks and the hardening in the cosmic ray spectrum. Monthly Notices of the Royal Astronomical Society: Letters, 2018, 474, pp.L42-L46. 10.1093/mnrasl/slx191. hal-01703765

HAL Id: hal-01703765

https://hal.science/hal-01703765

Submitted on 3 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Non-linear acceleration at supernova remnant shocks and the hardening in the cosmic ray spectrum

S. Recchia\* and S. Gabici\*

APC, Université Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne Paris Cité, France

Accepted 2017 November 22. Received 2017 November 22; in original form 2017 October 3

## ABSTRACT

In the last few years, several experiments have shown that the cosmic ray spectrum below the knee is not a perfect power law. In particular, the proton and helium spectra show a spectral hardening by ~0.1–0.2 in spectral index at particle energies of ~200–300 GeV nucleon⁻¹. Moreover, the helium spectrum is found to be harder than that of protons by ~0.1 and some evidence for a similar hardening was also found in the spectra of heavier elements. Here, we consider the possibility that the hardening may be the result of a dispersion in the slope of the spectrum of cosmic rays accelerated at supernova remnant shocks. Such a dispersion is indeed expected within the framework of non-linear theories of diffusive shock acceleration, which predict steeper (harder) particle spectra for larger (smaller) cosmic ray acceleration efficiencies.

**Key words:** acceleration of particles – cosmic rays – supernova remnants.

## 1 INTRODUCTION

In the standard picture of the origin of cosmic rays (CRs) the observed flux, at least below the energy of the 'knee' ($E_{\text{knee}} \approx 3$ PeV; see e.g. Höorandel 2006), is thought to be produced in the Galactic disc at supernova remnant (SNR) shocks through diffusive shock acceleration (DSA) (orgreenMalkov & Drury 2001). After leaving their sources, CRs are believed to propagate diffusively through the interstellar medium (ISM), and eventually escape from the Galaxy (see e.g. Blasi 2013). The observed CR spectrum below the knee resembles a power law in energy $\approx E^{-2.7}$ and can be roughly accounted for if one assumes that both the slope of the injection spectrum and the CR diffusion coefficient are power laws in energy, with slopes $-\gamma$ and $\delta$, respectively. Various observational constraints provide $\gamma + \delta \approx 2.7$ and $\delta \approx 0.3–0.6$ (see e.g. Strong & Moskalenko 1998; Evoli et al. 2008).

However, in the last few years much evidence points to a more complex scenario: most notably several experiments such as ATIC-2 (Panov et al. 2009), PAMELA (Adriani et al. 2011), CREAM (Yoon et al. 2011) and AMS-02 (Aguilar et al. 2015) found a spectral hardening in the proton and helium spectra at particle energies 200–300 GeV nucleon⁻¹. Moreover, the helium spectrum is found to be harder than that of protons by ~0.1 in spectral index. Some evidence for a similar hardening was also found in the spectra of heavier elements (Maestro et al. 2010). The PAMELA data (Adriani et al. 2011) suggest that the slope of protons changes at ~230 GeV from $\gamma_1 \sim 2.85$ to $\gamma_2 \sim 2.67$. Instead, the AMS-02 data found the break at ~335 GeV and a slope change of ~0.13.

This spectral feature is still not understood and several explanations have been put forward in which the spectral hardening is interpreted as: the result of a break in the CR diffusion coefficient (Tomassetti 2012; Aloisio & Blasi 2013; Genolini et al. 2017); the effect of a nearby source (Thoudam & Hörandel 2012); the consequence of a dispersion in spectral index at the sources (Yuan, Zhang & Bi 2011); the result of the spectral concavity due to the non-linear nature of DSA and of the contribution of particle acceleration at the reverse shock (Ptuskin, Zirakashvili & Seo 2013); the possible presence of distinct populations of CR sources (Zatsepin & Sokolskaya 2006); the consequence of a break in the energy loss rate (Krakau & Schlickeiser 2015).

In this paper, we suggest that the spectral hardening may be a natural prediction of the non-linear theory of DSA. It would result from the interplay of the efficient magnetic field amplification at SNR shocks and of the CR Alfvénic drift in the upstream region. Following Caprioli (2012), we will show that these two effects may result in a dispersion in the CR spectral slope at the sources which may lead to the observed hardening.

The reasons for considering such scenario are manifold: first of all, the presence of efficient magnetic field amplification has been detected in several SNRs (see e.g Ballet 2006; Vink 2012) and is widely considered as a crucial ingredient for the acceleration of CRs to the energy of the knee (see e.g. Bell 2004). Secondly, recent observations of γ-rays in SNRs both in the GeV (see e.g. Abdo et al. 2011) and in the TeV band (see e.g. Acciari et al. 2011) show that there may be a quite large dispersion in the slope of the CR spectrum at SNRs. In fact, in the cases in which the observed γ-ray emission is likely of hadronic origin (see e.g. Morlino & Caprioli 2012), the inferred CR spectrum shows a quite large dispersion in the spectral index, $\propto E^{-2.1} - E^{-2.5}$. Note that, based on this

<footnote>\* E-mail: recchia@apc.in2p3.fr (SR); gabici@apc.in2p3.fr (SG)</footnote>
observation, Yuan et al. (2011) made a fit of the CR spectrum with the assumption that the source spectrum has a slope with some dis- tribution between ≈2.1 and 2.5 and showed that the hardening can be reproduced. Here, we provide a possible physical explanation of the observed dispersion.

It is interesting to note that such spectra are significantly steeper than the universal spectrum ∝ $E^{-2.0}$ predicted by the *linear* (test particle) theory of DSA at SNR shocks. At first sight, the disagreement seems to be even larger if one considers *non-linear* theories of DSA (NLDSA), which account for the reaction of CRs on the shock dynamics (see e.g. Malkov & Drury 2001). In this context, the pressure exerted by CRs on to the upstream fluid induces the formation of a shock precursor, which in turn makes the CR spectrum at the shock concave, namely, steeper than ∝ $E^{-2.0}$ at low energies and harder at high energies. Moreover, the more efficient the CR acceleration the more evident becomes the concavity and for large efficiencies the spectrum above a few GeV becomes as flat as ∝ $E^{-1.5}$, clearly at odds with the γ-ray observations of SNRs reported above. Note that Ptuskin et al. (2013) invoked this concavity in order to explain the CR spectral hardening

For this reason, in the last few years a number of works have been devoted to the study of possible ways to reconcile the predictions of NLDSA theories with observations. In particular, it has been proposed that taking into account the velocity of the CR scattering centres (see e.g. Zirakashvili & Ptuskin 2008; Caprioli 2012) and the effect of the amplified magnetic field at the shock (see e.g. Vainio & Schlickeiser 1999; Caprioli 2012) the CR source spectrum may become significantly steeper than ∝ $E^{-2.0}$. In this paper, we focus on the work by Caprioli 2012, where a theory of NLDSA is developed, which includes both the magnetic field amplification, due to CR streaming instability, and the resulting enhanced velocity of the CR scattering centres (Alfvén waves propagate faster for larger values of the ambient magnetic field (see e.g. Zirakashvili & Ptuskin 2008, Caprioli 2012).

The main findings of Caprioli (2012) can be summarized as follows:

(i) the magnetic field amplification acts as a self-regulating mechanism of the acceleration process. The pressure of the amplified magnetic field makes the shock compression factor smaller than 4 (which is the strong shocks limit of the test particle regime of DSA, where the field is not amplified). As a consequence of that, the maximum acceleration efficiency for strong shocks turns out to be ~30 per cent and the shock modification induced by the CR pressure is quite modest. This is quite at odds with earlier formulations of NLDSA, in which the CR acceleration efficiency can reach values well above ~30 per cent and the compression factor can be well above 4 (see e.g. Blasi, Gabici & Vannoni 2005)

(ii) the spectrum at energies above a few GeV resembles quite well a power law whose spectral slope remains virtually constant, together with the acceleration efficiency, for a large part of the SNR lifetime (up to ~ 20 000-30 000 yr);

(iii) the combined effect of the magnetic field amplification and of the Alfvénic drift in the upstream region makes the spectrum steeper than $E^{-2.0}$, with slopes in the range 2.1–2.6 in agreement with the γ-ray observations of SNRs. In addition to that, the more efficient is the CR acceleration the steeper is the spectrum, contrary to the standard predictions of NLDSA, in which an efficient acceleration leads to harder spectra.

For a more extended discussion on this approach and its limits of validity, the reader is referred to Caprioli (2012).

Based on the results summarized above, in the following we treat the acceleration efficiency at SNR shocks as a free parameter in the range $\xi_{CR} \approx 0.03-0.3$, and we assume that the CR spectrum at SNR shock is a power law of slope $\gamma_{CR}$. We then compute the spectral slope as a function of the acceleration efficiency $\gamma_{CR}(\xi_{CR})$ by taking into account the magnetic field amplification by CR streaming instability and the effect of the velocity of the self-generated Alfvén waves, which act as scattering centres. We find that the dispersion in the acceleration efficiency induces a dispersion in the spectral slope, with steeper (harder) spectra corresponding to larger (lower) acceleration efficiencies. Finally, taking into account the dispersion in the spectral slope and the relation between acceleration efficiency and slope, the observed proton and helium spectral hardening at 200–300 GeV nucleon$^{-1}$ can be accounted for in a quite natural way.

The paper is organized as follows: in Section 2 we illustrate the calculation of the compression factor felt by CRs as a function of the CR acceleration efficiency and we show the resulting dispersion in slope. In Section 3, we use these results to estimate the observed proton and helium spectrum with the additional assumption that CRs propagate diffusively in the Galaxy with a power law diffusion coefficient (whose slope and normalization is chosen in order to fit the data) and we compare the obtained spectrum with the data. Finally, we conclude in Section 4.

## 2 MAGNETIC FIELD AMPLIFICATION AND THE SPECTRAL SLOPE

Keeping in mind the results of the previous section, here we illustrate a simple calculation which allows us to quickly estimate the slope of the CR spectral slope under the following assumptions:

(i) the CR acceleration efficiency $\xi_{CR}$ is an input parameter of the problem, and is in the range $\xi_{CR} \sim 0.03-0.3$. Thus, the shock modification is modest and the CR spectrum is nearly a perfect power law;

(ii) the magnetic field is amplified by the CR streaming instability and is assumed to be the same in the whole upstream region;

(iii) the Alfvén waves excited in the upstream propagate against the fluid at velocity $v_{A1}$, which is computed in the amplified magnetic field, while they are assumed to be isotropized in the downstream region, giving $v_{A2} \sim 0$ (the subscripts 1 and 2 refer to quantities calculated in the upstream and downstream region, respectively).

The slope of the CR spectrum depends on the effective compression factor felt by CRs, which in turn depends on the CR acceleration efficiency $\xi_{CR} \equiv P_{CR}/\rho_1 u_1^2$, on the fluid and Alfvénic Mach numbers upstream ($M_1^2 \equiv \rho_1 u_1^2 / \gamma P_g$ and $M_A^2 \equiv \rho_1 u_1^2 / 2P_w$, respectively) and on the jump conditions at the shock, where $\rho_1, u_1$, and $P_g$ are the upstream gas density, velocity and pressure, $P_{CR}$ is the CR pressure at the shock and $P_w$ is the pressure of the amplified (upstream) magnetic field. Here, we neglect the modest shock modification induced by the CR pressure in the upstream fluid, which implies that all the relevant physical quantities characterizing the fluid do not depend on the location upstream. Finally, we assume that the gas can be described by an adiabatic equation of state. Using equation (2.22) in Caprioli (2012), one can evaluate the effect of the magnetic field amplification due to CR streaming instability, and estimate the
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Effective compression factor felt by CRs at the shock as a function of the upstream fluid Mach number $M_1$ and for three different values of the CR acceleration efficiency $\xi_{CR}$.</figcaption></figure>

The Mach number is $M_1 = 100$.

$$M_A^{-2} = \frac{4}{25} \frac{\left[1 - (1 - \xi_{CR})^{\frac{5}{4}}\right]^2}{(1 - \xi_{CR})^{\frac{3}{4}}}. \qquad (1)$$

Another crucial parameter is the fluid compression factor $R = u_1/u_2$, which can be computed after taking into account the pressure of the amplified magnetic field. Following a procedure similar to that presented in Vainio & Schlickeiser (1999) and Caprioli (2012), we get:

$$\frac{M_1^2 \frac{\gamma+1}{R} - (\gamma-1)}{2(1 + \Lambda_B)} \approx 1, \quad \text{where} \qquad (2)$$

$$\Lambda_B = W \left[ 1 + R \left( \frac{2}{\gamma} - 1 \right) \right] \quad \text{and} \quad W = \frac{\gamma}{2} \frac{M_1^2}{M_A^2}.$$

The effective compression felt by CRs differs from $R$, because in the upstream region the Alfvén waves generated by the CR streaming instability propagate in the direction opposite to the fluid. Since CRs are coupled to waves through scattering, the effective advection velocity they experience in the upstream region is not $u_1$, but rather $u_1 - v_{A1}$. Thus, the effective compression factor felt by CRs is

$$R_{\text{eff}} = \frac{u_1 - v_{A1}}{u_2} = R \left( 1 - \frac{1}{M_A} \right). \qquad (3)$$

The CR spectral slope can then be estimated as (see e.g. Malkov & Drury (2001)):

$$\gamma_{CR} \sim \frac{3 R_{\text{eff}}}{R_{\text{eff}} - 1}. \qquad (4)$$

It is important to stress that only two physical parameters regulate the system: the fluid upstream Mach number $M_1$ and the CR acceleration efficiency $\xi_{CR}$.

In Fig. 1, we show the dependence of the effective compression ratio $R_{\text{eff}}$ on the Mach number $M_1$ for three different values of $\xi_{CR}$. For $M_1 \gtrsim 10$, the effective compression ratio (the same result also holds for the fluid compression ratio and for the spectral slope) is virtually independent on $M_1$. This implies that the slope of the spectrum of accelerated particles does not change for most of the SNR lifetime.

On the other hand, both the fluid and effective compression ratios, and thus also the spectral slope, strongly depend on the CR

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Fluid compression factor $R$, effective compression factor $R_{\text{eff}}$ and CR spectral slope $\gamma_{CR}$ as a function of the CR acceleration efficiency $\xi_{CR}$. The Mach number is $M_1 = 100$.</figcaption></figure>

acceleration efficiency. This is evident from Fig. 2, where one can see that, for $\xi_{CR}$ ranging in ~0.03–0.3, the slope ranges from ~4.1 to ~4.6. This result shows that the inclusion of the magnetic field amplification and of the Alfvénic drift in the calculation of the compression factor leads to quite steep source spectra, with larger $\xi_{CR}$ corresponding to steeper spectra. Note that within the present setup the fluid compression factor always remains $\lesssim 4$, while in the standard NLDSA theories compression ratios much larger than 4 are usually found and the CR acceleration efficiency can be well above ~30 per cent (see e.g. Malkov & Drury 2001).

### 3 COMPARISON WITH DATA

Here, we assume that SNR shocks accelerate CRs with an efficiency uniformly distributed in the range $\xi_{CR} \sim 0.03-0.3$, which implies, as shown in Section 2, a dispersion in the CR spectral slope, $\gamma_{CR}$, in the range ~4.1–4.6. Formally, this is the slope of the CR spectrum at the shock, and not that of the spectrum of particles escaping the SNR and injected in the ISM. However, under reasonable assumptions these two spectra are identical (see e.g. Gabici 2011).

The plausibility of assuming different acceleration efficiencies resides in the fact that the microphysics of particle injection in the acceleration process from the thermal plasma depends on several environmental parameters, such as the background magnetic field and the shock obliquity, the plasma temperature and the Mach number (see e.g. Blasi et al. 2005; Caprioli et al. 2015).

After escaping SNRs, CRs are believed to propagate diffusively in the Galaxy with a diffusion coefficient $D(R) = D_0(R/GV)^\delta$ ($R$ is the particle rigidity). The values of $D_0$ and $\delta$ are chosen in order to fit the observed proton spectrum in the energy range 40 GeV–10 TeV, namely around the spectral hardening at 200–300 GeV nucleon$^{-1}$. As for helium, we used the same injection spectral slopes and diffusion coefficient of protons, but we also took into account spallation. The proton and helium spectra below ~40 GeV nucleon$^{-1}$ are not considered since at these energies both the solar modulation and possible advection effects are important (see e.g. Aloisio, Blasi & Serpico 2015), which were not included in our calculation.

Under these assumptions the proton spectrum ($E$ is the particle energy) can be written as (see e.g. Blasi & Amato 2012):

$$f_{\text{p}}(E) = \int_{\xi_m}^{\xi_M} \frac{R_{\text{SN}}}{\pi R_d^2} \frac{H}{2D(E)} g_{\text{p}}(E) \frac{d\xi_{\text{CR}}}{\xi_M - \xi_m}, \qquad (5)$$
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Proton flux compared with the PAMELA, AMS-02 and CREAM data.</figcaption></figure>

where

$$g_p(E) = \frac{\xi_{CR} E_{SN}}{I(\gamma_{CR})(mc^2)^2} \left(\frac{E}{mc^2}\right)^{-\gamma_{CR}+2}. \quad (6)$$

$I(\gamma_{CR}) = \int_{x_0}^{\infty} dx \, x^{2-\gamma_{CR}} [\sqrt{1+x^2} - 1]$ is a normalization factor chosen in such a way that $\int_{E_0}^{\infty} g_p(E) E_k dE = \xi_{CR} E_{SN}$, where $x \equiv E/mc^2$ and $E_k$ is the particle kinetic energy.

The helium spectrum is given by

$$f_{He}(E) = \int_{\xi_{in}}^{\xi_{M}} \frac{R_{SN}}{\pi R_{d}^{2}} \frac{H}{2D(E)} g_{He}(E) \frac{1}{1 + \frac{h H n_{d} c \sigma_{sp}}{D(E)}} \frac{d\xi_{CR}}{\xi_{M} - \xi_{m}}, \quad (7)$$

where

$$g_{He}(E) = \eta_{He} \frac{\xi_{CR} E_{SN}}{I(\gamma_{CR})(mc^2)^2} \left(\frac{E}{mc^2}\right)^{-\gamma_{CR}+2}. \quad (8)$$

Here, $R_{SN} \approx 1/30$ yr is the SN explosion rate in the Galaxy, $R_d \approx 15$ kpc is the Galactic disc radius, $h \approx 250$ pc is the Galactic disc height, $H \approx 4$ kpc is the Galactic halo size, $n_d \approx 5 \text{ cm}^{-3}$ is the average gas density in the disc. $\eta_{He}$ is a factor chosen in such a way to reproduce the correct normalization of the helium spectrum. Finally, $\sigma_{sp}$ is the helium spallation cross-section (see e.g. Blasi & Amato 2012).

In Fig. 3, we show the proton flux as computed from equation (5) (red line) compared with the data by PAMELA (Adriani et al. 2011), by AMS-02 (Aguilar et al. 2015) and by CREAM (Yoon et al. 2011). The plot has been obtained with the diffusion coefficient parameters: $D_0 \sim 8 \times 10^{28} \text{ cm}^2 \text{ s}^{-1}$ and $\delta \sim 0.4$. The slope found for the diffusion coefficient is well within the observational constraints, namely $\delta \sim 0.3-0.6$. With this diffusion coefficient, the grammage traversed by CRs, namely $X = n_d h H c m_p / D(R)$, is $\approx 11 \text{ g cm}^{-2}$ at $10 \text{ GeV n}^{-1}$ (see e.g. Blasi 2013).

On the same figure, we also show the proton flux (green line) computed in the case of two distinct populations of CR sources, one with $\xi_{CR} \approx 0.03$ and the other with $\xi_{CR} \approx 0.3$. Also this plot has been obtained with a diffusion coefficient slope of $\delta = 0.4$, while the explosion rate of the population with the largest acceleration efficiency has been taken to be $\sim 3$ times smaller than that with the smallest efficiency. Note that taking into account such scenario could be motivated by a different behaviour of type I and II supernovae in the acceleration of CRs, due to the fact that the explosion happens in a different environment, denser and colder in the case of type I, in a hot rarefied bubble in the case of type II (see e.g. Ptuskin

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Helium flux compared with the PAMELA, AMS-02 and CREAM data.</figcaption></figure>

& Zirakashvili 2005; Zatsepin & Sokolskaya 2006). In Fig. 4, we show the same as in Fig. 3 for the helium flux.

Note that the dispersion in the CR acceleration efficiency, and the consequent dispersion in the CR spectral slope, naturally leads to a spectral hardening in the proton spectrum at $\lesssim \text{TeV energies and}$, overall, to a good agreement with the data. A similar hardening is found also in the helium spectrum. Moreover, in agreement with observations, this feature is less prominent in the helium spectrum compared to the proton spectrum and the helium spectrum is found to be harder than the one of protons. This is due to spallation, which hardens the spectrum, especially at lower energies (see also Blasi & Amato 2012). In the case of two distinct populations with different acceleration efficiencies the spectral hardening is also well reproduced, both in the proton and helium spectrum. However, in this case the spectral feature appears to be sharper (see e.g. Genolini et al. 2017) than in the case of uniformly distributed efficiency.

Finally, when comparing our results with data, one has to keep in mind that the AMS-02 and CREAM data for helium at $\sim 1$ TeV nucleon$^{-1}$ differs by $\sim 20-30$ per cent, making it impossible to obtain an equally accurate fit to both data sets.

## 4 DISCUSSION AND CONCLUSIONS

The magnetic field amplification at SNR shocks, which is thought to be necessary in order to accelerate CRs up to PeV energies, may also act as a feedback process which limits the maximum achievable CR acceleration efficiency to $\sim 30$ per cent, thus keeping the overall shock modification modest. Moreover, together with the Alfvénic drift, the magnetic field amplification leads to quite steep CR source spectra (spectral slope in momentum $\sim 4.1-4.6$), in agreement with the CR spectra in SNRs inferred from $\gamma$-ray observations (see e.g. Caprioli 2012).

In this paper, we studied the acceleration of CRs at SNR shocks under the following realistic assumptions:

(i) the CR acceleration efficiency it may vary within the range $\sim 0.03-0.3$;

(ii) the shock modification induced by the CR pressure at the shock is modest;

(iii) the magnetic field is significantly amplified by CR streaming instability and the Alfvén speed (computed in the amplified field) upstream of the shock is enhanced accordingly.

We showed that the dispersion in the CR acceleration efficiency produces a dispersion in the shock compression factor. This in turn
results in a dispersion in the CR spectral slope, with steeper spectra corresponding to larger acceleration efficiencies.

This result has then be used to demonstrate that, by assuming a diffusive propagation of CRs in the Galaxy with a spatial independent diffusion coefficient, the above mentioned dispersion in the slope of the injection spectrum can account in a quite natural way for the spectral hardening found in the proton and helium spectrum in the energy range ~200–300 GeV nucleon. Moreover, in agreement with observations, because of spallation the helium spectrum is found to be harder than the proton spectrum (even if their injection spectra are identical) and the helium spectral hardening is less prominent than that of protons.

In the scenario considered in this paper, the acceleration efficiency does not depend on the age of the remnant, as predicted by non-linear models of DSA (Caprioli 2012). A time dependent acceleration efficiency, due for instance to a change in the environment sampled by the shock, would somewhat affect our results. A possible observational test of the validity of this assumption might come from the observation of the γ-ray hadronic emission from an extended sample of individual SNRs. This probably goes beyond the capabilities of current γ-ray instruments, due to the limited number of detections (see e.g. Cristofari et al. 2013). However, the number of SNRs detected in TeV γ-rays will dramatically increase with the advent of next generation instruments such as the Cerenkov Telescope Array (see e.g. Cristofari et al. (2017)). The increase in statistics will most likely provide stringent constraints to our model.

## ACKNOWLEDGEMENTS

This work has been financially supported by the region Île-de-France under the Domaine d' Intérêt Majeur en Astrophysique et Conditions d' Apparition de la Vie (DIM-ACAV) programme, the Observatory of Paris (Action Fédératrice CTA), and the Programme National Hautes Energies (PNHE) funded by CNRS/INSU-IN2P3, CEA and CNES, France.

## REFERENCES

* Abdo A. A. et al., 2011, ApJ, 734, 28

* Acciari V. A. et al., 2011, ApJ, 730, L20

* Adriani O. et al., 2011, Science, 332, 69

* Aguilar M. et al., 2015, Phys. Rev. Lett., 114, 171103

* Aloisio R., Blasi P., 2013, J. Cosmol. Astropart. Phys., 7, 001

* Aloisio R., Blasi P., Serpico P. D., 2015, A&A, 583, A95

* Ballet J., 2006, Adv. Space Res., 37, 1902

* Bell A. R., 2004, MNRAS, 353, 550

* Blasi P., 2013, A&AR, 21, 70

* Blasi P., Amato E., 2012, J. Cosmol. Astropart. Phys., 1, 010

* Blasi P., Gabici S., Vannoni G., 2005, MNRAS, 361, 907

* Caprioli D., 2012, J. Cosmol. Astropart. Phys., 7, 038

* Caprioli D., Pop A.-R., Spitkovsky A., 2015, ApJ, 798, L28

* Cristofari P., Gabici S., Casanova S., Terrier R., Parizot E., 2013, MNRAS, 434, 2748

* Cristofari P., Gabici S., Humensky T. B., Santander M., Terrier R., Parizot E., Casanova S., 2017, MNRAS, 471, 201

* Evoli C., Gaggero D., Grasso D., Maccione L., 2008, J. Cosmol. Astropart. Phys., 10, 018

* Gabici S., 2011, Mem. Soc. Astron. Italiana, 82, 760

* Genolini Y. et al., 2017, preprint (arXiv:1706.09812)

* Höorandel J. R., 2006, J. Phys.: Conf. Ser., 47, 41

* Krakau S., Schlickeiser R., 2015, ApJ, 811, 11

* Maestro P. et al., 2010, preprint (arXiv:1003.5759)

* Malkov M. A., Drury L. O., 2001, Rep. Prog. Phys., 64, 429

* Morlino G., Caprioli D., 2012, A&A, 538, A81

* Panov A. D. et al., 2009, Bull. Russ. Acad. Sci.: Phys., 73, 564

* Ptuskin V. S., Zirakashvili V. N., 2005, A&A, 429, 755

* Ptuskin V., Zirakashvili V., Seo E.-S., 2013, ApJ, 763, 47

* Strong A. W., Moskalenko I. V., 1998, ApJ, 509, 212

* Thoudam S., Hörandel J. R., 2012, MNRAS, 421, 1209

* Tomassetti N., 2012, ApJ, 752, L13

* Vainio R., Schlickeiser R., 1999, A&A, 343, 303

* Vink J., 2012, A&AR, 20, 1

* Yoon Y. S. et al., 2011, ApJ, 728, 122

* Yuan Q., Zhang B., Bi X.-J., 2011, Phys. Rev. D, 84, 043002

* Zatsepin V. I., Sokolskaya N. V., 2006, A&A, 458, 1

* Zirakashvili V. N., Ptuskin V. S., 2008, in Aharonian F. A., Hofmann W., Rieger F., eds, AIP Conf. Series Vol. 1085, p. 336 preprint (arXiv:0807.2754)

This paper has been typeset from a *TEX/LATEX* file prepared by the author.