Fermi bubbles as sources of cosmic rays above 1 PeV

Fermi bubbles are giant gamma-ray structures extended north and south of the Galactic center with characteristic sizes of the order of 10 kpc discovered by the Fermi Large Area Telescope. Good correlation between radio and gamma-ray emission in the region covered by Fermi bubbles implies the presence of highenergy electrons in this area. Due to high energy losses it is rather problematic to transfer relativistic electrons from the Galactic disk toward the Fermi bubbles. Therefore it is natural to assume that these electrons are accelerated in-situ. Additionally this acceleration mechanism should also affect protons. In particular it may re-accelerate Galactic cosmic rays produced by supernova remnants. Unlike electrons, protons have huge lifetimes and therefore re-acceleration should not be a local effect but affect the whole Galaxy. The effect may even be observed near the Earth. In our model we propose that hadronic cosmic rays (CR) below the “knee” of the observed CR spectrum are produced by Galactic supernova remnants distributed in the Galactic disk. Re-acceleration of these particles in Fermi Bubbles produces CRs above the knee. This model provides a natural explanation of the observed CR flux, spectral indices, and matching of spectra at the knee.


Introduction
Fermi bubbles were discovered by [1] and [2] almost seven years ago and they still remain as one of the most attractive astrophysical events.The bubbles are presented by giant structures located above and below the Galactic plane visible in gamma-ray [3] and radio [4] energy bands.Despite their nature being still enigmatic, the location of these objects indicates their connection with past or present activity in the center of our Galaxy.Different models relate the bubbles to starburst activity [5], single [2,[6][7][8] or multiple [9] energy release events on a central black hole.
Such objects should not be unique phenomena observed only in our Galaxy.However, Fermi bubbles are too faint and it is impossible to detect exactly the same structures elsewhere.Yet some much more powerful objects with similar properties and probably similar nature are observed in some galaxies with active nuclei.For example even more giant structures are clearly seen in the direction of Cen-A in GHz radio [10,11], GeV [12] and TeV [13] gamma-ray ranges.Giant X-ray and radio lobes (bubbles) were found also in the galaxies NGC 3801 [14], Mrk 6 [15] and Circinus Galaxy [16].
Taking all facts together one can conclude that Fermi bubbles are indeed real structures.Since they enclose huge energies [2] and their volume is comparable to that of the Galactic disk, Fermi bubbles have a potential to affect the distribution of the cosmic rays (CRs) in the a e-mail: chernyshov@lpi.rub e-mail: hrspksc@hku.hkc e-mail: dogiel@lpi.rud e-mail: cmko@astro.ncu.edu.twGalaxy.In particular in our paper [17] we showed that Fermi bubbles may be responsible for the formation of the CR spectrum above the "knee".Below we will briefly recapitulate the major points of this model and also discuss some consequences.

Fermi bubbles as giant re-acceleration sites
The origin of the gamma-ray emission from Fermi bubbles is not yet established.However, is was shown by Cheng et al. [18] they should act like giant accelerators.Indeed even in the case of hadronic origin of the gamma-ray emission some additional primary leptonic component is necessary to reproduce microwave emission.The emission produced by secondary electrons have too steep a spectrum to match the experimental data [4].These electrons have relatively high energies and should be accelerated insitu, otherwise an unrealistic velocity of outflow from the central part of the Galaxy is required.
Such large accelerators which fill a substantial volume of the Galaxy should also modify the distribution of Galactic cosmic rays.It is generally accepted that supernova (SN) explosions in our Galaxy can provide enough energy to produce the observed total luminosity of the CR [see e.g., reference 19].Moreover diffusive shock acceleration naturally explains the power-law CR spectrum [20,21], and taking into account propagation effects, their spectral index.However, many fundamental questions related to the assumption of SNRs being the sources of the Galactic CR are still open.One of them is the maximum energy of CR which can be estimated from ISVHECRI 2016 the age of a typical SNR, T , expanding with a velocity of u sh E max ∼ Zeβ sh BT u sh , ( where β sh = u sh /c, B is the magnetic field strength at the shock and the term E = β sh B in this case can be interpreted as an effective electric field.
For the parameters of the standard Galactic SNR and for reasonable values of the Galactic magnetic field the maximum energy of CR protons cannot exceed 10 13 -10 14 eV.In more complex models outside quasilinear approximation it was shown that the magnetic field at the shock can be amplified.As a result for a conservative set of parameters the maximum energy reaches a value of about 10 15 eV.
The important point of the CR spectrum is a sudden steepening around 3 × 10 15 eV which indicates a change of the acceleration or propagation mechanism.Smooth attachment of the spectra above and below this energy and sharpness of the transition indicates that we are dealing with a single spectrum rather than with the sum of two distinct components.The review of different models suggested to explain this phenomenon can be found in the original paper [17].
In summary, it is generally agreed that SN shocks can only accelerate particles to energies E < 10 15 eV.However, Fermi bubbles with age and size exceeding that of typical SNR by 4 orders of magnitude can easily accelerate particles to much higher energies.Using Eq. ( 1) one can estimate that the bubbles have a potential to form the spectrum of CR in between 3 × 10 15 eV and 10 18 eV.We do not consider the possibility that Fermi bubbles can form a whole spectrum of CR both below and above 3 × 10 15 eV.The reason for this is that the total power needed for the luminosity of CRs in our Galaxy, L C R ∼ 10 41 erg/s [22], and the average energy release toward the Fermi bubbles is also 10 41 erg/s [9].To explain the origin of CR by acceleration in Fermi bubbles one requires the acceleration efficiency to be close to 100% which is unlikely.On the other hand the sudden change of the slope of the spectrum around the "knee" can be naturally explained due to a change of the acceleration properties.Thus our model can be described in the following way: SNRs in the disk accelerate particles with power-law distribution up to energies of 3 × 10 15 eV and Fermi bubbles further reaccelerate these particles up to 10 18 eV.

Structure of shocks in the Fermi bubble
In the model of Fermi bubbles origin we suggested that the activity responsible for the formation of the bubbles was due to stellar capture and tidal disruption by a central black hole [9,17].The average time between two successive captures in the Galaxy is between 10 4 years and 10 5 years [23].Thus the activity is periodic and the characteristic period between two events is shorter than the characteristic lifetime of the Fermi bubbles.Periodic energy releases in the Galactic center should form a series of shocks propagating through the halo.We note that from the numerical simulations it appears that only very powerful events expected in the case of capture of massive stars can form shocks. Thus the amount of shocks inside the Fermi bubbles should not be very large.
The separation between shocks can be estimated from their velocity, u, and the characteristic period of capture, τ cap , as l sh = τ cap u = 30(τ cap /3 × 10 4 yr)(u/10 8 cm/s) pc.(2) However, the exact separation between two consecutive shocks depends on the actual time separation of two consecutive capture events and their energy releases.There is another important spatial scale which characterizes processes of particle acceleration by a single shock, i.e., the acceleration length scale in a single shock, l D ∼ D/u, where D is the spatial diffusion coefficient near a shock whose value depends on particle interactions with small scale magnetic fluctuations and u is the shock velocity.

Particle acceleration by the bubble shocks
The problem of particle acceleration in conditions of supersonic turbulence (multi-shock structure) was extensively analyzed before.In a series of papers by [24][25][26] applied to acceleration processes in OB associations, which is quite similar to the structure of the Bubble, they introduce a nondimensional parameter characterizing the acceleration regime as The corresponding energy, E 1 , that separates different regimes can be estimated from the condition ψ ∼ 1 or l D (E 1 ) ∼ l sh which for the conditions of the Fermi bubble is In the case of ψ 1 or l sh l D analyzed in [24,25], there is a combined effect of a fast particle acceleration by a single shock, which generates the spectrum E −2 and relatively slow transformation of this spectrum due to interactions with other shocks (stochastic Fermi acceleration) into a hard E −1 spectrum in the intershock medium at relatively low energies.However, it is unclear if such slow transformation can be completed within the life time of the shocks in the Bubble.Furthermore, the hard E −1 spectrum requires a significantly higher power to be formed in comparison with a softer spectrum with a slope of E −2 .Thus the transformation of the spectrum requires a power which significantly exceeds 10 41 erg/s which the bubbles, as we mentioned above, can not supply.It is reasonable to assume that Fermi bubbles do not affect the spectrum of CR below E 1 .This statement is in agreement with our initial assumption that SNRs are the major contributors for CRs with energies E < 10 15 eV and the exact particle spectrum generated from the Bubble is unimportant in the energy range of E < 10 15 eV.
For ψ 1 or l sh l D the acceleration regime changes to a pure stochastic acceleration by a supersonic ISVHECRI 2016 turbulence.In the stationary case the equation for accelerated CRs can be presented in the form [25] (5) where ρ and z are the cylindrical spatial coordinates, p is the particle momentum, D(ρ) is the spatial diffusion coefficient and κ is the momentum diffusion coefficient.
Proton acceleration in the Bubble depends sensitively on the acceleration parameters and structure of the Bubble.In the following we present a detailed analysis.We present the bubble region as a cylinder extending above and below the Galactic plane from z = 0 to z = ±H with a radius ρ = ρ B .As boundary conditions we put the density of particles f equal to zero at the Galactic halo surface The diffusion coefficients inside and outside the bubble are supposed to be different where D B = Lc/3 is the coefficient inside the bubble due to interactions with a supersonic turbulence and D G is the average diffusion coefficient in the Galaxy defined for example in [19].The momentum diffusion coefficient is The momentum dependence of f can be presented by a power-law function, f ( p) ∝ p −γ , where γ should be determined from Eq. ( 5) and is approximately Thus a proper choice of D B allows to form a power-law distribution with a certain value of a spectral index.
To consolidate this idea, we work out a concrete numerical model.Essentially, we solve the stationary state of the CR transport Eq. ( 5) in our Galaxy with two Fermi Bubbles (one on each side of the Galactic plane).We modeled our Galactic halo as a cylinder of radius ρ G = 20 kpc, and the top and bottom at ±10 kpc from the midplane.Each Fermi Bubble is also a cylinder of the same height ±10 kpc, but with a radius ρ B = 3 kpc.The spatial diffusion coefficients are different inside and outside the bubble as described by Eq. ( 7).
Since we expect the average separation between shocks in the Fermi bubbles to be of the order of 100 pc [17] we consider a constant spatial diffusion coefficient and adopt D B = 2.08 × 10 30 cm 2 s −1 .Outside the bubble, we take into account the energy (or momentum) dependence of the spatial diffusion coefficient and adopt D G = D 0 ( pc/4 GeV) 0.6 , D 0 = 6.2 × 10 28 cm 2 s −1 [cf.27].
We assume that there is little or no stochastic acceleration outside the bubble (see Eq. ( 7)), and adopt κ B H 2 /D B = 1.9 (i.e., κ B = 4.4 × 10 −15 s −1 or the corresponding acceleration time scale is 7.6 Myr).
For Galactic SNRs we adopt the distribution suggested by [28] and modified it with a Gaussian thickness profile Figure 1.The CR spectrum at the Earth as a combination of the SNR contribution (in the Galactic disk) and the stochastic acceleration in Fermi Bubbles.The black solid line is the spectrum from our numerical model.For references to the experimental data see the original paper [17].In this model, D B = 2.08 × 10 30 cm 2 s −1 inside the bubbles and D G = 6.2 × 10 28 ( pc/4 GeV) 0.6 cm 2 s −1 outside, κ B H 2 /D B = 1.9, and the injection spectrum from SNR µ = 4.35 (see Eq. ( 10)).
here we take h = 100 pc, R = 8 kpc.We adopt the idea that SNRs inject energetic particles in the form of a power law with a high-energy cutoff at p max c ≈ 3 × 10 15 eV.Therefore, together with the SNR distribution (Eq.( 9)), the source function is Finally, the appropriate boundary conditions for the momentum coordinate are where the energy of the lower momentum boundary is p low c = 10 12 eV, and the upper momentum boundary is p up c = 3 × 10 18 eV.The condition at the lower momentum ensures that the spectral index matches that of low-energy CRs (say E < 10 12 eV).
The spatial boundary conditions are and The spectrum evaluated at the Earth's position is the solid line shown in Fig. 1.The model fits the data reasonably well and it is not coincident that the spectra join smoothly at the knee.The value of the parameter E 1 based on Eq. ( 4) is not reliable.To check how the result depends on this parameter we estimated the spectrum for different values of E 1 (see Eq. ( 7)).It turns out that the shape of the spectrum does not depend on E 1 as long as E 1 < 3 × 10 15 eV.We should note that in the estimations above we assumed that the CR spectrum consists only of protons.This assumption is for illustrative purposes only.To obtain the spectrum of different elements one should note that the propagation and acceleration parameters depend only on rigidity.Thus the spectra of all elements should have the same shape but should be shifted along the energy axis in accordance to their charge Z [see e.g.Fig. 5 in reference 29, for the details].Due to the same shape of the spectrum one should expect several "knees" in the spectrum corresponding to different elements at energies of 3Z × 10 15 eV.According to [30] there is a kneelike structure in the heavy component of CR spectrum at about 8 × 10 16 eV (apparently corresponding to iron) that strengthen this idea.
In this case the spectrum of particles above the knee should be slightly softer than the one estimated by us for the proton-only case.But as one can see from Eq. ( 4) the value of the spectral index can be easily adjusted by a proper choice of the acceleration parameters.
The important test of this model can be done via neutrino emission.Indeed we expect that highenergy CRs are concentrated in the area of Fermi bubbles [17].Therefore one should expect there should be an enhancement of the emission from this region.
It is interesting that IceCube has detected 4-5 out of 28 events above 30 TeV, which may originate from cosmic rays with energies > 10 15 eV from Fermi Bubbles, Lunardini et al. [31].This discovery seem to strengthen our theory.To test it we estimated expected intensities of the neutrino emission.
Since the statistics of IceCube data is very poor we cannot consider numerical values of the intensities as reliable data.Instead we estimated the neutrino flux from the Galactic disk and compared it to the neutrino flux from the Bubble area.The spatial distribution of relative intensity of the high-energy neutrino emission across the sky is shown in Fig. 2. Despite that there is a prominent excess at high latitudes, it turns out that the Bubble to disk ratio is very low: it does not exceed 10% in even the most promising scenarios.
Since in the IceCube data there is no prominent signal from the Galactic disk we have to conclude that the observed excess is not related to that expected in our model and is probably due to an external source.

Conclusion
We suggested a model of CR origin in the energy range 3 × 10 15 eV < E < 10 18 eV.In the model we assumed that CRs with energies below 3 × 10 15 eV are produced by SNRs in the Galactic disk.These particles diffuse out to the Galactic halo and are then re-accelerated in Fermi bubbles.For a proper choice of the acceleration parameters it is possible to reproduce the right value of the spectral index in the energy range above the "knee".The spectrum below the "knee" remains undisturbed.This model also naturally explains the smooth attachment of the spectra near the "knee".It also explains the change of chemical composition of CRs at high energies if the rigidity dependence of the acceleration and propagation is taken into account.
We also estimated the possible enhancement of the high-energy neutrino emission from the central area of the Galaxy.It turns out that the flux from Fermi bubbles should not exceed 10% of the flux expected from the Galactic disk.This value is very low and can not be detected at the present moment.

Figure 2 .
Figure 2. The relative intensity of high-energy neutrino emission across the sky.Here φ is the Galactic longitude and θ is the Galactic latitude.The picture is symmetrical across θ = 0 and across φ = 0 so only one corner is shown.