On particle acceleration and transport in plasmas in the Galaxy: theory and observations

Accelerated particles are ubiquitous in the Cosmos and play a fundamental role in many processes governing the evolution of the Universe at all scales, from the sub-AU scale relevant for the formation and evolution of stars and planets to the Mpc scale involved in Galaxy assembly. We reveal the presence of energetic particles in many classes of astrophysical sources thanks to their production of non-thermal radiation, and we detect them directly at the Earth as cosmic rays. In the last two decades both direct and indirect observations have provided us a wealth of new, high-quality data about cosmic rays and their interactions both in sources and during propagation, in the Galaxy and in the Solar System. Some of the new data have confirmed existing theories about particle acceleration and propagation and their interplay with the environment in which they occur. Some others have brought about interesting surprises, whose interpretation is not straightforward within the standard framework and may require a change of paradigm in terms of our ideas about the origin of cosmic rays of different species or in different energy ranges. In this article, we focus on cosmic rays of galactic origin, namely with energies below a few petaelectronvolts, where a steepening is observed in the spectrum of energetic particles detected at the Earth. We review the recent observational findings and the current status of the theory about the origin and propagation of galactic cosmic rays.


Introduction
Cosmic Rays have been known and studied for more than a century now (see e.g. Amato (2014); Blasi (2013) for recent reviews). They are highly energetic charged particles, mainly protons and He nuclei, with a minor fraction of heavier nuclei (1%), electrons (2%) and anti-matter particles (positrons and anti-protons, 1 ‰). Their origin was associated with Supernova (SN) explosions in the Galaxy already in the 1930s and the mechanism by which they would be accelerated up to very high energies in the blast waves emerging from SN explosions was proposed in the late 1970s (see § 2).
This association, however, has not yet found direct proof, and, as we discuss below, some recent developments, in both theory and high energy astrophysical observations, have actually cast doubts on whether Supernova Remnants (SNRs) should be considered as the primary CR accelerators in the Galaxy throughout the entire energy range, up to PeV (10 15 eV) energies. At the same time, direct CR measurements have shown that particle transport is more complex than traditionally considered, and that additional sources are likely needed at least for CR positrons (see e.g. Amato & Blasi (2018) for a recent review, and Sec. 6 for more details).
In recent years it has become increasingly clear that understanding the processes that govern the acceleration and transport of energetic particles is a necessary key to unveil a number of unsettled questions in Astrophysics. Energetic particles -which we will often refer to as Cosmic Rays, in the following -play a role in many astrophysical systems. Just to mention a few, from the smallest to the largest scales, CRs affect the behaviour of planetary magnetospheres (Grießmeier et al. 2015), they are believed to be the main ionizing agents penetrating dense gas clouds and regulating star formation (Padovani et al. 2020), they are suspected to be the primary drivers of galactic winds and the rulers of galactic feedback (Buck et al. 2020), and might even have a role in the generation of intergalactic magnetic fields (Donnert et al. 2018). In some of these situations the interaction between CRs and the ambient plasma is simply described in terms of momentum exchange, in some others the nature of CRs as charged particles, carrying an electric current, is instead fundamental. In all cases, the correct description of the interaction between CRs and the ambient medium is an extremely challenging plasma physics problem (see e.g. Zweibel (2013)).
All the information we can use to try to understand this problem comes from two kinds of observations: direct detection of CRs at Earth and detection of non-thermal radiation from astrophysical sources. In the last two decades both types of observations have experienced enormous developments.
In terms of indirect detection, most of the news have come from high energy telescopes observing the sky in the X-ray and gamma-ray bands. The launch of Chandra was quickly followed by two important results in terms of CR physics. First of all, the detection of non-thermal keV photons from SNRs highlighted the presence of multi-TeV electrons in these sources: since in most of the acceleration mechanisms considered in Astrophysics -and certainly in the one proposed for SNRs -the acceleration process only depends on particle rigidity, this was also taken as indirect evidence of the presence of multi-TeV hadrons in SNRs. In addition, the excellent spatial resolution of Chandra allowed to resolve the region of non-thermal X-ray emission, showing it to be very thin, with a typical size of few % of pc (see Vink (2012) for a review). Such a small size, when interpreted as due to synchrotron energy losses of the emitting electrons, indicates the presence of a magnetic field of several × 100 µG downstream of the SNR shock. Such a large value of the magnetic field cannot simply result from shock compression of the typical interstellar field (B ISM ≈ 3µG), but rather requires an extremely efficient amplification process to be at work.
At the same time of these discoveries, also theory was making a major step forward with the identification of a mechanism through which CRs could provide magnetic field amplification well beyond the quasi-linear theory limit of δB/B ≈ 1. The process invoked is the non-resonant streaming instability (Bell 2004) which we will discuss more in § 3: this is a current induced instability that operates in situations where the CR energy density dominates over that of the background magnetic field. This condition is easily realized in young, fast expanding SNRs, if CRs are accelerated with efficiency of order 10%, as required by the SNR-CR connection paradigm. The X-ray observations showing evidence for amplified magnetic fields were then taken as a proof of this connection: the amplified magnetic field implied that efficient CR acceleration was ongoing, and provided, at the same time, the ideal conditions for CRs to reach very high energies.
At higher photon energies, the gamma-ray satellites AGILE and Fermi, observing the sky in the 20 MeV -300 GeV energy range, and the Cherenkov detectors H.E.S.S., VERITAS and MAGIC operating between tens of GeV and ∼ 10 TeV, were providing, in the meantime, interesting but somewhat surprising data. For the first time, direct evidence of relativistic protons was found in two middle-aged (∼ 20 − 30 kyr old) SNRs interacting with molecular clouds, W44 (Abdo et al. 2010a;Giuliani et al. 2011) and IC443 (Tavani et al. 2010;Abdo et al. 2010b). These sources were not expected to be efficient CR accelerators and indeed the proton spectrum inferred from observations is cut-off at relatively low energies, ∼ 10 GeV in W44, ∼ 10 TeV in IC443. Follow-up theoretical studies showed that, in fact, at least in the case of W44, we are likely not witnessing fresh acceleration of CRs, with particles directly extracted from the interstellar plasma crossing the shock, but rather re-acceleration of particles from the galactic CR pool (e.g. Uchiyama et al. 2010;Cardillo et al. 2016). As far as young sources are concerned, gamma-ray observations have been partly disappointing, failing to provide indisputable evidence of hadronic acceleration (Funk 2015), and often showing steeper spectra than theoretically predicted in the case of efficient acceleration (Amato 2014). We will discuss gamma-ray observations of SNRs and the implications of these findings on CR acceleration in SNRs in § 2 and § 3, respectively.
In terms of direct detection, experiments such as PAMELA, AMS-02 and CREAM have provided us with very precise measurements of the CR spectrum and composition up to TeV energies. Among the many discoveries, the most important ones in terms of implications on our understanding of CR acceleration and transport are: 1) the detection of a hardening in the spectra of protons, He nuclei, and virtually all primary nuclei (Adriani et al. 2011;Aguilar et al. 2015a,b;Ahn et al. 2010a), at R ≈ 300 GV, where R = cp/Ze is the particle rigidity, with p the particle momentum, Z its atomic number, e the electron charge and c the speed of light; 2) the hardening found in the spectrum of secondary CRs at R ≈ 200 GV, with a change of spectral slope that is ∼ twice as large as that of primaries (Aguilar et al. 2018a); 3) the energy dependence of the ratio between secondary and primary cosmic rays, that has now been measured with excellent statistics up to TV rigidities (Aguilar et al. 2016c) and provides invaluable constraints on the energy dependence of particle transport in the Galaxy; 4) a rise in the fraction of positrons-to-electrons at energies larger than 30 GeV (Adriani et al. 2009a;Aguilar et al. 2013a), possibly suggesting the presence of an additional source of positrons in the Galaxy; 5) the spectrum of anti-p which is unexpectedly very close to that of protons and positrons (Aguilar et al. 2016a). We will discuss these discoveries and their implications for the origin of Cosmic Rays in § 4.

Testing the SNR paradigm
While already in the 1930s Baade & Zwicky (1934) suggested that CRs might originate from supernovae, the SNR paradigm was later formulated based on energetic arguments: SN explosions in the Galaxy can easily provide the power needed to sustain the CR population (Morrison 1957;Ginzburg & Syrovatskii 1964). Assuming that the locally measured CR energy density, w cr ∼ 1 eV/cm 3 , is representative of the CR energy density everywhere in the Galaxy, the CR production rate isẆ CR = V w cr /t conf with t conf and V the CR confinement time and volume respectively. Using t conf ≈ 15 Myr (Yanasak et al. 2001) and describing the Galaxy as a cylinder of radius R d ≈ 15kpc and height H d ≈ 5 kpc (Evoli et al. 2020b;Morlino & Amato 2020) we can estimateẆ CR ≈ 1−3×10 41 erg/s. SN events happen in the Galaxy about every 30 years and typically release 10 51 erg. Thus the power needed to sustain the galactic CR population turns out to be about 10 % of the power provided by SN events. The appeal of the SNR hypothesis was then increased by the formulation of the theory of diffusive shock acceleration (Bell 1978;Blandford & Ostriker 1978): this acceleration process, expected to be active at most shock waves, predicts the spectrum of accelerated particles to be a power-law, with an index close to −2 in the case of a strong shock, such as the blast wave of a SN explosion. Such a spectrum perfectly fits what inferred for the sources of Galactic CRs, lending support to the association. An equally remarkable, but much more recent finding, concerns the acceleration efficiency: kinetic simulations of diffusive shock acceleration (DSA) in a regime that is close to represent a SN blast wave have finally become available in the last decade, and they show that for parallel shock waves (magnetic field perpendicular to the shock surface) the acceleration efficiency is 10-20% (Caprioli & Spitkovsky 2014a), exactly as required for the SNR-CR paradigm to work.
Testing the SNR paradigm for the origin of CRs has long been one of the top priorities in High Energy Astrophysics. Since the early days, photons with energies in the range between GeVs and TeVs have provided a unique tracer of the CR population far from Earth: the interactions of CRs with the environment are expected, indeed, to produce radiation in this energy range. The process that most directly reveals the presence of hadrons is the decay of neutral pions produced when CR hadrons collide inelastically with ambient gas in the ISM: this process is commonly referred to as the hadronic production mechanism. As a rule of thumb, an inelastic nuclear collision will produce gamma rays with ∼ 10 % of the energy of the parent cosmic ray, so, for instance, photons of several tens to hundreds of TeV for CRs close to PeVs. Typically, the spectrum of the hadronic radiation at TeV mimics the parent CR spectrum shifted to lower energy by a factor 20-30 (Kelner et al. 2006;Kafexhiu et al. 2014). In the same energy range, radiation can originate from leptonic production mechanisms, mainly inverse Compton (IC) scattering of CR electrons off ambient radiation fields, and non thermal electron bremsstrahlung (Blumenthal & Gould 1970), which plays an important role in dense gas regions at sub GeV to GeV energies. A recurring difficulty when trying to infer the CR population from the GeV and TeV γ-ray emission from SNRs is to break the hadronic-leptonic degeneracy and unveil the dominant emission process: the spectral and morphological features of the emission are crucial to this task.
The gamma-ray spectrum is also affected by the particles' energy losses. For gamma-ray emitting particles, a number of loss mechanisms that are important at lower energies, like ionization and bremsstrahlung, have a negligible impact. The most relevant loss processes affect electrons, whose synchrotron and Inverse Compton emission must be properly taken into account when trying to disentangle the radiative contribution from leptons and hadrons. Electron synchrotron, particularly efficient in the strong magnetic fields believed to be associated with particle accelerators, has a fundamental role in cooling the CR electron population and thus shape the leptonic γ-ray spectrum. Therefore, observations of non-thermal radio and X-ray synchrotron emission can provide invaluable constraints on the electron population, and clues to solve the hadronic-leptonic degeneracy.
The SNR paradigm implies that young SNRs accelerate CRs up to the knee and that on average each SNR provides roughly 10 50 erg in accelerated particles. In order to test whether SNRs are the main contributors to the Galactic CR population, it is thus crucial to assess the energetics in electrons and protons and the spectra of accelerated particles in young SNRs.
Studies of the GeV-to-TeV spectra of the remnants, of their morphology and of the spatial correlation between TeV gamma-rays and the distribution of ambient gas and of X-ray radiation are the three most powerful instruments to extract the population of accelerated particles in the SNRs. GeV to TeV gamma-ray studies of both young SNRs , and dense molecular clouds close to middle-aged SNRs  were undertaken to assess the role of SNRs in the origin of Galactic CRs, by understanding and disentangling the dynamics of the coexisting and competing processes of acceleration and escape or release of particles in the interstellar medium (ISM). A useful repository of high energy observations of SNRs is the Manitoba Catalogue (http: //snrcat.physics.umanitoba.ca/).
Several young shell-type SNRs have been observed at TeV energies. Among these, the best studied ones are Cas A (Aharonian et al. 2001;Albert et al. 2007a;Ahnen et al. 2017), Tycho (Acciari et al. 2011;Archambault et al. 2017), and RX J1713.7-3946 H. E. S. S. Collaboration et al. 2018). We will discuss the latter as an example of how the SNR paradigm is tested with radiation from X-rays to gamma-rays.
As an example of multi-wavelength studies concerning middle-aged SNRs showing the so-called pion bump, a clear signature of the presence of relativistic protons, we will discuss the case of the brightest GeV source among these: W44.
2.1. Young supernova remnants: the case of RX J1713.7-3946 The supernova remnant RX J1713.7-3946, possibly associated with a star explosion seen by Chinese astronomers in the year AD393 (Pfeffermann & Aschenbach 1996;Wang et al. 1997), is one of the brightest sources of gamma rays in the TeV energy range. The sky map of RX J1713.7-3946, obtained with the H.E.S.S. telescope, is shown in the top left panel of Fig Fig. 1 show the radial profiles at keV and TeV photon energies from different regions of the shell: in four cases, the TeV supernova remnant extends, in radius, beyond the X-ray emitting shell Uchiyama et al. 2007;Tanaka et al. 2008;H. E. S. S. Collaboration et al. 2018). The main shock position and extent are visible in the Xray data and the γ-ray emission extending further is either due to accelerated particles escaping the acceleration (shock) region or particles in the shock precursor region. VHE (Very High Energy) electrons, up to hundreds of TeV, are accelerated in the shell of RX J1713-3946. These particles emit hard X-rays through synchrotron processes and TeV radiation through inverse Compton scattering off CMB photons very efficiently. A crucial piece of evidence in support of efficient particle acceleration comes from hard X-ray measurements with Suzaku (Uchiyama et al. 2007). Suzaku detected X-ray hotspots brightening and disappearing within one year timescale. If the X-ray variability is associated to the acceleration and immediate synchrotron cooling of the accelerated electrons, then the magnetic field in the shell must be roughly two orders of magnitude higher than the average magnetic field in the Galaxy, about 100 µGauss. Broadband X-ray spectrometric measurements of RXJ1713.7-3946 indicate also that electron acceleration proceeds in the most effective Bohm-diffusion regime. Finally, the presence of strongly amplified magnetic fields lends further support to the idea that not only electrons, but also protons up to 100 TeV are accelerated in the shell, as these amplified magnetic fields would be the result of CR induced instabilities (Uchiyama et al. 2007).
The high-energy part of the spectrum of RX J1713.7-3946, from X-rays to TeV photons, can be used to infer the population of particles producing this emission. The particle spectrum of the young SNR in the Left panel of Fig. 2 is characterised by a spectral index close to -2, as predicted by DSA, up to a few TeV (=10 12 eV). While the hard spectrum shows that particles are efficiently accelerated in the shell of this young SNR, the cut-off at few TeV suggests that RX J1713.7-3946 is currently accelerating electrons and protons up to 100 TeV but not up to PeV energies (H. E. S. S. Collaboration et al.

2018).
As mentioned, protons are likely accelerated in the supernova shell through the same In blue the contours of the shell as observed by XMM in X-rays (Cassam-Chenaï et al. 2004). The morphology of the TeV shell resembles closely the X-ray shell. Top right panel and following: radial profiles of the emission at selected places along the shell from H.E.S.S. and XMM-Newton. In four out of five regions of the shell, the TeV supernova remnant is more extended than the X-ray one as a result of possible escape of high energy particles from the shell. The coordinate "Radius", used for the one-dimensional profiles, refers to the mean radius "r" of the shell over which average is performed.  (Tanaka et al. 2008) are compared to leptonic (Left panel) and hadronic (Right panel) model predictions by Zabalza (2015). The total energetics in accelerated electrons and protons in the relevant leptonic and hadronic models of gamma-rays can be estimated by assuming a distance to the source of about 1 kpc (Pfeffermann & Aschenbach 1996;Wang et al. 1997). The required budget in electrons is determined only by the reported gamma-ray fluxes, if the target radiation is the CMB: one finds W e 1.2 × 10 47 erg. On the other hand, the total energy budget of protons in hadronic models depends on the highly uncertain ambient gas density: W p 5 × 10 49 (n H /1cm −3 ) −1 erg. The leptonic model in the Left panel corresponds to the electron IC scattering off CMB photons and an infrared radiation field with energy density 0.415 eV cm −3 and temperature T = 26.5 K.
The leptonic model requires a break in the electron spectrum at 2.5 TeV with the index of the electron energy distribution changing from Γ = 1.7 to Γ = 3 beyond the break. If due to synchrotron cooling, this spectral break implies a magnetic field of about 140 µG, which is much larger than what the X-ray flux allows within the same scenario, B ≈ 15µG. Similar difficulties arise if one tries to explain the break as a result of Inverse Compton cooling, in which case a photon field of energy density of 140 eV cm −3 , would be required, about 100 times larger than the average galactic value.
If the emission is hadronic in origin, a spectrum of protons with index Γ = 1.5, steepening to Γ = 1.9 at about 1 TeV, and with an exponential cut-off at 79 TeV is required. The break in the proton population spectrum can be explained if the hadronic emission is produced mostly in dense clumpy regions. In fact, low energy protons can be efficiently excluded from dense gas region during the timescale of ∼ 1000 years since the SNR explosion. The exclusion of low energy cosmic rays would also explain the hard γ-rays detected by the Fermi LAT telescope Gabici & Aharonian 2014;Inoue et al. 2012), and the lack of thermal X-ray emission from the shell. The latter has traditionally been one of the strongest argument against a possible hadronic origin of the emission from the shell of RXJ1713.7-3946, since it implies an ambient gas density as low as 0.1cm −3 (Ellison et al. 2010;Fukui et al. 2012;Gabici & Aharonian 2014). Finally hadronic gamma-ray production in gas condensations results in a narrow angular distribution of the radiation, which is beyond the reach of the current generation of Imaging Air Cherenkov Telescopes (IACTs), but could be tested with the upcoming Cherenkov Telescope Array (CTA) which will have an angular resolution of about 1-2 arcmin (Zirakashvili & Aharonian 2010;Ambrogi et al. 2018).
The conclusion one reaches, after investigating the spectrum and morphology of RX J1713.7−3946 as currently known, is that neither the hadronic nor the leptonic scenario is fully satisfactory. Each emission mechanism has strengths and weaknesses when compared with observations, suggesting that the ambient conditions might differ in different parts of the remnant, making one or the other process locally dominant. While for some specific SNRs one of the two scenarios might indeed be favored (see e.g. Funk (2015)), the general conclusions is that no known SNR has proved to accelerate particles beyond 100 TeV.

Molecular Clouds close to middle aged SNRs: the case of W44
Molecular clouds are regions of the Galaxy, typically a few tens of parsecs in radius, where the density of cold molecular gas is often orders of magnitude higher than elsewhere in the diffuse ISM. Stars are believed to be born in these clouds. Radio observations of the rotational 1 → 0 line emission of carbon monoxide are mainly used to trace the distribution of molecular gas (Heyer & Dame 2015). The cloud Galactocentric distance is usually estimated using a kinematical distance method, through each the radial velocity of the cloud is related to the rotation velocity of the Galaxy. Cross calibrations with the distance of spiral arms or known objects with precise parallax determination are also carried out (Roman-Duval et al. 2009;Reid et al. 2014).
Giant molecular clouds, which are typically 5 to 200 parsecs in diameter and have masses of 10 thousand to 10 million solar masses (Murray 2011), are excellent laboratories for CR physics. In these clouds, the hadronic channel of gamma-ray production is enhanced by the high target density, and easily dominates over the leptonic production mechanism. Contrary to the warmer atomic gas phase, which is homogeneously distributed in the Galaxy, one or a few giant molecular clouds are essentially the dominant contributions to the gas column density along a given direction in the sky. The emission enhancement associated with the clouds makes it possible to precisely locate along a given line of sight where the CR population produces the emission. Molecular clouds are thus used to perform a sort of tomography, and obtain a three dimensional view of the CR distribution in the Galaxy (Issa & Wolfendale 1981;Aharonian & Atoyan 1996;Aharonian 2001;Casanova et al. 2010a;Aharonian et al. 2020;Baghmanyan et al. 2020).
The plasma conditions in molecular clouds are generally different from those in the diffuse ISM. In addition to the gas density, also the magnetic field energy density and turbulence level are enhanced (see e.g. Crutcher (2012) for a comprehensive review). This will lead to a suppression of the diffusion coefficient and effective exclusion of lower energy CRs from the cloud (Gabici et al. 2007). If CRs can penetrate clouds, the γ-ray emission from π 0 -decay depends only upon the total mass of the cloud, M cl , its distance from the Earth, d, and the CR flux within the cloud, Φ CR . The latter is thus determined as Φ CR ∝ Φγ d 2 M cl , where Φ γ is the γ-ray flux from the cloud. Lower energy cosmic rays can be effectively excluded from penetrating clouds, which results in peculiar features in the gamma-ray spectrum from clouds, such as hardenings with respect to the average interstellar spectrum (Gabici et al. 2007).
The role of molecular clouds in testing the SNR paradigm is made particularly crucial by the time evolution of particle acceleration in SNRs. Within a DSA scenario, SNRs accelerate the highest energy CRs (in principle up to a few PeVs, but see § 3) at the transition between the free expansion and the Sedov phase, which typically happens a few tens to a few hundred years after the supernova explosion, depending on the explosion type and properties of the surrounding medium. During the Sedov phase the SN shock slows down and the magnetic field intensity decreases, so that the most energetic particles cannot be confined any longer and are free to escape. In practice, a SNR can accelerate the highest energy particles only for a short time. This fact, coupled to the low rate of PeV particles accelerating events (see § 3) makes the chances of observing a SNR when it still acts as a PeVatron very low. The runaway CRs can illuminate molecular clouds located close to the SNRs and this enhanced gamma-ray emission can thus provide crucial insights on the parental population of runaway CRs, which would otherwise escape the SNR without leaving a footprint (Montmerle & Cesarsky 1979;Issa & Wolfendale 1981;Gabici et al. 2007;Aharonian & Atoyan 1996;Aharonian 2001). The emission produced by these runaway CRs is essential to trace back acceleration up to PeV energies.
Gamma rays in association with dense molecular clouds located close to SNRs have been detected both at GeV and TeV energies. Depending on the location of massive clouds, on the acceleration history and on the timescales of the particle escape into the interstellar medium (which depend in turn on the diffusion coefficient), a broad variety of energy distributions of gamma rays is produced, from very hard spectra (much harder than the spectrum of the SNR itself) to very steep ones (Aharonian & Atoyan 1996;Gabici et al. 2009;Casanova et al. 2010b).
W44 is a 20,000 years old SNR located on the Galactic Plane at a distance of roughly 3 kpc from the Sun. W44, which is the brightest middle-aged SNR at GeV energies, is a perfect laboratory to test the coexisting processes of acceleration and escape of cosmic rays from SNRs. Due to the slow shock speed, high energy particles are expected to have already escaped from this source. A population of protons with spectral index close to -2.3 and a cutoff at about 80 GeV is likely responsible for the gamma-ray emission from the remnant measured by AGILE and Fermi-LAT (Abdo et al. 2010a;Giuliani et al. 2011). The presence of such a low energy cutoff might be the effect of the escape of the highest energy cosmic rays. Indeed, the remnant is thought to be currently a rather poor accelerator and most of its gamma-ray emission can be interpreted as due to reacceleration of ambient CRs, rather than to acceleration of fresh particles (Uchiyama et al. 2010;Cardillo et al. 2016).
Because of their short lives, SNRs are often found within the giant molecular gas complexes where they were born as massive stars. W44 resides within a giant gas complex of 10 6 solar masses, homogeneously distributed over an extended region surrounding the remnant (Dame et al. 2001). North west and south east of the remnant, Uchiyama et al. (2012) discovered two bright sources, which the authors associated to runaway cosmic rays colliding with the dense gas clouds. Peron et al. (2020) re-analysed both the Fermi-LAT data and the gas data from the clouds around W44 and noted that, despite the gas is homogeneously distributed, the GeV emissivity is enhanced only in the two regions which can be associated to regions of enhanced CR density, CR clouds, rather than gas clouds as previously thought. This phenomenon suggests the existence of a preferential path for CR escape, likely linked to the magnetic field structure in the vicinity of this source. The pion bump feature at about 1 GeV is evident in the spectra of both SNRs. For IC433 the data-points of MAGIC (Albert et al. 2007b) and VERITAS (Acciari et al. 2009) are also shown. Color-shaded areas denote the best-fit broadband smooth broken power law (60 MeV to 2 GeV); gray-shaded bands show Fermi LAT systematic errors below 2 GeV. The data points at the highest energy (4 ×10 10 -1×10 11 eV) suggest a hardening of the spectrum, which might be produced by runaway CRs illuminating a molecular cloud (MC) located in front of the shells.

Recent developments on CR acceleration in SNRs
In spite of all the recent observational developments, providing several hints of efficient CR acceleration in SNRs, a number of unsettled questions remain and challenge the association. In particular, one of the strongest arguments to look for alternative particle accelerators has to do with the difficulties SNRs have at reaching PeV energies. This aspect of the CR-SNR connection has long appeared as the most delicate.
As already mentioned, CR acceleration in SNRs is thought to be well described within the framework of diffusive shock acceleration (DSA). Acceleration mechanisms alternative to DSA have also been proposed and studied (see e.g. Lazarian et al. (2020)). In this article, however, we focus on DSA, which is still the best studied and the most promising, in this context, to reach the highest energies. The idea at the basis of the theory is that particles gain energy each time they cross the shock thanks to the discontinuity of the fluid velocity field at the shock, that leaves an unscreened electric field. The energy gain of the particle is a constant fraction (v sh /c with v sh the shock velocity) of its energy before the crossing. The particle moves diffusively between crossings, scattering on low frequency magnetic turbulence and continuously changing its pitch-angle. Reaching high energies requires a large number of crossings, which have to occur in a time shorter than the minimum between the duration of the system t life and the time-scale over which energy losses become important t loss . In the case of protons losses are negligible and the energy is limited by the time for which the SNR is active as an accelerator, t life < 10 4 yr. Then t life must be compared with the acceleration time, t acc , which depends on how quickly the particle is able to go back to the shock after each crossing: we can estimate t acc (E) ≈ D(E)/v 2 sh , where D(E) is the particle diffusion coefficient, typically an increasing function of the particle energy E. For relatively low turbulence levels one can estimate D(E) in quasi-linear theory, writing it as (see e.g. Amato (2014)): where v(E) is the particle velocity, r L (E) its Larmor radius, and W (k) is the spectrum of magnetic fluctuations causing the diffusion, with and k the wave-number resonant with the particle of energy E, namely k = 1/r L (E).
Assuming that the scattering is due to turbulence that is injected by SNRs on a typical scale L ≈ 100 pc with δB/B 0 ≈ 1, and then develops a Kolmogorov-type spectrum W (k) ∝ k −5/3 , from Eqs. 3.1 and 3.2 above †, one can estimate the diffusion coefficient as D(E) = 7 × 10 27 (E/GeV) 1/3 cm 2 /s, which, as we will discuss below, is not very different from what CR observations indicate (D(1GeV) ≈ 10 28 cm 2 /s). If this estimate of the diffusion coefficient were the appropriate one to describe particle transport in the vicinity of a SNR shock, the maximum achievable energy in these systems would only be E max ∼ few GeV, and could only improve to ∼ 100 TeV if particle transport were described by Bohm diffusion, namely δB/B 0 ≈ 1 at all relevant scales (Lagage & Cesarsky 1983). It is then clear that in order to reach PeV energies, not only efficient -Bohm like -scattering is needed, but also largely amplified magnetic fields.
As already mentioned above, this is exactly what X-ray observations of young SNRs were found to show: aside from the short variability time-scales already discussed in Sec. 2.1, X-ray synchrotron emission is seen to be confined to rims of thickness ∆ ∼ few× 0.01 pc (Vink 2012). The extremely small value of ∆ can be interpreted either as a result of radiative losses, that kill the electron population with increasing distance from the shock, or as a result of magnetic field damping (Rettig & Pohl 2012). If we interpret the thickness of the rims as the distance traveled by the emitting electrons in a synchrotron loss time, we can write with B −4 the magnetic field downstream of the shock in units of 10 −4 G. It is clear then that a magnetic field amplified by a factor of ∼100 with respect to the interstellar value (B ISM ≈ 3µG) is implied in the rims, if they are to be interpreted as the result of radiative losses. On the other hand, interpreting ∆ as a result of magnetic damping requires a value of the magnetic field which is still in the same range estimated above (Pohl et al. 2005), and hence largely amplified.
There are several mechanisms by which a magnetic field can be amplified at a shock (Bykov et al. 2012): some involve fluid instabilities (Giacalone & Jokipii 2007;Ohira 2016), some involve MHD instabilities and CR related effects (Beresnyak et al. 2009), some others CR induced instabilities, either of resonant type (Amato & Blasi (2006) and references therein) or of non-resonant type (Bell 2004;Reville et al. 2008;). Except for purely fluid instabilities (Giacalone & Jokipii 2007), all other classes require efficient CR acceleration, and in this sense the detection of amplified magnetic fields can be considered as indirect evidence of efficient CR acceleration in SNRs. However the different mechanisms proposed are not equivalent in terms of the spectrum of magnetic fluctuations they produce, and hence in terms of consequences on particle acceleration. Particle scattering is only efficient with resonant waves. Hence magnetic field amplification (MFA, hereafter) will increase scattering efficiency and lead to high maximum energies of the accelerated particles only if there is enough power on scales comparable with the particle Larmor radii. The most promising MFA mechanism in this sense turns out to be the non-resonant CR streaming instability (Bell's instability).
The super-Alfvénic streaming of energetic particles has long been known to induce the growth of magnetic fluctuations at wavelengths that are resonant with the gyroradii of the exciting particles (Skilling 1975): this is the so-called resonant streaming instability. While creating fluctuations on the right scales, this instability can only lead to MFA up to a level δB/B 0 1, and hence it is not powerful enough to explain the field strength deduced in SNRs, nor to guarantee particle acceleration up to the knee. In more recent times, however, it has been recognized that a more powerful instability arises if the CR current is large enough to twist the ambient magnetic field on the scale of the Larmor radius (r L,0 ) of the particles that carry the current (J CR > (c/4π)(B 0 /r L,0 )), or, equivalently, if the CR energy density, U CR , is larger than c/v D times the magnetic energy density (Bell 2004): where v D is the bulk velocity of CRs. If this condition is satisfied, the magnetic field grows very rapidly. The basic physical process can be described a follows: the CR current induces a compensating return current in the ISM plasma; the force J ret ∧ B induces transverse plasma motion; the current associated to this motion acts as a source of magnetic field; the result is that the magnetic field lines associated to right-hand polarized waves are stretched and for these modes the field is amplified.
In the vicinity of a shock that is accelerating particles, Eq. 3.3 can be turned into a condition for the acceleration efficiency ξ CR : Here ξ CR is the fraction of incoming flow energy (ρ ISM v 2 sh /2) that is converted into accelerated particles, ρ ISM is the mass density of the ISM plasma and v sh is the velocity of the blast wave, also coincident with the CR bulk velocity. Eq. 3.4 makes it clear that the possibility for the non-resonant streaming instability to operate strongly depends on the shock velocity: detailed calculations show that it is indeed a viable MFA mechanism in the vicinity of the fast shocks of young SNRs ), but then stops working early during the Sedov-Taylor phase of expansion of the blast wave, typically after few hundred to few thousand years, depending on the ambient medium density and magnetic field strength.
The growth initially occurs on very small scales (k ≈ (4π/c)(J CR /B 0 >> r L,0 , but quickly the power moves to larger and larger scales Riquelme & Spitkovsky (2009), possibly due to some mean-field dynamo process (Bykov et al. 2011;Rogachevskii et al. 2012), and numerical simulations show that the final outcome of the instability, when it develops at a shock, is a spectrum of fluctuations with W (k) ∝ k −1 , leading to Bohm diffusion of the particles (Caprioli & Spitkovsky 2014c). When saturation is reached, the magnetic energy density is a substantial fraction of the CR energy density (Caprioli & Spitkovsky 2014b).
The level of MFA that the non-resonant streaming instability can provide is in the correct range to explain the magnetic field strength inferred in SNRs (Schure et al. 2012), and as we just mentioned the associated scattering is efficient (Caprioli & Spitkovsky 2014c), and yet recent studies have cast doubt on the fact that these sources might be able to accelerate particles up to PeV energies (Cardillo et al. 2015). In fact, assuming that MFA is primarily due to the streaming of particles that leave the acceleration region, it is possible to build a description of the shock as a self-regulating system: when the level of turbulence in the upstream is low, a large fraction of particles can escape from the shock; when this happens, however, the large current in the upstream causes the growth of turbulence and escape is reduced; if the fraction of escaping particles becomes too small, then the turbulence level is reduced again favouring escape. The self-regulation mechanism qualitatively illustrated above translates into a quantitative prescription for the maximum energy of shock accelerated particles as a function of the system parameters. The current of escaping particles, which is what determines the level of MFA, is in turn determined by the shock velocity and the spectrum of accelerated particles, including its total energy content (which determines the amount of energy available), its slope and the maximum particle energy (which determine the current, once the energy content is fixed). When writing the equation describing E max as a function of time during the expansion of a SN blast wave, one finds that E max is an ever decreasing function of time during the SNR evolution, so that the relevant value of E max is that reached at the beginning of the Sedov-Taylor phase, namely the highest possible after a sufficient amount of mass has been processed by the blast wave. It is clear then that a higher maximum energy will be achieved in systems that enter the Sedov-Taylor phase earlier in time after the explosion. This occurs for type II explosions expanding in the dense and slow wind of a progenitor red super-giant star. In this case, assuming a ∝ E −2 spectrum, the maximum achievable energy reads: where ξ CR is the CR acceleration efficiency, E SN and M ej are the energy released and the mass of material ejected in the SN explosion, andṀ and v w are the mass loss rate and speed of the progenitor's wind. If one takes into account that the mass of the ejecta in a type II SN explosion is more likely around 10 M , it is immediately clear that reaching the knee is very challenging in this framework. One is forced to invoke extremely energetic events, or extreme acceleration efficiency, or finally extreme properties of the progenitor. In all cases these must be rare events. In fact, the product of ξ CR , E SN and the frequency in the Galaxy of the PeV producing SN explosions are constrained by the overall flux of CRs measured flux at Earth. When all of this is taken into account, one determines an event rate that cannot exceed a few in 10 4 yr (Cristofari et al. 2020).
In reality the possibility for SNRs to reach the knee becomes even more challenging when taking into account the fact that the accelerated CR spectrum is likely ∝ E −p with p > 2. Such a spectrum corresponds to a smaller current for a given E max and total energy in accelerated particles, hence requiring even more extreme parameters to reach E max ≈ PeV (e.g., Cardillo et al. (2015)). On the other hand a spectrum steeper than E −2 is exactly what gamma-ray observations of the majority of young SNRs require and what recent progress on propagation also seems to require : the CR spectrum released in the ISM will be steeper than E −2 only if the spectrum in the source is itself steeper than E −2 (Schure & Bell 2014;Cardillo et al. 2015) and a steep injection spectrum of CRs is exactly what the most recent CR data seem to suggest, as we discuss in § 4.

Cosmic ray transport in the Galaxy
As mentioned in § 1 the last decade has also been rich of observational progress providing interesting constraints on the properties of CR transport throughout the Galaxy. In this section we will review the main findings and try to put them in a coherent theoretical framework.
The first important discovery was that of a hardening in the spectrum of protons and He nuclei (Adriani et al. 2011;Aguilar et al. 2015c,d ) and also, though with somewhat lower statistics, of heavier nuclei (Ahn et al. 2010b). Within the framework of diffusive transport, the detection of a break in the spectrum of CRs can be interpreted as a signature of a change either in the injection spectrum or in the diffusion coefficient. During propagation through the Galaxy, CRs primarily loose energy due to adiabatic expansion and ionization. In addition, particles of a given species can be also lost due to spallation or decay. However, if one focuses on stable primary CRs, with energy above a few tens of GeV, energy losses during propagation can be neglected and a very simple expression for the diffuse steady state spectrum in the Galaxy can be found. If we assume that CRs are injected in the Galaxy at a constant rate and that particles then propagate diffusively with an average diffusion coefficient the steady state spectrum of stable CR nuclei in the Galaxy will be given by the product between injection and confinement time, Q(E) × τ esc . This will read, for primary nuclei: where the confinement time has been taken to be τ esc = H 2 /D(E), with H the size of the magnetized halo in which CRs are confined (see e.g. (Amato & Blasi 2018) for a more refined description). It is then clear then that the detected hardening implies a change in γ inj or δ. Several models have been proposed invoking either of the two (Tomassetti 2012;Thoudam & Hörandel 2012, and references therein). A possibility is that this feature is signaling the importance of non-linear effects in CR propagation. An early suggestion (Blasi et al. 2012) was that the break could point to the transition between scattering in self-generated and external turbulence. It had indeed been suggested since the '70s (Wentzel 1974;Cesarsky 1980) that at scales comparable with the Larmor radius of GeV particles, CRs could be an important source of turbulence in the galaxy through the resonant streaming instability. At larger scales, on the other hand, CRs become too few, given their steep spectrum, and the main source of scattering would become the large scale turbulence present in the Galaxy. The latter is usually assumed to be injected by SNRs on a scale of order tens to 100 pc and then cascade to smaller scales developing a Kolmogorov type spectrum k −5/3 . Such a spectrum translates in a diffusion coefficient with δ ≈ 1/3, flatter than the low energy value δ ≈ 0.7 that is appropriate to describe CR self-generated turbulence. A back of the envelope calculation places the transition between self-generated and external turbulence at a CR rigidity in the range 200-300 GV (Blasi et al. 2012), tantalizingly close to the value of 336 GV at which AMS-02 detects the hardening in the proton spectrum. This explanation of the hardening as due to a change in the properties of galactic transport entails a clear prediction for the spectrum of secondary CR nuclei. These are injected in the Galaxy as a result of spallation of primaries. If we approximate the spallation cross-section as independent of energy, their injection will be where n ISM is the target gas density. Their equilibrium spectrum in the Galaxy will then read: which implies that any hardening ∆ of the spectrum of primaries due to a change in the slope of the diffusion coefficient will reflect in a hardening of the spectrum of secondaries equal to 2∆: this expectation is perfectly consistent with the analysis of secondaries performed by AMS-02 (Aguilar et al. 2018b). Within this interpretation of the observed hardenings, the transport of CRs through the Galaxy becomes a complex non-linear problem, where the diffusion coefficient and the spectra of all nuclei need to be determined self-consistently. Once this complex problem is solved, however, and all available AMS-02 data are reproduced, the low energy Voyager data (Stone et al. 2013) are automatically reproduced (Aloisio et al. 2015). One important result that comes out of this analysis is directly related to the Boron-over-Carbon (B/C) measurements of AMS-02 at high energies (Aguilar et al. 2016d ). The B/C ratio has traditionally been considered the primary indicator of CR transport: B is the most abundant stable secondary and its mostly produced by the spallation of C, though the contribution by N and O is not negligible. In practice the B/C ratio provides a direct measurement of the grammage (mass per unit surface, or mass density integrated along the path-length) encountered by CRs during their propagation from their sources to Earth. If one compares Eqs. 4.5 and 4.3, for the spectrum of secondary and primary nuclei, respectively, one immediately sees that the ratio between the two is expected to scale with energy exactly as the diffusion coefficient. In fact this ratio, as measured by AMS-02 (Aguilar et al. 2016d ), well agrees with a high energy slope of the diffusion coefficient δ ≈ 0.4 (Aloisio et al. 2015;Evoli et al. 2019), but an additional, energy independent contribution seems to be required to well reproduce not only the highest, but also the lowest energy (Cummings et al. 2016) available data. The additional grammage needed, X s ≈ 0.15 g cm −2 , is independent of energy and of the order of what particles can accumulate within a SNR -or any source with an ambient density of order 1 cm −3 and a duration ≈ 10 4 yr -during acceleration (Aloisio et al. 2015;Bresci et al. 2019). Within this modeling, namely if hardenings are interpreted as a result of turbulence self-generation at low energies, the contribution X s turns out to be fundamental also for explaining another recent surprise found in CR data, namely the spectrum of antiprotons, which we will discuss in § 6.
Before concluding this section we would like to remark a most important consequence of a scenario in which CR propagation at high energy is described by a diffusion coefficient δ ≈ 0.4: CRs must be injected in the galaxy with a spectrum E −γinj with γ inj ≈ 2.3. This combination of injection and propagation parameters, and in particular the relatively weak energy dependence of the diffusion coefficient, helps to explain the low level of anisotropy detected at TeV energies, as shown by Blasi & Amato (2012a,b); Sveshnikova et al. (2013). On the other hand, as discussed in § 3, such a steep spectrum makes it very difficult for SNRs to be the primary sources of PeV CRs.
Within a DSA scenario SNRs accelerate the highest energy CRs (up to at least a few PeV) at the transition between the free expansion and the Sedov-Taylor phase, which typically happens a few hundred years after the supernova explosion. During the Sedov-Taylor phase the SN shock slows down and the magnetic field intensity decreases, so that the shock cannot confine any longer the most energetic particles, which escape the SNR. This means that a SNR can act as a PeVatron for a relatively short time. Considering the rate of SN explosions in the entire Galaxy (about 3 per century), the chances to observe a SNR when it is still a PeVatron are thus very low, and any proof of emission from PeV CRs even from young SNRs is challenging to find.

Alternative CR sources
While gamma-ray observations have proven that SNRs are efficient accelerators of cosmic ray electrons, and possibly protons, up to 100 TeV, the hypothesis that acceleration of cosmic rays proceeds up to PeV energies has been rejected in all known young SNRs because of the clear cut-offs detected at several TeVs in their spectra.
The GeV-to-TeV radiation from young shell-type SNRs, which had been long expected to provide final evidence to settle the question of the origin of cosmic rays, can be either of hadronic or leptonic origin. Typically, a 1 TeV γ-ray photon is emitted either by an electron or a proton of about 10 TeV. The cooling time for 10 TeV electrons is ≈ 5 × 10 4 yr (for IC scattering off the CMB photons) while for protons the cooling time is 5 × 10 7 (n/1cm −3 ) −1 yr (see e.g. Aharonian (2004)). Thus the ratio of the luminosity in IC gamma-rays to π 0 -decay gamma-rays is of the order of 10 3 × ( We The leptonic contribution to the emission is thus dominant over the hadronic one unless We Wp << 10 −3 or alternatively if the gas density inside the shell is n >> 1 cm −3 . This latter condition, however, if realized, would have the drawback of slowing down the shock wave very quickly and thus prevent efficient acceleration. On the other hand, a much larger energy density in protons than in electrons could result from rapid cooling of the electron population, especially in the presence of an amplified magnetic field, as inferred for young SNRs. If B > 10 µG, only a small fraction of their energy w M BR /w B ≈ 0.1(B/10µG) −2 , is released in IC gamma-rays, making the conditions for detection of π 0 decay more favourable.
To recap, no known SNR has proved to accelerate particles beyond 100 TeV. Additionally, the power in accelerated protons within SNRs derived from γ-ray observations depends on the highly uncertain local gas density n and on the magnetic and radiation fields within the remnant. As a result, it has so far been impossible to unequivocally prove that the population of SNRs in the Galaxy injects the necessary power (about 10 50 erg per supernova event) to sustain the CR population. In fact simulations of the SNR population show how the choice of acceleration and ISM parameters lead to remarkably different SNR populations as CR sources (Cristofari et al. 2017).
Some other classes of astrophysical sources, such as super-bubbles or star forming regions (see § 5.2) or remnants of GRBs in our Galaxy (Atoyan et al. 2006), and, in particular, the Centre of our Galaxy (see § 5.1), have long been proposed as alternative accelerators of particles and major contributors to the population of Galactic cosmic rays up to PeV energies (Casse & Paul 1980;Voelk & Forman 1982

Galactic Centre
A breakthrough in the quest for the origin of the highest energy CRs in the Galaxy was the discovery of a powerful PeVatron in the Centre of the Milky Way (HESS Collaboration et al. 2016).
The nucleus of the Milky Way is a very active region with numerous sources of nonthermal radiation and constitutes a unique laboratory for the study of very high energy astrophysical processes within the Galaxy and in external active galactic nuclei. The Galactic Centre (GC) is thought to host a super-massive black hole (SMBH) of 2.6 × 10 6 solar masses (Ghez et al. 2008;Gillessen et al. 2009) located very close to the dynamical centre of the Galaxy and coincident with the compact radio source Sagittarius A* (Sgr A*). Sgr A* emits radiation in X-rays and infrared through accretion of mass onto the BH. The region within 400 pc from the Centre of the Galaxy, called the Central Molecular Zone (CMZ), contains around 5% of the total Galactic molecular gas (Oka et al. 1998;Tsuboi et al. 1999). The molecular CO line towards the GC being somehow optically thick due to the high gas density, the gas distribution there is mapped additionally with other lines, such the CS radio line (Tsuboi et al. 1999). While the presence of this dense gas and dust prevents observations of this region in the optical and ultra-violet wavelengths, the Galactic Centre region is extremely bright at radio, infrared, X-ray and gamma-ray frequencies.
At very high energy (H. E. S. S. Collaboration et al. 2018), the GC hosts a bright pointlike source, HESS J1745-290, which, within errors, coincides spatially with Sgr A* and presents an energy spectrum with a steepening below 10 TeV. While the association of HESS J1745-290 with Sgr A* is well motivated in terms of spatial coincidence and of required energetics, the pulsar wind nebula (PWN) candidate G359.95-0.04 is also a viable counterpart of HESS J1745-290. HESS J1745-290 is surrounded by an extended component of VHE gamma-ray emission, correlated spatially with the CMZ. The spectrum of this radiation is a pure power law, with a spectral index -2.3 and without evidence of a spectral cutoff.
Because of their different energy distributions, the central point-like source and the extended emission cannot have the same origin. The diffuse VHE gamma-ray emission of the CMZ could be in principle produced through interactions of either relativistic protons with the ambient gas or of relativistic electrons with the radiation fields. However, a leptonic origin of the diffuse gamma-rays can be excluded. The power-law acceleration spectra of electrons should extend to about 100 TeV, which is extremely difficult because of the severe Inverse Compton and synchrotron radiative losses in the GC region. For the same reason, leptons hardly can escape the sites of their acceleration and propagate over the tens of parsec extended region of the HESS diffuse emission.
The spatial correlation between the gamma-ray emission from the Galactic Ridge and the ambient gas supports a hadronic origin of the emission. The spectrum of the parental protons -with a spectral index close to -2.4 -should extend to energies close to 1 PeV. Assuming a cutoff in the parent proton spectrum, the corresponding secondary gammaray spectrum deviates from the HESS data at 68%, 90% and 95% confidence levels for cutoffs at 2.9 PeV, 0.6 PeV and 0.4 PeV, respectively. This makes the discovery of the diffuse emission from the Galactic Centre the first robust claim of detection of a Galactic cosmic ray PeVatron.
The CR energy density in the CMZ, obtained combining the gamma-ray and gas distributions, is an order of magnitude higher than the sea of cosmic rays. The radial distribution of CRs follows a 1/r dependence, where r is the distance to the GC. This means that the PeVatron should be located within 10 pc from the GC and the total energy output in protons above 10 TeV over the whole region should amount to W CR = 10 49 erg. Such a modest energy output could be provided by a single supernova remnant event, such as Sgr A East, a Pulsar Wind Nebula. If the super-massive black hole, Sgr A*, is the GC PeVatron, then the particles producing the HESS extended emission are accelerated either in the vicinity of the SMBH, close to event horizon, or at the termination of a relativistic outflow, either a jet or a wind. This jet or wind should be injected close to the black hole and carry a substantial fraction of energy, extracted from the accretion disk. Indeed, the power required by the inferred CR distribution corresponds to 1% of the accretion power of the central black hole, and is 2-3 orders of magnitude larger than the bolometric luminosity of Sgr A*. The source, however, could have been more active in the past (HESS Collaboration et al. 2016).
Alternatively, the extended TeV emission from the Galactic center could be the combined result of acceleration within three powerful star clusters, the Arches, the Quintuplet and the Nuclear cluster (HESS Collaboration et al. 2016;Aharonian et al. 2019). In this case, the continuous injection over millions of years, required by the 1/r distribution of the diffuse emission, could result from the characteristic ages of massive clusters, roughly 10 6 years .
All the proposals discussed above assume that particle transport in the GC vicinity can be described through the same simple model of a spatially uniform diffusion coefficient adopted for the rest of the Galaxy. If this assumption is released, however, the GC excess can be interpreted as a result of a spatially dependent diffusion coefficient, as was proposed by Gaggero et al. (2017), based on Fermi-LAT data showing the spatial dependence of the gamma-ray galactic diffuse emission.

Star Clusters
Collective stellar winds and SNR shocks in clusters and associations of massive stars have long been suggested as possible, alternative or additional contributors to the Galactic cosmic ray flux (Casse & Paul 1980;Voelk & Forman 1982;Cesarsky & Montmerle 1983). Core-collapse SN progenitor stars and colliding wind binaries evolve in giant molecular clouds and mostly remain close to their birthplaces in groups of loosely bound associations or dense stellar clusters. The winds of multiple massive stars in such systems can collide and form collective cluster winds which drive a giant bubble, a so called superbubble, filled with a hot (T = 10 6 K) and tenuous (n < 0.01 cm −3 ) plasma. At the termination shock of the stellar cluster wind, turbulence can build up, in the form of MHD fluctuations and weak shocks (Torres et al. 2004;Reimer et al. 2006;Bykov & Toptygin 2001). Turbulence in the superbubble interiors can accelerate particles to very high energies, not only through the 1 st order Fermi process, but also via the 2 nd order mechanism (Bykov & Toptygin 2001). Supernova explosions of massive stars in thin and hot superbubbles can also produce efficient particle acceleration at the boundary of the superbubbles or at MHD turbulence and further amplify existing MHD turbulence (Ferrand & Marcowith 2010). It has also been recognized that multiple shocks can result in efficient acceleration beyond PeV energies (Klepach et al. 2000). The interaction of the accelerated particles with the ambient medium -often including dense molecular clouds -or with electromagnetic fields, leads to the efficient production of VHE gamma rays.
Recent observations at TeV energies of massive star-forming regions and stellar clusters, such as 30 Doradus in the LMC (H. E. S. S. Collaboration et al. 2015), Westerlund 1 (Abramowski et al. 2012) and the Cygnus region (Aharonian et al. 2002;Abdo et al. 2007;Konopelko et al. 2007;Abdo et al. 2009;Bartoli et al. 2012) in our own Galaxy, support the hypothesis that star forming regions are sites of high energy particle acceleration, and give new impulse to the γ-ray research in this field. The primary objectives of these gamma-ray observations are: 1) to constrain the fraction of mechanical energy of the stellar wind transferred to relativistic particles and hence gamma rays; 2) to unveil the physics of particle acceleration and propagation in Galactic stellar clusters and superbubbles. Furthermore, high-energy phenomena are receiving increasing attention also from the point of view of their impact on the life cycle of interstellar matter and starformation processes. The rate and efficiency of the star formation process depends, in fact, on the balance between the self-gravity of dense molecular cores and the countervailing forces which act to support the clouds. The most important of these are likely to be thermal pressure, turbulence, and magnetic fields. In order for magnetic support to be effective, a population of ionized particles must be present in the core. Since molecular clouds are opaque to ultraviolet radiation from stars, the main ionizing agent is thought to be low-energy CRs, and the magnetic support of the cloud is critically dependent on their abundance.
The Cygnus region hosts some of the most remarkable star-forming systems in the Milky Way, including Cygnus X, a star forming region at only 1.5 kpc from the Sun, with a total mass in molecular gas of a few million solar masses -at least 10 times the total mass in all other close-by star-forming regions, such as Carina or Orion -and a total mechanical stellar wind energy input of 10 39 erg s −1 , which corresponds to several per cent of the kinetic energy input by SNe in the entire Galaxy. Cygnus X hosts many young star clusters and several groups of O-and B-type stars, called OB associations. One of these associations, Cygnus OB2, contains 65 O stars and nearly 500 B stars. These super stars have created cavities filled with hot, thin gas surrounded by ridges of cool, dense gas where stars are now forming, which strongly emit at GeV energies, called the Fermi Cocoon (Ackermann et al. 2011). At TeV energies the Cygnus region shows two distinctive regions. One is possibly connected to the Cygnus X complex, the star association Cygnus OB2 and the Fermi cocoon observed at GeV energies. In the other region detected at TeV energies, the Cyg OB1 region, the Milagro collaboration discovered MGRO J2019+37 (Abdo et al. 2007), a very hard and extended source, possibly related to the massive star-forming region associated with the HII region Sharpless 104 (Sh 2-104) (Torres et al. 2004). Particle acceleration in shocks driven by the winds from the Wolf-Rayet stars in the young cluster Berk 87 in the Cyg OB1 association have also been proposed as a possible origin of the VHE gamma rays (Bednarek 2007;Bednarek & Sitarek 2007). VERITAS resolved the Milagro source into two sources, one of which, VER J2019+378, is a bright, 1 degree extended source, that likely accounts for the bulk of the Milagro emission, coincident with the star formation region Sh 2-104. Its spectrum in the range 1-30 TeV is well fitted with a power-law model of photon index 1.75, among the hardest values measured in the VHE band. The TeV counterpart of the Fermi-Cocoon has been studied with the ARGO detector up to about 10 TeV and with the HAWC detector up to 200 TeV (Bartoli et al. 2014;Hona et al. 2019). The spectral energy distribution shows a significant softening at a few TeV: this is revealed by the comparison between the ARGO and HAWC data and the Fermi-LAT data. This break in the γ-ray spectrum might hint at a cut-off in the injected CR spectrum (Hona et al. 2019), or possibly be explained as due to suppressed diffusion in the high turbulence environment of the Cygnus super-bubble, which would confine low energy particles while higher energy ones escape. Aharonian et al. (2019) argued that three ultra-compact clusters located in the Galactic Centre power the HESS diffuse emission from the CMZ. The 1/r dependence of the CR density on distance from the star cluster is, in fact, a distinct signature of continuous injection of CRs over the cluster lifetime and following diffusion through ISM. The efficiency of conversion of kinetic energy of powerful stellar winds can be as high as 10 percent. This implies that the population of young massive stars can provide production of CRs at a rate of up to 10 41 erg/s, which is sufficient to support the flux of Galactic CRs without invoking other source populations. This Galactic center PeVatrons together with the other massive star forming regions, such as Westerlund 1 and Cyg OB2, would represent the major factories of Galactic CRs . A broader review of this subject, including a discussion of the subtleties associated with acceleration and propagation of CRs in these environments, may be found in Bykov et al. (2020).

Implications of anti-matter data and CR leptons
One final subject we want to address in this paper concerns the implications of recent observational results on the spectrum of CR electrons, positrons and anti-protons. While electrons can be directly accelerated from the ISM plasma, positrons and anti-protons are extremely few in the ordinary ISM and hence they have long been considered to be purely secondaries, namely a byproduct of CR interactions during propagation in the Galaxy. If one assumes that e + andp are pure secondaries, the straightforward expectation is that the ratio between their fluxes and that of protons (their primary parent particles) should monotonically decrease with energy, as is true for the B/C ratio discussed in § 4. In fact, the ratio between e + and p should decrease even faster, because positrons are subject to radiation losses which further steepen the spectrum at high energies (see e.g. Amato & Blasi (2018) (2017)). An interesting possibility is that CRs accumulate sizeable grammage in the vicinity of their sources, where the diffusion coefficient is reduced with respect to the Galactic average. Here diffusion is energy dependent and this is where most of the Boron is produced. Then, in the rest of the Galaxy, diffusion must be energy independent and faster than usually assumed. Most e + andp would be produced by particles that accumulate most of their grammage in this second region: at a given rigidity their parent particles are ∼ 10 times more energetic than those that produce B nuclei, and hence spend a shorter time in the source vicinity.
While it is not fully clear that these assumptions would not violate other constraints, like those provided by anisotropy, a possible physical justification of such a scenario could stem from the effects of self-generated turbulence (see Sec. 3) which could considerably reduce the diffusion coefficient in the vicinity of CR sources. In terms of theory the viability of this scenario is not fully clear, especially due to the unknowns related to the abundance of neutrals, which can effectively damp the CR induced resonant streaming instability (D'Angelo et al. 2016;Nava et al. 2016). However, there are observational indications that enhanced confinement indeed occurs, at least in the vicinity of some sources. A striking example are the so-called TeV haloes (Abeysekara et al. 2017a) surrounding evolved pulsars. Many more cases are expected to be found with nextgeneration gamma-ray telescopes, such as CTA.
In reality, given the present uncertainties in many of the parameters involved, among which, most notably, thep production cross-section (Korsmeier et al. 2018), thep/p ratio is not in a statistically significant disagreement with the standard description of CR propagation through the Galaxy, nor with the B/C ratio, when all available information is taken into account and the subtle effects of re-acceleration are included in the description (Bresci et al. 2019). Indeed the high precision of the currently available data prompts theory to move forward and include second-order effects: one of these is the finite probability that during propagation particles might encounter an active CR accelerator (e.g. a SNR shock) and be re-accelerated. This has no effect on CR primaries, but flattens the spectrum of CR secondaries. Taking this phenomenon into account, together with the hardening of primaries, the energy dependence of thep production cross section, and, finally, the source grammage mentioned in § 4, allows one to reproduce the B/C ratio, thep/p ratio and all available data on stable nuclei (Bresci et al. 2019).
Differently fromp, the spectrum of e + is not easy to reproduce in the standard CR propagation scenario, unless additional sources of primary positrons are invoked.
The leptonic cosmic ray population is a small fraction of the total CR population. For instance at ∼1 GeV the ratio of the hadronic versus leptonic fluxes measured at Earth is roughly 100 : 1 (https://lpsc.in2p3.fr/crdb/). Differently from hadrons, leptons are heavily affected by energy losses due to synchrotron and inverse Compton processes in the interstellar magnetic and radiation fields. If positrons were purely secondary products of CR interactions, their spectrum would have to be steeper than that of protons at energies where losses are important. As recently shown by Diesing & Caprioli (2020), in order for this not to happen, the escape time from the Galaxy would have to be so short that protons would not be able to produce the observedp flux.
The anomalous positron spectrum has been at the origin of a huge number of papers, proposing either astrophysical or particle physics (Dark Matter annihilation related) solutions (Aguilar et al. 2019b). Limiting the discussion to astrophysical scenarios, important constraints on the plausibility of the different proposals come from the maximum age and distance of the highest energy lepton sources. As shown in Fig. 4, electrons and positrons are now detected up 20 TeV. The electron spectrum shows roughly a constant slope up to an energy of about 1 TeV, above which the spectrum becomes steeper. One possibility to interpret the break at 1 TeV is to relate it to the transition between a regime where a large number of lepton accelerators contribute to the spectrum, to a regime where only a few or even only one single electron source very close to the Sun contributes (Recchia et al. 2019).
The number of contributing sources depends, in turn, on the average intensity of the magnetic and radiation fields and on the diffusion coefficient in the Galaxy. The latter, in particular, seems to need some revision in light of the recent measurements of secondary nuclei, both stable and unstable ones, provided by AMS-02 (Aguilar et al. 2018b).
For the standard average values of the magnetic and radiation fields in the Galaxy, the cooling time of leptons with energies above few tens of GeV can be approximately estimated as: (6.1) Losses limit the distance λ from which leptons of a given energy can reach us to λ ≈ 2 D(E)t cool (E), with D(E) the galactic diffusion coefficient. Using for D(E) the estimate provided by Evoli et al. (2020b), which allows to reproduce all the available AMS-02 data on both primary and secondary nuclei (both stable and unstable), one can estimate that the highest energy leptons detected at the Earth cannot come from further than λ(1TeV) ≈ 3.5 kpc. This is a much larger distance than typically estimated in the past (see e.g. Atoyan et al. (1995)) and ensures that many sources contribute to the lepton flux even at the highest energies (Evoli et al. 2020a).
The most promising sources of primary positrons are likely Pulsar Wind Nebulae (PWNe). These are the nebulae formed by the highly relativistic magnetized wind produced by a fast spinning, highly magnetized neutron star (see e.g. Amato (2019) for a recent review). The magnetosphere of such a star is a very efficient anti-matter factory and the pulsar wind contains electrons and positrons in about equal amounts. At the wind termination shock, these particles are accelerated up to PeV energies with a spectrum in the form of a broken power-law, harder than E −2 at energies below ≈ 500 GeV and softer than E −2 at higher energies. As discussed by Bykov et al. (2017), the accelerated particles are released in the ISM after the pulsar leaves its associated Supernova Remnant to become a Bow Shock Pulsar Wind Nebula. Once taken into account, the contribution of CR leptons from these sources allows to very well explain both the overall flux of leptons and the positron fraction with very reasonable values of the few free parameters (Evoli et al. 2020a). telescopes, reveal a relatively smooth spectrum to approximately 1 TeV, where evidence of a cuto has been reported [24,26,28].
The PAMELA [30,31] and 38] satellite experiments measured the positron to electron ratio to increase above 10 GeV instead of the expected decrease [33] at higher energy, confirming earlier hints seen by the HEAT balloon-borne experiment [32]. The structure in the electron spectrum, as well as the increase in the positron fraction, may be related to contributions from individual nearby sources (supernova remnants or pulsars) emerging above a background suppressed at high energy by synchrotron losses [39]. Other explanations have invoked propagation e ects [35] or dark matter decay/annihilation processes (see, e.g., [34]). The significant disagreement in the ratio below ≥10 GeV is attributable to di erences in charge-sign dependent solar modulation e ects present near Earth at the times of measurement.
The ratio of antiprotons to protons is ≥ 2 ◊ 10 ≠4 [40] at around 10-20 GeV, and there is clear evidence [41]for the kinematic suppression at lower energy that is the signature of secondary antiprotons. The p/p ratio also shows a strong dependence on the phase and polarity of the solar cycle [42] in the opposite sense to that of the positron fraction. There is at this time no evidence for a significant primary component of antiprotons. No antihelium or antideuteron has been found in the cosmic radiation. The best measured upper limit on the ratio antihelium/helium is currently approximately 1 ◊ 10 ≠7 [43] The upper limit on the flux of antideuterons around 1 GeV/nucleon is approximately 2 ◊ 10 ≠4 (m 2 s sr GeV/nucleon) ≠1 [44].
A useful method for calculating the e ect of solar modulation including time, charge-sign, and rididity-dependent e ects is given in Ref. [45].  Most measurements are made at ground level or near the top of the atmosphere, but there are also measurements of muons and electrons from airplanes and balloons. Fig. 30.4 shows measurements of negative muons [46][47][48][49][50][51]. Since µ + (µ ≠ ) are produced in association with ‹µ(‹µ), the measurement of muons near the maximum of the intensity curve for the parent pions serves to calibrate the atmospheric ‹µ beam [52]. Because muons typically lose almost 2 GeV in passing through the atmosphere, the comparison near the production altitude is important for the sub-GeV range of ‹µ(‹µ) energies.

Cosmic Rays in the Atmosphere
The flux of cosmic rays through the atmosphere is described by a set of coupled cascade equations with boundary conditions at the top of the atmosphere to match the primary spectrum. Numerical or Monte Carlo calculations are needed to account accurately for decay and energy-loss processes, and for the energy-dependences of the cross sections and of the primary spectral index   Beatty et al. (2004) are reported, together with a number of model predictions: the black curve is for pure secondary production in the standard scenario (Moskalenko & Strong 1998); the blue one includes a more sophisticated propagation model (Gaggero et al. 2013); the green one invokes dark matter decay (Ibarra et al. 2013); the red one includes a contribution from pulsars (Yin et al. 2013

Conclusions and Outlook
In this article we have tried to review some recent theoretical and observational developments in the quest for the origin of galactic Cosmic Rays. High energy observations, in the X-rays, gamma-rays and VHE gamma-rays have opened new possibilities to test the long standing paradigm that wants these particles mainly produced in the blast waves of supernova explosions.
X-ray observations have provided undisputable proof that SNRs can accelerate electrons up to tens of TeV, and at the same time have shown the presence of amplified magnetic fields in their environment. These magnetic fields might be taken as a hint of efficient particle acceleration, although alternatives are possible. For sure their strengths are in the range needed to make SNRs act as PeVatrons, accelerating particles up to the knee energy, where the galactic cosmic ray proton spectrum is expected to end. In spite of this, theory points to a rare occurrence and an extremely short duration of the PeVatron phase for SNRs (Sec. 3).
Multi-TeV photons (10-100 TeV photons) are the most direct messengers of the presence of protons in the critical knee region. Secondary electrons, produced as decay products of charged pions, may provide complimentary information through their X-ray synchrotron emission, but this channel has a lower radiative efficiency and secondary leptons are often subdominant with respect to primaries. In the end, the spectrum of the highest energy gamma-rays qualifies as the best source of information about the extension of the parent proton spectrum, and hence is the only direct tool to identify the PeVatrons and obtain essential information about the physics of particle acceleration and the formation of the knee.
Detailed gamma-ray studies of young objects in this class have failed to show clear evidence of acceleration above ∼ 100 TeV (Sec. 2.1). Clear evidence of accelerated hadrons has only been found in middle-aged SNRs interacting with molecular clouds. The wealth of target for nuclear interactions makes these objects ideal candidates to be bright in gamma-rays of hadronic origin, without contamination from leptonic Inverse Compton emission. However, these sources are not very efficient as accelerators and the maximum energy they achieve is limited in all cases to 10 TeV (Sec. 2.2).
At the same time, observations in the TeV band have highlighted the existence of other kinds of sources showing relatively hard gamma-ray spectra extending above 100 TeV, in association with the galactic center region (Sec. 5.1) and with young stellar clusters (Sec. 5.2). How propagation in these complex environments, very likely more turbulent than the average ISM in the Galaxy, affects the spectrum finally released is yet to be thoroughly investigated.
On the other hand, direct measurements of CRs, and in particular high statistics spectra of primary and secondary nuclei, have provided new important insights on the properties of CR transport in the Galaxy, highlighting the presence of a change in propagation properties in the energy range 200-300 GeV. A possible interpretation of this finding is that below such energy CRs are mostly scattered by self-generated turbulence, another bit of information that stresses the importance of non-linearities in CR physics, and the importance of understanding how these particles affect not only the environment of their acceleration sites, but also the properties of the ISM at large (Sec. 4).
Finally, a fundamentally important piece of information that we discussed, has to do with the anti-matter component of CRs (Sec. 6). Direct detection experiments showed that the leptonic anti-matter component of CRs can hardly be interpreted as of pure secondary origin. Additional sources of primary positrons seem to be required. The most natural candidates are pulsars and their nebulae, well known to be excellent positron factories. These will release their positrons only at late stages of their evolution, after leaving the parent SNRs. Crucial information about the particle release by these sources is likely to come from deeper studies of the gamma-ray halos that have been found to surround some of them Abeysekara et al. (2017b).
In a time when a combination of new and surprising observations and theoretical difficulties are pushing the CR physics community to a deep re-evaluation of the SNR-CR connection, crucial tests are expected to come from high energy astrophysical observations. All the recent surprises and issues that we have discussed through this article suggest that different astrophysical sources might contribute at different energies and at different levels to the CR population. This prompts to the need of a census of all possible CR contributors, a purpose that is well served by an unbiased survey of the Galactic Plane in the crucial multi-TeV energy range. Such survey is currently being undertaken with the high altitude water Cherenkov detector, HAWC (Malone & HAWC Collaboration 2017). A PeVatron search conducted with detailed spectroscopy, higher sensitivity above 10 TeV and higher angular resolution is a key scientific goal of the upcoming Cherenkov Telescope Array (Cherenkov Telescope Array Consortium et al. 2019). In the search for PeVatrons in the Galaxy, CTA will benefit from a wider energy range and a much better angular resolution: while the latter will be essential to establish correlations with the gas maps and sources observed in other wavelengths, the former will finally give us access to the hundreds of TeV energy range, which is the realm where the hadronic mechanism is believed to be strongly dominant over the leptonic one.