A&A 481, 1-16 (2008)
DOI: 10.1051/0004-6361:20078645

Analytical view of diffusive and convective cosmic ray transport in elliptical galaxies

T. Hein - F. Spanier

Lehrstuhl für Astronomie, University of Würzburg, Am Hubland, 97074 Würzburg, Germany

Received 10 September 2007 / Accepted 7 November 2007

Abstract
Context. An analytical solution of the generalized diffusive and convective transport equation is derived to explain the transport of cosmic ray protons within elliptical galaxies.
Aims. Cosmic ray transport within elliptical galaxies is an interesting element in understanding the origin of high energetic particles measured on Earth. As probable sources of those high energetic particles, elliptical galaxies show a dense interstellar medium as a consequence of activity in the galactic nucleus or merging events between galaxies. Thus it is necessary for an appropriate description of cosmic ray transport to take the diffusive and convective processes in a dense interstellar environment into account. Here we show that the transport equations can be solved analytically with respect to the given geometry and boundary conditions in position space, as well as in momentum space.
Methods. From the relativistic Vlasov equation, which is the most fundamental equation for a kinetic description of charged particles within the interstellar medium in galaxies, one finds a generalized diffusion-convection equation in quasilinear theory. This has the form of a ``leaky box'' equation, meaning particles are able to escape the confinement region by diffusing out of the galaxy. We apply here the ``diffusion approximation'', meaning that diffusion in gyrophase and pitch angle are the fastest particle-wave interaction processes. An analytical solution can be obtained using the ``scattering time method'', i.e. separation of the spatial and momentum problems.
Results. The spatial solution is shown using a generalized source of cosmic rays. Additionally, the special case of a jet-like source is illustrated. We present the solution in momentum space with respect to an escape term for cosmic ray protons depending on the spatial shape of the galaxy. For a delta-shape injection function, the momentum solution is obtained analytically. We find that the spectral index measured on Earth can be obtained by appropriately choosing of the strength of Fermi I and Fermi II processes. From these results we calculate the gamma-ray flux from pion decay due to proton-proton interaction to give connection to observations. Additionally we determine the escape-spectrum of cosmic rays. The results show that both spectra are harder than the intrinsic power-law spectrum for cosmic rays in elliptical galaxies.

Key words: acceleration of particles - diffusion - galaxies: elliptical and lenticular, cD

1 Introduction

Since their discovery by Viktor Hess in 1912, cosmic rays have been one of the biggest fields of interest in astrophysics, and yet the origin of these particles is still an open question. Fully ionized atomic nuclei reach the Earth coming from outside the solar system with very high energies up to 10²⁰ eV. Most of them with energies <10¹⁷ eV seem to originate in the Milky Way, while the highest energetic ones are considered to have an extragalactic origin (for a review see Hoerandel 2007).

The most accepted model for the origin of ultrahigh energy cosmic rays (UHECRs) is acceleration in shock fronts due to Fermi processes (Fermi 1949); hence, the main sources are gamma ray bursts (GRBs), active galaxies like active galactic nuclei (AGN) (Tavecchio 2005), or colliding galaxies. The last two are specially interesting for two reasons. First UHECRs can be generated in elliptical galaxies. Second the increase in the interstellar medium density in objects of these types influences cosmic ray transport (see Bekki & Shioya 1998). But since all elliptical galaxies have an interstellar medium due to star winds, we conclude that the study of transport processes is interesting in general (see Knapp 1999).

In particular, as a result of the GZK-effect, very close AGN are the most probable candidates for UHECR sources (see Biermann 1995). One of these nearby AGN is the giant elliptical galaxy M 87. It is proposed that this galaxy is responsible for acceleration of cosmic ray protons due to Fermi I processes (Rieger et al. 2007; Blandford & Ostriker 1978) in shockfronts within the jet (Reimer et al. 2004). For an overview of proton acceleration in jets (see Mannheim 1993). Since particle acceleration in a jet is located within active galaxies surrounded by an interstellar medium, the high energetic protons undergo physical transport processes before they escape out of the galaxy and reach the detectors on Earth, making a model for cosmic ray transport in elliptical galaxies inevitable. This helps give an answer to the question about the origin of cosmic rays.

Progress has been made in the field of modeling cosmic ray transport with the numerical description of transport processes by Owens & Jokipii (1977) and Strong & Moskalenko (1998). Nevertheless we follow the basic ideas presented in the underlying papers of Lerche & Schlickeiser (1985), Wang & Schlickeiser (1987), and Lerche & Schlickeiser (1988), who used an analytical description of cosmic ray transport. In relation to our work such, a treatment has the following advantages: the model is adequate for cosmic ray transport within any kind of elliptical galaxy including arbitrary cosmic ray sources and the physical parameters involved in our model can be easily fitted to measurements. After all, our analytical model can serve as a test case for more profound numerical models.

In this paper we solve the cosmic ray transport equation analytically with respect to a kinetic description of the interstellar plasma in elliptical galaxies. Special attention is paid to the spatial transport of charged nuclei. In addition, the solution of the momentum equation is derived to explain general properties of this model. As a result, we present illustrative examples of spatial, as well as momentum, cosmic ray transport for given sources of charged nuclei. To show the connection to observations, we calculate the gamma-ray flux from neutral pion decay. These mesons are produced by inelastic scattering processes between cosmic ray protons. The resulting power-law spectrum is slightly harder than the intrinsic one for cosmic rays. Finally, we present the escape spectrum of charged particles leaving elliptical galaxies. Similar to the gamma-ray flux, this spectrum is flatter than the intrinsic one.

2 Basic equations

To describe the propagation of cosmic ray nuclei within elliptical galaxies, we follow Lerche & Schlickeiser (1985), Lerche & Schlickeiser (1988) and Schlickeiser (2002). For the description of transport processes they use the ``diffusion approximation'', which means that the fastest particle- plasma wave interaction processes are diffusion in gyrophase and pitch angle. Thus following Jokipii (1966), Hasselmann & Wibberenz (1968), and Skilling (1975), we take an isotrope particle distribution function in momentum space. Here we idealise the interstellar medium as a homogeneous volume containing primary cosmic rays being accelerated from the thermal background medium and secondaries resulting from fragmentation of primaries having a negligible abundance in the background medium (cf. Cowsik 1980; Hayakawa 1969).

The transport of these particles at large momenta (p > 10 GeV c^-1 nucleus^-1) is described by the steady-state transport equation (e.g. Schlickeiser 1983). Such a treatment is suitable to short timescales of diffusive and convective processes compared to the dynamical timescale of the galaxy ( $t_{{\rm dyn}}\approx 10^9$ years). This is true in the case of high energetic particles. We assume a spatial diffusion coefficient K(r) of 10²⁹ cm² s^-1 at p= 1 GeV, which is slightly larger than the value measured in the Milky Way ( $K(\vec{r})=10^{28}$ cm², (cf. Schlickeiser 2002) because of diffusive processes being less effective in elliptical galaxies. Consequently we get for protons with TeV-energy a timescale of $\approx$10⁸ years. Furthermore, the dynamical age of the galaxy has to be greater than the timescale of source variability to obtain an appropriate description. This is usually given, since the size of the accretion region onto the central black hole is of the order of a few light-days so that a maximal variability timescale of some days is assumed.

At large momenta, spatial diffusion in turbulent magnetic fields dominates convection in the galactic wind, so that we find a transport equation for the phase space density  $f(\vec{r},p)$ in spatial coordinates $\vec{r}$ and in the momentum coordinate p:

$$% \mathcal{L}_{\vec{r}}f+ \mathcal{L}_{p}f+S(\vec{r},p)=0.$$ (1)

The spatial operator $\mathcal{L}_{\vec{r}}$ is defined by

$$% \mathcal{L}_{\vec{r}}(\vec{r},p) \equiv \nabla \left[K(\vec{r},p)\nabla \right]$$ (2)

containing spatial diffusion with the spatial diffusion coefficient $K(\vec{r},p)=K(\vec{r})\kappa(p)$, where $\kappa(p)$ denotes the dimensionless dependence on the momentum variable p without loss of generality. The momentum operator

$$% \mathcal{L}_{p}(\vec{r},p) \equiv p^{-2} \frac{\partial}{\p... ...ain}}-p^2\dot{p}_{{\rm loss}}\right] -\frac{1}{\tau_{\rm c}}$$ (3)

describes momentum diffusion by second-order Fermi processes (D(p)), energy gain due to first order Fermi processes ( $\dot{p}_{\rm gain}$), as well as continuous ( $\dot{p}_{\rm loss}$) and catastrophic ( $\tau_{\rm c}$) momentum loss processes. As shown in Lerche & Schlickeiser (1985), fully-ionized particles heavier than protons have the same Fermi acceleration rates as protons, so hereinafter momentum means momentum per nucleon. We are interested primarily in the behaviour of the cosmic ray primary spectrum, so we do not take secondary particles due to fragmentation of primaries into account. On the other hand, we allow fragmentation of primaries as a general loss process.

A link between spatial and momentum diffusion processes can be seen in the relation between the two diffusion coefficients

$$% D(\vec{r},p)=\frac{C_1 p^2}{K(\vec{r},p)},$$ (4)

where C₁ stands for the proportionality factor being independent of $\vec{r}$ and p. This close connection arises from the same basic physical process behind spatial and momentum diffusion: protons are scattered in pitch angle due to the magnetic fields of MHD plasma waves causing spatial diffusion along ordered magnetic field lines, whereas cyclotron damping of the electric field associated with MHD waves affects diffusion in momentum space. For that reason it is necessary to solve this model in spatial coordinates as well as in momentum coordinates to get an acceptable description of transport processes in elliptical galaxies.

3 ``Scattering time'' method

We use the ``scattering time'' method proposed by Sunyaev & Titarchuk (1980) to get an important class of exact analytical solutions of Eq. (1) following Wang & Schlickeiser (1987). This implies, that the spatial and momentum operators can be separated as

$$% \mathcal{L}_{\vec{r}}(\vec{r},p)=h(p)~\mathcal{O}_{\vec{r}}... ...m} \mathcal{L}_{p}(\vec{r},p)=g(\vec{r})~\mathcal{O}_{p}(p).$$ (5)

For ease of exposition, the source function $S(\vec{r},p)$ is also a product of two separable functions, i.e.,

$$% S(\vec{r},p)=q(\vec{r})Q(p).$$ (6)

As an aside we note that the requirement Eq. (5) is trivially fulfilled, if $K(\vec{r})$ is constant ( $K(\vec{r})=K_0$) and $\dot{p}_{\rm gain}$, $\dot{p}_{\rm loss}$, as well as $\tau_{\rm e}$, are all independent of spatial variables. Thus we use the following model containing constant factors a_n with n=1,2,3,4 for the proportionality factors independent of spatial and momentum coordinates. First we take $D(p)=a_2 p^2/\kappa(p)$ to describe the momentum diffusion coefficient solely in momentum space. Furthermore we assume

$$% \dot{p}_{{\rm gain}}=a_1 p /\kappa(p)$$ (7)

telling that first-order Fermi acceleration is related to the (momentum dependend) spatial diffusion coefficient due to MHD plasma-wave-scattering interactions within the acceleration process. The continuous loss term  $\dot{p}_{\rm loss}$ is independent of spatial coordinates leading to

$$% \dot{p}_{{\rm loss}}(p)=-a_3 \rho (p).$$ (8)

Similarly, we set

$$% \tau_{c} (p)=a_{4}^{-1}\theta (p)$$ (9)

for the catastrophic loss time due to fragmentation.

Under these conditions, in addition to Eqs. (5) and (6), we can find the formal mathematical solution of Eq. (1) as a convolution of the spatial and momentum solution functions  $T(\vec{r})$ and M(p) following de Freitas Pacheco (1971):

$$% f(\vec{r},p)=\int_0^{\infty}{\rm d}u~T(\vec{r},u)M(p,u),$$ (10)

where $T(\vec{r},u)$ has to satisfy the given spatial boundary conditions, and

$$% \frac{\partial T}{\partial u}=\frac{1}{g(\vec{r})}\mathcal{O}_{\vec{r}} T$$ (11)

with

$$% \mathcal{O}_{\vec{r}}=\nabla \left[ K_0 \nabla \right] \equiv K_0 \Delta,\hspace{0.2cm} g(\vec{r})=1$$ (12)

and the conditions

  

$$% T(\vec{r},u=0) = q(\vec{r})/g(\vec{r})=q(\vec{r})$$ (13)

$$% T(\vec{r},u=\infty) = 0.$$ (14)

Here, M(p,u) has to satisfy the given spatial boundary conditions, and

$$% \frac{\partial M}{\partial u}=\frac{1}{h(p)}\mathcal{O}_{p} M$$ (15)

with

$$% \mathcal{O}_p=\frac{1}{p^2} \frac{\partial}{\partial p} \le... ...ght]- \frac{a_4}{\theta_4(p)};\hspace{0.2cm} h(p)=\kappa(p)$$ (16)

and the conditions

  

$$% M(p,u=0) = Q(p)/h(p)=Q(p)/\kappa(p)$$ (17)

$$% M(p,u=\infty) = 0.$$ (18)

As a consequence of the formal mathematical solution, we note that we have to solve two partial differential Eqs. (11) and (15) instead of the much more complicated differential Eq. (1).

4 Results

4.1 Spatial solution

The most convenient way to find the formal solution of Eqs. (11) and (15) is to start with the spatial problem. As can be seen from Eq. (12), the spatial operator  $\mathcal{O}_{\vec{r}}$ is of Sturm-Liouville type (cf. Arfken & Weber 2005) and therefore has a complete eigenfunction system  $E_i(\vec{r})$. As a consequence, the solution function  $T(\vec{r},u)$ can be expanded in this orthonormal system as

$$% T(\vec{r},u)=\sum_i A_i(\vec{r})~ {\rm e}^{-\lambda_i^2 u}.$$ (19)

Here the $A_i(\vec{r})$ are defined by

$$% A_i(\vec{r})=\alpha_i E_i(\vec{r}),$$ (20)

implying that the coefficients $\alpha _i$ weight each eigenfunction  $E_i(\vec{r})$. The $\lambda_i$ denote the eigenvalues of $\mathcal{O}_{\vec{r}}$ implying the special spatial geometry.

	Figure 1: Schematical view of an ellipse with its fundamental properties.
Open with DEXTER

	Figure 2: Left: gray-shaded plane cuts ellipse leading to a cut view like the graph on the left-hand side. Right: direction of unit vectors of the variables $\xi$, $\eta$, and $\phi$ in prolate spheroidal coordinates. From Weisstein (1999).
Open with DEXTER

The shape of elliptical galaxies is adjusted to the cosmic ray transport Eq. (12) using prolate spheroidal coordinates as they are defined by Abramowitz & Stegun (1972):

$$% \xi=\frac{r_1+r_2}{2f},\hspace{0.2cm} \eta=\frac{r_1-r_2}{2f}\cdot$$ (21)

As can be seen from Fig. 1 r₁ and r₂ are the distances to the foci of the confocal ellipse, where 2f denotes the distance between the two foci F₁ and F₂. Additionally, we use the variable $\phi$ for the usual azimuthal dependence like in spherical coordinates. The following relations give the relation between these coordinates and the semi-major axis a and the semi-minor axis b, respectively:

$$% a=f\xi,\hspace{0.2cm} b=f\sqrt{\xi^2-1}.$$ (22)

The numerical excentricity e has a direct relationship to the coordinate $\xi$ via

$$% e \equiv f/a=1/\xi.$$ (23)

The variable $\eta$ is defined as $\eta=\cos~\theta$ with $\theta$ the angle between the line on which the foci lie, and an arbitrary point on the ellipse. As a result the variables are defined in the range

$$% \xi~ \in~[1;\infty[,\hspace{0.2cm} \eta~ \in~[-1;1],\hspace{0.2cm} \phi~ \in~[0;2\pi].$$ (24)

Figure 2 shows an illustration of the definition of the three spatial variables $\xi$, $\eta$, and $\phi$. Because of these definitions, we can write Eq. (12) as

$$% \frac{\partial{T}}{\partial{u}}= \frac{K_0}{f^2(\xi^2-\eta^... ...1)(1-\eta^2)}\frac{\partial^2 T}{\partial \phi^2}\right\}\cdot$$ (25)

Here we used the Laplacian in prolate spheroidal coordinates. The general solution can be obtained by consecutive separation of variables (see Appendix A):

 

$$% T(\xi,\eta,\phi, u) \sum_{m,n} R_{mn}^{(1)}(c,\xi)\times S_{mn}^{(1)}(c,\eta)\times \cos(m\phi)\times \exp(-k^2u)$$ =

$$\sum_{m,n} \left\{\left[\sum_{r=0,1}^{\infty ^{\prime}} \frac{(2m... ...}{2}}(c\xi) \hat{d}_r^{mn}(c)P_{m+r}^{m}(\eta) \cos(m\phi)\exp(-k^2u) \right\}.$$ = (26)

Here the $P_{m+r}^{m}(\eta)$ denotes associated Legendre functions of the first kind, order m+r, and the $J_{n+\frac{1}{2}}(c\xi)$denote Bessel functions of the first kind and of order $n+\frac{1}{2}$. The sum is extended over even values of r as indicated by the mark $^{\prime}$. The factor c is given by $c=\frac{fk}{\sqrt{K_0}}$.

To define reasonable boundary conditions, we assume a ``leaky box'' model. Cosmic ray particles are trapped by disordered magnetic fields within the confinement region of an elliptical galaxy. In this they undergo diffusive and convective movements. At the edge of the box, leakage out of the confinement area is possible.

	Figure 3: Graphic demonstration of the cosmic ray particle density given by the solution Eq. (30) with the weighting factors (31). The cosmic ray particle density is normalised and given in arbitrary units. As an effect of chosen coordinates only ``one half'' of the galaxy is visible. For illustration we plotted $T(f,\eta ,u=0)$ with a constant $\xi \equiv 1.2$ specifying an E4 elliptical galaxy.
Open with DEXTER

As an illustrative example for spatial boundary conditions, we show the solution depending on a constant source function over the elliptical galaxy in Appendix B. To be more specific, we take a jet-like source function here. The jet points in the direction  $\eta_{{\rm inj}}$ (represented by a Dirac delta function) with a length scale chosen to be $f_{{\rm max}}$ for any choice of $\xi$being smaller than an arbitrary maximum value of the confinement region  $\xi_{\rm c}$. Particles leak out at the edge of this region. Such a boundary condition is known in the literature as a ``free-escape'' condition. For a realistic assumption we decide to let the jet end smoothly (see the ``Fermi'' function in Eq. (29)). We neglect any dependence on $\phi$ for an adequate illustration. These conditions are taken into account by

$$% T(\xi=\xi_{\rm c},\eta,\phi,u)=0,$$ (27)

by a periodical boundary condition in $\eta$

$$% T(\xi,\eta = \-1,\phi,u)=T(\xi,\eta = +1,\phi,u),$$ (28)

and by

$$% T(\xi,\eta,u=0)=q_0\frac{\delta(\xi-\xi_{\rm c})\delta(\eta_{{\rm inj}}-\eta)} {\left[\exp(4(f-f_{{\rm max}})+1\right]}\cdot$$ (29)

Finally after some calculations in which Eq. (26) has to match the boundary conditions Eqs. (27)-(29), we derive the general solution (cf. Appendix B):

$$% E_{i,r}(\vec{r}) \equiv T(\xi,u)=\sum_{i=1}^{\infty}\sum_{r... ...\exp \left[ -\frac{K_0 y_{ir}^2}{\xi_{\rm c}^2 f^2}u \right ].$$ (30)

The weighting factors are

$$% \alpha_{ir}=\frac{q_0 P_r^0(\eta_{{\rm inj}})}{\exp\left[4(... ...i_{\rm c}^2}{2}\left[ J_{r+\frac{3}{2}}(y_{ir})\right]^2}\cdot$$ (31)

Here the y_ir stands for zeros of $J_{r+\frac{1}{2}}(x)$. The sum over r is extended over even values of this parameter. We see that the general solution of Eq. (25) with respect to spatial boundary conditions indeed has the form of the eigenfunction expansion Eqs. (19) with (20). This solution is shown in Fig. 3 with respect to 100 spatial eigenfunctions ( $i_{{\rm max}}=10$, $r_{{\rm max}}=18$). The jet is responsible for the cosmic ray particles distributed over the whole galaxy due to diffusive processes. It is broadened with a larger distance to the centre. This effect is caused by the chosen geometry, whereas the loss of magnetic collimation in real astrophysical sources can provide this. To illustrate the shape of an elliptical galaxy we plotted $T(f,\eta ,u=0)$ with a constant $\xi \equiv 1.2$ specifying an E4-galaxy. Using a more complicated source function, a better physical description of particle distribution within elliptical galaxies can be obtained.

	Figure 4: Normalised weighting factors $\alpha _i$ dependent on the number of zeros yi in Eq. (32).
Open with DEXTER

4.2 Consistency checks of the formal spatial solution

To get a better understanding of our model, we prove the formal mathematical solution (Eq. (26)) of the spatial cosmic ray transport equation. Our discussion is related to the solution found in the illustrative example in Appendix B, but can similarly done with Eq. (26):

$$% T(\xi,u)=\sum_{i=1}^{\infty}\alpha_i \frac{1}{\sqrt{\xi}}J_... ...ht) \exp \left[ -\frac{K_0 y_i^2}{\xi_{\rm c}^2 f^2}u \right ]$$ (32)

with the weighting factors

$$% \alpha_i=\frac{q_0 \int_1^{\xi_g}\xi^{3/2}J_{\frac{1}{2}}\l... ...c{\xi_{\rm c}^2}{2}\left[ J_{\frac{3}{2}}(y_i)\right ]^2}\cdot$$ (33)

Here we used periodical boundary conditions for $\eta$ and $\phi$

$$% T(\xi,\eta = \-1,\phi,u)=T(\xi,\eta = +1,\phi,u),$$ (34)

$$% T(\xi,\eta,\phi=2\pi,u)=T(\xi,\eta,\phi=0,u),$$ (35)

and

$$% T(\xi,\eta,\phi,u=0)=q_0 \Theta(\xi_{\rm g}-\xi).$$ (36)

First, this solution has to accomplish the given spatial boundary conditions. As suggested in Eq. (36) we used a constant source function over the whole size of the galaxy in order to lose all angular dependencies. This behaviour is caused by the source function not having any dependence in $\eta$ and $\phi$, but also due to the periodic boundary conditions chosen for these variables. Such an effect can also be seen in a spherical or disc geometry. To the spherical geometry as performed by Schlickeiser et al. (1987), we added angular dependencies and a constant source function over the galaxy as an example. The weighting factors of the general solution functions turn out to disappear in any case except for the angular separation constants being zero (the angular solution functions are spherical harmonics in this case), which means that there is no dependence on these variables as it is expected.

Second, we show that the sum of solution functions (cf. Eq. (26)), together with the weighting factors in addition to the given spatial boundary conditions converges. For the illustrative example in Appendix B, the normalised expansion coefficients $\alpha _i$ are shown in Fig. 4. To obtain applicable results it is essential to include only the first few ones depending on the requested accuracy. The resulting solution function is shown in Fig. 5. In the special case of a constant source function, the solution reproduces the boundary condition  $\theta(\xi_{\rm g}-\xi)$ for the variable $\xi$. We used $\xi _{\rm g}=4$ for illustration. If we take many eigenfunctions into account, we see at the discontinuity points $\xi=0$ and $\xi=\xi_{\rm g}=4$ an oscillatory phenomenon (Gibbs phenomenon). But in this case a closer solution for the given boundary condition is obtained.

	Figure 5: Spatial solution Eq. (32) in addition with weighting factors $\alpha _i$ (Eq. (20)) as a function of $\xi$. As an illustration, we put $\xi _{\rm g}=4$. For the solution represented by the solid line we included 200 eigenvalues, whereas in the dashed solution only 10 eigenvalues are taken into account.
Open with DEXTER

In analytical calculations, it is a common assumption to include only the first spatial eigenvalue $\lambda_1$. This one is associated with the longest escape timescale being the most important one for modelling escape of particles out of the galaxy. For numerical purposes, the computing time gives an upper limit to the possible number of eigenvalues. Furthermore, the spatial solution as performed in this paper has to match the formal solution for a spherical geometry within the limit of small ellipticity ( $e \rightarrow 0$). Schlickeiser et al. (1987) found as the spatial solution

$$% T(R,u)=\sum_{m=1}^{\infty}c_m \frac{1}{\sqrt{R}}J_{\frac{1}... ...c{R}{R_1}\right) \exp\left(-\frac{K_0^2y_m^2}{R_1^2}u \right),$$ (37)

which already shows affinity to our solution Eq. (32). The y_m denotes the zeros of $J_{\frac{1}{2}}(x)$ (here we corrected Schlickeiser et al. 1987) and R₁ describes the edge of the galaxy. Generally we can write for both solutions:

$$% T(X,u)=\sum_{i=1}^{\infty} \alpha_i Z_i(c_iX)\exp(-\lambda_i^2u).$$ (38)

The asymptotic behaviour of the solution Eq. (32) is with $Z(c_iX)\equiv \xi^{-1/2} J_{\frac{1}{2}}(c_i\xi)$

$$% Z(c_i\xi)\,\includegraphics[width=1.2cm]{eq1}\, \frac{1}{c_i\xi}\cos \left(c_i\xi-\frac{1}{2}\pi \right).$$ (39)

The limit $c_i\xi\rightarrow \infty$ satisfies a spherical symmetry with a numerical excentricity equal to zero (cf. Eq. (23)). Furthermore with this limit, the semi-major axis $c_i\xi$ goes to infinity. Performing the same limit to Eq. (37) with $Z(c_iX)\equiv R^{-1/2}J_{\frac{1}{2}}(c_iR)$, we get

$$% Z(c_iR)\,\includegraphics[width=1.2cm]{eq2}\, \frac{1}{c_iR}\cos \left(c_iR-\frac{1}{2}\pi \right),$$ (40)

which is equal to Eq. (39). As a consequence we showed that our solution functions embody the spherical geometry within the limit of numerical excentricity equal to zero. It is also straightforward to show that the weighting factors have the same structure as those in Schlickeiser et al. (1987) so that the general solution is correct.

   
4.3 Momentum solution

For the formal momentum solution, it is necessary to have a closer look at the spatial solution. As noted above, the eigenfunction expansion Eq. (19), in addition to Eq. (20), is indeed the best way to solve the spatial transport Eq. (11). Inserting this expansion into the convolution Eq. (10), we can write

$$% f(\vec{r},p)=\sum_i ~A_i(\vec{r})R_i(p)$$ (41)

with

$$% R_i(p) \equiv \int_0^{\infty} {\rm d}u ~M(p,u){\rm e}^{-\lambda_i^2u}.$$ (42)

Consequently the momentum solution R_i(p) obeys the ordinary differential equation

$$% \mathcal{O}_p R_i(p) - \lambda_i^2 g(p)R_i(p)=-Q(p),$$ (43)

which can be seen after multiplying Eq. (15) by $T(\vec{r},u)$ as given in Eq. (19) and integrating over the convolution variable u from 0 to $\infty$. Each spatial eigenvalue  $\lambda_i^2$ enters this equation in the form of an inverse catastrophic loss time for particle escape out of the galaxy. Therefore Eq. (43) is called the ``leaky box'' equation in momentum space.

As the result of the eigenfunction Eq. (19), we have to solve one ordinary differential equation for each spatial eigenvalue instead of the partial differential Eq. (15). Following Lerche & Schlickeiser (1988), we introduce some simplifying assumptions in order to find analytic solutions of Eq. (43). We take in Eq. (16) $\theta_4=1$, $\rho(p)=p$, and $\kappa(p)=(p/p_1)^{s}$ where p₁ is a normalisation value. These assumptions imply that the most important continuous loss process in elliptical galaxies at energies >10 GeV is adiabatic energy loss due to a high galactic wind gradient, whereas pion production losses are neglected because of the low number density of HI and HII. The fragmentation lifetime ($\propto$$\theta_4$) is independent of momentum, and the momentum dependence of the diffusion coefficient is defined by the parameter s=2-q, where q is the spectral index of the magnetic turbulence power spectrum. Therefore we get

$$% \frac{p_{1}^{s}}{p^{2+s}} \frac{\partial}{\partial p} \left... ...}{p^{s}} +\lambda_{i}^{2} \right] R_{i}=-Q(p)p^{-s} p_{1}^{s}.$$ (44)

This ordinary differential equation can be solved by a standard technique taking the finiteness of R_i at $p\rightarrow 0$ and $p\rightarrow\infty$ into account. The formal mathematical solution is taken from Lerche & Schlickeiser (1988). However, we found after correction of some minor mistakes:

 

$$(a_2p_1^s)^{-1}\frac{\Gamma\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a... ...} \times p^a \exp\left[-\frac{\beta_1+(\beta_1^2+4\psi_i)^{1/2}}{2s}p^s \right]$$                                   R_i(p) =

$$\times\left\{U\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1}{... ..._i) ^{1/2}},\frac{a+3}{s},\frac{(\beta_1^2+4\psi_i)^{1/2}}{s}p^s\right] \right.$$

$$\left.\times\int_0^pdp_0~p_0^{s+1}Q(p_0)\exp\left[\frac{\beta_1-(\beta_1^2+4\psi_i)^{1/2}}{2s} p_0^s\right] \right.$$

$$\left. \times~ M\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1... ...) ^{1/2}},\frac{a+3}{s},\frac{(\beta_1^2+4\psi_i)^{1/2}}{s}p_0^s\right] \right.$$

$$\left.+~M\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1}{s(\be... ..._i) ^{1/2}},\frac{a+3}{s},\frac{(\beta_1^2+4\psi_i)^{1/2}}{s}p^s\right] \right.$$

$$\left.\times\int_p^{\infty}dp_0~p_0^{s+1}Q(p_0)\exp\left[\frac{\beta_1-(\beta_1^2+4\psi_i)^{1/2}}{2s} p_0^s\right] \right.$$

$$\left. \times~ U\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1... ... ^{1/2}},\frac{a+3}{s},\frac{(\beta_1^2+4\psi_i)^{1/2}}{s}p_0^s\right]\right\}.$$ (45)

Here, we took for brevity

$$% \beta_1 \equiv a_3/(a_2p_1^s),$$ (46)

$$% a \equiv a_1/a_2,$$ (47)

$$% \xi_1\equiv a_4/(a_2p_1^s),$$ (48)

and

$$% \psi_i\equiv \lambda_i^2/(a_2 p_1^{2s}).$$ (49)

The functions U and M denote confluent hypergeometric functions of first (Kummer) and second (Whittaker) order, respectively.

For the special case of a $\delta$-function injection of cosmic ray particles,

$$% Q(p)=\delta(p-p_{{\rm inj}}),$$ (50)

the solution Eq. (45) can be evaluated analytically. Therefore we get for $p<p_{{\rm inj}}$

 

$$(a_2p_1^s)^{-1}\frac{\Gamma\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a... ...} \times p^a \exp\left[-\frac{\beta_1+(\beta_1^2+4\psi_i)^{1/2}}{2s}p^s \right]$$                              R_i(p) =

$$\times\left\{M\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1}{... ...eft[\frac{\beta_1+(\beta_1^2+4\psi_i)^{1/2}}{2s} p_{{\rm inj}}^s\right] \right.$$

$$\left. \times ~U\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1... ...frac{a+3}{s},\frac{(\beta_1^2+4\psi_i)^{1/2}}{s}p_{{\rm inj}}^s\right]\right\},$$ (51)

and for $p>p_{{\rm inj}}$

 

$$(a_2p_1^s)^{-1}\frac{\Gamma\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a... ...} \times p^a \exp\left[-\frac{\beta_1+(\beta_1^2+4\psi_i)^{1/2}}{2s}p^s \right]$$                            R_i(p) =

$$\times\left\{U\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1}{... ...eft[\frac{\beta_1+(\beta_1^2+4\psi_i)^{1/2}}{2s} p_{{\rm inj}}^s\right] \right.$$

$$\left. \times ~M\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1... ...frac{a+3}{s},\frac{(\beta_1^2+4\psi_i)^{1/2}}{s}p_{{\rm inj}}^s\right]\right\}.$$ (52)

Because of the given structure of the solution Eq. (45), it is convenient to introduce the momentum value $p_{\star}$ via

$$% p_{\star}=\left(\frac{s}{\sqrt{\beta_1^2+4\psi_i}}\right)^{1/s}\cdot$$ (53)

This parameter indicates the momentum value above which the cosmic ray spectrum cuts off exponentially due to adiabatic losses ($\beta_1$) and escape losses ($\psi_i$). Now we can give some asymptotic behaviour for the solution:

In the case of low values of the momentum ( $p \ll p_{\star}$), we find, according to Abramowitz & Stegun (1972),

$$% R_i(p \ll p_{\star})\simeq p^{+s-3},$$ (54)

which is a flat power-law spectrum at very small momenta. For large arguments ( $p \gg p_{\star}$) of the Kummer function in Eq. (52), we get

$$% R_i(p \gg p_{\star}) \propto p^{\frac{a-3}{2}+\frac{a+3}{2}... ...[ -\frac{\beta_1+(\beta_1^2+4\psi_i)^{1/2}}{2s}p^s\right]\cdot$$ (55)

In this situation, that the exponential cutoff is not negligible. In the special case of dominating adiabatic losses, we take $\beta_1^2 \gg 4\psi_i$ to get

$$% R_i\left( p \gg (s/\beta_1)^{1/s}\right)=p^{+a-\xi_1/\beta_1} \exp\left[ -\frac{\beta_1}{s}p^s\right].$$ (56)

Also in this case the exponential cutoff dominates the spectral behavior. To summarize these results we expect a flat power-law spectrum at very low momenta, whereas the exponential cutoff dominates at very high momenta.

	Figure 6: Cosmic ray spectrum (grey line) for a delta shape injection at $p=p_{{\rm inj}}$ obtained by superposing the contributions from individual modes (i=1,2,3). The first mode R1(p) modelling the longest escape losstime is responsible for the high energy cutoff whereas the slope of the power law is determined by the most effective energy loss process (adiabatic loss).
Open with DEXTER

As an illustrative example a spectrum is modelled such that the result is a power law-spectrum with an exponent like the one observed from high energetic cosmic rays on Earth. This is demonstrated in Fig. 6. Here particles are injected at a momentum value of $p=p_{{\rm inj}}$ with a delta-shape injection function. Furthermore we assume an isotropic Kolmogorov turbulence model, i.e., the power law index of the turbulence is q=5/3. We found a steeper power-law spectrum as the predicted one from Eq. (54) over about two decades in momentum. In this regime the exponential cutoff already plays a nonvanishing role. As the result of the small momentum dependence of the exponential term in Eq. (51) ( $s\equiv 1/3$for Kolmogorov-like turbulence), a power-law spectrum over just two decades in momentum is obtained. Performing some delta-shape injections at increasing momenta with adequate normalisation values a power-law spectrum over a wide range in momentum is obviously possible.

The parameter that affects the final power-law index of the cosmic ray spectrum is the parameter a that describes the ratio of energy gains from the Fermi I process and the Fermi II momentum diffusion process. To match the -4.8 spectrum from observations, the value has been chosen as a=250. The adiabatic losses are assumed to be as strong as the Fermi II energy gains. In this context we need to remember that the solution function R_i(p) has to be multiplied by $4\pi p^2$ to obtain the number of particles at the momentum p:

$$% N_i(p)=4\pi p^2 R_i(p).$$ (57)

The spectrum shown in Fig. 6 comprises the first 10 eigenfunctions, and the cumulative solution (grey line) is given by $R_{\rm tot}(p)=R_1(p)+ R_2(p)+...+R_{10}(p)$. It can be seen that the spectrum is almost completely dominated by the first eigenfunction. Especially the cutoff is dominated by this first Eigenfunction. That result has to be compared with, e.g., Fig. 2 of Lerche & Schlickeiser (1988), where the situation in spiral galaxies is described. Both, elliptical and spiral galaxies have in common that the the cutoff is described by the first Eigenfunction. The difference between these two situations is that, at low energies, the higher Eigenfunctions describe the spectrum in spiral galaxies. The resulting power law is made by the sum of the single Eigenfunctions. In elliptical galaxies the situation is different, because the spectral shape is by oneself dominated by the first Eigenfunction.

The reason for this different behavior compared to spiral galaxies is based on the loss time scales. While in the latter the escape loss processes dominate because of the small galactic height compared to the radial size, this is different in elliptical galaxies, where the smallest ``edge'' of the confinement region is given by the semi-minor axis being larger than the galactic height in spirals. Therefore the dominating loss process in elliptical galaxies is adiabatic loss.

   
5 Connection to observations

The results of the solutions of the transport equations have been explained in the previous section. For astrophysical scenarios it is important to compare these with observations. Here we concentrate on the gamma-ray spectrum due to neutral pion decay and on the escaping cosmic ray spectrum from elliptical galaxies.

   
5.1 Gamma-ray flux from pion decay

The found solutions allow us to calculate the gamma-ray flux from pion decay. The pions are mainly produced due to proton-proton interactions, so we concentrate on the reaction $p+p\rightarrow p+p+x\times \pi^0$. The neutral pion decays after the mean lifetime $\tau=9$ $\times$ 10^-17 s in two gamma-photons. Other interaction channels can be treated analogously. Then the derivation is straightforward, but tedious.

We have to convolve the high-energy proton spectrum, Eq. (26) multiplied with Eq. (52), with the pion power of a single proton. The latter one is approximately given by (cf. Schlickeiser 2002; Mannheim & Schlickeiser 1994):

$$% P_{\pi}(\gamma_{\pi},\gamma_{p},\xi,\eta,\phi)= c\gamma_{\p... ...i}-\gamma_{p}^{3/4})H\left[\gamma_p-\gamma_{{\rm thr}}\right].$$ (58)

In this equation, $\Xi\sigma_{pp}^{\pi^0}(\gamma_p)$ gives the multiplicity $\Xi$ for neutral pion production in p-p interactions and $\sigma_{pp}^{\pi^0}(\gamma_p)$ the total cross-section. Within the high energy limit, this factor is given by (cf. Schlickeiser 2002):

$$% \Xi\sigma_{pp}^{\pi^0}(\gamma_p)=\sigma_{0,pp}^{\pi^0}(\gamma_p-1)^{0.53}.$$ (59)

The factor $\gamma_{{\rm thr}}$, together with the Heaviside function, gives the threshold value for the proton energy to produce neutral pions ( $\gamma_{{\rm thr}}m_pc^2=1.22$ GeV). Therefore we get, for the pion source function,

 

$$% Q_{\pi^0}(\gamma_{\pi},\xi,\eta,\phi) \frac{1}{\gamma_{\pi}m_{\pi}c^2}\int_1^{\infty}~{\rm d}\gamma_p T... ...eft[ -\frac{\beta_1+(\beta_1^2+4\psi_i)^{1/2}}{2s}(m_{p}c)^s\gamma_{p}^s\right]$$ =

$$\times U\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1}{s(\bet... ...,\frac{a+3}{s},\frac{(\beta_1^2+4\psi_i)^{1/2}}{s}(m_{p}c)^s\gamma_{p}^s\right]$$

$$\times~ c\gamma_{\pi}T(\xi,\eta,\phi,u)\sigma_{0,pp}^{\pi^0}(\gam... ...\delta(\gamma_{\pi}-\gamma_{p}^{3/4})H\left[\gamma_p-\gamma_{{\rm thr}}\right].$$ (60)

The power-law factor p^a+2 is due to Eq. (57) and we used

 

$$% C_{{\rm inj}} 4\pi(m_pc)^{a+2} \times(a_2p_1^s)^{-1}\frac{\Gamma\left[\frac{a+3... ...1}\exp\left[\frac{\beta_1+(\beta_1^2+4\psi_i)^{1/2}}{2s} p_{{\rm inj}}^s\right]$$ =

$$\times ~ M\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1}{s(\b... ...{1/2}},\frac{a+3}{s},\frac{(\beta_1^2+4\psi_i)^{1/2}}{s}p_{{\rm inj}}^s\right].$$ (61)

Performing the integration we get

 

$$% \hspace*{1cm} Q_{\pi^0}(\gamma_{\pi},\xi,\eta,\phi) \frac{4}{3}\frac{T^2(\xi,\eta,\phi)\sigma_{0,pp}^{\pi^0}}{m_{\pi ... ...-\frac{\beta_1+(\beta_1^2+4\psi_i)^{1/2}}{2s}(m_{p}c)^s\gamma_{p}^{4s/3}\right]$$ =

$$\times U\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1}{s(\bet... ...ma_{p}^{4s/3}\right]\times H\left[\gamma_{\pi}-\gamma_{{\rm thr}}^{3/4}\right].$$ (62)

It is obvious that the pion source function is proportional to the squared spatial distribution of the high-energy protons since pion production is a two-body process. Therefore the brightest luminosity is strongly correlated with the highest proton density in spatial coordinates because the mean lifetime of the neutral pion is extremely short so that the photons are produced very close to the p-p interaction region.

With this result, the differential gamma-ray source function is given by (cf. Schlickeiser 2002)

$$% Q_{\gamma}(E_{\gamma},\xi,\eta,\phi)=2\int_{\gamma_{{\rm mi... ...{\pi^0} (\gamma_{\pi},\xi,\eta,\phi)} {\sqrt{\gamma_{\pi}-1}},$$ (63)

where $\gamma_{{\rm min}}=\frac{E_{\gamma}}{m_{\pi}c^2}+\frac{m_{\pi}c^2}{4E_{\gamma}}$denotes the minimum energy of the gamma-ray photons after the decay of the neutral pion. We can state that the differential gamma-ray source function in the high-energy regime ( $E_{\gamma} \gg m_{\pi^0}c^2$), where the spectrum shows a power-law dependence, is approximately given by

 

$$% Q_{\gamma,\pi^0}(E_{\gamma},\xi,\eta,\phi) \propto \gamma_{\pi}^{4a/3+3}(\gamma_{\pi}^{4/3}-1)^{0.53}\exp\left[ -\frac{\beta_1+(\beta_1^2+4\psi_i)^{1/2}}{2s}(m_{p}c)^s\gamma_{p}^{4s/3}\right]$$

$$\times U\left[\frac{a+3}{2s}+\frac{\xi_1-[(3+a)/2]\beta_1}{s(\bet... ...{a+3}{s},\frac{(\beta_1^2+4\psi_i)^{1/2}}{s}(m_{p}c)^s\gamma_{p}^{4s/3}\right],$$ (64)

but the integral Eq. (63) cannot be performed analytically.

Consequently, as discussed before, we assume a power-law spectrum over a wide range in momentum until a maximum value of $\gamma_{{\rm p,~max}}m_pc^2$:

$$% N(\gamma_p)=N_0\gamma^{-z}H\left[\gamma_{{\rm p,~max}}-\gamma_p \right].$$ (65)

After performing the convolution (cf. Eq. (60)) the pion source function reads

$$% Q_{\pi^0}(\gamma_{\pi},\xi,\eta,\phi)=\frac{4}{3}\frac{T^2(... ...right] H\left[\gamma_{{\rm p,~max}}^{3/4}-\gamma_{\pi}\right].$$ (66)

From this calculation we get due to Eq. (63) in the high-energy limit ( $\sqrt{\gamma_{\pi}^2-1}\simeq \gamma_{\pi}$and $(\gamma_{\pi}^{4/3}-1)^{0.53}\simeq \gamma_{\pi}^{4/3\times 0.53}$) for gamma-ray energies above 130 GeV

$$% Q_{\gamma}(E_{\gamma},\xi,\eta,\phi)\simeq \frac{8}{3}\frac... ...max}}^{-z+0.78}-\gamma_{{\rm min}}^ {-4/3z+1.04}\right)\right]$$ (67)

as the gamma-ray source function. Thus the differential gamma-ray source function at high energies ( $E_{\gamma} \gg m_{\pi^0}c^2$) is given by the relation

$$% Q_{\gamma,\pi^0}(E_{\gamma},\xi,\eta,\phi)\propto E_{\gamma}^{-s_{\pi}}=E_{\gamma}^{-4/3z+1.04}$$ (68)

after Eq. (66). Taking an injection power-law index z=2.8, we find that $Q_{\gamma,\pi^0}(E_{\gamma},\xi,\eta,\phi)\propto E_{\gamma}^{-2.69}$, which is slightly harder than the injected proton spectrum. Therefore an observed gamma-ray flux from an elliptical galaxy gives constraints to the high-energy cosmic ray spectrum within that object since neutral pions are produced in the whole confinement volume. Therefore we expect an elliptical galaxy be an extended source of gamma rays.

5.2 Escape-spectrum of cosmic rays

We next derive the spectrum leaking out of elliptical galaxies. Therefore the mean free path  $\lambda _{{\rm mfp}}(p)$ for single scattering events between charged particles and plasma waves is needed. According to Schlickeiser (2002), this characteristic length is given by

$$% \lambda_{{\rm mfp}}(p)\equiv\frac{3K(\vec{r},p)}{v}=\lambda_{{\rm0,~mfp}}p^{2-q},$$ (69)

where $K(\vec{r},p)$ and v denote the spatial diffusion coefficient and the particle velocity respectively. Again the parameter q gives the momentum dependence of the turbulence spectrum. Here we use v=c for highly relativistic cosmic rays. Figure 7 gives the chosen geometry for calculating of the escape spectrum. We assume that cosmic rays will escape out of the galaxy, if their mean free path is larger than the distance to the edge of the galaxy. This border value is given by (cf. Sect. 4.3)

$$% r_{\rm g}(\xi,\eta,\phi)=f_{\rm e}\sqrt{(\xi\eta)^2+(\xi^2-1)(1-\eta^2)}.$$ (70)

	Figure 7: Schematic view of the geometry for calculating the escape-spectrum from elliptical galaxies. Particles within the gray-shaded area leave the galaxy by chance, if their mean free path  $\lambda _{{\rm mfp}}(p)$ is longer than the distance to the edge of the galaxy.
Open with DEXTER

All particles in the gray-shaded area will escape the galaxy, if $r_{\rm g}(\xi,\eta,\phi)-d_1\leq\lambda_{{\rm mfp}}(p)$, where we neglect a geometrical factor of order unity. The width of this spherical shell is constant in every direction so that $f_{\rm e}\sqrt{\xi^2-1}-d_1\equiv f_{\rm e}\xi-d_2$, as indicated in Fig. 7. The rate of the escaping high energetic particles is approximatively given by the inverse escape time

$$% \frac{K(\vec{r},p)}{r^2(\xi,\eta,\phi)}=\frac{c\lambda_{{\rm mfp}}(p)}{3f^2\left[(\xi\eta)^2+(\xi^2-1)(1-\eta^2)\right]},$$ (71)

where we neglect the geometrical factor for the direction of the trajectory of cosmic ray particles. Therefore the total number of particles escaping out from the gray-shaded area into the intergalactic medium per unit time and per unit momentum is given by

$$% N(p)\propto4\pi p^2\int_0^{2\pi}\int_{-1}^{1}\int_{d_1}^{r_... ...\xi^2-1)(1-\eta^2)\right]}f^2{\rm d}f {\rm d}\eta {\rm d}\phi.$$ (72)

The upper boundary in the f-integral gives the radius of the galaxy with respect to the variables $\xi$ and $\eta$. As an illustration we choose a constant cosmic-ray distribution like Eqs. (106) with (107) (see Appendix B). This one is independent of the variables $\eta$ and $\phi$. More realistic particle distribution functions can be treated analogously. Furthermore, we take a constant value of $\xi$ for a dedicated elliptical galaxy. The Heaviside function $H\left[ \lambda_{{\rm mfp}}(p)-(r_{\rm g}(\xi,\eta,\phi)-d_1)\right]$ implies that only high energetic particles with $\lambda_{{\rm mfp}}(p)\geq r_{\rm g}(\xi,\eta,\phi)-d_1$ can leave the galaxy.

Under these conditions we can convert Eq. (72) into

$$% N\propto 4\pi p^2\int_{0}^{2\pi}\int_{-1}^{1}\int_{r_{\rm g... ...^2+(\xi^2-1)(1-\eta^2)\right]}{\rm d}f{\rm d}\eta {\rm d}\phi.$$ (73)

The general solution of this integral is given by

$$% N\propto\frac{16\pi^2cp^2\lambda_{{\rm mfp}}^2(p)~\cot^{-1}... ...})}{3\sqrt{\xi^2-1}}\sum_i\alpha_iT_i(\xi={\rm const.})R_i(p).$$ (74)

From this calculation we recognise that the escaping spectrum we observe depends on the shape of the galaxy given by the choice of the parameter $\xi$ = const. Inserting Eqs. (69) into (74), we get

$$% N\propto \frac{16}{3}\pi^2c\lambda_{{\rm0,~mfp}}^2\sum_i\alpha_iT_i(\xi={\rm const.})R_i(p)p^{6-2q}.$$ (75)

Under the assumption that $R_{{\rm ges}}(p)$ is proportional to p^-4.8, we find for a Kolmogorov turbulence model (q=5/3) a power-law index of -2.13 for the escaping spectrum. This spectral behaviour is harder than the intrinsic spectrum, which has a spectral index of -2.8 (cf. Sect. 4.3).

Such an effect can be easily understood by the high energetic particles leaving the galactical confinement region more frequently than low energetic ones. As a consequence there are two possibilities explaining the observed high energy cosmic ray spectrum above the ``knee'' with spectral index of -3.1 assuming that elliptical galaxies provide a significant amount to the overall high-energy cosmic ray flux. First, the intrinsic spectrum in elliptical galaxies may be steeper then that one we have chosen here. In this context the measurement of the gamma-ray flux reaching the Earth from such an object would give us interesting constraints. Second, if the intrinsic spectrum in elliptical galaxies has nearly the same dependence as the one we measure here on Earth ( $N(p)\propto p^{-2.7}$), the transport of cosmic rays after escape from elliptical galaxies depends on energy. In these two scenarios it is possible to explain the high-energy component with our model.

Note that, due to the previously given arguments, the highest energetic particles cannot be confined within the galactical volume. In our examination the maximum energy of confined charged particles within elliptical galaxies is given by the relation

$$% f_{\rm e}\xi=\lambda_{{\rm mfp}}(p)=\lambda_{{\rm0,~mfp}}p^{2-q}.$$ (76)

For reasonable values of giant elliptical galaxies ( $\lambda_{{\rm mfp}}(p=10~{\rm GeV}~{c}^{-1})=10^{19}$ cm, $f_{\rm e}\xi=5$ $\times$ 10²² cm, cf. Schlickeiser 2002) we get from Eq. (76) a value of about 10²¹ eV for cosmic ray protons. Above this energy value, it is generally impossible that elliptical galaxies confine cosmic rays within their volume.

6 Conclusion

We showed that an analytical treatment of cosmic ray transport in elliptical galaxies based on the diffusion approximation is possible in general. The formal solution we found, combined with appropriate boundary conditions and source functions, can be used to study transport processes in elliptical galaxies. This model is valid for the complete physical parameter space as long as the diffusion approximation holds.

The first test case where we applied our model to is a jet-like injection shape like in M 87. As we can separate our problem into spatial and momentum problems, this test case probes the spatial problem. For this model we found that for very long timescales cosmic rays are distributed throughout the galaxy almost isotropically. We also provided a test case for the momentum problem. Under the assumption of a resulting power-law spectrum matching the observed power law-spectrum on Earth, we could identify the governing eigenfunctions and therefore the governing processes for cosmic ray acceleration and momentum diffusion in elliptical galaxies. It turned out that in elliptical galaxies adiabatic losses are responsible for the high energy cutoff, whereas the slope of the spectrum is given by the ratio of the strength of Fermi I and Fermi II processes. As a result of the basic possibility of explaining this power law-with physical parameters that might be found in M 87, our calculation gives rise to the theory of M 87 as a source of ultra-high energy cosmic rays. In this context the gamma-ray flux from pion decay and the escaping cosmic ray spectrum are essential for a more adequate view of cosmic ray astrophysics. We find that the gamma-ray spectrum with a power law index of -2.69 is a bit harder than the intrinsic cosmic ray spectrum (power law index -2.8). In the context of understanding the sources of high-energy cosmic rays, the flux of charged particles escaping from the confinement region of elliptical galaxies is also very interesting. For the same intrinsic spectrum as above, we find -2.13 as the spectral index for escaping particles. Again the spectrum is harder than the value within elliptical galaxies.

Our model may also serve as a testbed for some more advanced questions:

• Are elliptical galaxies, and especially AGN within these objects, reasonable sources of the ultrahigh energy cosmic rays?

• What is the total flux leaving a single elliptical galaxy with respect to a measured source of cosmic rays in elliptical galaxies?

• What is the total cosmic ray flux we can expect to measure on Earth in the special case of the nearby active elliptical galaxy M 87 or Centaurus A?

Clearly, to answer these questions, a detailed quantitative analysis of the properties of cosmic ray transport within elliptical galaxies is required. Therefore we need to determine from observations the exact physical parameters of the interstellar medium and of the emission processes within the source. With our analytical model at hand, the task of constraining the possible parameter space gets easier. Especially by assuming that UHECR are coming mostly from sources like M 87, we may be able to derive the total cosmic ray content of this type of elliptical galaxies.

Nevertheless this work is considered as a starting point for more sophisticated (numerical) models. Another interesting point is the expansion of this model to cosmic ray electrons. While the spatial description may be derived analogously, we need to add the physical processes involved in the transport of electrons in momentum space. Synchrotron losses and the inverse Compton effect will then play a major role (see, e.g., Casadei & Bindi 2004). The output of this model can be easily tested with the radio data of elliptical galaxies since the electrons provide the main contribution to this radiation.

Acknowledgements

T.H. acknowledges support by Graduiertenkolleg 1147 and F.S. acknowledges support by the Deutsche Forschungsgemeinschaft, DFG project number Sp 1124/1-1. We would like to thank R. Schlickeiser for his useful comments.

Appendix A: General solution of Eq. (25)

The general solution of Eq. (25) can be obtained by consecutive separation of variables. First we use

$$% T(\xi,\eta,\phi,u)=X(\xi,\eta,\phi)U(u),$$ (77)

leading to the following differential equations with -k² as the separation constant:

$$% \frac{\partial U}{\partial u}=-k^2 U$$ (78)

$$% \frac{\partial}{\partial \xi}\left[(\xi^2-1)\frac{\partial ... ...^2)}\frac{\partial^2 X}{\partial \phi^2}+c^2(\xi^2-\eta^2)X=0.$$ (79)

The factor c is given by

$$% c=\frac{fk}{\sqrt{K_0}}\cdot$$ (80)

Equation (78) is easily solved by

$$% U(u)=\exp(-k^2 u),$$ (81)

where we put the integration constant to equality without loss of generality. This result serves an a consistency check to the eigenfunction expansion Eq. (19). The separation constant k will be calculated by use of spatial boundary conditions.

Equation (79) can be solved through another separation of variables. With

$$% X(\xi,\eta,\phi)=A(\xi,\eta)B(\phi),$$ (82)

one gets

$$% -\frac{1}{B}\frac{\partial^2 B}{\partial \phi^2}= \frac{(\x... ...tial A}{\partial \eta}\right]\right\} +(\xi^2-1)(1-\eta^2)c^2.$$ (83)

The separation constant is set to be +m². Therefore we get

$$% \frac{\partial^2 B}{\partial \phi^2}+m^2 B=0,$$ (84)

and

$$% \frac{(\xi^2-1)(1-\eta^2)}{(\xi^2-\eta^2) \times A} \left\{... ...right]\right\}+\left\{(\xi^2-1)(1-\eta^2)c^2-m^2\right\}A = 0.$$ (85)

The solution of Eq. (84) is given by

$$% B(\phi)=\sin(m\phi)+ \cos(m\phi),$$ (86)

where we again used integrating constants equal to 1. The last separation

$$% A(\xi,\eta)=R(\xi)S(\eta)$$ (87)

solves Eq. (85), which can be written by

$$% \frac{\partial}{\partial \xi}\left[(\xi^2-1)\frac{\partial ... ...2 \left(\frac{1}{\xi^2-1}+\frac{1}{1-\eta^2}\right)\right]A=0.$$ (88)

With the separation constant $+\beta_{mn}$ we get

$$% \frac{\partial}{\partial \xi}\left[(\xi^2-1)\frac{\partial ... ...ight]-\left(\beta_{mn}-\xi^2c^2+\frac{m^2}{\xi^2-1}\right)R=0,$$ (89)

and

$$% \frac{\partial}{\partial \eta}\left[(1-\eta^2)\frac{\partia... ...ht]+\left(\beta_{mn}-\eta^2c^2-\frac{m^2}{1-\eta^2}\right)S=0.$$ (90)

It is notable that the $\xi$ and $\eta$ solution functions satisfy similar differential equations. From Abramowitz & Stegun (1972) we find

$$% R_{mn}^{(p)}(c,\xi)= \left[\sum_{r=0,1}^{\infty ^{\prime}} ... ...ime}} i^{r+m-n} \frac{(2m+r)!}{r!}d_r^{mn} Z_{m+r}^{(p)}(c\xi)$$ (91)

as the solution for the $\xi$ dependence of Eq. (25). It is $Z_{n}^{(p)}(c\xi)=\sqrt{\frac{\pi}{2c\xi}}J_{n+\frac{1}{2}}(c\xi)$for the first kind (p=1) and $Z_{n}^{(p)}(c\xi)=\sqrt{\frac{\pi}{2c\xi}}Y_{n+\frac{1}{2}}(c\xi)$for the second kind (p=2). J and Y are Bessel functions of the first and second kind respectively. For our model we can neglect the solution of second kind because of physical motivated spatial boundary conditions (see Appendix B). It is clear that n=m+r. So i^r+m-n in Eq. (91) is equal to 1. The mark $^{\prime}$ connected with the sum means a summation either over even or over odd values of r. The latter ones as well as the weighting factors d_r^mn have also been adjusted to spatial boundary conditions.

Equation (90) has the formal mathematical solution given by

$$% S_{mn}^{(1)}(c,\eta)=\sum_{r=0,1}^{\infty ^{\prime}}\hat{d}_r^{mn}(c)P_{m+r}^m(\eta),$$ (92)

and

$$% S_{mn}^{(2)}(c,\eta)=\sum_{r=0,1}^{\infty ^{\prime}}\hat{d}_r^{mn}(c)Q_{m+r}^m(\eta).$$ (93)

Here $P_{m+r}^{m}(\eta)$ and $Q_{m+r}^m(\eta)$ denote associated Legendre polynomials of first and second kind respectively and $\hat{d}_r^{mn}(c)$ are weighting factors. Again $^{\prime}$ means a summation over even or odd values of r. We can neglect the solution function of second kind with respect to spatial boundary conditions (see the consistency checks).

Under these conditions we can write the formal mathematical solution as

 

$$% T(\xi,\eta,\phi, u) \sum_{m,n} R_{mn}^{(1)}(c,\xi)\times S_{mn}^{(1)}(c,\eta)\times \cos(m\phi)\times \exp(-k^2u)$$ =

$$\sum_{m,n} \left\{\left[\sum_{r=0,1}^{\infty ^{\prime}} \frac{(2m... ...}{2}}(c\xi) \hat{d}_r^{mn}(c)P_{m+r}^{m}(\eta) \cos(m\phi)\exp(-k^2u) \right\}.$$ = (94)

Appendix B: Effect of spatial boundary conditions to the general solution Eq. (26)

As an illustrative example we calculate the solution depending on the following boundary conditions. First we define one for the variable $\xi$. Remembering, that $f\xi$ is equal to the semi-major axis a we assume ``free escape'' boundary conditions in the form

$$% T(\xi=\xi_{\rm c},\eta,\phi,u)=0.$$ (95)

Hence the particle distribution function $T(\xi,\eta,\phi,u)$ is set to be zero at the maximum size $f\xi_{\rm c}$ of the confinement region of cosmic ray particles being a convenient way to model leakage out of the galaxy. This region is assumed to be larger than the (optical) size of the galaxy ( $f\xi_{\rm g}$). Furthermore we assume periodical boundary conditions for $\eta$ and $\phi$:

$$% T(\xi,\eta = \-1,\phi,u)=T(\xi,\eta = +1,\phi,u)$$ (96)

$$% T(\xi,\eta,\phi=2\pi,u)=T(\xi,\eta,\phi=0,u).$$ (97)

Finally we take

$$% T(\xi,\eta,\phi,u=0)=q_0 \Theta(\xi_{\rm g}-\xi)$$ (98)

for a given distance of the foci (2f) as a simple model to explain a constant source function over the whole size of the galaxy. Here we neglect the $\phi$-dependence (taking a $\phi$-dependence into account is straightforward). As a consequence we set m=0 in Eq. (26) leading to n=r (cf. Appendix A).

In order to match Eq. (96) only even solution functions are practical. This is done by $P_{r}^0(\eta)$. Therefore we neglect $Q_{r}^0(\eta)$ having singularities at the points $\eta=1$and $\eta=-1$. From Eq. (95) we recognise that $J_{r+\frac{1}{2}}(c\xi)=0$ if $\xi=\xi_{\rm c}$. Let y_ri be the zeros of $J_{r+\frac{1}{2}}(\xi)$, then we find

$$% c=\frac{y_{ri}}{\xi_{\rm c}},\hspace{0.2cm} k_{ri}=\frac{y_{ri}\sqrt{K_0}}{\xi_{\rm c}f}\cdot$$ (99)

As a result we write

$$% k^2 \equiv k_{ri}^2=\frac{K_0 y_{ri}^2}{\xi_{\rm c}^2 f^2}\cdot$$ (100)

This means that we have to sum over all zeros y_ri aditionally. The last boundary condition Eq. (98) is useful to determine the weighting factors d_ri and $\hat{d}_{ri}$ respectively. From

$$% \sum_{i=1}^{\infty} \left\{ \frac{\sum_{r=0,1}^{\infty^{\pr... ...}}\hat{d_{ri}}P_r^0(\eta)\right\} =q_0 \Theta(\xi_{\rm g}-\xi)$$ (101)

we see that we can cancel the factors d_ri.

Adjusting Eq. (26) to the boundary condition Eq. (96) it becomes clear that only even values of rmatch the required periodicity. Using the orthonormality relation for Legendre polynomials

$$% \int_{-1}^1P_r^0(\eta)P_{r^{\prime}}^0(\eta){\rm d}\eta=\frac{2}{2r+1}\delta_{rr^{^{\prime}}},$$ (102)

the RHS of Eq. (101) evolves to be

$$% q_0 \Theta(\xi_{\rm g}-\xi) \int_{-1}^1 P_{r^{^{\prime}}}^0(\eta) {\rm d}\eta.$$ (103)

Performing the integral it is obvious that only $r^{^{\prime}}=0$gives a non-vanishing value for $\int_{-1}^1 P_{0}^0(\eta) {\rm d}\eta=2$. Consequently we set $r=r^{^{\prime}}=0$. As a result all summations over r disappear so there is only

$$% \sum_{i=1}^{\infty}\hat{d_{i}} \sqrt{\frac{\pi \xi_{\rm c}}... ..._i\frac{\xi}{\xi_{\rm c}}\right) = q_0 \Theta(\xi_{\rm g}-\xi)$$ (104)

left. Here we use the orthonormality relation for Bessel functions, i.e.

$$% \int_0^a J_{\nu}\left( \alpha_{\nu u} \frac{\rho}{a}\right)... ...\rho = \frac{a^2}{2}[J_{\nu +1}^2(\alpha_{\nu u})]\delta_{uw},$$ (105)

to calculate the weighting factors. The Bessel function of second kind diverges at $\xi=0$ and is therefore not appropriate for the general mathematical solution (cf. Appendix A). Finally we find (we substitute $\alpha_i \equiv \sqrt{\frac{\pi\xi_{c}}{2y_i}} \times \hat{d_{i}}$)

$$% E_i(\vec{r}) \equiv T(\xi,u)=\sum_{i=1}^{\infty}\alpha_i \f... ...ht) \exp \left[ -\frac{K_0 y_i^2}{\xi_{\rm c}^2 f^2}u \right ]$$ (106)

with the weighting factors

$$% \alpha_i=\frac{q_0 \int_1^{\xi_{\rm g}}\xi^{3/2}J_{\frac{1}... ...c{\xi_{\rm c}^2}{2}\left[ J_{\frac{3}{2}}(y_i)\right ]^2}\cdot$$ (107)

Again we see that the general solution of Eq. (25) with respect to spatial boundary conditions is indeed of the form of the eigenfunction expansion Eqs. (19) with (20).

References

• Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions (New York: Dover) (In the text)
• Arfken, G. B., & Weber, H. J. 2005, Materials and Manufacturing Processes (In the text)
• Bekki, K., & Shioya, Y. 1998, ApJ, 497, 108 [NASA ADS] [CrossRef] (In the text)
• Biermann, P. L. 1995, Nucl. Phys. B Proc. Suppl., 43, 221 [NASA ADS] [CrossRef] (In the text)
• Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29 [NASA ADS] [CrossRef]
• Casadei, D., & Bindi, V. 2004, ApJ, 612, 262 [NASA ADS] [CrossRef] (In the text)
• Cowsik, R. 1980, ApJ, 241, 1195 [NASA ADS] [CrossRef]
• de Freitas Pacheco, J. A. 1971, A&A, 13, 58 [NASA ADS] (In the text)
• Fermi, E. 1949, Phys. Rev., 75, 1169 [NASA ADS] [CrossRef] (In the text)
• Hasselmann, K., & Wibberenz, G. 1968, Z. Geophys., 34, 353 (In the text)
• Hayakawa, S. 1969, Cosmic ray physics. Nuclear and astrophysical aspects, Interscience Monographs and Texts in Physics and Astronomy (New York: Wiley-Interscience)
• Hoerandel, J. R. 2007 [arXiv:astro-ph/0702370v1] (In the text)
• Jokipii, J. R. 1966, ApJ, 146, 480 [NASA ADS] [CrossRef] (In the text)
• Knapp, G. R. 1999, in Star Formation in Early Type Galaxies, ed. P. Carral, & J. Cepa, ASP Conf. Ser., 163, 119 (In the text)
• Lerche, I., & Schlickeiser, R. 1985, A&A, 151, 408 [NASA ADS] (In the text)
• Lerche, I., & Schlickeiser, R. 1988, Ap&SS, 145, 319 [NASA ADS] [CrossRef] (In the text)
• Mannheim, K. 1993, A&A, 269, 67 [NASA ADS] (In the text)
• Mannheim, K., & Schlickeiser, R. 1994, A&A, 286, 983 [NASA ADS]
• Owens, A. J., & Jokipii, J. R. 1977, ApJ, 215, 677 [NASA ADS] [CrossRef] (In the text)
• Reimer, A., Protheroe, R. J., & Donea, A.-C. 2004, A&A, 419, 89 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Rieger, F. M., Bosch-Ramon, V., & Duffy, P. 2007, Ap&SS, 309, 119 [NASA ADS] [CrossRef]
• Schlickeiser, R. 1983, in International Cosmic Ray Conference, 12, 193 (In the text)
• Schlickeiser, R. 2002, Cosmic Ray Astrophysics/Reinhard Schlickeiser, Astronomy and Astrophysics Library, Physics and Astronomy Online Library (Berlin: Springer) (In the text)
• Schlickeiser, R., Sievers, A., & Thiemann, H. 1987, A&A, 182, 21 [NASA ADS] (In the text)
• Skilling, J. 1975, MNRAS, 172, 557 [NASA ADS] (In the text)
• Strong, A. W., & Moskalenko, I. V. 1998, ApJ, 509, 212 [NASA ADS] [CrossRef] (In the text)
• Sunyaev, R. A., & Titarchuk, L. G. 1980, A&A, 86, 121 [NASA ADS] (In the text)
• Tavecchio, F. 2005, in The Tenth Marcel Grossmann Meeting, On recent developments in theoretical and experimental general relativity, gravitation and relativistic field theories, ed. M. Novello, S. Perez Bergliaffa, & R. Ruffini, 512 (In the text)
• Wang, Y.-M., & Schlickeiser, R. 1987, ApJ, 313, 200 [NASA ADS] [CrossRef] (In the text)
• Weisstein, E. W. 1999, MathWorld - A Wolfram Web Resource, http://mathworld.wolfram.com/ProlateSpheroidalCoordinates.html (In the text)