The gradient of diffuse ray emission in the Galaxy
Key Words.:
cosmic rays – MHD – Gamma rays: observations ISM: supernova remnantsWe show that the wellknown discrepancy, known for about two decades, between the radial dependence of the Galactic cosmic ray nucleon distribution, as inferred most recently from EGRET observations of diffuse rays above 100 MeV, and of the most likely cosmic ray source distribution (supernova remnants, superbubbles, pulsars) can be explained purely by propagation effects. Contrary to previous claims, we demonstrate that this is possible, if the dynamical coupling between the escaping cosmic rays and the thermal plasma is taken into account, and thus a selfconsistent calculation of a Galactic Wind is carried out. Given a dependence of the cosmic ray source distribution on Galactocentric radius , our numerical wind solutions show that the cosmic ray outflow velocity, , also depends both on , as well as on vertical distance , with and denoting the thermal gas and the Alfvén velocities, respectively, at a reference level . The latter is by definition the transition boundary from diffusion to advection dominated cosmic ray transport and is therefore also a function of . In fact, the cosmic ray escape time averaged over particle energies decreases with increasing cosmic ray source strength. Thus an increase in cosmic ray source strength is counteracted by a reduced average cosmic ray residence time in the gas disk. This means that pronounced peaks in the radial distribution of the source strength result in mild radial ray gradients at GeV energies, as it has been observed. The effect might be enhanced by anisotropic diffusion, assuming different radial and vertical diffusion coefficients. In order to better understand the mechanism described, we have calculated analytic solutions of the stationary diffusionadvection equation, including anisotropic diffusion in an axisymmetric geometry, for a given cosmic ray source distribution and a realistic outflow velocity field , as inferred from the selfconsistent numerical Galactic Wind simulations performed simultaneously. At TeV energies the rays from the sources themselves are expected to dominate the observed “diffuse” flux from the disk. Its observation should therefore allow an empirical test of the theory presented.
1 Introduction
Our information on the spatial distribution of cosmic rays (CRs) in the Galaxy stems largely from measurements of nonthermal emission, generated by the energetic charged particles interacting with matter and electromagnetic fields. For ray energies above 100 MeV, the main production process is probably via decay, resulting from nuclear collisions between high energy particles and interstellar matter. Past and recent observations in the GeV range have shown a roughly uniform distribution of diffuse ray emissivity in the Galactic plane, exhibiting only a shallow radial gradient (in a cylindrical Coordinate system). Hence, if rays were to map the spatial CR distribution, we would expect it to be uniform as well. However, associating CR production regions with star formation regions, all possible Galactic CR source distributions are strongly peaked towards a Galactocentric distance at which a ring of molecular gas resides. It is commonly believed that the bulk of the CR nucleons below about is produced in supernova remnants (SNRs)  the majority being core collapse SNRs  most likely by diffusive shock acceleration (Krymsky 1977; Axford et al. 1977; Bell 1978a,b; Blandford & Ostriker 1978). It has been argued that up to some tenths of the hydrodynamic explosion energy might be converted into CR energy (e.g., Berezhko et al. 1994). Since highmass star formation mostly occurs in a spatially nonuniform manner, i.e. in OB associations predominantly located in the spiral arms of late type galaxies, we are confronted with the problem of reproducing a mild radial gradient in the diffuse Galactic ray emission as it has been observed for the first time by the COSB satellite (Strong et al. 1988) and, more recently, with higher angular and energy resolution, higher sensitivity and lower background by the EGRET instrument of the CGRO satellite (Strong & Mattox 1996, Digel et al. 1996). If the SNRs are the sources of the CR nucleon component and if the source distribution is inhomogeneous, this discrepancy must arise during the propagation of CRs from their sources through the interstellar medium. Unlike e.g. the interpretation of radio synchrotron emission, generated by relativistic electrons, ray data open the possibility of studying the nucleonic component of the CRs, in which almost all of the energy is stored (see e.g., Dogiel & Schönfelder 1997). The distribution of the ray emissivity in the Galactic disk therefore bears important information on the origin of CRs and on the conditions of CR propagation in the Galaxy.
2 Gradient of the Ray Emissivity in the Galactic Disk
The first data on the radial distribution (i.e. gradient) of the ray emissivity in the disk were obtained with the SASII satellite (energy range 30 – 200 MeV). The data showed rather strong variations of the emissivity along the Galactic plane which dropped rapidly with radius (see e.g., Stecker & Jones 1977), in rough agreement with the distribution of candidate CR sources in the Galaxy such as SNRs or pulsars. A noticeable discrepancy however emerged from the COSB data (energy range 70 – 5000 MeV), in which the emissivity gradient was found to be rather small, when compared to the SNR distribution, especially at high energies (300 MeV), with a maximum variation by a factor of only 2 (see Strong et al. 1988). An energy dependence of the gradient was mentioned by almost all groups analyzing the COSB data and was usually interpreted as a steeper gradient for CR electrons, producing the soft part of the ray spectrum, compared to nuclei. Detailed analyses of different models of CR propagation based on the gradient value were usually performed with the COSB data because of their much better statistics compared with the SASII data.
The recent measurements with the EGRET instrument onboard CGRO, obtained by different methods at energies 100 – 10000 MeV, showed that the emissivity drops at the edge of the disk. Digel et al. (1996), performed a study of the outer part of the Galaxy towards molecular complexes, the Cepheus and Polaris Flares in the local arm and the larger molecular complex in the Perseus arm. Since the total masses of these complexes are known, it was possible to infer the CR density from the measured ray fluxes in the direction of these objects. It was discovered that the apparent emissivity decreases by a factor 1.7, which is somewhat smaller than that for COSB. The analysis of Strong & Mattox (1996), based on a model of the average gas density distribution in the disk, shows a smaller intensity gradient which does not differ significantly from that of COSB in the Galactic disk.
One of the important conclusions which follows from all these data, and which we will discuss below, is that the distribution of the emissivity in the disk in the GeV range is rather uniform compared with the most probable distribution of CR sources.
3 Basic Ideas for a New Approach of CR Propagation in the Galaxy
A natural explanation of a uniform CR distribution would be effective radial mixing due to the diffusion of CRs produced in different parts of the Galactic disk. It is then straightforward to infer the mixing volume, which usually includes the Galactic disk plus a large Galactic halo; the details of such a model were discussed in the book of Ginzburg and Syrovatskii (1964). This 3dimensional cylindrical model, which includes diffusiveadvective transport, a free escape of cosmic rays from the halo boundary into intergalactic space, and the observed supernova shells as sources of cosmic rays in the Galactic disk, explained very well characteristics of the Galactic radio and ray emission as well as the data on cosmic ray spectra and their chemical composition including stable and radioactive secondary nuclei, intensities of positrons and antiprotons etc. (a summary of this analysis can be found in Berezinskii et al. 1990).
In general the propagation equation for CRs is described by
(1) 
where is the diffusion tensor and is the advection velocity, which in general are functions of the coordinates, and is the term which describes particle energy losses including adiabatic losses if the velocities vary spatially. The right hand side of the equation describes the spatial and energy dependence of the source density, derived from radio observations of the SNR distribution and the injection spectrum of CRs.
The main conclusion was that almost all observations can be reasonably explained if the halo extension is about several kpc, the injection spectrum of electrons and protons is a powerlaw ( with ), and the radially uniform velocity of advection does not exceed the value of 20 km/s. Compton scattering of relativistic electrons was found to play a significant rôle in the halo ray emission. However, this model failed to explain a smooth emissivity distribution of rays in the Galactic disk (see for details Dogiel & Uryson 1988). Even in the case of a very extended halo with a radius larger than 10 kpc the derived emissivity gradient (calculated for the observed SN distribution) was larger than observed. Only for a hypothetical uniform distribution of the sources in the disk the calculations can reproduce the data.
This model was further developed by Bloemen et al. (1991, 1993). Extensive investigations of the 3dimensional diffusionadvection transport equation for nucleons, lowenergy electrons and the rays showed that even in the most favorable case of an extended halo with a vertical height as large as 20 kpc, their model, although reproducing the COS B data marginally, is not able to remove the signature of the observationally inferred SNR source distribution, i.e. the distribution of the calculated emissivity was still steeper than permitted by the data. The vertical gradient of the advection velocity derived from the fluxes of stable and radioactive nuclei near Earth had to be smaller than 15 km/s/kpc. The other unexpected conclusion was that the halo extension obtained from the nuclear data was significantly less than estimated from the ray data (see also Webber et al. 1992).
Recently, Strong and collaborators in a series of papers (see Strong et al. 2000, and references therein) have developed a numerical method for this model and made a new attempt to analyze the ray emission and the CR data, based on the latest data from COMPTEL, EGRET and OSSE. They limited the analysis to radially uniform diffusionadvection and to diffusionreacceleration transport models, and concluded that no such diffusionadvection model can adequately describe the data, in particular the B/C ratio and the energy dependence. These authors claimed, however, that by including reacceleration one can account for all the observational data. A peculiar consequence of their analysis was that they did not use the observed supernova source distribution (see e.g., Kodeira 1974, Leahy & Xinji 1989, Case & Bhattacharya 1996, 1998) as the input for Eq. (1), but rather derived it from the ray data to reproduce the observed spatial variations of the emissivity in the disk. The result was a source distribution that is flat in the radial direction.
Thus we see that the conventional model of unifrom diffusionadvection has serious problems in spite of its evident achievements. One solution is to assume that some of the observational data are not significant like the SN distribution derived from the radio observations (Strong et al. 2000). The alternative is to conclude that it is time to abandon the standard model, which is what we do in this paper. We shall demonstrate, that strong radial source gradients will be removed by a strong advection velocity in the halo (due to a Galactic wind driven by the CRs themselves, see below) that varies with radius and height . In addition anisotropic diffusion with different diffusion coefficients and in the disk and the halo, respectively, might also play a rôle. It should be emphasized that a radially varying advection velocity occurs naturally in spiral galaxies, even for a uniform source distribution, because the gravitational potential increases towards the centre, thus inducing stronger velocity gradients in this direction (see Fig. 2), as has been shown by Breitschwerdt et al. (1991).
The existence of strong advective CR transport in the Galactic halo has been shown on dynamical grounds in a number of papers in the past (Ipavich 1975; Breitschwerdt et al. 1987, 1991, 1993; Fichtner et al. 1991; Zirakashvili et al. 1996; Ptuskin et al. 1997). The key element of halo transport theory is that CRs, which by observations are known to escape from the Galaxy, resonantly generate waves by the socalled streaming instability (Kulsrud & Pearce 1969) leading to strong scattering of CRs. Therefore, even in the case of strong nonlinear wave damping, advection is at least as important a CR transport mechanism out to large distances in the halo as diffusion, provided that the level of MHDturbulence is high enough for coupling between CRs and MHD waves (Dogiel et al. 1994, Ptuskin et al. 1997). There is also growing indirect observational evidence of outflows from the interpretation of soft Xray data of galactic halos in edgeon galaxies like NGC 4631 (Wang et al. 1995), and also of the soft Xray background in our own Galaxy (Breitschwerdt & Schmutzler 1994, 1999). Furthermore, the near constancy of the spectral index of nonthermal radio continuum emission over large distances along the minor axis in the halo of edgeon galaxies is most naturally explained by an advective transport velocity of relativistic electrons (along with the nucleons) that is ever increasing with distance from the Galactic plane (Breitschwerdt 1994).
In the disk of spiral galaxies, the regular magnetic field is following roughly the spiral arms and is therefore mostly parallel to the disk, with noticeable deviations in some regions where outflow is expected. Here, also a regular vertical component seems to be present (eg. Hummel et al. 1988), which has been detected in a number of galaxies like NGC 4631, NGC 5775 (Tüllman et al. 2000), and NGC 4217. Multiwavelength observations of the galaxy NGC 253 show a local correlation between nonthermal radio continuum, H and Xray emission near the diskhalo interface in offnuclear regions (Dettmar 1992; M. Ehle, private communication). This also spatially coincides with enhanced star formation activity in the disk as can be seen from FIR data.
Since the disk is not fully ionized in contrast to the halo and since waves are efficiently dissipated there by ionneutral damping, the most important contribution to the random field in the disk is by turbulent mass motions, induced by supernova explosions and other stellar mass loss activity. Thus the wave spectrum will be very different from the one in the halo, where selfexcited waves subject to nonlinear wave damping (Dogiel et al. 1994, Ptuskin et al. 1997), satisfying the gyroresonance condition, dominate. Consequently, the averaged diffusion coefficient will be different in the Galactic disk and the halo (see Sect. 4). Therefore, we believe that anisotropic diffusion, together with radially varying advection, is the most general and most probable mode of CR transport to occur.
We investigate the CR transport processes of diffusion and advection and discuss the possibility of flattening radial CR source gradients of a given SNR distribution (despite observational uncertainties) by particle propagation. In our view it is using the “natural” boundary condition when calculating the response of CR transport to a given CR source distribution (whatever its observational limitations may be at the present time). This is in contrast to adjusting the source distribution a posteriori. To implement the transport processes properly we shall allow for radially varying advection, anisotropic diffusion, (different values of the diffusion coefficient parallel and perpendicular to the disk) and the appropriate boundary conditions to employ a cosmic ray distribution and then to calculate from these the resulting spatial distribution of CRs.
In Sect. 4 we discuss the question of anisotropic CR diffusion and in Sect. 5 we introduce a simple model that describes the observed CR source distribution. In the following Sect. 6 we show heuristically (and in Sect. 8 also analytically) how a radially varying CR source distribution induces variations in the CR energy density, which in turn leads to a radial variation of the diffusionadvection boundary, , and the outflow velocity, respectively, and thus to a tendency to flatten the radial CR distribution and hence the ray emissivity gradient. In Sect. 7 we demonstrate by numerical calculations how such a source distribution will naturally generate a radially varying outflow velocity. In Sect. 8 we discuss in detail models of different complexity for advectivediffusive transport with radially varying outflow velocity, and show analytically how in each case for a given CR source distribution its radial signature on the energetic particle and ray distribution is reduced. The most advanced of these models include a dependence of galactic wind velocity on the CR source strength. In addition a full analytic solution of the 3dimensional diffusionadvection equation in axisymmetry for a given realistic velocity dependence , parallel and perpendicular to the Galactic disk similar to the one derived in Sect. 7 is calculated. In Sect. 9 the source contributions to the “diffuse” ray background at TeV energies are taken into account and shown to be highly significant. This should allow an empirical test of our theoretical picture. In Sect. 10 we discuss and summarize our results. A number of detailed calculations can be found in the Appendices A, B and C.
4 Anisotropic Diffusion Coefficient
In most models of diffusive CR propagation, the diffusion tensor is approximated by a scalar quantity , representing spatially uniform transport, , where describes the energy dependence. However, there are good reasons why the diffusion coefficient may be anisotropic.
The propagation of CRs in the interstellar medium is mainly determined by their interaction with electrical and magnetic fields. CRs interact strongly with fluctuations of the magnetic field (MHDwaves) and are scattered by them in pitch angle. In the simplest case the total magnetic field consists of two components and can be written as
(2) 
where is the average large scale magnetic field and is the field of small scale fluctuations.
Effective scattering of particles by these fluctuations occurs when the interaction is resonant, i.e. when the scale of the fluctuations in the magnetic field is of the order of the particle gyroradius. This leads to a stochastic motion of particles through space; the associated diffusion coefficient along the magnetic field is of the order of
(3) 
where is the particle velocity. The value is the transport length of the particle along the magnetic field lines, and is the rate of resonant particle scattering. The main assumption of this description is that the fluctuation amplitude is much less than the average magnetic field
(4) 
The displacement of a charged particle perpendicular to the regular magnetic field due to scattering is therefore small:
(5) 
i.e.,
(6) 
where is the particle gyroradius, and .
The situation is more complex if there exists also a large scale random magnetic field whose scale is much larger than the particle gyroradius:
(7) 
In this case the diffusion transverse to the magnetic field occurs due to the divergence of the particle trajectories, which is much faster than the resonant diffusion from Eq. 5.
The procedure to derive the transport equation for CRs in this case was described , e.g. by Toptygin (1985), who showed that the maximum value of is:
(8) 
Thus, diffusion in the Galaxy is quasiisotropic only in the case .
A more precise analysis (see Berezinsky et al. 1990, chap. 9) shows, however, that in the interstellar medium the correlation between the components of the diffusion tensor leads to an effective perpendicular diffusion coefficient:
(9) 
Observations show that in the Galactic disk the large scale magnetic field is typically parallel to the disk, whereas in the halo also a significant perpendicular component may exist, as it has been observed in the case of the edgeon Sc galaxy NGC 4631. Therefore one can expect in the disk, and in the halo^{1}^{1}1Note that in the disk, , and in the halo, , are both parallel to the large scale regular magnetic field direction., with and being the components of the diffusion tensor perpendicular and parallel to the Galactic plane. Thus, in general, CRs spend a considerable amount of their transport time in the halo, subject to anisotropic diffusion.
5 A simple model for the Galactic CR source distribution
The inference on the spatial distribution of CR sources from direct observations is plagued by a number of problems. From energy requirements for the bulk of CRs below eV, it is known that the only nonhypothetical Galactic candidates for the sources are SNRs and pulsars, with global energetic requirements favouring the former. SNRs can be best studied in the radio continuum and in soft Xrays, but as low surface brightness objects larger and older remnants are systematically missed. Since samples are usually flux limited, the more distant objects will be lost as well. Pulsars on the other hand should be present in the Galaxy in large numbers. However they are only detectable if their narrow beams happen to cross the line of sight of the observer or, in Xrays, as isolated old neutron stars; so far only four candidates are known from the ROSAT AllSky Survey (Neuhäuser & Trümper 1999). Therefore selection effects will bias samples heavily in both cases.
Based on the observations of SNRs by Kodeira (1974) and pulsars by Seiradakis (1976), Stecker & Jones (1977) have given a simple radial distribution of the form
(10) 
where, is a SNR surface density in , with , and and are matched to the observations. Here we adopt and . Since both SNRs and pulsars are associated with Population I objects, it is reassuring that and are very similar for both classes of objects.
To determine the proportionality constant for the distribution given in Eq. (10), we use the results obtained by Leahy & Xinji (1989) who used the catalogue of Li (1985) derived from radio observations. Leahy & Xinji (1989) considered shelltype SNRs and applied empirical correction factors due to the incompleteness of the flux limited sample. The spatial distribution thus obtained shows a peak at 46 kpc from the Galactic center (see Fig. 1), assuming an overall rotational symmetry. A more systematic study has recently been undertaken by Case & Bhattacharya (1998), who have made a careful analysis of an enhanced Galactic SNR sample, using an improved relation. These authors find a peak in the distribution at around 5 kpc^{2}^{2}2The difference is due to the fact that with the older value of kpc in Eq. (10) the peak of the distribution appears at kpc and a scale length of kpc, thus confirming the gross features of the older work.
The number of SNRs in an annulus of width is given by
(11)  
and thus the total number of SNRs in the disk of radius kpc is
(13) 
for , and . Numerical integration of Eq. (5) yields .
On the other hand, the total number of SNRs that should be observable at present is their production rate times their average life time before they merge with the hot intercloud medium and lose their individual appearance, or observationally escape the detection limit. The rate at which SNe occur in the disk is per century for randomly exploding stars (van den Bergh 1990, 1991) and 0.45 per century for explosions occurring in OB associations (Evans et al. 1989), giving a total SN rate in the Milky Way disk of per century. The total number of SNRs in the disk is then given by
(14) 
for years we reproduce the value given by Leahy & Xinji (1989), who estimate for Galactic SNRs with a limiting surface brightness of . From Eq. (13) we obtain: .
Dragicevich et al. (1999) have analyzed also the radial distribution of SNe in a sample of 218 external galaxies of different Hubble types, corrected for inclination angle. The data were sorted into radial bins and the numbers converted into SNe surface densities. They showed that the radial SN surface distribution can be well fitted by an exponential radial decrease of the form
(15) 
with . Taking a subsample of 36 SbSbc galaxies, they were able to derive a value of kpc for the exponential scale length of the Milky Way as a representative SbSbc galaxy. The number of SNe derived is consistent with historical records of Galactic SNe within 4 kpc from the Sun. In Fig. 1 it is shown that both the exponential decrease and, interestingly, but rather fortuitously, also the absolute numbers in the SN surface density agree fairly well for kpc with the earlier derived relation (10) for the SNR surface density.
Therefore the value calculated for the proportionality constant, , provides a reliable estimate for the number of Galactic SNRs in a circular ring at a distance .
6 Flattening of the radial ray emissivity gradient
In the following we will argue that the CR propagation picture needs to be changed for two physical reasons. Firstly, as we have already pointed out, diffusive transport in the disk and in the halo is not the same due to a different magnetic field geometry and to different wave excitation and damping mechanisms. Thus anisotropic diffusion seems to be rather the rule than the exception. Secondly, advection of CRs is the dominant transport mode above a certain transition boundary located at (), with varying with Galactic radius.
More specifically, we consider SNRs as the primary sources for the Galactic CRs. Each remnant results from an energy deposition of , which is converted into relativistic particles with a certain efficiency (Drury et al. 1989; Berezkho & Völk 2000); the individual value of is poorly known from theory due to uncertainties in the injection process and depends e.g., on the value and orientation of the circumstellar magnetic field. SNe exploding inside a superbubble, i.e. a hot tenuous medium, initially generate low Mach number shocks, which are less efficient in accelerating particles. After a time of the order of the sound crossing time, however, the shock impinges on the much colder and denser surrounding shell and becomes progressively stronger thereby accelerating particles more efficiently; this should to lowest order compensate for the initially decreased efficiency in the hot ambient medium. We would expect that due to the continuous energy input by successive SN explosions also a long lived shock would be able to accelerate particles to energies in excess of eV, albeit adiabatic energy losses would become more and more severe with time. Since the diffusion coefficient of CRs increases with energy, advective transport of particles significantly above 1 GeV will be the dominant mode of transport only at larger distances from the Galactic plane. Therefore these particles will quickly fill an extended halo and not generate many rays in the disk via decay.
Bykov & Fleishman (1992) have argued that successive explosions inside a bubble can generate strong turbulence, which should transform a significant amount of the total free energy to cosmic rays. However, at the same time the injection rate at the shocks may be reduced as a result of shock modification due to previously generated CR particles. With the details of the acceleration mechanism in superbubbles being still debatable, we believe that to lowest order there is no difference in the overall energy transfer from thermal plasma to CRs, if a SN explodes inside a superbubble or just forms a single remnant. Thus the energy production rate of CRs should be roughly proportional to the number of SN explosions, regardless whether they occur in the general ISM or inside a superbubble. The numerical difference in the derived CR energy density (and CR pressure) (cf. Eqs. (18) and (17)), however, between treating particle acceleration in superbubbles as equal to single remnants, and disregarding acceleration in superbubbles altogether, is small. According to the previous Sect. it amounts to a factor , and is therefore well below the uncertainty in the acceleration efficiency . In the following, we tend to be conservative and retain a low value of .
Using the results from Sect. 5 we can estimate the local (in radius) number of SNRs within a circular ring of the Galaxy with a width of, say, 2 kpc. Relating the SN rate directly to the number of observable SNRs by Eq. (14) and writing (cf. Eq. (11))
(16) 
we can now define the SNR rate per unit radial confinement volume, , by
(17)  
The local CR energy density, in the disk at , can then be estimated to be
(18) 
with , and being the local diffusive CR storage volume and time, respectively; the factor accounts for one halfspace, e.g. . The CR pressure in the disk is denoted by and is the specific heat ratio for the CRs. Working in the framework of an average transition boundary, , which divides space below and above in purely diffusion and advection regions, respectively, the probability for particles at to return and hence contribute to the CR pressure in the disk becomes exponentially small. We have therefore used as the vertical extension for the CR storage volume. As we can see from Eq. (17), , depends on the local and total number of sources and the global SN rate, but not on the SNR life time. The vertical extent of the diffusionadvection boundary is determined by the balance of diffusive and advective CR flux at this vertical distance, i.e. to lowest order by
(19) 
and the local CR escape time is
(20) 
Thus the radial CR pressure variation is given by
(21)  
It is easy to show that Eq. (21) follows directly from the kinetic equation for CRs, as can be seen from Eq. (47).
The diffuse ray intensity resulting from decay photons is
(22) 
with and being the lineofsight averaged CR and hydrogen number density, respectively, and is the line of sight. From Eq. (21) we have
(23) 
disregarding the changes in the CR distribution function due to energydependent diffusion. is an observationally fixed quantity, for which we have to use its line of sight averaged value (see Eq. (22)) if we are interested in quantitative calculations of the ray intensity; here we are primarily concerned with the radial variations of .
In order to obtain for a weak ray gradient, we should assume that
(24) 
if ; if , we have and therefore should require
(25) 
7 Radially Dependent Cosmic Ray Driven Outflows
As we have deduced in the last Sect., the CR pressure in the disk is a radially dependent quantity, and therefore we expect the outflow velocity, , and the mass loss rate, , to be also radially dependent. In an earlier paper (Breitschwerdt et al. 1991) we have shown that such a behaviour already exists as a consequence of the radial dependence of the gravitational potential. The net result was a monotonic decrease (increase) of terminal velocity (mass loss rate) with increasing Galactic radius for a radially constant mass density, . Now we have superposed the radial variation of the CR source density and investigate in the following how this changes the outflow.
However, Eq. (21) is an implicit equation, since depends on , which in turn depends on the energy density available in CRs, i.e. and thus itself, to drive the outflow. To that end we have performed selfconsistent galactic wind calculations of the fully nonlinear equations, in which for a given gravitational potential of the Milky Way, and a relativistic CR gas (), together with a spatially averaged mass density , an average thermal pressure , an averaged halo magnetic field , and a small average level of wave amplitude , the advectiondiffusion boundary , and are calculated selfconsistently, using , the value derived from the numerical integration of Eq. (5) (see Sect. 5). The form of the potential (including, disk, bulge and dark matter halo) and the opening of the flux tube due to geometrical divergence are the same as used by Breitschwerdt et al. (1991).
We have used a radial grid for our calculations in the Galactic radius interval from kpc, with spacings of kpc. The net radial dependence of the outflow velocity we obtained from our numerical calculations is shown in Fig. 2, and its  and dependence in Fig. 7. The outflow velocity in Fig. 2 (solid line) for a fixed Galactic frame of reference at the radially varying transition boundary , shows a similar convex shape as the source distribution (short dashed curve), although much less pronounced. This is due to the fact, that the decrease of the gravitational potential, as can be seen from the escape speed, , with for bulge, disk and halo contributions, respectively, (dotted curve in Fig. 2), keeps the outflow velocity from dropping off too rapidly with increasing Galactocentric distance. Therefore, even a uniform SNR distribution would enforce a radially varying outflow velocity, as has been shown in Fig. 7 of Breitschwerdt et al. (1991). The radial variation of the outflow velocity at a distance of kpc from the disk, , is shown by the long dashed curve, which both shows the transition from the flow velocity at reference level to the escape speed, and also a roughly dependence, which we will use in our analytical treatment later (cf. Eq. (98)).
It can be seen from comparison of and , as obtained from the fully nonlinear calculations, that the simple ansatz of Eq. (19) is indeed fulfilled. The functional dependence of with radius is straightforward to understand. Close to the Galactic centre, the gravitational pull is strongest, as can be seen from in Fig. 2; since we chose all other quantities being the same across the disk (constant density, thermal pressure, magnetic field strength), the outflow velocity, and hence mass loss rate, are smallest here. Eq. (19) then tells us, that must be largest. At the outer parts of the Galaxy, the gravitational field becomes weaker, but also the source distribution, and hence the CR pressure, decrease exponentially, and so the outflow velocity (cf. Fig. 2) also decreases, and must increase again (see Fig. 3). It is noteworthy that the maximum of and the minimum of at kpc do not coincide with the maximum of and , respectively, at kpc. This must be a consequence of the interplay between the gravitational field and the source distribution in the fully nonlinear equations (cf. Breitschwerdt et al. 1991).
Finally we mention that the choice of constant boundary conditions across the Galactic disk is a conservative one. In reality we should also enhance the thermal temperature and pressure in regions of higher supernova activity. The net effect would be a more pronounced peak in outflow velocity and a deeper minimum in , respectively, and therefore an even “better” quantitative proportionality between and according to Eq. (24).
8 Models of advective and anisotropic diffusive CR transport
In the following we discuss several advection models (some including also anisotropic spatial diffusion), in which we examine in detail the ideas presented in Sect. 6. We show that for a given radial CR source distribution function , we can find a function , for which the effect of the radially varying Galactic wind velocity leads to an almost uniform CR distribution in the disk. We refer to relations as “compensation equations”. It is thus possible to explain the observational data by pure CR propagation effects.
To fix ideas, we start out with a functional relation for the Galactic CR source distribution of the form
(27) 
neglecting the thickness of the Galactic disk; here denotes the Dirac delta function. Then for a given source distribution function, , we try to find a wind structure of the form
(28) 
which makes the CR distribution almost uniform. In the following we describe models of increasing complexity with respect to geometry and outflow velocity.
8.1 Onedimensional model with a constant galactic wind velocity
We start from the simplest onedimensional CR transport model, which is an extreme case of anisotropic diffusion, since only propagation in the direction is allowed. In this case the galactocentric radius is simply a parameter in the model.
The onedimensional diffusionadvection equation for CR nucleons can be written as
(29) 
where for simplicity we take the advection velocity to be constant and directed away from the Galactic plane
(30) 
Here is the Heavyside step function, and the sign ”” corresponds to the velocity above, and the sign ”” to the one below the Galactic plane.
The boundary conditions for CRs are determined either by free escape into intergalactic space, if the density of electromagnetic fluctuations generated by the CR flux decreases fast enough (see Dogiel et al. 1994), or by CR advection in a galactic wind to infinity (Breitschwerdt et al. 1991; Ptuskin et al. 1997), if the level of fluctuations is high enough. Which of these cases is relevant for the Galaxy, is the subject of a separate investigation. Fortunately, CR spectra and densities in the disk are independent of the boundary conditions far away from the Galactic plane, if there is an outer region of advective transport. In the case we discuss here, it is assumed that the CR propagation region can be formally divided into a diffusion halo wrapping around the Galactic plane and an adjacent advection region, reaching out to intergalactic space. These two regions are separated by a boundary surface, , at which both the CR density and the flux have to be continuous.
The sources are concentrated in the disk and are supposed to emit a powerlaw spectrum of particles
(31) 
where is the thickness of the Galactic disk and the coefficient determines the CR production by SNRs in the disk.
The location of the transition boundary for a constant diffusion coefficient and advection velocity is determined in this approximation by
(32) 
Then for the transport equation reads
(33) 
and for
(34) 
Since adiabatic energy losses do not change the particle spectral index, . The solution of Eq. (33) is
(35) 
and the CR density is constant in the advection region,
(36) 
The constants can be determined from the boundary conditions at
(37) 
and at by continuity of the diffusive and advective CR fluxes
(38) 
Since at also , we obtain
(39) 
and therefore
(40) 
We see that at
(41) 
and the “compensation” equation must have the form
(42) 
for the CR particle density to be uniform, in other words, the radial dependence of the outflow velocity must be the same as that of the source distribution.
From Eq. (40) we can also estimate the total pressure, , of CRs in the disk, which is
(43)  
(44) 
The physicals units of are , and equals the CR power in the Galactic disk, hence (neglecting constants of order unity)
(45) 
On the other hand the normalized source distribution (cf. Eq. (16)) is equivalent to
(46) 
where and are the area and volume of the Galactic disk, respectively, and identifying with , Eq. (44) becomes
(47)  
which corresponds to Eq. (21), thus verifying the validity of the intuitively derived relation between CR pressure, source distribution and Galactic SN rate in a Galactic wind flow.
8.2 Onedimensional model with dependence of the wind
Based on the conclusions of the previous Sections, we can generalize the onedimensional solution obtained by Bloemen et al. (1993), taking into account radial variations of the sources and the wind velocity. Let us suppose that the advection velocity has the form
(48) 
the solution of Eq. (29) is then given by
(49) 
where is a constant; note that is a typical time constant (efolding time) of the flow.
We find that a vanishing radial dependence of the CR distribution, , can be obtained if
(50) 
In these simple models the size of the halo can be derived only locally, because global or largescale characteristics, like e.g. the CR distribution along the Galactic plane, are independent of the CR propagation itself.
These simple analytical solutions of the onedimensional transport equation show that the main effect, which formally leads to the ”compensation”, is a curved transition boundary between diffusion and advection regions, i.e.
(51) 
note that such a behaviour is indeed apparent from Figs. 2 and 3. The main point here is that the compensation effect is a consequence of the dependence of the CR lifetime, .
In the framework of the pure onedimensional model it is unimportant how far from the Galactic plane the boundary is, but as we shall see, this figure will be essential for the threedimensional case.
8.3 Threedimensional axisymmetric model with radial velocity dependence
To demonstrate the effect of radial changes of the boundary between regions of diffusive and advective propagation of CRs, we investigate the diffusionadvection equation with a given radial dependence of the wind velocity and the source density in a more realistic cylindrical geometry (axial symmetry, i.e. ) with the velocity varying as
(52) 
and the source distribution as
(53) 
Then the equation for the CR distribution function we have to solve reads
(54) 
where we have stressed the anisotropic diffusion by two different constant diffusion coefficients in the radial and perpendicular direction, and , respectively. Note that has the same dimensions as the diffusion coefficients, and that is the particle distribution function in phase space. The boundary conditions at are implicitly included in Eq. (54). Here we treat a special case of radial variation of the boundary :
(55) 
corresponding to the velocity dependence given in Eq. (52). Comparison with Figs. 2 and 3 shows that such a behaviour may indeed represent the outer Galaxy, where the outflow velocity decreases with increasing radius and the distance of the boundary from the Galactic disk increases accordingly.
Eq. (54) can be solved analytically if we restrict ourselves to selfsimilar solutions, in which the distribution function does not vary with and independently. The price we have to pay is that there is no unique solution that covers both and . Instead a similarity solution is found for and , and one for , which have to match in overlapping regions. We start out with the latter one, applying a transformation of independent variables in the form
(56) 
and introducing the function
(57) 
Equation (54) in terms of and then reads:
(58)  
where the coefficients are given by
(59)  
This equation can be cast into a selfsimilar form for two sets of the parameters and : , , and , .
We can derive an analytical solution of this equation in the form (see Bakhareva & Smirnova 1980)
(60) 
where is a constant determined by the CR source power.
The unknown coefficients , and are determined from the system of algebraic equations which one can obtain from Eq. (58) by equating the coefficients for the different powers of to zero. The resulting system of algebraic equations is given by
(61)  
(62)  
(63)  
(64) 
We note that this system for the unknowns , and is overdetermined; therefore we have to check for consistency.
From the system of equations (61)(64) it can be deduced that
(65)  
(66)  
(67)  
(68) 
It can be verified that these solutions are all consistent. However, implies and , and hence does not yield any physical solution. Taking gives . Therefore the solution is given by
(69) 
where the constant is determined from the conservation of particle flux at .
The limit corresponds to so that
(70)  
because at large distances from the disk (), the particle transport is completely determined by advective transport, and the number of particles ejected by the sources is conserved.
The full solution of Eq. (54) reads
(72)  
From the solution Eq. (72) we can obtain the particle density behaviour at large and , which are:

for :
(73) 
and for :
(74)
From Eq. (54) it is clear that we are unable to describe the particle distribution near with the selfsimilar variable (Eq. (56)) and the function from Eq. (57), since it is impossible to satisfy the boundary conditions at .
In order to study the radial behaviour of the distribution function near a different selfsimilar ansatz is used, which can be extended down to the spatial region near the sources. The similarity variable has the form
(75) 
and the transformed distribution function reads
(76) 
Then Eq. (54) becomes
(77) 
and if
(78) 
the RHS of Eq. (77) can also be cast into selfsimilar form
(79) 
The boundary condition at immediately follows from integration of Eq. (77) for small , in which case the equation becomes
(80) 
We integrate Eq. (80) along a box, encircling the Galactic plane with , from to , and subsequently let . The formal integration near gives the boundary condition at ,
(81) 
The function is obtained by solving the homogeneous equation
(82) 
The solution of Eq. (82) is found to be
(83)  
where
(84) 
Our goal is to study the dependence of the particle density near , in other words we have to determine the value of .
From the asymptotic expansion of hypergeometrical functions at , we find that the solution Eq. (83) depends on as
(85) 
This also follows directly from Eq. (77), which for large has the form
(86) 
It is easy to verify that behaves like a power law in , i.e.
(87) 
and the characteristic equation
(88) 
has just one degenerate solution:
(89) 
We have already shown that at large the distribution function has the form (cf. Eq. (73)). On the other hand it follows from Eq. (76) for
(90) 
and therefore , thus matching the two solutions obtained for different regions in space. It is clear that by the choice of the similarity variable (Eq. (75)), it is impossible to describe the distribution function near the Galactic centre (), but everywhere else in the disk.
Since close to the sources, , , we derive from Eq. (76)
(91) 
whereas the source distribution is required to have a radial exponent as follows from Eq. (78), and therefore
(92) 
Thus radially dependent advection flattens the distribution of CRs in the disk as compared to the source distribution
and the “compensation” condition, which in this case means a partial compensation since rather than , is realized in this case if
(93) 