THE FERMI BUBBLE AS A SOURCE OF COSMIC RAYS IN THE ENERGY RANGE >10¹⁵ eV

The recent discovery of a pair of giant Fermi Bubbles in the Galactic center (GC) is one of the most remarkable events in astrophysics. The first indication of the structure appeared in a paper by Dobler et al. (2010), which they called the Fermi Haze. Using the special procedure of background subtraction of the Fermi Large Area Telescope (LAT) data, Su et al. (2010) discovered a pair of symmetric structures above and below the Galactic plane in the GC direction. The origin of the bubble, if its existence can be proved, is still enigmatic, and up to now a few models have been presented in the literature. The team that subtracted this structured gamma-ray emission from the total diffuse Galactic emission presented different explanations of the phenomenon, but they seemed to favor the model that described a single massive release of energy in the GC when a huge cloud of gas or a star cluster was captured by the central supermassive black hole about 10 Myr ago.

A similar explanation was suggested by Guo & Mathews (2011) and Guo et al. (2011). They assumed that the Fermi Bubbles were created by a recent active galactic nucleus jet activity about 1–2 Myr ago, which was active for a duration of ∼0.1–0.5 Myr, releasing energy totaling ∼(1–8) × 10⁵⁷ erg. The bipolar jets were ejected into the Galactic halo along the symmetric axis perpendicular to the Galactic plane.

It is important to note that the existence of the bubbles was first evidenced in X-rays detected by ROSAT as a narrow envelope with very sharp edges (Bland-Hawthorn & Cohen 2003) and later the Wilkinson Microwave Anisotropy Probe (WMAP) detected an excess of radio signals at the location of the gamma-ray bubbles (Finkbeiner 2004; Dobler & Finkbeiner 2008; Dobler 2011). The ROSAT structure is explained as a fast wind with a velocity u_w ∼ 10⁸ cm s⁻¹ driving a shock into the halo gas. This phenomenon requires an energy release of about 10⁵⁵ erg at the GC and this activity should be periodic on a timescale of the order of 10 Myr. This requirement of energy release in the GC is consistent with the observations of the existence of hot plasma with a temperature of about 10 keV in the GC region but with a radius of 30–50 pc only. This cannot be confined and will escape from the GC with a speed of u_w ∼ 10⁸ cm s⁻¹ (Koyama et al. 2007). Therefore, sources with a power of about 10⁴¹ erg s⁻¹ must have heated the plasma or released ∼10⁵⁵ erg from the GC in the past 10 Myr. However, Crocker & Aharonian (2011) have proposed a relatively slow energy release (∼10³⁹ erg s⁻¹) from supernova (SN) explosions as a source of proton production in the GC. The observed gamma rays come from hadronic processes of the protons in the halo. The plasma in the halo is extremely turbulent and the protons are trapped for a time comparable to the Hubble time. But this model requires a separate origin of electrons, which have a much shorter lifetime than protons, to explain the WMAP data.

Intensive energy release has been observed, indeed, at the center of normal galaxies as strong variations of X-ray radiation. There are common characteristics of these X-ray sources. First, all of them have been bright sources, and their X-ray luminosity could go up to about 10⁴⁴ erg s⁻¹. Second, they have shown a high level of variability in their X-ray light curve within years. In the "high state," the luminosity of one source could be at least 100 times higher than the luminosity in its "low state." Third, most of them have a super-soft spectrum during the flare, with effective blackbody temperatures of only about 10–100 eV (Komossa & Bade 1999; Halpern et al. 2004). The classic examples that satisfy these characteristics are RX J1624.9+7554, RX J1242.6−1119A, RX J1420+5334, RX J1331−3243, and NGC 5905. Many scenarios have been proposed to explain these phenomena, but most of them fail to explain some of the observed results. A detailed discussion of these scenarios can be found in Komossa & Bade (1999). Among all of the listed models, the tidal disruption model is the most commonly accepted, and it gives the most satisfactory explanation to the observations by considering the radiation from the disk. In this model, when a star passes by a black hole within a capture radius, where the black hole tidal force becomes stronger than the self-gravity of the star, the star can be captured. The detailed capture and disruption process of a main-sequence star has been studied by several authors (e.g., Rees 1988; Cannizzo et al. 1990). The capture rate of main-sequence stars in our Galaxy and in other galaxies is about 10⁻⁴ yr⁻¹ to 10⁻⁵ yr⁻¹ (see Syer & Ulmer 1999; Alexander 2005). Recently, more stellar capture events have been observed (e.g., Esquej et al. 2008; Gezari et al. 2008, 2009; Komossa & Bade 1999; Cappelluti et al. 2009). Dynamical studies of nearby galaxies suggest that most, if not all, galaxies with a bulge component host a central supermassive black hole and that the bulge and black hole mass are tightly correlated (Magorrian et al. 1998; Tremaine et al. 2002; Greene & Ho 2007).

The Burst Alert Telescope on board Swift has identified a transient X-ray source called GRB110328A (Cummings et al. 2011) with later optical identification (Cenko et al. 2011; Leloudas et al. 2011), which is located in the direction of the constellation Draco when it erupted in a series of X-ray flares. The distance to this source is determined to be z ∼ 0.35 by using H_β and O iii emission lines by the Gemini telescope (Levan et al. 2011b). The characteristics of GRB 110328A appear inconsistent with those of a gamma-ray burst (Barthelmy et al. 2011). In fact, its time-dependent characteristics, including various timescales in light curve, multi-wavelengths, etc., seem to be better explained in terms of the tidal disruption of a star by a supermassive black hole (Barres de Almeida & De Angelis 2011; Bloom et al. 2011; Levan et al. 2011a; Burrows et al. 2011; Zauderer et al. 2011). All of these recent observations suggest that stellar capture processes are quite common in other normal galaxies.

Observations have also revealed much evidence of unusual processes occurring in the central region of our Galaxy, for instance, the enigmatic 511 keV annihilation emission discovered by the International Gamma-Ray Astrophysics Laboratory (see, e.g., Knoedlseder et al. 2005) whose origin is still debated. The hot plasma has a temperature of about 10 keV, which cannot be confined in the GC and, therefore, sources with a power of about 10⁴¹ erg s⁻¹ are required to heat the plasma (see Koyama et al. 2007 and references therein). In fact, plasma outflows with velocities of ≳ 10⁷ cm s⁻¹ are observed from the nuclear region of our Galaxy (see Crocker et al. 2010) and from the nucleus of Andromeda (Bogdan & Gilfanov 2010). Time variations of the 6.4 keV line and X-ray continuum emission observed in the direction of molecular clouds in the GC which are supposed to be a reflection of a giant X-ray flare occurred several hundred years ago (see Inui et al. 2009; Ponti et al. 2010; Terrier et al. 2010 and references therein). HESS observations of the GC in the TeV energy range indicated an explosive injection of cosmic rays (CRs) there, which might be associated with the supermassive black hole Sgr A* (e.g., Aharonian et al. 2006).

In a series of papers (Cheng et al. 2006, 2007; Dogiel et al. 2009a, 2009b, 2009c, 2011), we developed a model of energy release in the GC due to star accretion onto the central black hole for the interpretation of X-ray and gamma-ray emission from the GC. Our goal was to explain these observational data in the framework of a single model. Basic assumptions in these models are (1) the Galactic supermassive black hole Sgr A* can capture a star at a rate of ν_s ∼ 10⁻⁴–10⁻⁵ yr⁻¹, and (2) the energy release from each capture in the form of a flux of subrelativistic protons is W ∼ 4 × 10⁵² M²_*R_*⁻¹M^1/3_bh(b/0.1)⁻² erg, where M_* (in units of M_☉) and R_* (in units of R_☉) are the mass and the radius of the captured star, M_bh (in units of 10⁶ M_☉) is the mass of the supermassive black hole, and b is the ratio of the periapse distance r_p to the tidal radius R_T (see the review of Alexander 2005). In a timescale much longer than the capture timescale, this model can be treated to have an average power injection $\dot{W}\sim 3 \times 10^{40}$ erg s⁻¹. These protons heat the surrounding plasma by Coulomb losses to 10 KeV.

Based on this model, Cheng et al. (2011; CCDKI model) argued that up to several hundred capture events might have occurred in the past 10 Myr, which may have generated a series of shocks propagating through the central part of the Galactic halo and thus produced accelerated relativistic electrons responsible for the bubble emission. Processes of charged particle acceleration by the bubble shocks in terms of sizes of the envelope, maximum energy of accelerated particles, etc., may differ significantly from those obtained for SNe. In this paper, we examine whether a "signal" from charged particles accelerated in the bubble region can be seen in the spectrum of CRs observed at the Earth. We present simple estimations of hadronic CR acceleration by the Fermi Bubble shocks up to energies above 10¹⁵ eV. The paper is organized as follows. In Section 2, we review the current understanding of CR acceleration by supernova remnants (SNRs) and conclude that this process can only produce CRs with energies less than 10¹⁵ eV. We present a simple solution of the multiple-shock structure in the halo in Section 3. In Section 4 we discuss the protons accelerated by the bubble shocks. We emphasize that a broken power law of particle distribution must be formed because of the finite spacing between consecutive shocks and the spectral break naturally occurs at 10¹⁵(u/10⁸ cm s⁻¹)(l_sh/30 pc)(B/5 μG) eV. Charged particles below and above this critical energy are accelerated by two different acceleration mechanisms. In Section 5, we calculate the total particle spectrum by summing up the contribution from all shocks in the bubble and compare it with the observed hadronic CR spectrum with energies larger than 10¹⁵ eV. In Section 6, we suggest a model that can produce the CR spectrum within and beyond the "knee" (around 3 × 10¹⁵ eV). Summary and discussion is presented in Section 7.

From a general point of view, SN explosions are enough to supply the power needed for the luminosity of CRs in our Galaxy, L_CR ∼ 10⁴¹ erg s⁻¹ (for a general review see Berezhko et al. 1994; Reynolds 2008). Diffusive shock acceleration (see Krymskii 1977; Bell 1978) is considered to be a viable and natural mechanism for CR accelerated by SNRs. The mechanism produces a power-law spectrum with the necessary spectral index that is observed experimentally. The simplest one-dimensional kinetic equation describing this process has the form

$$\begin{equation} \frac{\partial f}{\partial t}+\frac{\partial }{\partial z} \left(u(z)f-D\frac{\partial f}{\partial z}\right) -\frac{1}{3p^2}\frac{du(z)}{dz}\frac{\partial }{\partial p}\left(p^3f\right)=0, \end{equation} \tag{ 1 }$$

where z is the coordinate perpendicular to the shock front, p is the particle momentum, u(z) is the velocity distribution which describes a velocity jump at the shock, and D is the spatial diffusion coefficient. The solution of this equation is a power-law function, f(p)∝p^−γ, in which the spectral index γ is a function of the velocity jump at the shock. For strong shocks with a Mach number much larger than unity, γ = 4. The corresponding energy spectral index is ν = 2 (N(E)∝E^−ν).

The current status of the observations of middle-aged SNRs by LAT on board Fermi with an energy range from 0.2 to 100 GeV has provided some insight into the shock-acceleration theory of SNRs (Castro & Slane 2010; Uchiyama 2011). Assuming that the gamma rays are produced by hadronic processes, Castro & Slane (2010) deduced the spectral index of four SNRs ranging between 2.1 and 2.4. However, whether the observed GeV gamma rays are produced by hadronic processes or leptonic processes is very difficult to differentiate. On the other hand the ambiguity can be removed if broadband emissions are observed. In particular, if GeV and TeV spectra can be described by a single power law, which is steeper than E⁻², the hadronic processes could be the more favorable mechanisms. Currently, about 10 SNRs have been detected in both GeV and TeV bands, including Tycho and CTB37A, whose GeV–TeV gamma-ray emission shows uniformly steep spectral indices of about 2.3 and 2.2, respectively (see Table 1 of Caprioli 2011). All of these recent observations are consistent with conventional SN shock-acceleration theories, which suggest that the spectrum of CRs is roughly described by E⁻².

Many fundamental questions related to the assumption that SNRs are the sources of Galactic CRs are still open. The maximum energy of the accelerated particles is the main concern for this scenario, which can be roughly estimated from a very simple relation. The acceleration time at the shock is τ_acc(E) ∼ D(E)/u²_sh, where u_sh is the shock velocity (∼a few 10⁸ cm s⁻¹). The minimum value of the diffusion coefficient at the shock follows from the Bohm diffusion scenario, i.e., D_Bohm(E) = (c/3)r_L(E), where r_L is the Larmor radius of the particle. Equating the acceleration time with the lifetime of the shock T we get an estimate for the maximum energy of accelerated particles after time T,

$$\begin{equation} E_{\rm max} \sim Ze\beta _{{\rm sh}}u_{{\rm sh}} B T , \end{equation} \tag{ 2 }$$

where β_sh = u_sh/c, B is the magnetic field strength at the shock. The combination u_shB/c can be interpreted as an effective electric field.

For an SNR of typical age τ_SNR ∼ 1000 yr, the maximum energy of protons is easily estimated by requiring that the acceleration time remains smaller than τ_SNR. Lagage & Cesarsky (1983) and Berezhko & Völk (2000) demonstrated that the maximum energy of protons within the scenario of Bohm diffusion is as large as E_max ∼ 10¹³–10¹⁴ eV for standard galactic SNRs. Berezhko et al. (1994) estimated the efficiency of acceleration when a feedback reaction of accelerated particles on the front structure was included and they showed that in the Bohm limit CRs absorb about 20% of the explosion energy. The acceleration process acts as an effective viscosity in widening the region of the shock velocity jump and eventually the acceleration process stops.

However, outside the quasi-linear model the acceleration of CRs at the shock fronts of SNRs may make the acceleration of particles more effective. As Bell (2004; see also Bykov et al. 2009) showed that during acceleration at shocks of SNR, the magnetic non-resonant fluctuations were strongly driven. A nonlinear MHD simulation indicated that CR-excited turbulence could amplify the magnetic field. It appears that acceleration to the spectral break at 10¹⁵ eV normally ceases as an SNR enters the Sedov phase. Thus, CR acceleration by SN shocks can only provide particles with energies less than 10¹⁵ eV.

The spectral index of the observed CR flux changes from 2.7 to 3.1 around energy 10¹⁵ eV, and this is known as the knee. The standard model of CR acceleration by SN shocks cannot explain CR energies above the knee because it only produces a single power-law spectrum up to the energy around 10¹⁵ eV. In addition, the CR spectrum flattens again for energies above 10¹⁸ eV, and this is known as the "ankle." Large size is required to accelerate and to confine charged particles above the ankle (the Larmor radius at these energies is comparable with the halo height). The origin of CRs above the ankle is generally attributed to an extragalactic origin because those particles could not be confined inside the Galaxy and known potential Galactic accelerators could hardly accelerate particles to such high energies.

The origin of the steepening for E ≳ 10¹⁵ eV is still an open question and different mechanisms of CR acceleration in the range 10¹⁵–10¹⁸ eV in the Galaxy have been proposed. Ptuskin et al. (2010) assumed that the CR flux with energies above 10¹⁵ eV is produced by very young galactic SNRs. They modeled the particle acceleration by spherical shocks with back-reaction of CR pressure on the shock structure. The significant magnetic field amplification in young SNRs produced by CR streaming instability may lead to a flux of CRs with the maximum energy of accelerated particles about 5 × 10¹⁸ eV. In this model, the steepening of the CR spectrum at the knee position is due to distortion of the spectrum ejected from young SNRs by the propagation process.

Another interpretation was suggested by Erlykin & Wolfendale (2006), Erlykin et al. (2011), and Lagutin et al. (2008) who assumed that CRs at the knee were produced by a single, recent local SN. Recently Butt (2009) summarized problems (including energies within the knee) of the conception that isolated SNRs are the main accelerators of CRs and discussed alternative scenarios of CR acceleration.

The steepening of the CR spectrum at high energies may also be the result of a change in properties of diffusion in the interstellar medium. This effect of propagation was mentioned first by Syrovatskii (1971) who noticed that the standard diffusion of CRs in the interstellar medium might be changed by convection due to a drift of these particles in the large-scale Galactic magnetic fields. This model was developed in Ptuskin et al. (1993) who assumed two types of CR diffusion in the Galaxy: the usual diffusion due to particle scattering on fluctuations of random magnetic fields and the Hall anisotropic diffusion (drift motion) due to the large-scale Galactic magnetic field whose effect might become important just above the knee energy.

As an alternative model, Jokipii & Morfill (1987) suggested a mechanism of acceleration in the Galaxy of ultrahigh energy CRs in a Galactic wind and its hypothetical termination shock. In this scenario, SNRs accelerate the bulk of CRs up to 10¹⁵ eV. These particles are further accelerated up to 10¹⁹–10²⁰ eV at a termination shock which is at a distance of a few hundred kiloparsecs from the disk. Ip & Axford (1992) analyzed multiple interactions of particles with SNRs in the Galactic disk as a source for CR acceleration above the knee. However, too many shocks are required in the disk in order to produce CR flux at the knee. Bykov & Toptygin (1993) showed that regions of CR acceleration to energies above 10¹⁵ eV might be OB associations where the concentration of shock fronts is very high. We will discuss this model in Section 4.

In summary, it is generally agreed that SNR shocks can only accelerate particles to energies less than 10¹⁵ eV. On the other hand, accretion processes in the GC may generate giant shocks which are effective for particle acceleration above 10¹⁵ eV. It is interesting to note that the notion of recurrent activities in the Galactic center as a source of cosmic rays was proposed some 30 years ago (Khazan & Ptuskin 1977; Said et al. 1981; Ptuskin & Khazan 1981), although not in the same perspective as the present paper.

As was assumed in CCDKI, the central massive black hole captures a star every τ_cap ∼ 3 × 10⁴ yr; as a result, about W ∼ 3 × 10⁵² erg of energy in the form of subrelativistic particles is released. This heats up the surrounding gas in the central region of our Galaxy. The hot gas expands into the halo and forms a propagating upward shock. The situation is very similar to that of the stellar wind of a massive star blowing into its surrounding medium (see, e.g., Weaver et al. 1977; Bisnovatyi-Kogan & Silich 1995).

The gas distribution in the disk and in the halo was derived in Cordes et al. (1991) and can be presented as a double exponential distribution

$$\begin{eqnarray} n(\rho,z) = &&0.025 \exp \left(- \frac{z}{1\,{\rm kpc}}\right) \times \exp \left[- \left(\frac{\rho }{20\,{\rm kpc}}\right)^2\right] \nonumber \\ && +\,0.2 \exp \left(- \frac{z}{0.15\,{\rm kpc}}\right)\nonumber \\ &&\times \exp \left[{-}\!\left(\frac{\rho -4\,{\rm kpc}}{2\,{\rm kpc}}\right)^2\right] \ {\rm cm}^{-3}. \end{eqnarray} \tag{ 3 }$$

The energy release in the GC as a result of star capture can be either impulsive or continuous depending on the characteristic times of star capture and energy dissipation of subrelativistic protons (plasma heating by Coulomb losses). The capture time is roughly τ_cap ∼ 3 × 10⁴ yr and the dissipation time τ_diss in the CCDKI model is determined by the rate of ionization losses of protons injected with energy E_p, which is given by

$$\begin{equation} \tau _{\rm diss}\simeq 10^6 \left(\frac{n}{1\,{\rm cm}^{-3}}\right)^{-1}\sqrt{\frac{E_p}{100\,{\rm MeV}}}\ {\rm yr}, \end{equation} \tag{ 4 }$$

where n is the gas density in the vicinity of the GC, which can be quite high (see the discussion in Cheng et al. 2007), e.g., in nearby molecular clouds n > 10⁴ cm⁻³.

If τ_cap ≪ τ_diss, we have the case of a stationary energy injection in the GC. In this case the region of heated gas is bounded by a single shock (see Weaver et al. 1977). For τ_cap ≫ τ_diss a multi-shock structure is formed in the halo with shocks of different ages. A similar multi-shock structure can also be created if there are epochs of high-frequency star captures in the GC. Thus, the number of shocks is determined by the injection and dissipation parameters.

For a highly simplified case of the exponential atmosphere with the scale height z₀, i.e., the plasma density n(z) = n₀exp (− z/z₀), an analytic solution of shock propagation was obtained by Kompaneets (1960; see also the review of Bisnovatyi-Kogan & Silich 1995). This solution gives a qualitative picture of the shock propagation and parameters of the medium bounded by the shock that roughly described the situation expected in the Fermi Bubble. If the rate of energy injection is L then the radius of the shock as a function of the height z and the time t is

$$\begin{eqnarray} \rho (z,t)&=&2z_0\arccos \left\lbrace \frac{1}{2}\exp \left(\frac{z}{2z_0}\right)\right.\nonumber\\ &&\times\, \left.\left[1-\left(\frac{y}{2z_0}\right)^2+ \exp \left(-\frac{z}{z_0}\right)\right]\right\rbrace . \end{eqnarray} \tag{ 5 }$$

Here z is the coordinate perpendicular to the Galactic plane,

$$\begin{equation} y=\int \limits ^t_0\left[\frac{\left(\gamma _p^2-1\right)}{2}\frac{\lambda Lt}{V(t)mn_0}\right]^{1/2}dt, \end{equation} \tag{ 6 }$$

V(t) is the current volume enveloped by the shock

$$\begin{equation} V(t)=2\pi \int \limits ^{a(t)}_0\rho ^2(z,t)\,dz, \end{equation} \tag{ 7 }$$

a(t) is the position of the top of the shock

$$\begin{equation} a(t)=-2z_0\ln \left(1-\frac{y}{2z_0}\right), \end{equation} \tag{ 8 }$$

L = W/τ_cap is the average luminosity of the central source, γ_p is the polytropic coefficient, and λ describes the fraction of explosion energy converted into the thermal energy of gas (see Bisnovatyi-Kogan & Silich 1995). As follows from Equations (5) and (8), for the finite time t₁ determined from the condition y(t₁) = 2z₀, the shock breaks through the exponential atmosphere and the bubble top a(t₁) tends to infinity while the bubble radius in the Galactic plane (z = 0) tends asymptotically to ρ = 2z₀cos ⁻¹(1/2) ≃ 2z₀ and for z ≫ z₀ to ρ ≃ πz₀.

A numerical solution of the system (5)–(8) in dimensionless coordinates for shock waves at different ages is shown in Figure 1. This figure is meant to be illustrative only. In reality the distribution of shocks should be far more complicated. Neglecting the thin 4 kpc torus component of Equation (3), the scale height of the "atmosphere" is z₀ = 1 kpc. The shock distribution in Figure 1 is suggestive of a cylindrical bubble with an edge at a radius ρ_B ≃ 3 kpc.

Zoom In Zoom Out Reset image size

Based on this, we put forward a quasi-stationary model of the Fermi Bubble which we regard as the source of hadronic CRs with energy larger than 10¹⁵ eV. The essence of the CCDKI model is that energy is quasi-periodically injected into the halo when the stellar capture processes take place and may exist over a timescale comparable with the age of Milky Way. Consequently, the Fermi Bubble should have a stationary structure. The idealized Kompaneets solution above shows that there is a stationary sideway boundary for shocks. For quasi-periodic star capture, the bubble interior is filled with shocks propagating in series and eventually stopping at ρ_B ≃ 3 kpc. However, in the Kompaneets solution dissipation processes are ignored. The realistic situation has to be described by a set of dissipative hydrodynamic equations, which takes account the shock propagation in a non-uniform medium and various dissipation processes, including shock heating, energy transfer into CRs, slowing down due to accumulation of material, etc. Fitting with the observed gamma-ray spectrum, Cheng et al. (2011) concluded that electrons should have an escape timescale of 15 Myr in order to explain the spectral break position. The characteristic dissipation timescale of the shocks should be of the same order if these electrons are transported away by the shocks. The shocks will be mostly dissipated when they arrive at the sideway boundary (ρ_B) and there will be no pile-up of shocks at ρ_B. As the speed of shocks along the bubble axis is progressively larger than the sideway speed, the upper and lower boundaries of the bubble will be the same as the halo boundary (z = ±H = ±10 kpc from the midplane). Thus, the bubble has a stationary structure.

The same result for the dimension of the bubble can be obtained if we use the swept-up mass model proposed by Cheng et al. (2011). In this model the radius of the sideway shock front or the swept-up front is given by $\rho _s=\sqrt{\vphantom{A^A}\smash{\hbox{${2\lambda W/\pi m n \Delta z u_{s\rho }^2}$}}}$, where u_sρ is the speed of the sideway shock front (or the swept-up front) and Δz = uτ_cap. With u ∼ 10⁸ cm s⁻¹, τ_cap ∼ 3 × 10⁴ yr, Δz ∼ 30 pc. We argue that the sideway shock front or the swept-up front will disappear when its speed is smaller than the local sound speed. In Cheng et al. (2011), we have estimated that the shocks heat up the halo to ∼1 keV and the characteristic sound speed is of the order of v_s ∼ 3 × 10⁷ cm s⁻¹. Thus putting u_sρ ∼ v_s, the sideway boundary ρ_B ≈ 3.2 kpc(λW/2 × 10⁵² erg)^1/2(Δz/30 pc)^−1/2(n/10⁻³ cm⁻³)^−1/2(v_s/3 × 10⁷ cm s⁻¹)⁻¹ (note that λW is the fraction of injected energy converted into thermal energy of gas).

Correct analysis of shock acceleration in the bubble requires sophisticated calculations in each stage of this process which we perform later. Now we present simple estimates of the characteristics of the spectra of the accelerated particle in the framework of the CCDKI model.

Below we analyze the spectrum of protons accelerated in the bubble and discuss whether the bubble's contribution to the total flux of CRs in the Galaxy may explain the knee steepening. We remind the reader that the generally accepted point of view is that the flux of relatively low energy CRs (<10¹⁵ eV) is generated by SNRs, which eject a power-law spectrum E⁻² into the interstellar medium. This spectrum is steepened by propagation (escape) processes in the Galaxy in accordance with the spectrum observed near Earth (for details, see Berezinskii et al. 1990). However, these sources can hardly produce CRs with energies >10¹⁵ eV, at which a steepening (the knee) in the CR spectrum is observed. For characteristics of the knee spectrum and models of its origin, see the review of Kotera & Olinto (2011). We suggest that the bubble could generate the flux of CRs at energies >10¹⁵ eV because the shocks in the bubble have much larger length scales and longer lifetimes in comparison with those in SNRs.

In the framework of CCDKI, the bubble may fill with hundreds of shocks propagating in series one after another, though a single shock structure cannot be excluded. The average separation between two shocks is given by

$$\begin{equation} l_{{\rm sh}}=\tau _{\rm cap}u=30\,\left({\tau _{\rm cap}\over 3\times 10^4\,{\rm yr}}\right) \left({u\over 10^8\,{\rm cm\,s^{-1}}}\right)\ {\rm pc}. \end{equation} \tag{ 9 }$$

However, the exact amount of time between two consecutive shocks depends on the actual time between two consecutive capture events and their energy releases. There is another important spatial scale which characterizes processes of particle acceleration by a single shock: the diffusion length scale at a single shock l_D ∼ D/u. Here u is the shock velocity and D is the spatial diffusion coefficient of the energetic particles near a shock which depends on particle interaction with small-scale magnetic fluctuations. In the Bohm limit, D ∼ cr_L(E)/3, where r_L(E) = E/ZeB is the particle Larmor radius. In this case

$$\begin{equation} l_D\sim {cr_L\over u}={cE\over ZeB u}. \end{equation} \tag{ 10 }$$

The problem of particle acceleration in conditions of supersonic turbulence (multiple-shock structure) has been extensively analyzed (e.g., Spruit 1988; Achterberg 1990; Schneider 1993; Melrose & Pope 1993) as well as quasi-periodic flows (e.g., Webb et al. 2003). In a series of papers by Bykov & Toptygin (1993, 2001) and Bykov & Fleishman (1992) the idea was applied to acceleration processes in OB associations, which is quite similar to the structure of the bubble. They introduced a dimensionless parameter characterizing the acceleration regimes:

$$\begin{equation} \psi =\frac{l_{{\rm sh}}}{l_D}\sim \frac{ul_{{\rm sh}}}{D}\sim \frac{ul_{{\rm sh}}}{cr_L}. \end{equation} \tag{ 11 }$$

The critical energy E₁ that separates two regimes of acceleration can be estimated from the condition ψ ∼ 1 or l_D(E₁) ∼ l_sh. For the conditions of the Fermi Bubble the critical energy is

$$\begin{equation} E_1\approx {ZeBul_{{\rm sh}}\over c} = 10^{15}\,Z\left({B\over 5\,\mu {\rm G}}\right) \left({l_{{\rm sh}}\over 30\,{\rm pc}}\right) \left({u\over 10^8\,{\rm cm}}\right)\ {\rm eV}. \end{equation} \tag{ 12 }$$

In the case of ψ ≫ 1 or l_D ≪ l_sh, the analysis in Bykov & Toptygin (1993) and Bykov & Fleishman (1992) showed that there is a combined effect of a fast particle acceleration by a single shock, which generates the spectrum E⁻² and relatively slow transformation of this spectrum due to interaction with other shocks (stochastic Fermi acceleration) into a hard E⁻¹ spectrum in the intershock medium at relatively low energies. However, it is unclear if such slow transformation can be completed within the lifetime of the shocks in the bubble. A detailed numerical analysis is needed. Furthermore, from the general point of view the characteristic acceleration time is quite short in the range E ≲ E₁, which is roughly given by the shock-acceleration time cr_L/u². In the range E ≳ E₁, the acceleration timescale increases to the time of stochastic acceleration, cl_sh/u². With an average Galactic spatial diffusion coefficient D_G outside the bubble and a Galactic halo of height H, the characteristic escape time is τ_esc ∼ H²/D_G. We expect that escape processes, which play a crucial role in determining the particle spectrum shown in our next analysis, are insignificant in the range E ≲ E₁. Therefore the particle spectrum produced by the bubble should be ∼E^−ν for E < 10¹⁵ eV, where 2 > ν > 1. As discussed in Section 2, SNRs are the major contributors for CRs with energies E ≲ 10¹⁵ eV; the exact particle spectrum generated from the bubble is unimportant in the energy range of E ≲ 10¹⁵ eV.

In the case of ψ ≪ 1 or l_D ≫ l_sh, inside the region of supersonic turbulence the acceleration regime shifts to a pure stochastic acceleration by the supersonic turbulence. We extend the equation derived by Bykov & Toptygin (1993) for CR acceleration in supersonic turbulence in a stationary state and axisymmetric geometry to include spatial dependent diffusion coefficient and external source,

$$\begin{eqnarray} &&\frac{\partial }{\partial z}\left(D(\rho,p)\frac{\partial f}{\partial z}\right) +\frac{1}{\rho }\frac{\partial }{\partial \rho } \left(D(\rho,p)\rho \frac{\partial f}{\partial \rho }\right)\nonumber\\ &&\quad +\,\frac{1}{p^2}\frac{\partial }{\partial p}\left(\kappa (\rho,p) p^2\frac{\partial f}{\partial p}\right)=-Q(\rho,z,p), \end{eqnarray} \tag{ 13 }$$

where ρ and z are the cylindrical spatial coordinates and p is the particle momentum. D(ρ, p) is the spatial diffusion coefficient and κ(ρ, p) is the momentum diffusion coefficient. Their spatial and momentum dependence in our model is described below. Q(ρ, z, p) is the possible CR source which will be useful in our numerical example in Section 6.

As we mentioned above, it is reasonable to assume that CRs with energies E ≲ 10¹⁵ eV are supplied by SNRs. Therefore, in this section, we concentrate on the analysis of the acceleration of CRs in the energy range E ≳ E₁ in the bubble by supersonic turbulence. In Section 6, we will treat the case with SNRs and bubble and deal with energy from less than 10¹² eV to larger than 10¹⁸ eV.

We set the boundary condition of the distribution function at the Galactic halo outer boundary

$$\begin{equation} f|_\Sigma =0,\quad {\rm at}\quad \rho =\rho _G\quad {\rm and}\quad z=\pm H. \end{equation} \tag{ 14 }$$

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 model the bubble region as a cylinder extending above and below the Galactic plane from z = 0 to z = ±H with a radius of ρ = ρ_B and assuming there is no CR source inside the bubble (i.e., Q = 0). The diffusion coefficients inside and outside the bubble are supposed to be different:

$$\begin{eqnarray} &&D(\rho)=D_B\,\theta (\rho _B-\rho)+D_G\,\theta (\rho -\rho _B), \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&\kappa (\rho,p) =\kappa _B\,p^2\theta (\rho _B-\rho), \end{eqnarray} \tag{ 16 }$$

where D_B ∼ cl_sh/3 is the coefficient inside the bubble as a result of interactions with a supersonic turbulence and D_G is the average diffusion coefficient in the Galaxy, e.g., defined in Berezinskii et al. (1990). The momentum diffusion coefficient is κ_B ∼ u²/D_B. The momentum dependence of f is represented by a power-law function, f(p)∝p^−γ, where γ should be determined from Equation (13).

To understand the dependence of γ on other parameters, we make two simplifications of Equation (13), which do not affect the value of γ significantly. First, for H < ρ_G, as expected from Strong & Moskalenko (1998), particles that escape through the radial boundaries at ρ = ρ_G are insignificant (see Berezinskii et al. 1990), and we can shift the halo boundary to infinity, i.e., ρ_G = ∞. Second, we model the axisymmetric geometry of the problem as planar geometry (i.e., we assume ∂f/∂ρ ≫ f/ρ). We go back to the axisymmetric geometry afterward.

As in Bulanov et al. (1972) and Bulanov & Dogel (1974), we search for solutions to Equation (13) by the method of separation variables, f = R(ρ)Z(z)p^−γ. The solution for Z(z) has a very simple form

$$\begin{equation} Z_n(z)=\cos (k_n z/H), \end{equation} \tag{ 17 }$$

where k_n = π(n + 1/2). We should point out that to make sure f is non-negative we must take n = 0 for physical solutions.

To illustrate ideas, we consider the case D_G = D_B and approximate the axisymmetric geometry as planar (i.e., d²R/dρ² ≫ 1/ρ dR/dρ). Using the dimensionless variable ϱ = ρ/H, Equation (13) can be simplified as

$$\begin{equation} {{d}^2 R\over {d}\varrho ^2}- \left[k_n^2+\gamma (3-\gamma){\kappa _B\theta (\varrho _B-\varrho)H^2\over D_B}\right]R=0, \end{equation} \tag{ 18 }$$

which has the exact form of the Schrödinger equation for a rectangular potential well:

$$\begin{equation} {{d}^2\Psi \over {d}\varrho ^2} +{2m\over \hbar ^2} \left({\cal E}-U_0\right)\Psi =0,\quad {\rm for}\quad 0<\rho <\rho _B \end{equation} \tag{ 19 }$$

and

$$\begin{equation} {{d}^2\Psi \over {d}\varrho ^2} +{2m\over \hbar ^2}{\cal E}\Psi =0,\quad {\rm for}\quad \rho >\rho _B \end{equation} \tag{ 20 }$$

when we define ${\cal E}$ and U₀ as

$$\begin{equation} {\cal E}=-\,{\hbar ^2k_n^2\over 2m}\,\quad {\rm and}\quad U_0={\hbar ^2\over 2m}\gamma (3-\gamma){\kappa _B H^2\over D_B}. \end{equation} \tag{ 21 }$$

The solution of this equation is well known (e.g., Landau & Lifshitz 1991, chapter 3, Section 22):

$$\begin{eqnarray} \Psi (\varrho)&=&C_1\exp (k_n\varrho ),\quad {\rm for}\ \varrho <0,\,{\rm where}\ k_n=\sqrt{-2m{\cal E}}/\hbar \nonumber \\ \Psi (\varrho)&=&C_2\exp (-k_n\varrho ),\quad {\rm for}\ \varrho >\varrho _B, \nonumber \\ \Psi (\varrho)&=&C\sin (\varsigma \varrho +\delta ),\quad {\rm for}\ 0<\varrho <\varrho _B,\nonumber\\ && {\rm where}\ \varsigma =\sqrt{2m({\cal E}-U_0)}/\hbar. \end{eqnarray} \tag{ 22 }$$

From the continuity of the logarithmic derivative of Ψ at the well boundaries we have the condition

$$\begin{equation} \arcsin \left(\frac{\hbar k_n}{\sqrt{-2mU_0}}\right)=\frac{(j\pi -k_n\rho _B)}{2} \end{equation} \tag{ 23 }$$

and

$$\begin{equation} \cos \xi =\pm \chi \xi ,\quad {\rm for\ odd}\ j, \end{equation} \tag{ 24 }$$

$$\begin{equation} \sin \xi =\pm \chi \xi ,\quad {\rm for\ even}\ j, \end{equation} \tag{ 25 }$$

where

$$\begin{equation} \xi =\frac{k_n\rho _B}{2}\,\quad {\rm and}\quad \chi =\frac{\hbar }{\rho _B}\sqrt{\frac{-2}{mU_0}}. \end{equation} \tag{ 26 }$$

For strong acceleration, u²H²/D²_B ≫ 1, Equations (24) and (25) determine a finite number of levels. From these equations, the first term of the series gives

$$\begin{equation} \gamma =\frac{3}{2}+\sqrt{\frac{9}{4}+\frac{\pi ^2D_B}{\rho _B^2\kappa _B}}\simeq 3. \end{equation} \tag{ 27 }$$

For weak acceleration, u²H²/D²_B ≪ 1, we can use the shallow-well solution presented in Landau & Lifshitz (1991). In this case there is only one level at E₀ ≃ U₀ that gives

$$\begin{equation} \gamma \simeq \frac{3}{2}+\sqrt{\frac{9}{4}+\frac{\pi ^2D_B}{H^2\kappa _B}}\gg 1. \end{equation} \tag{ 28 }$$

A rough estimate of the power of CR production by the bubble can be done in the same way as presented in Berezinskii et al. (1990) for GeV CRs. The energy density of CRs at E = 3 × 10¹⁵ eV is n_CR ≃ 6.7 × 10⁻¹⁷ erg cm⁻³ (see Kotera & Olinto 2011). Then the power required for the bubble to produce the knee at the Earth is

$$\begin{equation} W_B\sim \frac{cn_{{\rm CR}}M_H}{x}, \end{equation} \tag{ 29 }$$

where M_H is the total mass of hydrogen in our Galaxy, which is about 10⁴³ g, and x is an extrapolation of energies >10¹⁵ eV of the CR grammage derived from the chemical composition by Jones et al. (2001) up to energies of about several hundred GeV, x(E) ∼ 11.8 × (4.9 GeV/E)^0.54 g cm⁻². Then we obtain the required power of CR sources at the knee energy range W_B ∼ 2 × 10³⁹ erg s⁻¹, which can easily be supplied by star capture processes.

More accurate values of γ can be derived from numerical calculation of the axisymmetric case:

$$\begin{equation} {1\over \varrho }{\!{d}\over {d}\varrho } \left(\varrho {{d} R_i\over {d}\varrho }\right)- \left[k_n^2+\gamma (3-\gamma){\kappa _{Bi}H^2\over D_B}\right]R_i=0. \end{equation} \tag{ 30 }$$

Here the index i = 1, 2 denotes the regions inside (ϱ < ϱ_B) and outside (ϱ ⩾ ϱ_B) the bubble, respectively. Note that κ_B1 = κ_B and κ_B2 = 0 (see Equation (16)). A solution inside and outside the bubble is searched as series of the Bessel functions (J_ν, inside the bubble) and the McDonald functions (K_ν, outside the bubble).

The boundary conditions at the bubble radius, ϱ = ϱ_B, are

$$\begin{equation} R_1(\varrho _B) = R_2(\varrho _B)\,\quad {\rm and}\quad D_B\left.{{d}R_1\over {d}\varrho }\right|_{\varrho _B} = D_G\left.{{d}R_2\over {d}\varrho }\right|_{\varrho _B}. \end{equation} \tag{ 31 }$$

These relations can be satisfied if

$$\begin{equation} \alpha _1 = k_n^2 + \gamma (3-\gamma)\kappa _B H^2/D_B < 0, \end{equation} \tag{ 32 }$$

and the above requirement implies γ > 3. For n = 0 we have

$$\begin{equation} R_1(\varrho) = C_1 J_0(\sqrt{-\alpha _1}\varrho)\,\quad {\rm and}\quad R_2(\varrho) = C_2 K_0(\pi \varrho /2), \end{equation} \tag{ 33 }$$

or

$$\begin{equation} \frac{D_B\sqrt{-\alpha _1}J_1(\sqrt{-\alpha _1} \varrho _B)}{J_0(\sqrt{-\alpha _1 } \varrho _B)} = \frac{\pi D_G K_1(\pi \varrho _B/2)}{2K_0(\pi \varrho _B/2)}. \end{equation} \tag{ 34 }$$

Numerical results of different ratios of D_G/D_B as a function of the ratio between the times of particle escape from the Galaxy H²/D_B and the time of Fermi acceleration 1/κ_B ∼ D_B/u² are shown in Figure 2. The dashed curve shows the approximate estimate given by Equation (28), which indicates that Equation (28) is a good approximation and we will use it in the next section of data fitting for simplicity.

Zoom In Zoom Out Reset image size

The spatial distribution of accelerated protons for different ratios D_G/D_B is shown in Figure 3. We can see that these two distributions do not have a qualitative difference even if the ratios differ by a factor of 10.

Zoom In Zoom Out Reset image size

In this section, we use the general models developed in the last section to fit the CR spectrum for E > 10¹⁵ eV but we ignore any possible spectral modulation effect as a result of propagation from the bubble to the Earth. We have shown that the charged particles in the bubble can be described by a power-law spectrum,

$$\begin{equation} \frac{d\dot{N}}{dE}=\frac{\dot{N}_{0}}{E_1}\left({E\over E_1}\right)^{-\nu }, \end{equation} \tag{ 35 }$$

where E₁ is given by Equation (12) and ν = γ − 2, which can be expressed as

$$\begin{equation} \nu =-\frac{1}{2}+\sqrt{\frac{9}{4}+\frac{\pi ^2D_B^2}{u^2H^2}}, \end{equation} \tag{ 36 }$$

for E₁ < E < E₂ where E₂ is the high-energy cutoff and will be discussed below. We have used the momentum diffusion coefficient κ_B = u²/D_B. Here $\dot{N}_{0}$ is a normalization constant which is chosen to fit the observed CR spectrum by assuming that the CR spectrum with energy E > 10¹⁵ eV is entirely contributed by particles from the bubble. E₂ can be estimated by Equation (2), which gives

$$\begin{equation} (E_2)_{\rm min} \approx 3\times 10^{17}\left({B\over 5\, \mu {\rm G}}\right) \left({T\over 10\,{\rm Myr}}\right)\left({u\over 10^8\,{\rm cm\,s^{-1}}}\right)^2 \ {\rm eV}, \end{equation} \tag{ 37 }$$

where T is the time taken by a single shock propagating from the disk to the top of the bubble. However, this estimate assumes that particles are only accelerated by a single shock. As described in the previous section, the high-energy particles are accelerated by multiple shocks in the bubble; in other words, particles can diffuse downward and continue to be accelerated by younger shocks in the lower part of the bubble. Therefore T ∼ 10 Myr is only a minimum lifetime, which means that the above estimate can be considered as the lower limit of E₂. Another possible way to restrict the maximum value of E₂ is when the Larmor radius of the particles (r_L = E/eB) is larger than the radius of the bubble (H/2),

$$\begin{equation} (E_2)_{\rm max}\approx 10^{19}\left({B\over 5\,\mu {\rm G}}\right) \left({H\over 10\,{\rm kpc}}\right)\ {\rm eV}. \end{equation} \tag{ 38 }$$

In fact, this estimate should be more appropriate for E₂ and we will use it in our model fitting process.

In our model particle spectrum there are four parameters: the lower cutoff E₁ (Equation (12)), the upper cutoff E₂ (Equation (38)), the spectral index ν (Equation (36)), and the normalization ${\dot{N}}_0$ (Equation (35)). We consider l_sh = 30 pc, B = 5 μG, u = 10⁸ cm s⁻¹, and H = 10 kpc to be very reasonable mean values in the bubble, and therefore in our model fitting we fix E₁ = 10¹⁵ eV and E₂ = 10¹⁹ eV. On the other hand, the conversion efficiency from the shock energy into particle energy and the mean diffusion coefficient D_B in the bubble are the most uncertain parameters. Therefore, in our model fitting we treat these two parameters as fitting parameters. The solid line in Figure 4 indicates the best fit. The best fit gives ν = 3.12, which corresponds to D_B ∼ 3 × 10³⁰ cm² s⁻¹. This is in agreement with the estimation D_B ∼ cl_sh/3 if we take the average separation between shocks as l_sh ∼ 100 pc (see Equation (9)). This seemingly large diffusion coefficient is in fact at least an order of magnitude smaller than the coefficient in the halo for E > 5 × 10¹⁵ eV (e.g., Jones et al. 2001; suggested D_G ≈ 2.0 × 10²⁸(E/4.9 GeV)^0.54 cm² s⁻¹).

Zoom In Zoom Out Reset image size

If the CRs in the energy range 10¹⁵ eV <E < 10¹⁹ eV arriving at the Earth come from the bubble, then the power provided by the bubble for CR in this energy range is given by $\dot{W}_{{\rm CR}} \approx \int _{10^{15}{\rm \,eV}}^{10^{19}{\rm \,eV}}4\pi R^2 F_{{\rm CR}}(E)dE \sim 3\times 10^{39}(R/10\,{\rm kpc})^2$ erg s⁻¹, where F_CR(E) is the observed CR energy flux and R is the mean distance to the bubble. We find that the conversion efficient from shock power, $\dot{W}\sim W/\tau _{\rm cap}\sim 3 \times 10^{40}$ erg s⁻¹, is about 10%, which is consistent with recent estimation by using the Fermi-LAT data (Abdo et al. 2010). Realistically, particles could escape from the bubble through various locations of the bubble's surface, where the local strength of the magnetic field may be different. In addition, we have pointed out that there are two possible ways to estimate the value of E₂ (Equations (37) and (38)). Therefore, it is likely that E₁ and E₂ should have some distribution coming out from the bubble. In the dotted line in Figure 4 we keep the spectral index as that of the solid line but assume that E₁ and E₂ have a uniform distribution between 2 × 10¹⁴ eV <E₁ < 2 × 10¹⁵ eV (which corresponds to 1 μG <B < 10 μG) and (E₂)_min < E₂ < (E₂)_max, respectively. We can see that the model curve starts to drop around 3 × 10¹⁸ eV and matches the data better. Furthermore, we have argued that the injected spectrum from the shock is E⁻² for E < E₁ and slowly becomes harder for E⁻¹ due to interactions with multiple shocks after a sufficiently long time. However, it is unclear if this spectral hardening process can be completed within the finite lifetime of the particles in the bubble (see Section 4 for a more detailed discussion). Therefore we cannot predict the exact spectral index for E < E₁. For reference purposes, in Figure 4 we show these two possibilities: the dashed line for E⁻² and the dash-dotted line for E⁻¹, respectively. Both E₁ and E₂ of these two lines also assume uniform distributions.

The observed high-energy CR spectrum is a broken power law. The energy spectral index for CR energy smaller than the knee (around 3 × 10¹⁵ eV) is 2.7. The index increases to 3.1 for larger energy (until around 10¹⁹ eV where the spectrum becomes harder again). The general argument is that the sources (i.e., acceleration mechanisms and/or sites) responsible for the energy range within and beyond the knee should be different. The intriguing fact is that the two power laws match quite well at the knee. This coincident problem is difficult to solve if the two sources are totally unrelated.

CRs within (i.e., energy less than) the knee are generally attributed to SNRs in our Galaxy (see Section 2). In the previous two sections, we alluded to the acceleration site of CRs beyond the knee (i.e., energy larger than the knee) to the Fermi Bubble. In Section 4, we studied Equation (13) without a source in the bubble. The solution gave the characteristic spectrum of the system. In reality we need a source or a seed population for the bubble. We deem that we do not need a new seed population for bubble acceleration. Instead we propose the following model: some of the CRs produced by SNRs in the Galactic disk are re-accelerated in the bubble to energy beyond the knee. This model has the potential to solve the coincident problem naturally because the source of CRs beyond the knee is seeded by the source within the knee.

To consolidate our idea, we work out a concrete numerical model. Essentially, we solve the stationary state CR transport equation (13) 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 of ρ_B = 3 kpc.

The Fermi Bubbles are filled with shocks as described in Section 3 (see Figure 1). The spatial diffusion coefficient are different inside and outside the bubble as described by Equation (15). Due to the very turbulent conditions inside the bubble, we consider a constant spatial diffusion coefficient and adopt D_B = 2.08 × 10³⁰ cm² s⁻¹ (cf. estimate value by the fitting process in Section 5). Outside the bubble, we take into account the energy (or momentum) dependence of the spatial diffusion coefficient and adopt D_G = D₀(pc/4 GeV)^0.6, D₀ = 6.2 × 10²⁸ cm² s⁻¹ (cf. Jones et al. 2001).

According to the analysis in Section 4, the acceleration of energetic particles in the bubble is facilitated by stochastic acceleration (second-order Fermi acceleration). Assuming that there is little or no stochastic acceleration outside the bubble, we model the momentum diffusion coefficient as a step function as in Equation (16) and adopt κ_BH²/D_B = 1.9 (i.e., κ_B = 4.4 × 10⁻¹⁵ s⁻¹ or the corresponding acceleration timescale is 7.6 Myr).

The Galactic disk contains SNRs. We adopt the distribution suggested by Stecker & Jones (1977) and modified it with a Gaussian thickness profile

$$\begin{equation} Q_{\rm SNR}(\rho,z)\propto \left({\rho \over R_\odot }\right)^{1.2} \exp \left({-}{3.22\rho \over R_\odot }\right)\exp \left({-}{z^2\over h^2}\right), \end{equation} \tag{ 39 }$$

where ρ is the galactocentric radius and z is the distance perpendicular to the midplane. 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¹⁵ eV. Therefore, together with the SNR distribution (Equation (39)), the source function is

$$\begin{eqnarray} Q(\rho,z,p)&=&Q_0 \left({p\over p_{\rm max}}\right)^{-\mu }\exp \left({-}{p\over p_{\rm max}}\right) \left({\rho \over R_\odot }\right)^{1.2}\nonumber\\ &&\times\, \exp \left({-}{3.22\rho \over R_\odot }\right)\exp \left({-}{z^2\over h^2}\right). \end{eqnarray} \tag{ 40 }$$

As mentioned in Section 2, we take the SNR injection spectrum to be μ = 4.35 (see also Biermann & Strom 1993). The normalization Q₀ is obtained by fitting the simulation result to the observed spectrum and the value is 1.5 × 10¹⁴ particles s⁻¹ kpc⁻³ (GeV/c)⁻³. Integrating Q(ρ, z, p) over the Galaxy and momentum (from 1 GeV c⁻¹ to p_max) gives the total luminosity of CRs 4 × 10⁴⁰ erg s⁻¹, which is consistent with the value in the literature (e.g., Berezinskii et al. 1990).

Finally, the appropriate boundary conditions for the momentum coordinate are

$$\begin{equation} \left.{p\over f}{\partial f\over \partial p}\right|_{p=p_{\rm low}}= -4.7,\quad f|_{p=p_{\rm up}} = 0, \end{equation} \tag{ 41 }$$

where the energy of the lower momentum boundary is p_lowc = 10¹² eV, and the upper momentum boundary is p_upc = 3 × 10¹⁸ eV. The condition at the lower momentum ensures that the spectral index matches that of low-energy CRs (say E < 10¹² eV).

The spatial boundary conditions are

$$\begin{equation} \left.{\partial f\over \partial \rho }\right|_{\rho =0} = 0,\quad f|_{\rho =\rho _G} = 0, \end{equation} \tag{ 42 }$$

where the radius of the Galactic disk was taken to be ρ_G = 20 kpc:

$$\begin{equation} \left.{\partial f\over \partial z}\right|_{z=0} = 0,\quad f|_{z=\pm H} = 0, \end{equation} \tag{ 43 }$$

where the height of the halo H = 10 kpc.

The spectrum evaluated at Earth's position is the solid line shown in Figure 5. The model fits the data reasonably well and it is not coincident that the spectra join smoothly at the knee. For reference, in Figures 6 and 7 we show the spatial distribution of CRs from our simulation. The labels "Low" and "High" refer to the number density of particles in the energy range 1 × 10¹³ ∼ 3 × 10¹⁵ eV and 3 × 10¹⁵ ∼ 3 × 10¹⁸ eV, respectively. Figure 8 is a contour plot of the number density distribution of re-accelerated CRs (E > 3 × 10¹⁵ eV) (thick lines). In the figure we also plot the distribution of seed particles from SNRs (thin lines). We point out that the spatial distributions of seed and re-accelerated CRs in the disk are quite different. In principle, CR distribution can be derived from gamma-ray data. For instance, Breitschwerdt et al. (2002) used the gradient of gamma-ray emissivity in the disk derived from the EGRET data for the analysis of CR propagation in the Galaxy. If the diffuse gamma-ray data at E > 10¹⁵ eV were available, the gradient test would be a nice tool to investigate possible proton sources in this energy range and might lend support to our model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have summarized our current understanding of the origin of CRs. It is generally believed that most CR power can be provided by SNRs. However, the CRs with energies E > 10¹⁵ eV are quite difficult to achieve in SNe due to the limited acceleration time and energy content in SN shocks. On the other hand, we argue that shocks in the Fermi Bubble produced by stellar capture events can have a much longer lifetime >10⁷ yr and larger energy content ∼3 × 10⁵² erg, which allow them to produce CRs with energies E > 10¹⁵ eV. If processes of CRs which escape from the Galaxy are taken into account, the predicted CR spectrum contributed by the bubble is E^−ν, where

$$\begin{eqnarray*} \nu &=&-\frac{1}{2}+\sqrt{\frac{9}{4}+10\frac{(D/3\times 10^{30}\,{\rm cm}^2\,{\rm s^{-1}})^2}{(u/10^8\,{\rm cm\,s^{-1}})^2(H/10\,{\rm kpc})^2}}\sim 3 \end{eqnarray*}$$

for 10¹⁵ eV < E < 10¹⁹ eV.

However, it is very difficult to predict the exact value of the spectral index ν due to the poorly measured value of the diffusion coefficient in the bubble. So we fit the observed CR spectrum between 10¹⁵ eV and ∼10¹⁹ eV and find that the diffusion coefficient is about 3 × 10³⁰ cm² s⁻¹. By matching with the observed flux in the knee region we find that the conversion efficiency from shocks in the bubble to CRs is about 10%, which is quite consistent with numerical simulations. Other input parameters in this model, such as the capture timescale τ_cap ∼ 3 × 10⁴ yr, the mean rate of energy release $\dot{W}\sim 3 \times 10^{40}$ erg s⁻¹ per capture, and the injected wind speed u ∼ 10⁸ cm s⁻¹, have been estimated and used in other observed phenomena in the GC (e.g., Cheng et al. 2006, 2007; Dogiel et al. 2009a, 2009b, 2009c, 2011).

We put forth the idea that the Fermi Bubble acts as the re-acceleration site for the CRs produced by SNRs in the Galactic disk. The re-accelerated particles form the part of the observed spectrum that is beyond the knee (about 3 × 10¹⁵ eV). The part within the knee is formed by the CRs produced by the SNRs. We demonstrated this idea by a numerical model. We solve the stationary transport equation in our Galaxy with two Fermi Bubbles (see Section 3 for our model of the bubbles). The re-acceleration process in the bubble is facilitated by stochastic acceleration. Our model simulated the observed spectrum nicely. Therefore we consider this model to provide a natural explanation of the flux, spectral index, and matching at the knee of CRs in this energy range. In a related issue, Mertsch & Sarkar (2011) showed that the gamma ray from the bubbles can be produced by the stochastic acceleration of electrons throughout the bubbles.

As described in Section 3 there are many shocks propagating in the Fermi Bubbles. After being re-accelerated inside the bubbles by the multiple shocks, protons (and nuclei) escape the bubbles. The lifetime of protons by the pp collision is of the order of 10⁹ ∼ 10¹⁰ yr (e.g., Crocker & Aharonian 2011). The diffusion time for protons at these energies to escape the Galaxy is of the order of 10⁷ yr. Thus, after leaving the bubbles, protons (and nuclei as well) diffuse throughout the whole Galaxy (including the Earth) without any attenuation of energy. Not only can protons (and nuclei) be accelerated by the multiple shocks, but electrons can as well. Nevertheless, energetic electrons lose energy efficiently. The lifetime of electrons can be estimated by τ_e = 1/β_eE (e.g., inverse Compton and synchrotron). Taking β_e = 3 × 10⁻²⁵ eV⁻¹ s⁻¹, the lifetime of electrons of energy 0.1 ∼ 1 TeV is 1 ∼ 0.1 Myr. Once they leave the acceleration site (i.e., Fermi Bubbles), they lose most of their energy within a short distance (less than 1 kpc). Hence, the electrons are mostly confined in the bubbles.

Su et al. (2010) suggested that the electron spectrum must be E^−α for E < TeV and the spectral index α ∼ 2.4–2.8. It is clear that the energy loss processes for protons and electrons are very different; therefore, protons and electrons in the bubble can have different spectra. In the CCDKI model we have assumed that the shocks can produce an injected electron spectrum ∼E⁻², in which the spectrum is modified by processes of energy losses and escape. In our subsequent work we will take into account the stochastic acceleration processes (multiple-shock) to see how the electron spectrum is affected. We will also further analyze the spatial spectral distribution of electrons. Preliminary results can be found in Chernyshov et al. (2011), in which the spatial distribution of gamma-ray emission is reproduced nicely by the multi-shock model.

We are very grateful to A. M. Bykov, R. M. Crocker, A. D. Erlykin, Y. Uchiyama, and V. N. Zirakashvili for useful discussions, and to K. F. Sinkov who performed graphic and numerical illustrations of the Kompaneets solution presented in Section 3. D.O.C. and V.A.D. are partly supported by the NSC-RFBR Joint Research Project RP09N04 and 09-02-92000-HHC-a. K.S.C. is supported by the GRF Grants of the Government of the Hong Kong SAR under HKU 7011/10P. C.M.K. is supported by the Taiwan National Science Council Grants NSC 98-2923-M-008-01-MY3 and NSC 99-2112-M-008-015-MY3. W.H.I. is supported by the Taiwan National Science Council Grant NSC 97-2112-M-008-011-MY3 and Taiwan Ministry of Education under the Aim for Top University Program National Central University.