Escapelimited model of cosmicray acceleration revisited
Y. Ohira^{1}  K. Murase^{2}  R. Yamazaki^{3}
1  Department of Earth and Space Science, Osaka University, Toyonaka 5600043, Japan
2 
Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 6068502, Japan
3 
Department of Physical Science, Hiroshima University, HigashiHiroshima 7398526, 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 gammaray spectrum of young supernova remnants (SNRs) such as
RX J1713.73946, but also in the spectrum of CRs showering the
Earth.
Aims. The escapelimited model of cosmicray 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 longstanding problems. The number spectrum of nuclear CRs observed at the Earth, , shows a break at the ``knee'' energy ( eV), below which the spectral index is about (Cronin 1999). Because of the energydependent 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 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 gammaray 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 highenergy gammaray observations have revealed the existence of highenergy 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 backreactions 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 (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 inferred from the CR spectrum on Earth. Even in the testparticle limit of DSA, such a soft source spectrum requires a shock with the small Mach number (), 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 agelimited 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 escapelimited acceleration, to explain the observed CR spectrum on Earth (Sect. 3). This model is preferable to the agelimited acceleration when we consider observational results for young SNR RX J1713.73946, of which TeV 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 ( eV) may be of Galactic origin, the highest energy CRs above 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), gammaray 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 eV to 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 highenergy CRs, because they can explain the ultra high energy cosmicray (UHECR) spectrum above eV assuming the proton composition. In such a proton dip model, the source spectrum of UHECRs is with s=2.42.7, depending on the models of the source evolution (Berezinsky et al. 2006). The required source spectral index of s=2.42.7 can be explained by several possibilities. First, it can be attributed to the acceleration mechanism itself. One can consider nonFermi acceleration mechanisms (Berezinsky et al. 2006) or the twostep 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 agelimited model of the CR acceleration (Sect. 2), we study the escapelimited model in Sect. 3. For a simple understanding, the general argument in a stationary, testparticle 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 agelimited acceleration
For comparison with the escapelimited acceleration, we briefly
summarize the case of the agelimited acceleration. In this case,
the maximum momentum of accelerated particles
is determined by
,
where
and
are the age of the shock and the acceleration time scale, respectively.
When we consider DSA,
is given by (Drury 1983)
(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 Bohmtype diffusion, i.e.,
where B, and are the magnetic field strength, the gyrofactor and the proton mass, respectively. Taking into account that the fluid velocities are related to the shock velocity, , as , where is the compression ratio of the shock, we derive (e.g., Aharonian & Atoyan 1999)
3 Escapelimited 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 of the turbulence spectrum (Drury et al. 2009; Reynolds 1998). Particles with the Lorentz factor above satisfying the resonance condition, , where 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 highenergy 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 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 . As long as , 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 , 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 . Below we consider the escapelimited model, where the maximum energy is essentially determined by .
3.1 A simple case of stationary, testparticle approximation
In order to take an essential feature of the escapelimited
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
with the boundary condition , where () is the upstream escape boundary. The fluid velocity u(x) is given by
where u_{1} and u_{2} are constants. The solution to the transport equation in the testparticle approximation is derived as (Caprioli et al. 2009)
(6) 
where f_{0}(p)=f(x=0,p) is given by
and q=3u_{1}/(u_{1}u_{2}). The escape flux at is
=  
=  (8) 
Let us introduce a new variable and a new function . Then we expand around its maximum value at . To do this, we calculate the first and the second derivatives as
(9) 
=  
(10) 
Below we consider the case of Bohm diffusion, . Then one can find
(where ) so that we obtain
=  
=  (12) 
where . Note that the quantity given by Eq. (11) also plays the role of maximum momentum of the accelerated particles at the shock (x=0) because one can see for while for (Caprioli et al. 2009).
Finally, going back to the function ,
we obtain
One can clearly see from Eq. (13) that particles with momentum around 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, testparticle case the quantity
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 escapelimited case the maximum momentum
is determined by
Given that
(15) 
which is the same as in the agelimited case (Eq. (2) in Sect. 2), we obtain
Since , we find from Eqs. (3) and (16)
Hence, as long as , which is expected in the early phase of the shock, we do not need the escapelimited acceleration. However, if the shock evolves so that , we have to consider the effect of the particle escape and cannot simply apply the wellknown result of agelimited acceleration, Eq. (3).
3.3 The spectrum of CRs dispersed from an accelerator
In this subsection, we derive the timeintegrated 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
,
at a certain epoch labeled by a parameter ,
is defined as the number of protons with momentum between p and
,
which is produced in the interval between
and
.
Here
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
contains the term of exponential cutoff at the momentum
which depends on
(see, for example, Eq. (25)). The number of protons with momentum between p and
which is escaping from the accelerator at the epoch between
and
,
is denoted by
,
and we assume
where
and . This is because we expect that particles with momentum around are escaping most efficiently from the source. Indeed, as shown in Sect. 3.1, Eqs. (18) and (19) are a good approximation in the testparticle, stationary case. Here, we simply expect these assumptions are also correct in general.
The timeintegrated spectrum of protons which have escaped at the source
is obtained by
In order to derive a simple analytical form, we approximate Eq. (19) as
If we use a general mathematical formula for functions for an arbitrary function
(22) 
where , then Eq. (21) is rewritten as
where is the inverse function of , that is, and . Using Eqs. (18) and (23), we calculate as
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
is assumed to have
which is a power law with the index s and the exponential cut off at . Substituting Eqs. (25) into (24), we obtain
In particular, if and are written by the power law forms like and , then, , so that Eq. (26) becomes
where
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, timedependent 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 as , the final form of Eqs. (27) and (28) is correct as long as (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 is different from Eq. (28) (see Appendix).
4 Application to young supernova remnants
4.1 Inconsistency of agelimited acceleration with observed results of RX J1713.73946
RX J1713.73946 is a representative SNR from which bright TeV 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 ( mG) are derived from the width of synchrotron Xray filaments (Parizot et al. 2006; see also Bamba et al. 2005b,a,2003; Vink & Laming 2003) and from time variation of synchrotron Xray hot spots (Uchiyama et al. 2007). This also supports the hadronic origin of TeV rays, because the leptonic, onezone emission model (e.g., Aharonian & Atoyan 1999) cannot explain the TeVtoXray flux ratio. Hence it is natural to assume that the TeV rays are produced by the hadronic process, although there are several arguments against this interpretation (Katz & Waxman 2008; Plaga 2008; Butt 2008).
In the agelimited case, Eq. (3) reads
where we adopt , and is the magnetic field strength in units of mG.
HESS observation revealed that the cutoff energy of TeV ray spectrum is low (Aharonian et al. 2006,2007), so that in the onezone hadronic scenario the maximum energy of protons,
is estimated as 30100 TeV (Villante & Vissani 2007). If
TeV and mG, then Eq. (29) tells us
,
implying far from the ``Bohm limit'' (
)
which is inferred from the Xray 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 Xray
observations. The precise Xray spectrum of RX J1713.73946 is
revealed, which gives
cm s^{1} (Tanaka et al. 2008).
Then Eq. (29) can be rewritten as (Yamazaki et al. 2009)
(30) 
Hence, in order to obtain TeV, we need G in the context of the hadronic scenario of TeV rays, which contradicts current estimates of the magnetic field.
A possible solution is to consider the escapelimited acceleration. One can find from Eq. (17) that if we take
pc, the maximum energy becomes
=  (31) 
which is consistent with the observed gammaray spectrum. Below we consider the model of escapelimited 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 of .
4.2 Evolution of p_{m}
Time evolution of the maximum momentum of accelerated particle has been so far discussed in many contexts (e.g., Ptuskin & Zirakashvili 2003). One way to estimate 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 is also unknown. This prevents us from predicting 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
is expected to increase up to the knee energy (10^{15.5} eV) until the end of the free expansion phase
and decreases from that epoch. As seen in Sect. 4.1, is limited by the escape at
,
that is
.
Then, to reproduce the observed CR spectrum from GeV to the knee, we may assume a functional form of
(32) 
where , so that at . For later convenience, we change variables from t to the SNR radius in order to take , where is a variable which appeared in Sect. 3.3. We further assume the dynamics of as
(33) 
where is the shock radius at . Then we obtain
so that we have (see the last paragraph of Sect. 3.3). If we adopt ad hoc, in addition to which is expected in the Sedov phase, we find as a phenomenologically required value in the escapelimited 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 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
.
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; BisnovatyiKogan & 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
where and are the ejecta mass and the density of ambient gas, respectively. The equation of motion of the thin shell is given by (BisnovatyiKogan & Silich 1995; Ostriker & McKee 1988)
where the gas velocity u is related to the shock velocity and the adiabatic index as . Quantities and are the pressures of the postshock gas and the ambient gas, respectively. For strong shocks, one can neglect . The explosion energy consists of the internal energy , the kinetic energy, and the energy , which is carried by the escaping CRs. Then Eq. (36) is rewritten as
Since , the lefthand side of Eq. (37) becomes
(38) 
where . Hence we obtain
Once and are given, we can integrate Eqs. (35) and (39) to derive the dynamics of the SNR, . We mainly consider the evolution of in the Sedov phase. Hence, below we simply assume that is constant with r because even for the corecollapse 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 is written as
where is the Heaviside step function, i.e., for x>0 while for x<0. Substituting Eqs. (27), (28), and (34) into Eq. (40), we obtain
where is the ratio of the total energy of escaping CRs to the explosion energy. Note that 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 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,
.
The CR pressure at the shock
is given by
(42) 
where we neglect the contribution of nonrelativistic particles. Then one can find that is related to as
(43) 
where we used the fact that the distribution function of CRs at the shock front is essentially for . Note that is approximately proportional to when s>2, because in this case only weakly depends on . Hereafter, the cases of s>2 and s<2 are called PS and PH, respectively. Since , the normalization factor of the CR spectrum inside the SNR shock is calculated as
(44) 
where the mechanical energy of the ejecta is (up to the numerical coefficient)
(45) 
If the explosion is adiabatic, is constant with during the Sedov phase because . However, if the modification of SNR dynamics via CR cooling is taken into account (Sect. 4.3), then is no longer proportional to , and it is obtained by solving Eq. (39). In Figs. 1 and 2 we show as a function of time. Figure 1 is for the case PS in which we adopt , s=2.23, , and , while Fig. 2 is for the case PH with parameters , , and . One can find that after , is around 0.2, so that is approximately given by with constant ^{}. For PS, the effect of CR cooling appears at the late time because lowenergy 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.
Figure 1: The index 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 
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 f_{0}(p) to the downstream Maxwelian at the injection momentum
, namely
(46) 
With this we derive
(47) 
where is the spectral index in the nonrelativistic regime, that is for . We further assume the constant and . Then we obtain , so that . Although the value of is uncertain, one can expect . In the test particle approximation, one derives . If the nonlinear model of shock acceleration is considered, the spectrum f_{0}(p) is softer than p^{4} in the nonrelativistic regime (Berezhko & Ellison 1999). In particular, if , then , so that the value of becomes similar to those of PS case.
In summary, the evolution of the normalization factor of the accelerated particles spectrum is given by
, with
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
(49) 
Now, we discuss which injection model is suitable to reproduce the galactic CR spectrum observed on Earth. We adopt as a typical value (see Sect. 4.2) and to depict the Galactic CR spectrum. We assume , 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 , our results are qualitatively unchanged.
For PS is smaller than s because , so that the model predicts the harder spectrum of escaping particles rather than that of the source. However, since , the difference is small. In order to reproduce the observed Galactic CR spectrum , the source spectrum should be . 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 is small, is always near 2, which is the value predicted by the diffusive shock acceleration theory in the strong shock, testparticle limit. In particular, if  indeed, even in the testparticle limit where the cooling via CR escape can be neglected, the value of is not exactly zero unless (see Fig. 2)  then, one can find . This is what Berezhko & Krymsky (1998) and Ptuskin & Zirakashvili (2005) showed. Note that from Fig. 2 is negative for a long time, so that . Therefore, model PH cannot reproduce the observed Galactic CR spectrum on Earth.
For TL, neglecting ,
we obtain
(50) 
Hence if , then . For the test particle acceleration, we have , so that the source spectrum should be in order to obtain . For the nonlinear acceleration, one typically expects and s<2, so that . Therefore, the nonlinear model seems to explain the inferred source spectrum marginally. Possible rigorous determination of the value of 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 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 FanaroffRiley (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 nonrelativistic DSA theory can be applied.
In this scenario, extragalactic CRs with energies higher than the second knee ( 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 , the heads of the cocoon advances into the IGM with a velocity . At the same time, the cocoon expands sideways with a velocity . Since the typical cocoon shock is nonrelativistic, we apply the escapelimited 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 nonthermal 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 and to derive the spectral index . Here we adopt , where is the variable which appeared in Sect. 3.3.
First, let us consider the evolution of the maximum momentum
in a phenomenological way. For young SNR (Sect. 4.2), we phenomenologically expect
.
Then, by using the SedovTaylor solution (
and
), we can easily obtain
where c is a phenomenologically introduced parameter since one can expect in general. Note that , however, B does not depend on but on as long as B is generated by plasma instabilities like CR streaming instabilities. In the cocoon shock scenario, and are replaced with and , respectively.
In order to obtain the value of , the dynamics of the AGN cocoon are necessary. A simple consideration of the cocoon dynamics for the constant density IGM tells us that is almost timeindependent and evolves as (Begelman & Cioffi 1989), so that the cocoon radius evolves as and the jet radius evolves as . Then, we obtain and for c=0 and c=1, respectively.
Next let us consider the time dependence of the normalization factor of the spectrum of accelerated CRs,
.
The volume of the acceleration region swept by the cocoon shock is
,
where
is the total area of the shock surface. If we assume the elliptical shape of the cocoon^{}, then
,
so that
.
For PS and PH,
is related to
.
The dependence of
is written as
,
where we neglect for simplicity the evolution of the acceleration efficiency which was considered in Berezhko (2008). Accordingly, for PS
leads to
(52) 
while for PH, leads to
(53) 
Finally, for the model TL (see Sect. 4.4), leads to
(54) 
where we assume the constant IGM density^{}. Therefore we obtain
With the above results we can obtain the spectrum of escaping particles. Here, we assume expected for c=1 ( ). Then we derive from Eq. (28)
(56) 
Interestingly, all three cases (PS, PH and TL) lead to the source spectral index , which is required in the proton dip model and in other extragalactic scenarios. For the model PS, s=2.2 leads to . For the model PH, s < 2 leads to . Note that in this case the value of does not depend on s but on generally. Finally for the model TL, (based on the test particle acceleration) leads to , 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 can give harder indices of .
6 Summary and discussion
We investigated the escapelimited 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 escapelimited 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 escapelimited model, momentum or spatial boundary. The reality is somewhere between these two extremes in any case. However, we note that once a deltafunction 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 ( ). On the other hand, can be achieved for PS and TL.
We also applied our escapelimited 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 protondip model inferring . Whether the proton dip model is real or not can also be tested by future UHECR and highenergy 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 escapelimited 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 escapelimited model. For young SNRs, we determined the evolution of the maximum energy in the phenomenological way, and adopted . With Eqs. (16) and (51), we obtained , where . 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 , which means that B rapidly decreases with radius (or time). In principle, both the dependence of B on 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 escapelimited 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 while . For example, when we consider the dynamics of an old SNR, we have (e.g., Yamazaki et al. 2006), so that we have , i.e., . On the other hand, the value of 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, (see Appendix A.2), we find when and where we assume . 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 would be suppressed. The spectrum of highenergy gamma rays could give us important information on both the acceleration and escape processes of CRs with energies much lower than the knee energy.
AcknowledgementsWe 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 GrantinAid from JSPS. This work was supported in part by grantinaid 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
is obtained from Eq. (23) of Ptuskin & Zirakashvili (2005),
where and 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
(A.2) 
(A.3) 
where A is a normalization factor of the distribution function, and p and are normalized by while is normalized by . The distribution function of CRs at can be found by the same method as Ptuskin & Zirakashvili (2005). The fluid velocity at is
and we solve the following equation
(A.5) 
Then one can get the following solution
(A.6) 
Therefore the total spectrum of escaping CRs is
(A.7) 
where is the first term of Eq. (A.1) and is the second term of Eq. (A.1). These are given by
(A.8) 
The integral of Eq. (A.9) is approximately given by
where R_{0} is the radius at which ( ) changes to ( ). Here, we assume that changed when , that is, . Then from Eq. (A.4), .
A.1
In this case, the
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.
is calculated as
=  
=  (A.11) 
Hence, as long as , the energy spectrum index of escaping CRs is .
A.2
In this case, the
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
is important, and then
is calculated as
=  
=  (A.12) 
Because is softer than (at ), is larger than . Hence the spectral index of escaping CRs is and does not depend on . The largest possible hardening from the spectral index of the acceleration site is
(A.13) 
where we assume the relation between the spectrum index and the compression ratio, .
References
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009a, ApJS, 183, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009b, ApJ, 706, L1 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Aharonian, F. A., & Atoyan, A. M. 1999, A&A, 351, 330 [NASA ADS] [Google Scholar]
 Aharonian, F. A., Akhperjanian, A. G., Aye, K.M., et al. 2004, Nature, 432, 75 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2005, A&A, 437, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2006, A&A, 449, 223 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aharonian, F., Akhperjanian, A. G., BazerBachi, A. R., et al. 2007, A&A, 464, 235 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Aloisio, R., Berezinsky, V., Blasi, P., et al. 2007, Astropart. Phys., 27, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Arons, J. 2003, ApJ, 589, 871 [NASA ADS] [CrossRef] [Google Scholar]
 Axford, W. I., Leer, E., & Skadron, G. 1977, Proc. 15th Int. Cosmic Ray Conf., Plovdiv, 11, 132 [Google Scholar]
 Bamba, A., Yamazaki, R., Ueno, M., & Koyama, K. 2003, ApJ, 589, 827 [NASA ADS] [CrossRef] [Google Scholar]
 Bamba, A., Yamazaki, R., Yoshida, T., Terasawa, T., & Koyama, K. 2005a, ApJ, 621, 793 [NASA ADS] [CrossRef] [Google Scholar]
 Bamba, A., Yamazaki, R., & Hiraga, J. S. 2005b, ApJ, 632, 294 [NASA ADS] [CrossRef] [Google Scholar]
 Begelman, M. C., & Cioffi, D. F. 1989, ApJ, 345, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 1978, MNRAS, 182, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 2004, MNRAS, 353, 550 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhko, E. G. 2008, ApJ, 684, L69 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhko, E. G., & Ellison, D. C. 1999, ApJ, 526, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhko, E. G., & Krymsky, G. F. 1988, Soviet Phys.Uspekhi, 12, 155 [Google Scholar]
 Berezinsky, V., Gazizov, A., & Grigorieva, S. 2006, Phys. Rev. D, 74, 043005 [NASA ADS] [CrossRef] [Google Scholar]
 Biermann, P. L., & Strittmatter, P. A. 1987, ApJ, 322, 643 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29 [NASA ADS] [CrossRef] [Google Scholar]
 BisnovatyiKogan, G. S., & Silich, S. A. 1995, Rev. Mod. Phys., 67, 661 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Eichler, D. 1987, Phys. Rep., 154,1 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Butt, Y. 2008, MNRAS, 386, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Caprioli, D., Blasi, P., & Amato, E. 2009, MNRAS, 396, 2065 [NASA ADS] [CrossRef] [Google Scholar]
 Cronin, J. W. 1999, Rev. Mod. Phys., 71, S165 [CrossRef] [Google Scholar]
 Drury, L. O'C. 1983, Rep. Prog. Phys., 46, 973 [NASA ADS] [CrossRef] [Google Scholar]
 Drury, L. O'C., & Völk, H. J. 1981, ApJ, 248, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Drury, L. O'C., Markiewicz, W. J., & Voelk, H. J. 1989, A&A, 225, 179 [NASA ADS] [Google Scholar]
 Drury, L. O'C., Aharonian, F. A., Malyshev, D., & Gabici, S. 2009, A&A, 496, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Enomoto, R., Tanimori, T., Naito, T., et al. 2002, Nature, 416, 823 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Fujita, Y., Kohri, K., Yamazaki, R., & Kino, M. 2007, ApJ, 663, L61 [NASA ADS] [CrossRef] [Google Scholar]
 Fujita, Y., Ohira, Y., Tanaka, S. J., & Takahara, F. 2009, ApJ, 707, L179 [NASA ADS] [CrossRef] [Google Scholar]
 Gabici, S., Aharonian, F. A., & Casanova, S. 2009, MNRAS, 396, 1629 [NASA ADS] [CrossRef] [Google Scholar]
 Giacalone, J., & Jokipii, J. A. 2007, ApJ, 663, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Helder, E. A., Vink, J., Bassa, C. G., et al. 2009, Science, 325, 719 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Inoue, S., Sigl, G., Miniati, F., & Armengaud, E. 2007, Proceedings of the 30th ICRC, Merida, Mexico [arXiv:0711.1027] [Google Scholar]
 Inoue, T., Yamazaki, R., & Inutsuka, S. 2009, ApJ, 695, 825 [NASA ADS] [CrossRef] [Google Scholar]
 Kang, H., Rachen, J. P., & Biermann, P. L. 1997, MNRAS, 286, 257 [NASA ADS] [Google Scholar]
 Kang, H., Jones, T. W., LeVeque, R. J., & Shyue, K. M. 2001, ApJ, 550, 737 [NASA ADS] [CrossRef] [Google Scholar]
 Kachelriess, M., & Semikoz, D. V. 2006, Phys. Lett. B, 634, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Katagiri, H., Enomoto, R., Ksenofontov, L. T., et al. 2005, ApJ, 619, L163 [NASA ADS] [CrossRef] [Google Scholar]
 Katz, B., & Waxman, E. 2008, J. Cosmology Astropart. Phys., 01, 018 [NASA ADS] [CrossRef] [Google Scholar]
 Krymsky, G. F. 1977, Doki. Akad. Nauk SSSR, 234, 1306 [NASA ADS] [Google Scholar]
 Lucek, S. G., & Bell, A. R. 2000, MNRAS, 314, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Luo, Q., & Melrose, D. 2009, MNRAS, 397, 1402 [NASA ADS] [CrossRef] [Google Scholar]
 Malkov, M. A. 1997, ApJ, 485, 638 [NASA ADS] [CrossRef] [Google Scholar]
 Malkov, M. A., & Völk, H. J. 1995, A&A, 300, 605 [NASA ADS] [Google Scholar]
 Malkov, M. A., & Drury, L. O'C. 2001, Rep. Prog. Phys., 64, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Murase, K., Ioka, K., Nagataki, S., & Nakamura, T. 2006, ApJ, 651, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Murase, K., Inoue, S., & Nagataki, S. 2008, ApJ, 689, L105 [NASA ADS] [CrossRef] [Google Scholar]
 Murase, K., Mészáros, P., & Zhang, B. 2009, Phys. Rev. D, 79, 103001 [NASA ADS] [CrossRef] [Google Scholar]
 Niemiec, J., Pohl, M., Stroman, T., & Nishikawa, K. 2008, ApJ684, 1174 [NASA ADS] [CrossRef] [Google Scholar]
 Ohira, Y., & Takahara, F. 2009 [arXiv:0912.2859] [Google Scholar]
 Ohira, Y., Reville, B., Kirk, J. G., & Takahara, F. 2009a, ApJ, 698, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Ohira, Y., Terasawa, T., & Takahara, F. 2009b, ApJ, 703, L59 [NASA ADS] [CrossRef] [Google Scholar]
 Ostriker, J. P., & McKee, C. F. 1988, Rev. Mod. Phys., 60, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Pelltier, G., Lemoine, M., & Marcowith, A. 2006, A&A, 453, 181 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Parizot, E., Marcowith, A., Ballet, J., & Gallant, Y. A. 2006, A&A, 453, 387 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pe'er, A., Murase, K., & Mészáros, P. 2009, Phys. Rev. D., 80, 123018 [NASA ADS] [CrossRef] [Google Scholar]
 Plaga, R. 2008, New A. 13, 73 [Google Scholar]
 Ptuskin, V. S., & Zirakashvili, V. N. 2003, A&A, 403, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ptuskin, V. S., & Zirakashvili, V. N. 2005, A&A, 429, 755 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Putze, A., Derome, L., Maurin, D., Perotto, L., & Taillet, R. 2009, A&A, 497, 991 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rachen, J. P., & Biermann, P. L. 1993, A&A, 272, 161 [NASA ADS] [Google Scholar]
 Reville, B., Kirk, J. G., Duffy, P., & O'Sullivan, S. 2007, A&A, 475, 435 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reville, B., Kirk, J. G., & Duffy, P. 2009, ApJ, 694, 951 [NASA ADS] [CrossRef] [Google Scholar]
 Reynolds, S. P. 1998, ApJ, 493, 375 [NASA ADS] [CrossRef] [Google Scholar]
 Riquelme, M. A., & Spitkovsky, A. 2009, ApJ, 694, 626 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Strong, A. W., & Moskalenko, I. V. 1998, ApJ, 509, 212 [NASA ADS] [CrossRef] [Google Scholar]
 Strong, A. W., Moskalenko, I. V., & Reimer, O. 2000, ApJ, 537, 763 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Takahara, F. 1990, Prog. of Theo. Phys., 83, 1071 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, T., Uchiyama, Y., Aharonian, F. A., et al. 2008, ApJ, 685, 988 [NASA ADS] [CrossRef] [Google Scholar]
 Uchiyama, Y., Aharonian, F. A., Tanaka, T., Takahashi, T., & Maeda, Y. 2007, Nature, 449, 576 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Vietri, M. 1995, ApJ, 453, 883 [NASA ADS] [CrossRef] [Google Scholar]
 Villante, F. L., & Vissani, F. 2007, Phys. Rev. D, 76, 125019 [NASA ADS] [CrossRef] [Google Scholar]
 Vink, J., & Laming, J. M. 2003, ApJ, 584, 758 [NASA ADS] [CrossRef] [Google Scholar]
 Völk, H. J., Berezhko, E. G., & Ksenofontov, L. T. 2008, A&A, 483, 529 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wang, X.Y., Razzaque, S., Mészáros, & Dai, Z.G. 2007, Phys. Rev. D, 76, 083009 [NASA ADS] [CrossRef] [Google Scholar]
 Warren, J. S., Hughes, J. P., Badenes, C., et al. 2005, ApJ, 634, 376 [NASA ADS] [CrossRef] [Google Scholar]
 Waxman, E. 1995, Phys. Rev. Lett., 75, 386 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Yamazaki, R., Yoshida, T., Terasawa, T., Bamba, A., & Koyama, K. 2004, A&A, 416, 595 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yamazaki, R., Kohri, K., Bamba, A., et al. 2006, MNRAS, 371, 1975 [NASA ADS] [CrossRef] [Google Scholar]
 Yamazaki, R., Kohri, K., & Katagiri, H. 2009, A&A, 495, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Footnotes
 ...^{}
 We search for approximate solutions to Eqs. (39), (41) and (48), which determine (and ) for given parameters , s, and . The procedure is as follows. For a trial value of which is assumed to be constant, we solve Eqs. (39) and (41) to obtain so as to calculate E and as functions of . It is found that is almost timeindependent after , hence we can derive its average value in the epoch . Then we update the value of with Eq. (48) with the average value of . We repeat this procedure until the iteration converges.
 ... cocoon^{}
 Berezhko (2008) assumed a kind of spherical cocoon, i.e., , and used which is similar to the relation obtained for the model PS. However, when the cocoon becomes spherical, and evolve according to the adiabatic solution as in the SedovTaylor 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 and .
All Figures
Figure 1: The index 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 
Figure 2: The same as Fig. 1, but for the case PH. 

Open with DEXTER  
In the text 
Copyright ESO 2010