EDP Sciences
Home All issues Volume 513 (April 2010) A&A, 513 (2010) A17 Full HTML
Free Access
Issue
A&A
Volume 513, April 2010
Article Number A17
Number of page(s) 10
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/200913495
Published online 15 April 2010
A&A 513, A17 (2010)

Escape-limited model of cosmic-ray acceleration revisited

Y. Ohira1 - K. Murase2 - R. Yamazaki3

1 - Department of Earth and Space Science, Osaka University, Toyonaka 560-0043, Japan
2 - Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan
3 - Department of Physical Science, Hiroshima University, Higashi-Hiroshima 739-8526, Japan

Received 19 October 2009 / Accepted 10 January 2010

Abstract
Context. The spectrum of cosmic rays (CRs) is affected by their escape from an acceleration site. This may be observed not only in the gamma-ray spectrum of young supernova remnants (SNRs) such as RX J1713.7-3946, but also in the spectrum of CRs showering the Earth.
Aims. The escape-limited model of cosmic-ray acceleration is studied in general. We discuss the spectrum of runaway CRs from the acceleration site. The model will also be able to constrain the spectral index at the acceleration site and the ansatz with respect to the unknown injection process into the particle acceleration.
Methods. Our methods are analytical derivations. We apply our model to CR acceleration in SNRs and in active galactic nuclei (AGN), which are plausible candidates of Galactic and extragalactic CRs, respectively. In particular we take account into the shock evolution with cooling of escaping CRs in the Sedov phase for young SNRs.
Results. The spectrum of escaping CRs generally depends on the physical quantities at the acceleration site like the spectral index, the evolution of the maximum energy of CRs and the evolution of the normalization factor of the spectrum. It is found that the spectrum of runaway particles can be both softer and harder than that of the acceleration site.
Conclusions. The model explains spectral indices of both Galactic and extragalactic CRs produced by SNRs and AGNs, respectively, suggesting the unified picture of CR acceleration.

Key words: acceleration of particles - cosmic rays - ISM: supernova remnants - galaxies: jets

1 Introduction

The origin of cosmic rays (CRs) has been one of the long-standing problems. The number spectrum of nuclear CRs observed at the Earth, $\mathcal{N} (E) \propto E^{-s}$, shows a break at the ``knee'' energy ( ${\approx}10^{15.5}$ eV), below which the spectral index is about $s \approx -2.7$ (Cronin 1999). Because of the energy-dependent propagation of CRs, the spectral shape at the source is different from that observed on Earth. Taking into account the propagation effect, the source spectral index has been well constrained as $s \approx 2.2{-}2.4$ in various models (e.g., Strong & Moskalenko 1998; Putze et al. 2009). This value of s has been also inferred to explain the Galactic diffuse gamma-ray emission (e.g., Strong et al. 2000). This may give us valuable insights into the acceleration mechanism of CRs.

Mechanisms of CR acceleration have also been studied for a long time, and the most plausible process is a diffusive shock acceleration (DSA) (Blandford & Ostriker 1978; Krymsky 1977; Axford et al. 1977; Bell 1978). Very high-energy gamma-ray observations have revealed the existence of high-energy particles at the shock of young supernova remnants (SNRs), which supports the DSA mechanism as well as the paradigm that the Galactic CRs are produced by young SNRs (e.g., Enomoto et al. 2002; Katagiri et al. 2005; Aharonian et al. 2005,2004). Recent progress of the theory of DSA revealed that the back-reactions of accelerated CRs are important if a large number of nuclear particles are accelerated (Malkov & Drury 2001; Drury & Völk 1981). There are several observational facts which are consistent with the predictions of such a nonlinear model (Bamba et al. 2005b; Warren et al. 2005; Bamba et al. 2005a; Uchiyama et al. 2007; Bamba et al. 2003; Helder et al. 2009; Vink & Laming 2003). The model predicts, however, the harder spectrum of accelerated particles at the shock than $f (p) \propto p^{-4}$ (corresponding to s=2) where p is the momentum of CRs, in particular, near the knee energy (Malkov 1997; Kang et al. 2001; Berezhko & Ellison 1999). This apparently contradicts the source spectral index of $s \approx 2.2{-}2.4$ inferred from the CR spectrum on Earth. Even in the test-particle limit of DSA, such a soft source spectrum requires a shock with the small Mach number ($M \la10$), which is unexpected for young SNRs.

There are several models of DSA, depending on the boundary conditions imposed. Different models predict different spectra of CRs dispersed from the shock region. So far, the age-limited acceleration has been frequently considered as a representative case (Sect. 2). In this model, all the particles are stored around the shock while accelerated. When the confinement becomes inefficient, all the particles run away from the region at a time. Then the source spectrum of CRs which has just escaped from the acceleration region is expected to be the same as that at the shock front. Therefore, this model predicts that the source spectrum is the same as that of accelerated particles, which is typically harder than the observed one. In this paper, we consider an alternative model, the escape-limited acceleration, to explain the observed CR spectrum on Earth (Sect. 3). This model is preferable to the age-limited acceleration when we consider observational results for young SNR RX J1713.7-3946, of which TeV $\gamma$-ray emission is more precisely measured than any others (Sect. 4.1).

The nature of CRs with energies much higher than the knee energy is also still uncertain. While CRs below the second knee ( ${\sim}{10}^{17.5}$ eV) may be of Galactic origin, the highest energy CRs above ${\sim}{10}^{18.5}$ eV are believed to be extragalactic. Possible candidates are active galactic nuclei (AGNs) (e.g., Rachen & Biermann 1993; Takahara 1990; Pe'er et al. 2009; Biermann & Strittmatter 1987), gamma-ray bursts (Murase et al. 2006; Waxman 1995; Vietri 1995), magnetars (Arons 2003; Murase et al. 2009) and clusters of galaxies (Inoue et al. 2007; Kang et al. 1997). The intermediate energy range from ${\sim}{10}^{17.5}$ eV to ${\sim}{10}^{18.5}$ eV is more uncertain. Both the Galactic and extragalactic origins are possible and it may just be a transition between the two. As the extragalactic origin, AGNs (Berezinsky et al. 2006), clusters of galaxies (Murase et al. 2008) and hypernovae (Wang et al. 2007) have been proposed so far.

Among these possibilities, AGN are one of the most plausible candidates for accelerators of high-energy CRs, because they can explain the ultra high energy cosmic-ray (UHECR) spectrum above ${\sim}{10}^{17.5}$ eV assuming the proton composition. In such a proton dip model, the source spectrum of UHECRs is $\mathcal{N} \propto E^{-s}$ with s=2.4-2.7, depending on the models of the source evolution (Berezinsky et al. 2006). The required source spectral index of s=2.4-2.7 can be explained by several possibilities. First, it can be attributed to the acceleration mechanism itself. One can consider non-Fermi acceleration mechanisms (Berezinsky et al. 2006) or the two-step diffusive shock acceleration in two different shocks (Aloisio et al. 2007). Second, the index can be attributed to a superposition of many AGNs with different maximum energies, and one can suppose that AGNs with different luminosities may have different maximum energies (Kachelriess & Semikoz 2006). Recently, Berezhko (2008) proposed another possibility with the cocoon shock model. In this cocoon shock scenario, different maximum energies can be interpreted as maximum energies of escaping particles at different ages of AGN jets. Although it is very uncertain whether efficient CR acceleration occurs there, this scenario would also be one of the possibilities to be investigated in detail.

The organization of the paper is as follows. After the brief introduction of the age-limited model of the CR acceleration (Sect. 2), we study the escape-limited model in Sect. 3. For a simple understanding, the general argument in a stationary, test-particle approximation is given in Sect. 3.1. Then we derive the formulae of the maximum energy of accelerated particles in Sect. 3.2 and of the spectrum of escaping particles in Sect. 3.3. We consider the applications to young SNR and AGN in Sects. 4 and 5, respectively. Section 6 is devoted to a discussion.

2 Maximum attainable momentum in the age-limited acceleration

For comparison with the escape-limited acceleration, we briefly summarize the case of the age-limited acceleration. In this case, the maximum momentum of accelerated particles $p_{\rm m}^{\rm (age)}$ is determined by $t_{\rm acc}=t_{\rm age}$, where $t_{\rm age}$ and $t_{\rm acc}$ are the age of the shock and the acceleration time scale, respectively. When we consider DSA, $t_{\rm acc}$ is given by (Drury 1983)

\begin{displaymath}% t_{\rm acc} = \frac{3}{u_1-u_2} \left(\frac{D_1(p)}{u_1}+\frac{D_2(p)}{u_2}\right), \end{displaymath} (1)

where D(p) and u are the diffusion coefficient as a function of the momentum of accelerated particles and the velocity of the background fluid, respectively. Subscripts 1 and 2 represent upstream and downstream regions, respectively. For simplicity, we assume the Bohm-type diffusion, i.e.,

\begin{displaymath}% D_1(p)=D_2(p)= \frac{\eta_{\rm g}m_{\rm p}c^3}{3eB}\left(\frac{p}{m_{\rm p}c}\right), \end{displaymath} (2)

where B, $\eta_{\rm g}$ and $m_{\rm p}$are the magnetic field strength, the gyro-factor and the proton mass, respectively. Taking into account that the fluid velocities are related to the shock velocity,  $u_{\rm sh}$, as  $u_1 = \sigma u_2=u_{\rm sh}$, where $\sigma$ is the compression ratio of the shock, we derive (e.g., Aharonian & Atoyan 1999)

\begin{displaymath}% p_{\rm m}^{\rm (age)} = \frac{\sigma-1}{\sigma(\sigma+1)}\frac{eB}{\eta_{\rm g}c^2}u_{\rm sh}^2t_{\rm age}. \end{displaymath} (3)

3 Escape-limited acceleration

In the framework of DSA, accelerated particles are scattered by the turbulent magnetic field and go back and forth across the shock front. Upstream turbulence may be excited by the accelerated particles themselves (Bell 1978), and the magnetic field strength of this turbulence is theoretically expected to be strong (e.g., Lucek & Bell 2000). There are observational evidences suggesting that CRs are responsible for substantial amplification of the ambient magnetic field in the precursors of shock fronts in SNRs and that such magnetic turbulence well confines the particles around the shock front (Yamazaki et al. 2004; Bamba et al. 2005b,a; Uchiyama et al. 2007; Bamba et al. 2003; Vink & Laming 2003; Parizot et al. 2006), leading to the efficient CR acceleration.

The spectrum of accelerated particles is affected by the spatial and spectral structures of the magnetic turbulence through the process in which the particles escape from the shock toward far upstream regions. There are mainly two scenarios of the escape model considered so far; one causes the effect on the boundary in the momentum space, and the other causes the effect on the spatial boundary. The former comes from significant decay of the wave amplitude below the wave number $k_{\min}$ of the turbulence spectrum (Drury et al. 2009; Reynolds 1998). Particles with the Lorentz factor above  $\gamma_{\rm c}$ satisfying the resonance condition, $\omega-k_{\min}v-\Omega_{\rm c}(\gamma_{\rm c})=0$, where $\Omega_{\rm c}(\gamma)=eB/\gamma m_{{\rm p}}c$ is the cyclotron frequency, are not confined around the shock front and escape into far upstream regions. In this context, the escape flux was calculated previously (e.g., Ptuskin & Zirakashvili 2005; Drury et al. 2009). The latter effect has been recently discussed by several authors (Ptuskin & Zirakashvili 2005; Reville et al. 2009; Caprioli et al. 2009). The turbulence generation may be connected with the flux of accelerated particles themselves. Hence, in the region far from the shock front, the flux of high-energy particles is small and wave excitation is less significant. If the accelerated particles reach the region, they are dispersed into the far upstream region. Let $\ell$ be the distance from the shock beyond which the amplitude of the upstream turbulence becomes negligible. Characteristic spatial length of particles penetrating into the upstream region is given by  $D(p)/u_{\rm sh}$. As long as $D(p)/u_{\rm sh}\ll\ell$, the particles are confined without the significant escape loss, and they are accelerated to higher energies. On the other hand, when their momentum increases up to sufficiently high energies satisfying $D(p)/u_{\rm sh}\ga \ell$, their acceleration ceases and they escape into the far upstream. Therefore, the maximum momentum of accelerated particles in this scenario is given by the condition  $D(p)/u_{\rm sh}\sim\ell$. Below we consider the escape-limited model, where the maximum energy is essentially determined by  $D(p)/u_{\rm sh} \la \ell$.

3.1 A simple case of stationary, test-particle approximation

In order to take an essential feature of the escape-limited acceleration, we calculated the escape flux and the maximum attainable energy of accelerated particles for the simplest case (see also Caprioli et al. 2009). Let us consider the stationary transport equation

\begin{displaymath}% u(x)\frac{\partial f}{\partial x}= \frac{\partial}{\partial... ...3}\frac{{\rm d} u}{{\rm d} x}\frac{\partial f}{\partial p} +Q, \end{displaymath} (4)

with the boundary condition $f(x=-\ell)=0$, where $x=-\ell$ ($\ell>0$) is the upstream escape boundary. The fluid velocity u(x) is given by

\begin{displaymath}% u(x) = \left\{ \begin{array} {c@{} l@{}} u_1 &\ \ (x<0) \\ u_2 &\ \ (x>0), \end{array} \right. \end{displaymath} (5)

where u1 and u2 are constants. The solution to the transport equation in the test-particle approximation is derived as (Caprioli et al. 2009)

\begin{displaymath}% f(x,p)=f_0(p) \frac{\exp\left[\frac{u_1x}{D(p)}\right]-\exp... ...ell}{D(p)}\right]} {1-\exp\left[-\frac{u_1\ell}{D(p)}\right]}, \end{displaymath} (6)

where f0(p)=f(x=0,p) is given by

\begin{displaymath}% f_0(p)=K\exp\left\{ -q\int_{p_{\rm inj}}^p\frac{{\rm d}\log p'} {1-\exp\left[-\frac{u_1\ell}{D(p')}\right]} \right\}, \end{displaymath} (7)

and q=3u1/(u1-u2). The escape flux at $x=-\ell$ is
                   $\displaystyle % \phi(p)$ = $\displaystyle u_1 \left. f \right\vert _{x=-\ell}-D(p) \left. \frac{\partial f}{\partial x} \right\vert _{x=-\ell}$  
  = $\displaystyle \frac{u_1f_0(p)}{1-\exp\left[\frac{u_1\ell}{D(p)}\right]}\cdot$ (8)

Let us introduce a new variable $y=\ln p$ and a new function $\Phi(y)=\ln\phi(p)$. Then we expand $\Phi(y)$ around its maximum value at  $y=y_{\rm m}$. To do this, we calculate the first and the second derivatives as

\begin{displaymath}% \frac{{\rm d}\Phi}{{\rm d} y} = -\frac{1}{1-\exp\left[-\fra... ...t[q-\frac{u_1\ell}{D(y)}\frac{{\rm d}\ln D}{{\rm d} y}\right], \end{displaymath} (9)


                              $\displaystyle % \frac{{\rm d}^2\Phi}{{\rm d} y^2}$ = $\displaystyle \frac{\frac{u_1\ell}{D(y)}} {\left\{1-\exp\left[-\frac{u_1\ell}{D... ...)}\frac{{\rm d}\ln D}{{\rm d} y} \right] \frac{{\rm d}\ln D}{{\rm d} y} \right.$  
    $\displaystyle \left. + ~ \left[1-\exp\left[-\frac{u_1\ell}{D(y)}\right]\right] ... ...}\ln D}{{\rm d} y}\right)^2- \frac{\partial^2\ln D}{\partial y^2} \right\}\cdot$ (10)

Below we consider the case of Bohm diffusion, $D(p)=D_0(p/m_{\rm p}c)$. Then one can find

\begin{displaymath}% p_{\rm m} = e^{y_{\rm m}} = \frac{u_1\ell m_{\rm p}c}{qD_0}, \end{displaymath} (11)

(where ${\rm d} \Phi/{\rm d} y=0$) so that we obtain
                       $\displaystyle % \Phi(y)$ = $\displaystyle \Phi(y_{\rm m}) + \frac{1}{2}\frac{\partial^2\Phi}{\partial y^2}(y_{\rm m})(y-y_{\rm m})^2 + \cdots,$  
  = $\displaystyle \ln\left[ \frac{u_1f_0(p_{\rm m})}{1-e^{-q}}\right] -\left(\frac{y-y_{\rm m}}{\xi}\right)^2 + \cdots,$ (12)

where $\xi^2 = 2(1-e^{-q})/q<1$. Note that the quantity $p_{\rm m}$ given by Eq. (11) also plays the role of maximum momentum of the accelerated particles at the shock (x=0) because one can see $f_0(p)\propto p^{-q}$ for $p\ll p_{\rm m}$ while $f_0(p)\propto\exp(-p/p_{\rm m})$ for $p_{\rm m}\ll p$ (Caprioli et al. 2009).

Finally, going back to the function $\phi(p)$, we obtain

\begin{displaymath}% \phi(p) = \frac{u_1f_0(p_{\rm m})}{1-e^{-q}} \exp\left[-\left(\frac{\ln p-\ln p_{\rm m}}{\xi}\right)^2 +\cdots\right]. \end{displaymath} (13)

One can clearly see from Eq. (13) that particles with momentum around $p_{\rm m}$ escape the shock region most efficiently.

3.2 The maximum energy of accelerated particles

We have seen in Sect. 3.1 that in the stationary, test-particle case the quantity $p_{\rm m}$ given by Eq. (11) plays the role of maximum momentum of the accelerated particles at the shock. Taking this into account, we assume that in the more general escape-limited case the maximum momentum  $p_{\rm m}^{\rm (esc)}$ is determined by

\begin{displaymath}% q \frac{D(p_{\rm m}^{\rm (esc)})}{u_1}=\ell. \end{displaymath} (14)

Given that

\begin{displaymath}% D(p)= \frac{\eta_{\rm g}m_{\rm p}c^3}{3eB}\left(\frac{p}{m_{\rm p}c}\right), \end{displaymath} (15)

which is the same as in the age-limited case (Eq. (2) in Sect. 2), we obtain

\begin{displaymath}% p_{\rm m}^{\rm (esc)} = \frac{3}{q}\frac{eB}{\eta_{\rm g}c^2}u_{\rm sh} \ell. \end{displaymath} (16)

Since $q=3\sigma/(\sigma-1)$, we find from Eqs. (3) and (16)

\begin{displaymath}% \frac{p_{\rm m}^{\rm (esc)}}{p_{\rm m}^{\rm (age)}} =(\sigma+1) \frac{\ell}{u_{\rm sh}t_{\rm age}}\cdot \end{displaymath} (17)

Hence, as long as $\ell\gg u_{\rm sh}t_{\rm age}$, which is expected in the early phase of the shock, we do not need the escape-limited acceleration. However, if the shock evolves so that $\ell\ll u_{\rm sh}t_{\rm age}$, we have to consider the effect of the particle escape and cannot simply apply the well-known result of age-limited acceleration, Eq. (3).

3.3 The spectrum of CRs dispersed from an accelerator

In this subsection, we derive the time-integrated spectrum of CR particles which is dispersed from an accelerator. The derivation is essentially identical to that of Ptuskin & Zirakashvili (2005). However, our argument is simpler and more general, so the final form of the spectrum (Eqs. (27) and (28)) is more general. Note that our formalism is applicable not only to DSA but also to arbitrary acceleration processes.

The proton production rate $N(p,\chi){\rm d}\chi {\rm d} p$, at a certain epoch labeled by a parameter $\chi$, is defined as the number of protons with momentum between p and $p+{\rm d}p$, which is produced in the interval between $\chi$ and  $\chi+{\rm d}\chi$. Here $\chi$ is the parameter which describes the dynamical evolution of the accelerator - it can be either simply the age or the position of the shock front. It is expected that $N(p,\chi)$ contains the term of exponential cutoff at the momentum  $p_{\rm m}(\chi)$ which depends on $\chi$ (see, for example, Eq. (25)). The number of protons with momentum between p and  $p+{\rm d}p$ which is escaping from the accelerator at the epoch between $\chi$ and  ${\rm d}\chi$, is denoted by $\phi(p,\chi){\rm d}\chi {\rm d} p$, and we assume

\begin{displaymath}% 4 \pi p^2 \phi(p,\chi) {\rm d}\chi {\rm d} p \propto N(p_{\rm m}(\chi),\chi)G(p,\chi){\rm d}\chi {\rm d} p, \end{displaymath} (18)

where

\begin{displaymath}% G(p,\chi)= \exp\left[-\left(\frac{\ln p-\ln p_{\rm m}(\chi)}{\xi}\right)^2\right], \end{displaymath} (19)

and $\xi<1$. This is because we expect that particles with momentum around $p_{\rm m}(\chi)$ are escaping most efficiently from the source. Indeed, as shown in Sect. 3.1, Eqs. (18) and (19) are a good approximation in the test-particle, stationary case. Here, we simply expect these assumptions are also correct in general.

The time-integrated spectrum of protons which have escaped at the source  $N_{\rm esc}(p)$ is obtained by

\begin{displaymath}% N_{\rm esc}(p) =4 \pi p^2 \int \phi(p,\chi){\rm d}\chi. \end{displaymath} (20)

In order to derive a simple analytical form, we approximate Eq. (19) as

\begin{displaymath}% G(p,\chi)\approx\sqrt{\pi}\xi\cdot \delta(\ln p-\ln p_{\rm m}(\chi)). \end{displaymath} (21)

If we use a general mathematical formula for $\delta$-functions for an arbitrary function $g(\chi)$

\begin{displaymath}% \delta(g(\chi)) = \frac{\delta(\chi-\chi_0)} {\left[\frac{{\rm d} g}{{\rm d}\chi}\right]_{\chi=\chi_0}}, \end{displaymath} (22)

where $g(\chi_0)=0$, then Eq. (21) is rewritten as

\begin{displaymath}% G(p,\chi)\propto \frac{p_{\rm m}(\chi)} {\left[\frac{{\rm ... ...ight]_{\chi=p_{\rm m}^{-1}(p)}}\delta(\chi-p_{\rm m}^{-1}(p)), \end{displaymath} (23)

where $p_{\rm m}^{-1}(p)$ is the inverse function of $p_{\rm m}(\chi)$, that is, $p_{\rm m}(p_{\rm m}^{-1}(p))=p$ and $p_{\rm m}^{-1}(p_{\rm m}(\chi))=\chi$. Using Eqs. (18) and (23), we calculate $N_{\rm esc}(p)$ as

\begin{displaymath}% N_{\rm esc}(p) \propto \frac{p N(p,p_{\rm m}^{-1}(p))}{p_{\... ... p_{\rm m}}{{\rm d}\chi}\right]_{\chi=p_{\rm m}^{-1}(p)}}\cdot \end{displaymath} (24)

This is the most general analytical formula of the spectrum of protons dispersed from an accelerator.

For the remainder of the paper, the form of $N(p,\chi)$ is assumed to have

                $\displaystyle % N(p,\chi){\rm d}\chi {\rm d} p$ = $\displaystyle K(\chi)\left(\frac{p}{m_{\rm p}c}\right)^{-s}$  
    $\displaystyle \times ~ \exp\left[-\left(\frac{p}{p_{\rm m}(\chi)}\right)\right] {\rm d}\log\chi~{\rm d} p,$ (25)

which is a power law with the index s and the exponential cut off at  $p_{\rm m}(\chi)$. Substituting Eqs. (25) into (24), we obtain

\begin{displaymath}% N_{\rm esc}(p)\propto \frac{p^{1-s}K\left(p_{\rm m}^{-1}(p)... ... p_{\rm m}}{{\rm d}\chi}\right]_{\chi=p_{\rm m}^{-1}(p)}}\cdot \end{displaymath} (26)

In particular, if $K(\chi)$ and $p_{\rm m}(\chi)$ are written by the power law forms like $K(\chi)\propto\chi^\beta$ and $p_{\rm m}(\chi)\propto\chi^{-\alpha}$, then, $p_{\rm m}^{-1}(p)\propto p^{-1/\alpha}$, so that Eq. (26) becomes

\begin{displaymath}% N_{\rm esc}(p) \propto p^{-s_{\rm esc}}, \end{displaymath} (27)

where

\begin{displaymath}% s_{\rm esc} = s + \frac{\beta}{\alpha}\cdot \end{displaymath} (28)

This is the simplest form of the spectrum of CRs which are dispersed from the acceleration region.

Generally speaking, to obtain the energy spectrum of accelerated particles, time-dependent kinetic equation should be solved. Instead, we have assumed that at an arbitrary epoch the spectral form is given by Eq. (25). This assumption is justified if the spectrum at the given epoch is dominated by those which are being accelerated at that time, in other words, if the particle spectrum does not so much depend on the past acceleration history. For example, in the case of the spherical expansion, accelerated particles suffer adiabatic expansion after they are transported downstream of the shock and lose their energy (e.g., Yamazaki et al. 2006), so that the contribution of the previously accelerated particles is negligible. Strictly speaking, even if we consider the energy loss via adiabatic expansion, the energy spectrum of accelerated particles does depend on the past acceleration history in some cases. When we use the shock radius  $R_{\rm sh}$ as $\chi$, the final form of Eqs. (27) and (28) is correct as long as $\beta>\left(s-1\right)\left(\sigma^{-1}-1\right)$ (see Appendix), which is satisfied in the cases considered in Sects. 4 and 5. Otherwise, the form of Eq. (25) is no longer a good approximation, and the final form of  $s_{\rm esc}$ is different from Eq. (28) (see Appendix).

4 Application to young supernova remnants

4.1 Inconsistency of age-limited acceleration with observed results of RX J1713.7-3946

RX J1713.7-3946 is a representative SNR from which bright TeV $\gamma$-rays have been detected. The HESS experiment measured the TeV spectrum and claimed that its shape was better explained by the hadronic model (Aharonian et al. 2006,2007). Furthermore, evidences of amplified magnetic field ($B\ga0.1$ mG) are derived from the width of synchrotron X-ray filaments (Parizot et al. 2006; see also Bamba et al. 2005b,a,2003; Vink & Laming 2003) and from time variation of synchrotron X-ray hot spots (Uchiyama et al. 2007). This also supports the hadronic origin of TeV $\gamma$-rays, because the leptonic, one-zone emission model (e.g., Aharonian & Atoyan 1999) cannot explain the TeV-to-X-ray flux ratio. Hence it is natural to assume that the TeV $\gamma$-rays are produced by the hadronic process, although there are several arguments against this interpretation (Katz & Waxman 2008; Plaga 2008; Butt 2008).

In the age-limited case, Eq. (3) reads

\begin{displaymath}% cp_{\rm m}^{\rm (age)}=8\times10^3 \frac{B_{\rm mG}}{\eta_{... ...ht)^2 \left(\frac{t_{\rm age}}{10^3~{\rm yr}}\right){\rm TeV}, \end{displaymath} (29)

where we adopt $\sigma=4$, and $B_{\rm mG}$ is the magnetic field strength in units of mG.

HESS observation revealed that the cutoff energy of TeV $\gamma$-ray spectrum is low (Aharonian et al. 2006,2007), so that in the one-zone hadronic scenario the maximum energy of protons, $E_{\max,p}=cp_{\rm m}^{\rm (age)}$ is estimated as 30-100 TeV (Villante & Vissani 2007). If  $E_{\max,p}<100$ TeV and $B\approx1$ mG, then Eq. (29) tells us $\eta_{\rm g}\ga80$, implying far from the ``Bohm limit'' ( $\eta_{\rm g}\approx1$) which is inferred from the X-ray observation (Yamazaki et al. 2004; Parizot et al. 2006) or expected theoretically (Lucek & Bell 2000; Giacalone & Jokipii 2007; Reville et al. 2007; Ohira et al. 2009b; Bell 2004; Inoue et al. 2009; Ohira & Takahara 2009). This statement is recast if we involve recent results of X-ray observations. The precise X-ray spectrum of RX J1713.7-3946 is revealed, which gives $u_{\rm sh}=3.3$ $\times$ $10^8\eta_{\rm g}^{1/2}$ cm s-1 (Tanaka et al. 2008). Then Eq. (29) can be rewritten as (Yamazaki et al. 2009)

\begin{displaymath}% cp_{\rm m}^{\rm (age)}=5\times10^3 B_{\rm mG} \left(\frac{t_{\rm age}}{10^3{\rm yr}}\right){\rm TeV}. \end{displaymath} (30)

Hence, in order to obtain $cp_{\rm m}^{\rm (age)}<100$ TeV, we need $B\la20~\mu$G in the context of the hadronic scenario of TeV $\gamma$-rays, which contradicts current estimates of the magnetic field.

A possible solution is to consider the escape-limited acceleration. One can find from Eq. (17) that if we take $\ell\sim10^{-2}u_{\rm sh}t_{\rm age}\sim0.04$ pc, the maximum energy becomes

                          $\displaystyle % E_{\max,p}$ = $\displaystyle \min\left\{ cp_{\rm m}^{\rm (age)}, cp_{\rm m}^{\rm (esc)} \right\} = cp_{\rm m}^{\rm (esc)} \sim 100~{\rm TeV},$ (31)

which is consistent with the observed gamma-ray spectrum. Below we consider the model of escape-limited acceleration under simple assumptions, estimating the evolution of the number density and the maximum momentum of accelerated particles so as to discuss the spectral index  $s_{\rm esc}$ of  $N_{\rm esc}(p)$.

4.2 Evolution of pm

Time evolution of the maximum momentum of accelerated particle $p_{\rm m}$ has been so far discussed in many contexts (e.g., Ptuskin & Zirakashvili 2003). One way to estimate $p_{\rm m}$ is to use Eq. (16). In this approach, a key parameter is the magnetic field, which is likely amplified around the shock front (Bamba et al. 2005b; Yamazaki et al. 2004; Bamba et al. 2005a; Uchiyama et al. 2007; Bamba et al. 2003; Vink & Laming 2003) and may depend on various physical quantities like the shock velocity, the ambient density, and so on. At present, the evolution of the magnetic field is not well understood despite many works on the subject (e.g., Riquelme & Spitkovsky 2009; Niemiec et al. 2008; Ohira et al. 2009a; Luo & Melrose 2009). In addition, the evolution of another parameter  $\eta_{\rm g}$ is also unknown. This prevents us from predicting $p_{\rm m}$ rigorously.

Here we adopt a different phenomenological approach based on the assumption that young SNRs are responsible for observed CRs below the knee (Gabici et al. 2009). The maximum energy  $cp_{\rm m}$ is expected to increase up to the knee energy (1015.5 eV) until the end of the free expansion phase  $t_{\rm Sedov}$ and decreases from that epoch. As seen in Sect. 4.1, $p_{\rm m}$ is limited by the escape at $t>t_{\rm Sedov}$, that is $p_{\rm m}=p_{\rm m}^{\rm (esc)}$. Then, to reproduce the observed CR spectrum from $\sim$GeV to the knee, we may assume a functional form of

\begin{displaymath}% p_{\rm m}^{\rm (esc)}(t)=p_{\rm knee} \left(\frac{t}{t_{\rm Sedov}}\right)^{-6.5/a}, \end{displaymath} (32)

where $cp_{\rm knee}=10^{15.5}~{\rm eV}$, so that $cp_{\rm m}^{\rm (esc)}=1~{\rm GeV}$ at $t=10^at_{\rm Sedov}$. For later convenience, we change variables from t to the SNR radius  $R_{\rm sh}$ in order to take $\chi=R_{\rm sh}$, where $\chi$ is a variable which appeared in Sect. 3.3. We further assume the dynamics of  $R_{\rm sh}$ as

\begin{displaymath}% R_{\rm sh} = R_{\rm Sedov} \left(\frac{t}{t_{\rm Sedov}}\right)^b, \end{displaymath} (33)

where $R_{\rm Sedov}$ is the shock radius at $t=t_{\rm Sedov}$. Then we obtain

\begin{displaymath}% p_{\rm m}^{\rm (esc)}(R_{\rm sh}) = p_{\rm knee} \left(\frac{R_{\rm sh}}{R_{\rm Sedov}}\right)^{-6.5/ab}, \end{displaymath} (34)

so that we have $\alpha = 6.5/ab$ (see the last paragraph of Sect. 3.3). If we adopt $a\approx2.5$ ad hoc, in addition to $b\approx2/5$ which is expected in the Sedov phase, we find $\alpha\approx6.5$ as a phenomenologically required value in the escape-limited model. Note that even if the SNR dynamics are modified by the CR escape, b is almost the same as one in the adiabatic case and its effect on $\alpha$ is expected to be so small that our conclusion is not affected qualitatively.

4.3 Dynamics of SNR shock waves

In this subsection, we consider the dynamics of the SNR shock to estimate the evolution of the normalization factor of the spectrum  $K(R_{\rm sh})$. A simple treatment of the dynamics of the SNR shock from the free expansion to the adiabatic expansion (Sedov) phase has been given by several authors (Drury et al. 1989; Bisnovatyi-Kogan & Silich 1995; Ostriker & McKee 1988). Here we extend their method taking into account the cooling by CR escape. The total mass of the SNR shock shell is calculated as

\begin{displaymath}% M(R_{\rm sh}) = M_{\rm ej} + 4\pi \int_0^{R_{\rm sh}} {\rm d} r r^2 \rho_{\rm am}(r), \end{displaymath} (35)

where $M_{\rm ej}$ and $\rho_{\rm am}(r)$ are the ejecta mass and the density of ambient gas, respectively. The equation of motion of the thin shell is given by (Bisnovatyi-Kogan & Silich 1995; Ostriker & McKee 1988)

\begin{displaymath}% \frac{{\rm d}(Mu)}{{\rm d}t} = 4\pi R_{\rm sh}^2 (P_{\rm in}-P_{\rm am}), \end{displaymath} (36)

where the gas velocity u is related to the shock velocity  $u_{\rm sh}$ and the adiabatic index  $\gamma_{\rm ad}$ as $u=2u_{\rm sh}/(\gamma_{\rm ad}+1)$. Quantities  $P_{\rm in}$ and  $P_{\rm am}$ are the pressures of the post-shock gas and the ambient gas, respectively. For strong shocks, one can neglect  $P_{\rm am}$. The explosion energy ${\cal E}={\cal E}_{\rm th}+Mu^2/2 +Q_{\rm esc}$ consists of the internal energy ${\cal E}_{\rm th}=4\pi R_{\rm sh}^3P_{\rm in}/(3(\gamma_{\rm ad}-1))$, the kinetic energy, and the energy  $Q_{\rm esc}$, which is carried by the escaping CRs. Then Eq. (36) is rewritten as

\begin{displaymath}% \frac{{\rm d}(Mu)}{{\rm d}t} + \frac{3(\gamma_{\rm ad}-1)}{... ...\rm sh}} \left\{{\cal E} - Q_{\rm esc}(R_{\rm sh}) \right\}. \end{displaymath} (37)

Since ${\rm d} t=2{\rm d} R_{\rm sh}/(\gamma_{\rm ad}+1)u$, the left-hand side of Eq. (37) becomes

\begin{displaymath}% \frac{\gamma_{\rm ad}+1}{4M}R_{\rm sh}^{-\omega} \frac{{\rm d}}{{\rm d} R_{\rm sh}}\left[ R_{\rm sh}^\omega(Mu)^2\right], \end{displaymath} (38)

where $\omega=6(\gamma_{\rm ad}-1)/(\gamma_{\rm ad}+1)$. Hence we obtain
                       $\displaystyle % u_{\rm sh}(R_{\rm sh})$ = $\displaystyle \frac{\gamma_{\rm ad}+1}{2} \left[ \frac{2\omega}{M^2(R_{\rm sh})R_{\rm sh}^{\omega}} \right.$  
    $\displaystyle \times ~ \left. \int_0^{R_{\rm sh}} {\rm d} r ({\cal E} -Q_{\rm esc}(r)) r^{\omega -1} M(r) \right]^{1/2}.$ (39)

Once $Q_{\rm esc}(r)$ and $\rho_{\rm am}(r)$ are given, we can integrate Eqs. (35) and (39) to derive the dynamics of the SNR,  $u_{\rm sh}(R_{\rm sh})$. We mainly consider the evolution of  $u_{\rm sh}$ in the Sedov phase. Hence, below we simply assume that $\rho_{\rm am}$ is constant with r because even for the core-collapse supernova, the wind region is about a few pc and has been already passed by the SNR shock when the Sedov phase begins. The energy carried by escaping particles $Q_{\rm esc}$ is written as

\begin{displaymath}% Q_{\rm esc}(R_{\rm sh}) = \Theta (R_{\rm sh}-R_{\rm Sedov})... ...rm m}(R_{\rm sh})}^{p_{\rm knee}} cp N_{\rm esc}(p) {\rm d} p, \end{displaymath} (40)

where $\Theta(x)$ is the Heaviside step function, i.e.,  $\Theta(x)=1$ for x>0 while $\Theta(x)=0$ for x<0. Substituting Eqs. (27), (28), and (34) into Eq. (40), we obtain
                       $\displaystyle % Q_{\rm esc}(R_{\rm sh})$ = $\displaystyle \eta {\cal E}\Theta (R_{\rm sh}-R_{\rm Sedov})$  
    $\displaystyle \times \left\{ \begin{array}{ll} \frac{1-(R_{\rm sh}/R_{\rm Sedov... ...log(p_{\rm knee}/m_{{\rm p}}c)} & (s+\beta/\alpha = 2), \\ \end{array} \right.$ (41)

where $\eta$ is the ratio of the total energy of escaping CRs to the explosion energy. Note that  $\gamma_{\rm ad}$ is assumed to be constant. Integrating Eq. (39) with the aid of Eqs. (35) and (41), the dynamics of the SNR shock are determined. To make a more realistic case, we should perform fluid simulations with CR back reaction, which is beyond the scope of this paper (e.g., Drury et al. 1989; Völk et al. 2008).

4.4 Evolution of K

In this subsection we discuss the evolution of the normalization factor  $K(R_{\rm sh})$ of the spectrum of accelerated particles. At present, the injection process for CR acceleration at the shock is not well understood. Hence we consider two representative scenarios of the injection process to model the amount of the accelerated particles.

At first, we consider the same injection model as that of Ptuskin & Zirakashvili (2005). The model requires that the CR pressure at the shock is proportional to the fluid ram pressure, that is, $P_{\rm CR}\propto \rho_{\rm am}u_{\rm sh}^2$. The CR pressure at the shock  $P_{\rm CR}$ is given by

\begin{displaymath}% P_{\rm CR} \approx \frac{4\pi}{3} \int_{m_{{\rm p}}c}^{p_{\rm m}(R_{\rm sh})} {\rm d} pvp^3f_0(p), \end{displaymath} (42)

where we neglect the contribution of non-relativistic particles. Then one can find that  $P_{\rm CR}$ is related to $f_0(m_{{\rm p}}c)$ as
$\displaystyle % P_{\rm CR}\propto \left\{ \begin{array}{ll} f_0(m_{{\rm p}}c) &... ..._{{\rm p}}c)R_{\rm sh}^{\alpha(s-2)} & (s<2, {\rm PH}), \\ \end{array} \right.$   (43)

where we used the fact that the distribution function of CRs at the shock front is essentially $f_0(p)\propto p^{-s-2}$ for $p>m_{{\rm p}}c$. Note that $P_{\rm CR}$ is approximately proportional to $f_0(m_{{\rm p}}c)$ when s>2, because in this case  $P_{\rm CR}$ only weakly depends on $p_{\rm m}$. Hereafter, the cases of s>2 and s<2 are called PS and PH, respectively. Since $P_{\rm CR}\propto \rho_{\rm am}u_{\rm sh}^2$, the normalization factor of the CR spectrum inside the SNR shock $K(R_{\rm sh})$ is calculated as
$\displaystyle % K(R_{\rm sh}) \propto R_{\rm sh}^3 f_0(m_{{\rm p}}c) \propto \l... ...\ E(R_{\rm sh}) R_{\rm sh}^{\alpha(2-s)} & ({\rm PH}), \\ \end{array} \right.$   (44)

where the mechanical energy of the ejecta is (up to the numerical coefficient)

\begin{displaymath}% E(R_{\rm sh}) \approx \frac{4}{3} \pi \rho_{\rm am}u_{\rm sh}^2R_{\rm sh}^3 \propto u_{\rm sh}^2R_{\rm sh}^3. \end{displaymath} (45)

If the explosion is adiabatic, $E(R_{\rm sh})$ is constant with $R_{\rm sh}$ during the Sedov phase because $u_{\rm sh}\propto R_{\rm sh}^{-3/2}$. However, if the modification of SNR dynamics via CR cooling is taken into account (Sect. 4.3), then  $u_{\rm sh}$ is no longer proportional to $R_{\rm sh}^{-3/2}$, and it is obtained by solving Eq. (39). In Figs. 1 and 2 we show $\delta={\rm d} \ln E/{\rm d}\ln R_{\rm sh}$ as a function of time. Figure 1 is for the case PS in which we adopt $\eta=0.2$, s=2.23, $\gamma_{\rm ad}=5/3$, and $\beta=-0.2$, while Fig. 2 is for the case PH with parameters $\eta=0.2$, $\gamma_{\rm ad}=4/3$, and $\beta=\alpha(2-s)-0.2$. One can find that after $t\sim t_{\rm Sedov}$, $\delta$ is around -0.2, so that  $E(R_{\rm sh})$ is approximately given by $E(R_{\rm sh})\propto R_{\rm sh}^{\delta}$ with constant $\delta\approx-0.2$[*]. For PS, the effect of CR cooling appears at the late time because low-energy CRs have a large fraction of the total energy of the escaping CRs. On the other hand, for PH, the effect emerges at an earlier epoch.

\begin{figure} \par\includegraphics[width=8.5cm,clip]{13495fg1.eps} \end{figure} Figure 1:

The index $\delta={\rm d} \ln E/{\rm d}\ln R_{\rm sh}$as a function of time for the case PS. The solid and dashed lines are for the adiabatic expansion and the expansion with CR cooling, respectively. See the text for details.

Open with DEXTER

\begin{figure} \par\includegraphics[width=8.5cm,clip]{13495fg2.eps} \end{figure} Figure 2:

The same as Fig. 1, but for the case PH.

Open with DEXTER

Next we consider the thermal leakage (TL) model (Malkov & Völk 1995). This model requires the continuity of the distribution function f0(p) to the downstream Maxwelian at the injection momentum  $p_{{\rm inj}}$, namely

\begin{displaymath}% f_0(p_{{\rm inj}}) \propto \rho_{\rm am} /p_{{\rm inj}}^3. \end{displaymath} (46)

With this we derive
$\displaystyle % K(R_{\rm sh}) \propto R_{\rm sh}^3 f_0(m_{{\rm p}}c) \propto \rho_{\rm am} R_{\rm sh}^3 p_{\rm inj}^{s_{\rm low}-1},$   (47)

where $s_{\rm low}$ is the spectral index in the non-relativistic regime, that is $f_0(p)\propto p^{-s_{\rm low}-2}$ for $p_{\rm inj}<p<m_{{\rm p}}c$. We further assume the constant  $\rho_{\rm am}$ and $p_{\rm inj} \propto u_{\rm sh}$. Then we obtain $K(R_{\rm sh}) \propto R_{\rm sh}^{3+(s_{\rm low}-1)(\delta-3)/2}$, so that $\beta=3+(s_{\rm low}-1)(\delta-3)/2$. Although the value of $s_{\rm low}$ is uncertain, one can expect $s_{\rm low}\geq2$. In the test particle approximation, one derives $s_{\rm low}=s\geq2$. If the nonlinear model of shock acceleration is considered, the spectrum f0(p) is softer than p-4 in the non-relativistic regime (Berezhko & Ellison 1999). In particular, if  $s_{\rm low}\approx3$, then $\beta\approx\delta$, so that the value of $\delta$ becomes similar to those of PS case.

In summary, the evolution of the normalization factor of the accelerated particles spectrum is given by $K(R_{\rm sh})\propto R_{\rm sh}^\beta$, with

$\displaystyle % \beta=\left\{ \begin{array}{ll} \delta & ({\rm PS}) \\ \delta+... ...PH}) \\ 3 + (s_{\rm low}-1)(\delta -3)/2 & ({\rm TL}). \\ \end{array} \right.$   (48)

4.5 The spectrum of escaping particles

We obtain from Eqs. (27), (28), and (48) the index of the momentum spectrum of escaping particles as

$\displaystyle % s_{\rm esc}=\left\{ \begin{array}{ll} s+\delta/\alpha & ({\rm P... ...+\{3+(s_{\rm low}-1)(\delta-3)/2\}/\alpha & ({\rm TL}). \\ \end{array} \right.$   (49)

Now, we discuss which injection model is suitable to reproduce the galactic CR spectrum observed on Earth. We adopt $\alpha=6.5$ as a typical value (see Sect. 4.2) and $s_{\rm esc}=2.2$ to depict the Galactic CR spectrum. We assume $s_{\rm esc} \approx 2.2{-}2.4$, where we believe that the deviation from 2.0 is significant. This can be tested not only by Galactic CR observations but also by observations of extragalactic galaxies. Note that even if we adopt $s_{\rm esc}=2.4$, our results are qualitatively unchanged.

For PS $s_{\rm esc}$ is smaller than s because $\delta\la0$, so that the model predicts the harder spectrum of escaping particles rather than that of the source. However, since $\delta/\alpha\approx0.03$, the difference is small. In order to reproduce the observed Galactic CR spectrum  $s_{\rm esc}$, the source spectrum should be $s\approx s_{\rm esc}\approx2.2$. Hence the PS model requires s>2 at the source. This condition is satisfied if we consider the diffusive shock acceleration at a moderate Mach number (Fujita et al. 2009). It is also possible to derive s>2 if we consider the effects of neutral particles (Ohira et al. 2009b; Ohira & Takahara 2009).

For PH, because the value of $\delta/\alpha$ is small, $s_{\rm esc}$ is always near 2, which is the value predicted by the diffusive shock acceleration theory in the strong shock, test-particle limit. In particular, if $\delta=0$ - indeed, even in the test-particle limit where the cooling via CR escape can be neglected, the value of $\delta$ is not exactly zero unless $t\gg t_{\rm Sedov}$ (see Fig. 2) - then, one can find $s_{\rm esc}=2$. This is what Berezhko & Krymsky (1998) and Ptuskin & Zirakashvili (2005) showed. Note that from Fig. 2 $\delta$ is negative for a long time, so that $s_{\rm esc}<2$. Therefore, model PH cannot reproduce the observed Galactic CR spectrum on Earth.

For TL, neglecting $\delta$, we obtain

\begin{displaymath}% s_{\rm esc}\approx s+\frac{3}{13}(3-s_{\rm low}). \end{displaymath} (50)

Hence if $s_{\rm low}<3$, then $s_{\rm esc}>s$. For the test particle acceleration, we have $s=s_{\rm low}$, so that the source spectrum should be $s\approx2.2$ in order to obtain $s_{\rm esc}=2.4$. For the nonlinear acceleration, one typically expects $s_{\rm low}>2$ and s<2, so that $s_{\rm esc}\leq 29/13 \approx 2.23$. Therefore, the nonlinear model seems to explain the inferred source spectrum marginally. Possible rigorous determination of the value of  $s_{\rm esc}$ may give a constraint on the theory of nonlinear acceleration.

5 Application to AGN cocoon shocks

In this section, taking account of a constraint derived from the spectrum on Earth, we study the origin of CRs with energies higher than ${\sim}{10}^{17.5}$ eV and their acceleration mechanism at AGNs. There are many works which discuss UHECR production in AGNs (e.g., Berezinsky et al. 2006; Rachen & Biermann 1993; Takahara 1990; Pe'er et al. 2009; Biermann & Strittmatter 1987; Berezhko 2008). Many of them focus on UHECR acceleration in radio galaxies including Fanaroff-Riley (FR) I and II galaxies, which typically have powerful jets. In the context of DSA, one can basically suppose three acceleration zones; internal shocks in jets, hot spots, and cocoon shocks. The former two are the most widely discussed scenarios, but the detailed study of DSA at such mildly relativistic shocks has not yet been achieved. We concentrate on the cocoon shock scenario proposed by Berezhko (2008), where the non-relativistic DSA theory can be applied.

In this scenario, extragalactic CRs with energies higher than the second knee ( ${\sim}{10}^{17.5}$ eV) may be accelerated at the outer cocoon shock running into the intergalactic medium (IGM). As powerful jets penetrating into a uniform ambient medium with a density  $\rho_{\rm am}$, the heads of the cocoon advances into the IGM with a velocity $u_{\rm h}$. At the same time, the cocoon expands sideways with a velocity $u_{\rm c}$. Since the typical cocoon shock is non-relativistic, we apply the escape-limited model considered in previous sections. Although we hereafter focus on this scenario, note that it is very uncertain whether the efficient acceleration occurs there since the observed non-thermal emission is much weaker than that from hot spots and lobes.

Below we investigate whether the CR spectrum above the second knee can be explained by the AGN cocoon shock scenario with the same parameters for young SNRs explaining the CR spectrum below the knee, which were discussed in Sect. 4. Similar to the previous calculations in Sect. 4, we hereafter calculate the values of $\tilde{\alpha}={\rm d}\ln p_m^{(\rm esc)}/{\rm d}\ln t$ and $\tilde{\beta}={\rm d}\ln K/{\rm d}\ln t$ to derive the spectral index  $s_{\rm esc}$. Here we adopt $\chi=t$, where $\chi$ is the variable which appeared in Sect. 3.3.

First, let us consider the evolution of the maximum momentum  $p_m^{\rm (esc)}(t)$ in a phenomenological way. For young SNR (Sect. 4.2), we phenomenologically expect $p_m \propto t^{-\frac{13}{5}}$. Then, by using the Sedov-Taylor solution ( $u_{\rm sh} \propto t^{-3/5}$ and $R_{\rm sh}\propto t^{2/5}$), we can easily obtain

\begin{displaymath}% p_m^{(\rm esc)} (t) \propto u_{\rm sh}^{\frac{13+2c}{3}} R_{\rm sh}^c, \end{displaymath} (51)

where c is a phenomenologically introduced parameter since one can expect $\ell \propto R_{\rm sh}^c$ in general. Note that $p_m^{\rm (esc)} \propto B \ell u_{\rm sh}$, however, B does not depend on $R_{\rm sh}$ but on $u_{\rm sh}$ as long as B is generated by plasma instabilities like CR streaming instabilities. In the cocoon shock scenario, $u_{\rm sh}$ and $R_{\rm sh}$ are replaced with $u_{\rm c}$ and $R_{\rm c}$, respectively.

In order to obtain the value of $\tilde{\alpha}={\rm d}\ln p_m^{(\rm esc)}/{\rm d}\ln t$, the dynamics of the AGN cocoon are necessary. A simple consideration of the cocoon dynamics for the constant density IGM tells us that $u_{\rm h}$ is almost time-independent and $u_{\rm c}$ evolves as $u_{\rm c} \propto t^{-1/2}$ (Begelman & Cioffi 1989), so that the cocoon radius evolves as $R_{\rm c} \propto u_{\rm c}t\propto t^{1/2}$ and the jet radius evolves as $R_{\rm h} \propto u_{\rm h} t\propto t$. Then, we obtain $\tilde{\alpha}=\frac{13}{6}$ and $\tilde{\alpha}=2$ for c=0 and c=1, respectively.

Next let us consider the time dependence of the normalization factor of the spectrum of accelerated CRs, $K(t) \propto t^{\tilde{\beta}}$. The volume of the acceleration region swept by the cocoon shock is ${\sim}\mathcal{A} u_{\rm c} t$, where  $\mathcal{A}$ is the total area of the shock surface. If we assume the elliptical shape of the cocoon[*], then $\mathcal{A} \propto R_{\rm c} u_{\rm h}t\propto t^{3/2}$, so that $\mathcal{A} u_{\rm c} t\propto t^{2}$. For PS and PH,  $f_0(m_{{\rm p}}c)$ is related to  $P_{\rm CR}$. The dependence of $P_{\rm CR}$ is written as $P_{\rm CR} \propto \rho_{\rm am} u_{\rm c}^2\propto t^{-1}$, where we neglect for simplicity the evolution of the acceleration efficiency which was considered in Berezhko (2008). Accordingly, for PS $f_0(m_{\rm p}c) \propto P_{\rm CR}$ leads to

$\displaystyle % K(t) \propto (\mathcal{A} u_{\rm c} t) \times f_0(m_{\rm p}c) \propto t^{1},$   (52)

while for PH, $f_0(m_{\rm p}c) \propto P_{\rm CR} t^{\tilde{\alpha}(2-s)}$ leads to
$\displaystyle % K(t) \propto (\mathcal{A} u_{\rm c} t) \times f_0(m_{\rm p}c) \propto t^{1+\tilde{\alpha}(2-s)}.$   (53)

Finally, for the model TL (see Sect. 4.4), $f_0(m_{\rm p}c) \propto \rho_{\rm am} p_{\rm inj}^{s_{\rm low}-1}$ leads to
$\displaystyle % K(t) \propto (\mathcal{A} u_{\rm c} t) \times f_0(m_{\rm p}c) \propto t^{(5-s_{\rm low})/2},$   (54)

where we assume the constant IGM density[*]. Therefore we obtain
$\displaystyle % \tilde{\beta}=\left\{ \begin{array}{ll} 1 & ({\rm PS}) \\ 1+\t... ...}(2-s) & ({\rm PH}) \\ (5-s_{\rm low})/2 & ({\rm TL}). \\ \end{array} \right.$   (55)

With the above results we can obtain the spectrum of escaping particles. Here, we assume $\tilde{\alpha}=2$ expected for c=1 ( $\ell \propto R_{\rm c}$). Then we derive from Eq. (28)
$\displaystyle % s_{\rm esc}= \left\{ \begin{array}{ll} s+1/2 & ({\rm PS}) \\ 5/2 & ({\rm PH}) \\ s+(5-s_{\rm low})/4 & ({\rm TL}). \\ \end{array} \right.$   (56)

Interestingly, all three cases (PS, PH and TL) lead to the source spectral index $s_{\rm esc} \approx 2.4{-}2.7$, which is required in the proton dip model and in other extragalactic scenarios. For the model PS, s=2.2 leads to $s_{\rm esc} = s + 1/2 = 2.7$. For the model PH, s < 2 leads to $s_{\rm esc} = 2 + 1/2 =2.5$. Note that in this case the value of $s_{\rm esc}$ does not depend on s but on $\tilde{\alpha}$ generally. Finally for the model TL, $s=s_{\rm low}=2$ (based on the test particle acceleration) leads to $s_{\rm esc} = 2+3/4=2.75$, which is also consistent with the observed UHECR spectrum in the proton dip model unless the redshift evolution of UHECR sources are fast. If the effect of nonlinear acceleration is prominent, s < 2 and $s_{\rm low}>2$ can give harder indices of  $s_{\rm esc} < 2.75$.

6 Summary and discussion

We investigated the escape-limited model of CR acceleration, in which the maximum energy of CRs of an accelerator is limited by the escape from the acceleration site. The typical energy of escaping CRs decreases as the shock decelerates because the diffusion length becomes longer. After revisiting the escape-limited model and reconsidering its details more generally, we derived a simple relation between the spectrum of escaping particles and that in the accelerator. Then, using the obtained relation, we discussed which model of injection is potentially suitable to make the Galactic and extragalactic CRs observed at the Earth. As discussed in the beginning of Sect. 3, there are two approaches to the maximum momentum of accelerated particles in the escape-limited model, momentum or spatial boundary. The reality is somewhere between these two extremes in any case. However, we note that once a delta-function approximation is made as in Eq. (21), the two are essentially identical.

For young SNRs we considered the shock evolution with cooling by escaping CRs and those spectra for the three injection models. As a result, we found that for PH, it is difficult to satisfy the condition for the source spectrum of Galactic CRs ( $s_{\rm esc} \approx 2.2{-}2.4$). On the other hand, $s_{\rm esc} \approx 2.2{-}2.4$ can be achieved for PS and TL.

We also applied our escape-limited model to AGN cocoon shocks as well as young SNRs. This model is just one of the various candidates proposed so far, even if AGNs are UHECR accelerators. Nevertheless, it is interesting that the young SNR and the AGN cocoon shock scenarios can explain the Galactic and extragalactic cosmic rays observed on Earth in the same picture for all the three injection models if we accept the proton-dip model inferring $s_{\rm esc} \approx 2.4{-}2.7$. Whether the proton dip model is real or not can also be tested by future UHECR and high-energy neutrino observations. We focused on the proton case. Obviously, heavier nuclei become important above the knee so that we need to take them into account to explain the CR spectrum over the whole energy range. We can also apply the escape-limited model to heavy nuclei CRs for this purpose, although this is beyond the scope of this paper.

We point out a potential problem for the magnetic field amplification in the escape-limited model. For young SNRs, we determined the evolution of the maximum energy in the phenomenological way, and adopted $\alpha\approx6.5$. With Eqs. (16) and (51), we obtained $B \propto u_{\rm sh}^{\frac{10+2c}{3}}$, where $\ell \propto R_{\rm sh}^c$. The same result is obtained for AGN cocoon shocks because we considered the case in which the same parameters describe both the young SNR shocks and the AGN cocoon shocks. In particular, for c =1 as used in Eq. (31), we obtained $B^2 \propto u_{\rm sh}^8$, which means that B rapidly decreases with radius (or time). In principle, both the dependence of B on $u_{\rm sh}$ and the value of c can be determined theoretically, and then the evolution of the maximum energy should be predicted. Some previous works are based on theoretical arguments on the magnetic field evolution (e.g., Bell 2004; Caprioli et al. 2009; Berezhko 2008; Pelltier et al. 2006), which seem to be different from our phenomenological one. At present, the mechanisms of the particle acceleration and the magnetic field amplification are still highly uncertain despite many theoretical efforts (e.g., Riquelme & Spitkovsky 2009; Niemiec et al. 2008; Ohira et al. 2009a; Luo & Melrose 2009). Hence we expect that further theoretical and observational studies can resolve this discrepancy or exclude the possibility of escape-limited acceleration in the future.

We mainly considered spectra of dispersed CRs around young SNRs and AGN cocoon shocks. However, applications to other astrophysical objects are, of course, possible. For example, the old SNRs detected by Fermi LAT, like W28, W44, W51 and IC 443 (Abdo et al. 2009b,a) have been of great interest because they likely generate escaping CRs. In fact, the number of CRs around these old SNRs is likely to decrease with time or the shock radius, that is $\beta<0$ while $\alpha>0$. For example, when we consider the dynamics of an old SNR, we have $u_{\rm sh} \propto t^{-2/3} \propto R_{\rm sh}^{-2}$ (e.g., Yamazaki et al. 2006), so that we have $E \propto u_{\rm s}^2R_{\rm s}^3\propto R_{\rm sh}^{-1}$, i.e., $\delta = -1 < 0$. On the other hand, the value of $\alpha$ may be different from 6.5, which could be attributed to various complications like the interaction with the dense molecular cloud, and so on. For example, for the maximum hardening case, that is, $s_{\rm esc}=s + (s-1)(\sigma^{-1}-1)/\alpha$ (see Appendix A.2), we find $s_{\rm esc} \approx 1.5$ when $\alpha \approx 1$ and $s \approx 2.5$ where we assume $s=(\sigma +2)/(\sigma -1)$. This might be the case for old SNRs such as W51C (Abdo et al. 2009b). In addition, the maximum energy may be rather small for the old SNRs, so that the spectrum above  $p_m^{(\rm esc)}$ would be suppressed. The spectrum of high-energy gamma rays could give us important information on both the acceleration and escape processes of CRs with energies much lower than the knee energy.

Acknowledgements
We thank Akira Okumura and Yutaka Fujita for useful comments. We also thank the referee, Luke Drury, for valuable comments to improve the paper. Y.O. and K.M. acknowledge Grant-in-Aid from JSPS. This work was supported in part by grant-in-aid from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, No. 19047004, No. 21740184, No. 21540259 (R. Y.).

Appendix A: Effect of adiabatic loss on the spectrum of escaping CRs

The instantaneous spectrum of CRs escaping from the acceleration site when the shock radius  $R_{\rm sh}$ is obtained from Eq. (23) of Ptuskin & Zirakashvili (2005),

                               $\displaystyle % \frac{{\rm d}f_{\rm esc}}{{\rm d}R_{\rm sh}}$ = $\displaystyle \frac{1}{u_{\rm sh}}\frac{{\rm d}f_{\rm esc}}{{\rm d}t} = 4\pi\delta(p-p_{\rm m}(R_{\rm sh}))$  
    $\displaystyle \times~ \left\{ \frac{\left(1-\frac{1}{\sigma}-\frac{\xi_{\rm cr}}{2} \right)}{3}R_{\rm sh}^2p_{\rm m}f(p_{\rm m},R_{\rm sh}) \right.$  
    $\displaystyle \left. - ~\int_{0}^{R_{\rm sh}}{\rm d}rr^2f(p_{\rm m},r)\left(\fr... ...}}+\left(1-\frac{1}{\sigma}\right)\frac{p_{\rm m}}{R_{\rm sh}}\right) \right\},$ (A.1)

where $\sigma$ and $\xi_{\rm cr}$ are the total shock compression ratio and the ratio of the CR pressure to the shock ram pressure. The spectrum and the maximum momentum at the shock are assumed as

\begin{displaymath}% f(p,R_{\rm sh}) = Ap^{-s-2}R_{\rm sh}^{\beta-3}, \end{displaymath} (A.2)

\begin{displaymath}% p_{\rm m}(R_{\rm sh}) = R_{\rm sh}^{-\alpha}, \end{displaymath} (A.3)

where A is a normalization factor of the distribution function, and p and $p_{\rm m}$ are normalized by $p_{\rm knee}$ while $R_{\rm sh}$ is normalized by $R_{\rm Sedov}$. The distribution function of CRs at $r<R_{\rm sh}$ can be found by the same method as Ptuskin & Zirakashvili (2005). The fluid velocity at $r<R_{\rm sh}$ is

\begin{displaymath}% u(r,t) = \left(1-\frac{1}{\sigma}\right)\frac{r}{R_{\rm sh}(t)}u_{\rm sh}(t), \end{displaymath} (A.4)

and we solve the following equation

\begin{displaymath}% \frac{\partial f}{\partial t}+u(r)\frac{\partial f}{\partia... ...ft(r^2u\left(r\right)\right)p \frac{\partial f}{\partial p}=0. \end{displaymath} (A.5)

Then one can get the following solution

\begin{displaymath}% f(p,r)=Ap^{-s-2}r^{\sigma (s+\beta-1)-s-2}R_{\rm sh}^{-(s+\beta-1)(\sigma-1)}. \end{displaymath} (A.6)

Therefore the total spectrum of escaping CRs is
$\displaystyle % f_{\rm esc}=\int \frac{{\rm d}f_{\rm esc}}{{\rm d}R_{\rm sh}}{\rm d}R_{\rm sh} =f_{\rm esc,surface}+f_{\rm esc,inside},$   (A.7)

where $f_{\rm esc,surface}$ is the first term of Eq. (A.1) and $f_{\rm esc,inside}$ is the second term of Eq. (A.1). These are given by

\begin{displaymath}% f_{\rm esc,surface}=\frac{4\pi A}{\alpha}\frac{1-\frac{1}{\sigma}-\frac{\xi_{\rm cr}}{2}}{3}p^{-s-2-\beta/\alpha}, \end{displaymath} (A.8)


                             $\displaystyle % f_{\rm esc,inside}$ = $\displaystyle 4\pi A\left(\alpha-1+\frac{1}{\sigma}\right)$  
    $\displaystyle \times ~\int_{0}^{R_{\rm sh}}{\rm d}R \delta(p-p_{\rm m}(R))p_{\rm m}^{-1-s}R^{-(s+\beta-1)(\sigma-1)-1}$  
    $\displaystyle \times ~ \int_{0}^{R}{\rm d}rr^{\sigma(s+\beta-1)-s}.$ (A.9)

The ${\rm d}r$ integral of Eq. (A.9) is approximately given by

  $\displaystyle \int_{0}^{R}{\rm d}rr^{\sigma(s+\beta-1)-s} \propto$  
  $\displaystyle \qquad \qquad \left\{ \begin{array}{ll} R^{\sigma(s+\beta-1)-s+1}... ...\left(s-1\right)\left(\frac{1}{\sigma}-1\right)\right], \\ \end{array} \right.$ (A.10)

where R0 is the radius at which $\beta_{\rm free}$ ( ${>}\left(s-1\right)\left(\sigma^{-1}-1\right)$) changes to $\beta$ ( ${<}\left(s-1\right)\left(\sigma^{-1}-1\right)$). Here, we assume that $\beta$ changed when $R_{\rm sh}=1$, that is, $t=t_{\rm sedov}$. Then from Eq. (A.4), $R_0 = R_{\rm sh}^{1-\sigma^{-1}}$.

A.1  ${\beta>\left(s-1\right)\left(\frac{1}{\sigma}-1\right)}$

In this case, the ${\rm d}r$ integral of Eq. (A.9) is dominated by the outer region, that is, the spectrum of escaping particles does not depend on the past acceleration history. $f_{\rm esc,inside}$ is calculated as

                         $\displaystyle % f_{\rm esc,inside}$ = $\displaystyle \frac{4\pi A(\alpha-1+\frac{1}{\sigma})}{1-s+\sigma (s+\beta-1)}$  
    $\displaystyle \times ~ \int_{0}^{R_{\rm sh}}{\rm d}R \delta(p-p_{\rm m}(R))p_{\rm m}^{-1-s}R^{\beta-1}$  
  = $\displaystyle \frac{4\pi A}{\alpha}\frac{\alpha-1+\frac{1}{\sigma}}{1-s+\sigma (s+\beta-1)}p^{-s-2-\beta/\alpha}.$ (A.11)

Hence, as long as $\beta>\left(s-1\right)\left(\sigma^{-1}-1\right)$, the energy spectrum index of escaping CRs is $s_{\rm esc}=s+\beta/\alpha$.

A.2 ${\beta<\left(s-1\right)\left(\frac{1}{\sigma}-1\right)}$

In this case, the ${\rm d}r$ integral of Eq. (A.9) is dominated by the inner region, that is, the spectrum of escaping particles depends on the past acceleration history. Especially, only the acceleration at $t=t_{\rm sedov}$ is important, and then $f_{\rm esc,inside}$ is calculated as

                            $\displaystyle % f_{\rm esc,inside}$ = $\displaystyle \frac{4\pi A(\alpha-1+\frac{1}{\sigma})}{1-s+\sigma (s+\beta_{\rm free}-1)}$  
    $\displaystyle \times ~ \int_{R_0}^{R_{\rm sh}}{\rm d}R \delta(p-p_{\rm m}(R))p_{\rm m}^{-1-s}R^{(s-1)(\frac{1}{\sigma}-1)-1}$  
  = $\displaystyle \frac{4\pi A}{\alpha}\frac{\alpha-1+\frac{1}{\sigma}}{1\!-\!s\!+\!\sigma (s\!+\!\beta_{\rm free}\!-\!1)}p^{-s-2-(s-1)(\frac{1}{\sigma}-1)/\alpha}.$ (A.12)

Because $f_{\rm esc,inside}$ is softer than $f_{\rm esc,surface}$ (at  $p < p_{\rm m}$), $f_{\rm esc,inside}$ is larger than $f_{\rm esc,surface}$. Hence the spectral index of escaping CRs is $s_{\rm esc}=s + (s-1)(\sigma^{-1}-1)/\alpha$ and does not depend on $\beta$. The largest possible hardening from the spectral index of the acceleration site $\Delta s = s-s_{\rm esc}$ is
$\displaystyle % \Delta s =\frac{s-1}{\alpha}\left(1-\frac{1}{\sigma}\right) =\frac{3(s-1)}{\alpha(s+2)},$   (A.13)

where we assume the relation between the spectrum index and the compression ratio, $s=(\sigma +2)/(\sigma -1)$.



References

Footnotes

...$\delta\approx-0.2$[*]
We search for approximate solutions to Eqs. (39), (41) and (48), which determine $\delta$ (and $\beta$) for given parameters $\alpha$, s, $\eta$ and  $\gamma_{\rm ad}$. The procedure is as follows. For a trial value of $\beta$ which is assumed to be constant, we solve Eqs. (39) and (41) to obtain $u_{\rm sh}(R_{\rm sh})$ so as to calculate E and $\delta$ as functions of  $R_{\rm sh}$. It is found that $\delta$ is almost time-independent after  $t_{\rm Sedov}$, hence we can derive its average value in the epoch $t>t_{\rm Sedov}$. Then we update the value of $\beta$ with Eq. (48) with the average value of $\delta$. We repeat this procedure until the iteration converges.
... cocoon[*]
Berezhko (2008) assumed a kind of spherical cocoon, i.e.,  $\mathcal{A} \propto R_{\rm c} u_{\rm c} t$, and used $K(t) \propto \rho_{\rm am} R_{\rm c}^3$ which is similar to the relation obtained for the model PS. However, when the cocoon becomes spherical, $u_{\rm c}$ and $u_{\rm h}$ evolve according to the adiabatic solution as in the Sedov-Taylor solution for SNRs (Fujita et al. 2007; Begelman & Cioffi 1989).
... density[*]
However, it may not be the case since the IGM density may drop with the distance from the nucleus. Then we need to use different adiabatic solutions for $u_{\rm c}$ and $R_{\rm c}$.
Copyright ESO 2010

All Figures

  \begin{figure} \par\includegraphics[width=8.5cm,clip]{13495fg1.eps} \end{figure} Figure 1:

The index $\delta={\rm d} \ln E/{\rm d}\ln R_{\rm sh}$as a function of time for the case PS. The solid and dashed lines are for the adiabatic expansion and the expansion with CR cooling, respectively. See the text for details.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=8.5cm,clip]{13495fg2.eps} \end{figure} Figure 2:

The same as Fig. 1, but for the case PH.

Open with DEXTER
In the text

Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.