text:
# Relic Radio ”Bubbles” and Cluster Cooling Flows

## 1 Introduction

The interaction of extended radio sources in clusters of galaxies with the intracluster medium (ICM) has been inferred since the first observations of ”head-tail” radio sources over two decades ago. However, direct observation of this interaction required imaging of the ICM itself, and early hints of this came from the Einstein Observatory x-ray observations of Cygnus A (Harris et al. 1994). Unmistakable evidence for interaction came from the $`ROSAT`$ high resolution images of NGC 1275 in the Perseus cluster (Böhringer et al. 1993), and more recent high resolution x-ray observations with the Chandra Observatory have provided clear and dramatic evidence for the interaction of radio sources with the hot intracluster medium (e.g., Fabian et al. 2000). These images, when coupled with high resolution radio observations, clearly show the extended radio source lobes displacing the x-ray emitting gas and forming a cavity in the ICM. In general these cavities show a rather regular and spheroidal morphology, suggesting the presence of a ”relaxed” structure that is in approximate equilibrium with its surroundings. There is no evidence for a high temperature layer surrounding the cavities that would be produced by a shock front; in some cases the region of bright x-ray emission surrounding cavity has a lower temperature than the ambient ICM (Nulsen et al. 2002). In addition, the x-ray observations of the cavity in A2597 show that the pressure in the limb brightened region surrounding the cavity is not higher than the ambient ICM pressure, but that it is comparable to or somewhat lower than the ambient value (McNamara 2002a, McNamara et al. 2001). The densities in the limb brightened region are, however, higher than that in the undisturbed ICM. This is not what would be expected if the surface of the cavity were expanding supersonically into the ICM, since the resulting shock in front of the cavity surface would yield a pressure significantly higher than ambient: $`p_1/p_0=(2\gamma M_0^2\gamma +1)/(\gamma +1)`$, where $`M_0>1`$ is the Mach number of the shock relative to the ambient gas, the subscripts $`0`$ and $`1`$ refer to pre- and post- shock conditions, and $`\gamma `$ is the usual ratio of specific heats. The above relation is valid for arbitrary shock strengths (e.g., De Young 2002). Moreover, the geometry of the cavities, and especially the ”relic” cavities discussed below, is not suggestive of the supersonic expansion that would be expected for a radio lobe that is still being supplied by collimated jets emanating from the nucleus of the parent galaxy. As is well known, these energetic lobes have a much more elongated geometry, such as that seen in the ”classical” FR-II radio sources. The nearly spherical geometry of the x-ray and radio morphologies under discussion here is consistent with a radio lobe that has come into rough pressure equilibrium with the ambient medium. Again, this is particularly true of the relic cavities discussed below, where there is no clear connection seen in radio emission between the relic and the active nucleus. In addition, the inner radio source in A2597 shows a complex and contorted geometry that suggests that strong deceleration and ”frustration” may be taking place in the inner regions where the extended source is still being fed by outflow from the nucleus. This again would lead to prompt deceleration of the radio source relative to the ICM, and when the outflow from the nucleus ceases, such objects would rapidly come into pressure equilibrium with the ambient gas. All of this is suggestive of rather slow, subsonic inflation rates for the cavities, and if they are in pressure equilibrium with the ambient ICM, then their internal energy densities exceed the minimum equipartition energy densities by roughly an order of magnitude (McNamara et al. 2001). It is important to note that all the above arguments about the expansion speed of the relics into the ICM are indirect. Since no expansion speeds can be measured, it is possible that the expansion of these objects is still supersonic, in excess of $`10^8`$ cm s-1, in which case their internal energies are even higher than if they were in pressure equilibrium. In addressing this issue, Soker, Blanton & Sarazin (2002) argue that certain combinations of effects can result in a mildly supersonic expansion producing x-ray signatures that are consistent with the data, though strong shocks appear to be ruled out. The consequences of possible supersonic expansion will be discussed subsequently.

### 1.1 The Influence Upon Cooling Flows

The geometries and energy densities of the radio cavities naturally suggest the occurrence of buoyant motion through the hot ICM in the presence of the gravitational potential of the central galaxy and cluster core. The evolution of such buoyant cavities has been suggested to be an important factor in the evolution of cluster cooling flows. It has been known for many years that the total energies present in extended radio sources in clusters ($`10^{57}10^{59}`$ erg) are enough to significantly influence the overall energy budget of the ICM. The problem has been, and may still be, in finding an effective mechanism for mixing this energy uniformly throughout the inner regions of the ICM. Several authors have recently suggested that buoyant radio cavities can accomplish this, either through advective mixing of differing regions of the ICM at different temperatures or through the dissolution of the radio cavities and the dispersal of their energetic and hot radio plasma into the ICM. The two dimensional axisymmetric hydrodynamic calculations of Reynolds et al. (2002) show some advection of cooler intracluster material into hotter regions of the ICM, though it is less clear that a truly buoyant cavity with a geometry similar to that seen in the x-ray data is produced by these simulations. This may be due in part to the effects of the axisymmetry constraint. The three dimensional hydrodynamic calculations of Brüggen et al. (2002) show more clearly the development of a buoyant cavity, and these simulations also show the development of Rayleigh-Taylor instabilities that lead to mixing with the ICM and eventual destruction of the cavity as a separate entity. This process could in principle inject a significant amount of energy into the ICM and can be more effective for reheating than the mixing of different regions of the ICM via advection, though this is also seen in the simulation of Brüggen et al. The two dimensional very high resolution hydrodynamic simulations of Brüggen and Kaiser (2002) show very clearly the onset and non-linear development of surface instabilities in the rising cavity, and these elegant calculations show that the cavity is destroyed and mixing well underway after the cavity has risen only a few of its own diameters. Presumably in three dimensions this process could be even more effective due to the larger number of modes available for the development of the instability, (assuming that the resolution of the numerical simulation is unchanged). Thus hydrodynamic consideration of buoyant cavities in the ICM show that very significant mixing of the cavity material with the ambient ICM can occur on short timescales, of order $`10^7`$ years, and this can transfer a significant amount of the energy from the radio source plasma to the intracluster medium. The final mixing of this material throughout the ICM has yet to be calculated; Brüggen & Kaiser (2002) provide an estimate of what might be the final state of the ICM after such mixing by averaging the energy input from the buoyant plumes over azimuth. This estimate shows that, if complete mixing can occur on a time short compared with inflow times, the injection of energy from radio sources can be a significant factor in reheating cooling flows and hence may be one solution to the long standing cooling flow problem.

## 2  The Evolution of Relic Radio Bubbles

Recently, new radio and x-ray observations have revealed an additional feature in the intracluster medium of some rich clusters, and this is the presence of pairs of what appear to be relic radio ”bubbles” that lie well outside the more luminous radio cavities in the inner core of the ICM. Spectacular examples of this phenomenon are found in the Perseus (A426) and A2597 clusters, though other examples have also been found (McNamara 2002b). These relic bubbles are coherent objects that are nearly spherical in appearance and have nearly featureless radio emission at low frequencies but are not easily seen at high frequencies (Fabian et al. 2002). They also appear to be in equipartition with the ambient ICM at pressures of $`14\times 10^{10}`$ dyne cm-2, and if so then they have internal energy densities that are again in excess of the equipartition values by factors of ten. The buoyant rise times to their current positions from the central galaxy are of order $`10^8`$ years, which exceeds the radiative lifetimes of the electrons in the inner lobes by factors of ten (McNamara et al. 2001, Churazov et al. 2001). This may imply the need for electron re-energization (see also Brüggen et al. 2002).

### 2.1 The Role of Intracluster Magnetic Fields

In light of the results from the numerical simulations, one of the most surprising aspects of the relic radio bubbles is that they are intact. If such objects were to follow the evolutionary path described by hydrodynamic simulations, one would expect that at distances $`30`$ kpc from the cluster center and at times $`10^8`$ years the radio remnants would have become fragmented and assimilated into the ambient intracluster medium. The fact that they have not done so may have important implications for the overall energy balance in the intracluster medium. The reason for the preservation of the relic radio bubbles may possibly be found in the magnetic fields that permeate the ICM. It has been known for some time that the hot gas in clusters of galaxies often has within it a significant magnetic field (e.g., Carilli & Taylor 2002, Taylor et al. 2002), with typical magnetic field strengths of order $`5\times 10^6`$ G. The origins of such fields remain somewhat obscure, but their effects on the evolution of radio sources in clusters are significant. As a radio source near the cluster center begins to inflate a cavity in the hot ICM, the ambient magnetic field is excluded from this cavity along with the hot intracluster gas. This will result in the external field forming a sheath around the cavity in which the field is primarily tangential to the cavity surface and has a higher value than the ambient field due to the effects of compression, as can be seen from simple flux conservation arguments. As mentioned above, the x-ray data imply that the cavity expands subsonically or transonically, and thus this region of compressed and largely azimuthal field will expand into the ambient ICM at the local signal speed, slightly ahead of the advancing boundary of the inflating cavity.

This layer of tangential magnetic field around the buoyant cavity will suppress the short wavelength and fastest growing modes of the Rayleigh-Taylor (R-T) and Kelvin-Helmholtz (K-H) instabilities. Because the lifting and mixing of different layers of the ICM is mediated primarily by these instabilities, as is large-scale mixing of the energetic radio emitting plasma with the ICM, the presence of this external layer of magnetic field will have an important influence on the evolution of the buoyant radio cavities. A magnetic field parallel to the interface between two fluids of differing density will suppress the R-T instability because the field acts as a source of surface tension at the interface and thus stabilizes small wavelength perturbations. This stabilization occurs for all wavenumbers larger than (e.g., Chandrasekhar 1961; Shore 1992)

$$k_c=\frac{2\pi g(\rho _2\rho _1)}{B^2Cos^2\theta },$$
(1)
where $`g`$ is the acceleration due to gravity, $`\rho _1`$ is the density of the lighter fluid (the radio bubble) and $`\theta `$ is the angle between the magnetic field vector and the wave vector of the perturbation. For a tangled azimuthal field an average value of $`<Cos^2\theta >=1/2`$ can be used. The value of $`k_c`$ provides the wavenumber of marginal stability, and formally the growth rate of this mode is zero. Perturbations of longer wavelength grow initially at a rate given by

$$n^2=\frac{gk}{(\rho _2+\rho _1)}\left[(\rho _2\rho _1)\frac{kB^2cos^2\theta }{2\pi g}\right],$$
(2)
where $`n`$ is the coefficient in the initial growth rate given by $`\mathrm{exp}nt`$. The wavenumber with the most rapid growth rate is found by differentiating Eq. 2 with respect to $`k`$. Note that the instability treated here is that of the surface of a buoyant bubble rising in a gravitational field and not that of the surface of a decelerating expansion. It is assumed here that the radio relics are in approximate pressure equilibrium with the ambient gas. If the relics are expanding supersonically then this will not be the case, but the age of these relics and their disconnection from the active nucleus, together with their geometry, strongly suggests that they have come into approximate equilibrium with their surroundings. The non-relic radio bubbles in the inner regions may be expanding at a mildly supersonic rate, as mentioned in Introduction, and such supersonic propagation speeds can act to suppress the R-T instability in those objects.

A conservative estimate that produces the largest value of $`k_c`$ and the fastest growth rate is to assume that $`\rho _2\rho _1`$. For the relic bubbles in A2597 and Perseus the ambient number densities are of order $`10^2`$ cm-3 at distances of $``$ 30 kpc from the cluster center, which is the appropriate distance for these relic radio bubbles. Using the gravitational potential of a central galaxy of mass $`10^{12}`$ solar masses and an average value of the magnetic field appropriate for cooling flow clusters of $`5\times 10^6`$ G (e.g., Carilli & Taylor 2002), which assumes no significant amplification of the field in the compressed sheath surrounding the bubble and which will maximize the value of $`k_c`$, one finds $`k_c=1.37\times 10^{22}`$, or

$$\lambda _c=2\pi /k_c=15.2kpc.$$

This is the shortest wavelength for the onset of the R-T instability under these conditions. The wavenumber of the fastest growing mode in the linear regime, obtained via differentiation of Eq. 2, is $`k_{}=k_c/\sqrt{3}`$, and substitution of this into Eq. 2 gives the maximum growth rate as

$$\mathrm{\Gamma }_{}=1/n_{}=4.2\times 10^7yr.$$

This is the time required for an initial perturbation of wavelength $`\lambda _{}`$ to grow by a factor of $`e`$; subsequent growth into the non-linear regime that will result in disruption of the radio cavity will require times greater than this by at least factors of two.

The interface between the buoyant cavity and the ambient ICM is also subject to the Kelvin-Helmholtz (K-H) instability between two fluids in relative motion. Again, the effective surface tension of the tangential ICM field that surrounds the cavity can suppress the onset of this instability. In the absence of a gravitational restoring force perpendicular to the fluid interface, which is the case here, the flow is stable against the K-H instability if (e.g., Chandrasekhar 1961)

$$(U_1U_2)^2\frac{B^2(\rho _1+\rho _2)}{2\pi \rho _1\rho _2},$$
(3)
where $`U_1`$ and $`U_2`$ are the velocities of the two fluids (in this case $`U_2=0`$ and $`U_1=U_{rel}`$), $`B`$ is the average value of the tangential magnetic field, and the densities $`\rho `$ have the same meaning as previously. The Kelvin-Helmholtz instability is of less importance as an agent for mixing the intracluster medium with itself and with the hot radio source plasma than is the Rayleigh-Taylor instability. This is because the K-H instability will, in its fully developed non-linear form, lead to a turbulent mixing layer along the surface of the radio source cavity (e.g., De Young 2002). This will influence a much smaller volume of the ICM than will the R-T instability that leads to the destruction and complete mixing of the radio source bubble with the ICM. The ratio of these two volumes is of order $`V_{KH}/V_{RT}\mathrm{\Delta }R/R`$, where $`R`$ is the radius of the bubble and $`\mathrm{\Delta }RR`$ is the thickness of the turbulent mixing layer.

The relative speed of the buoyantly rising bubble is clearly subsonic; the three-dimensional simulations of Brüggen et al. (2002) show a relative speed of $`2.5\times 10^7`$ cm s-1, and their analytic estimate of the terminal speed of a rising bubble yields a similar value. (See also Churazov et al. 2001.) Equation 3 actually applies to the case where the same magnetic field exists on both sides of the interface. The stability criterion basically states that if the relative speed of the two flows is less than the Alfvén speed, then the Kelvin-Helmholtz modes are stabilized, since perturbations are then damped by the fast moving Alfvén waves moving along the interface. For the case of the relic radio bubbles the magnetic field is different inside and outside the interface, and hence the relevant Alfvén speed for the two different regions must be used. In the exterior, the same intracluster conditions as used in calculating the onset of the Rayleigh-Taylor instability gives stability against the K-H modes if $`U_{rel}1.5\times 10^7`$ cm s-1. For the interior the conditions are less well known; pressure equilibrium plus the lack of x-ray emission in the Chandra band pass suggests the presence of a hot and rarefied gas, and if the interior number density is of order $`10^4`$ (cf. Sect. 3) and the magnetic fields are near the equipartition values of a few $`\mu `$G, then again stability against K-H modes occurs for $`U_{rel}3\times 10^7`$ cm s-1. Thus these speeds are comparable to or greater than the relative speeds of the buoyantly rising bubbles obtained from the numerical simulations, and the interface is marginally stable against the Kelvin-Helmholtz instability. The growth rates of this instability are of order $`nkU_{rel}`$, and for any appreciable mixing to occur (and for any observable deformation of the bubble interface) the wavelength of the instability should be of order $`1`$ kpc or more. This then gives growth times of order $`10^7`$ years or more.

### 2.2 Implications for Re-energization of Cooling Flows

Use of the buoyant speeds obtained from the simulations (Brüggen et al. 2002, Churazov et al. 2001) gives lifetimes for the relic radio bubbles in A2597 of $`10^8`$ years and $`5\times 10^7`$ years for those in A426. The three dimensional hydrodynamic simulations of Brüggen et al. (2002) and the similar high resolution two dimensional simulations of Brüggen & Kaiser (2002), which are of the appropriate scale for these objects, show that the rising cavities are clearly being disrupted at times of $`34\times 10^7`$ years at a distance from the nucleus of $`15`$ kpc for the three dimensional case and at $`6\times 10^7`$ yr and $`20`$ kpc for the two dimensional calculation. (The two dimensional case has a higher energy input than is appropriate for these cluster radio sources, and it also may be more stable due to suppression of some 3-D modes.) However, the relic radio bubbles are still intact and show no signs of disruption at distances of $`30`$ kpc and ages of $`10^8`$ years. Hence some additional processes other than purely hydrodynamical ones must be at work, and the above calculations show that the displaced intracluster magnetic field may provide the stabilizing influence that keeps the relic bubbles intact. In this case the Rayleigh-Taylor instability, which is the most disruptive, does not even commence its growth in the small amplitude linear phase until times of about $`5\times 10^7`$ years, and the *shortest* wavelength of the initial instability is $`15`$ kpc, which is comparable to the radii of the relic bubbles in A2597 and is comparable to the overall size of the bubbles in A426. This is to be compared with the purely hydrodynamic simulations, which show that the R-T instability has proceeded far beyond the linear scale at times of $`5\times 10^7`$ years, and at wavelengths that are much less than the overall radius of the rising bubble. Hence both the spatial and temporal scales for disruption of these relics in the presence of external magnetic fields are such that this instability has not caused disruption by their current age. This then implies a significant reduction of the mixing of energetic radio source material with the ambient ICM, and in addition it may imply a lowered efficiency in ”lifting” one part of the ICM into another. This is because, in the purely hydrodynamic case, the ”lifting” is due to boundary layer (as opposed to mixing layer) effects, and the unperturbed surface of the bubble has both a thinner boundary layer than a surface roughened by small wavelength instabilities and a smaller total surface area than that of a fully perturbed surface that is deformed by large wavelength instabilities. However, in the MHD case the presence of external magnetic fields may actually increase the drag on the bubble, and hence its lifting effect, because of the long range nature of the fields. A self consistent MHD calculation is required to determine the magnitude of this effect.

Additional support for magnetic stabilization comes from the two-dimensional MHD calculations of Brüggen and Kaiser (2001). These calculations consider a scale much larger than that appropriate for the relic radio bubbles, with initial configurations extending from $`200`$ to $`400`$ kpc and with source energies more appropriate for FR-II radio sources than for the FR-I objects considered here. In addition these calculations do not consider the effects of the ambient magnetic field in the ICM, but when a circumferential field is placed around the radio bubble ”by hand”, the resulting geometry is similar to that seen in the relic radio bubbles. Moreover, this configuration shows no signs of disruption or mixing with the ICM. It is important to note that this calculation of Brüggen and Kaiser employs a magnetic field strength that has an energy density of 10 percent that of the thermal energy in the ambient ICM (their ”weak field” case). While this field is dynamically unimportant, it is stronger than the cluster magnetic fields used here. However, the results are clearly indicative of the effects of tangential magnetic fields around a radio source bubble. The obvious next step is a self-consistent three dimensional MHD calculation that includes the ambient magnetic field in the intracluster medium, and this is underway. However, the above calculations for the onset of disruptive instabilities in the presence of the magnetic field in the ICM show that the transfer of energy from the radio bubbles to the intracluster medium may not be prompt, and that the effectiveness of such radio sources in reheating cooling flows may be less than originally suggested.

The final stage of transfer of energy from a radio relic to the ICM takes place through the turbulent dissolution of the radio source and its ultimate transfer into heat. The time scale for this very non-linear process is calculated next.

### 2.3 Turbulent Mixing of Radio source Debris and ICM Heating

Once a radio bubble has been disrupted by the Rayleigh-Taylor instability, it develops large substructures which then break down into ever smaller eddies and cells. This process can be seen in the high resolution 2-D hydrodynamic simulations of Brüggen & Kaiser (2002), and it ultimately results in turbulent flow. The current numerical simulations do not follow the evolution of the flow into this regime (Brüggen et al. 2002). Though such flows contain very fine structure, they do not actually heat the ambient ICM until the turbulent cascade has proceeded down to scales corresponding to the dissipation region. A question of relevance to the reheating of the ICM is the time scale for this process to occur. Using the wavelengths of the initial instability found above, it is possible to calculate the development of fully non-linear MHD turbulence in three dimensions (e.g., Orszag 1970, De Young 1992). The method employed solves the time dependent MHD equations in their non-linear form, and the use of Fourier transformed second order moments of velocity and magnetic field allows a calculation of the time evolution of both kinetic and magnetic energy over many orders of magnitude in spatial resolution. A wide variety of initial and boundary conditions can be accommodated, and details of the method of calculation, including the energy transfer among different wavenumbers and the closure relations and dissipation scales, can be found in the above references.

The calculation here begins with continuous injection of energy at a single wavelength, which is the wavelength of the most unstable and fastest growing mode calculated in Sect 2.1. Energy injected at this scale quickly causes perturbations of the flow at smaller scales, and this causes the development of a turbulent cascade of energy to ever smaller scales or larger wavenumbers. The multiple interactions among turbulent eddies on different scales results in the creation of a steady state turbulent flow with a loss free propagation of energy through an inertial range into a dissipation scale. The calculation assumes that the magnetic field is homogeneous and weakly isotropic (allows helicity), since the non-linear stage of the R-T instability, once established, will disrupt the original tangential field and advect it with the large scale eddies as they rotate. Equipartition of the turbulent magnetic field with the kinetic energy will occur only on the smallest scales. If for any reason the large scale tangential field persists, it will inhibit the turbulent cascade and increase the time required for the turbulent energy to be injected as heat into the intracluster medium. Once an equilibrium state has been established, heat is being injected into the ICM at the same rate as energy is being extracted from the dissolving radio relic. Figure 1 shows the results of such a calculation in dimensionless form. The time development of the turbulence spectrum from an initial delta function injection of energy at $`k=1`$ is clearly seen, and it is also clear that an equilibrium spectrum is obtained after about $`10`$ large scale eddy turnover times. (Figure 1 shows the spectrum of turbulent kinetic energy; the spectrum of the turbulent magnetic energy is very similar for wavenumbers above $`k=10`$ where the magnetic and kinetic energy quickly come into equipartition.)

The connection to the relic radio bubbles is made by noting that the shortest wavelength for onset of the Rayleigh-Taylor instability is $`15`$ kpc, and that the wavelength for the fastest growing mode is $`\sqrt{3}`$ times this. Using this as the scale of the large eddy injection region, and using an eddy rotation speed comparable to that of the speed of the rising radio bubbles from the hydrodynamic simulations ($`300`$ km s-1) gives

$$\tau _{edd}=\lambda _{edd}/v_{edd}=4.4\times 10^7yr.$$

This provides the scaling factor for the scale free result shown in Figure 1, and it then implies that once the instability begins, the actual heating of the ICM commences about $`4\times 10^8`$ years later. The calculation shown in Figure 1 does not extend to the actual dissipation range, which lies at much smaller scales than considered in the figure. However, extension of the turbulent spectrum to the dissipation range will occur on times much smaller than the time required to establish equilibrium over the range shown due to the much more rapid turnover time of the small scale turbulent cells. Hence this calculation, when coupled with the calculation for the onset of the Rayleigh-Taylor instability in the presence of the displaced ICM magnetic field, shows that the transformation of radio source energy into heat input to the ICM may not take place until about $`5\times 10^8`$ years after the initial formation of the buoyant radio source bubble. This is a time comparable to the cooling time in the central regions of rich clusters (e.g., Fabian 1994), and hence the overall energy balance between cooling of the ICM and possible reheating due to the energy injected from radio sources may need to be reconsidered.

As mentioned at the end of Section 1, once the energy is placed into the ICM as heat, it must be distributed throughout the ICM if it is to have a global effect on the cooling flows in clusters, and Brüggen & Kaiser (2002) estimated this effect via an averaging calculation. This global dissipation and mixing will take additional time, presumably comparable to the orbital periods of the galaxies in the cluster in the absence of large scale ICM motions due to subcluster mergers. However, it is possible that local injection of heat can have a significant effect on the cooling flow phenomenon because local injection of heat can result in convective turnover being induced in the ICM on scales comparable to that of the relic and now virtually destroyed radio bubble. This convective activity can serve to mix hotter and cooler regions of the ICM, and if multiple events such as this were to occur, an overall effect on the cooling flow could result. The timescale for this to take place, together with its spatial extent, must be determined by a numerical calculation, but $`5\times 10^8`$ years seems an appropriate lower limit, since the convection begins as soon as heating commences.

### 2.4 Relation to Other Radio ”Relics” in Clusters

Diffuse radio emission from clusters of galaxies has been observed in many clusters (e.g., Carilli & Taylor 2002), and in some cases the radio morphology of this emission is suggestive of that seen in FR-I or FR-II radio sources, even though a parent galaxy cannot be clearly identified. Recent high resolution observations of four of these possible ”relics” (Slee et al. 2001) show very extended, diffuse and filamentary morphologies. This appearance is very different from the relic radio bubbles considered here, and the sources observed by Slee et al. do not have clear and unambiguous identifications of the parent galaxies. The radio morphologies are suggestive of much later stages of evolution of radio relics, possibly under the influence of large-scale shearing motions in the intracluster medium. Kaiser & Cotter (2002) have modeled these objects as relics of FR-II radio sources. Though FR-II objects are not commonly seen in rich clusters, if these are in fact FR-II remnants then their evolutionary tracks may be different than those of the lower powered relic radio bubbles considered here. In any case, it is interesting that even these sources do not appear to have been completely assimilated into the ambient ICM, even at the ages of greater than $`10^8`$ years considered by Kaiser & Cotter. This suggests again that the magnetic fields in both the ICM and in the radio relics naturally inhibit the complete mixing of the radio source with its environment, at least at some level. Calculation of the transfer of energy from the radio source to the ambient ICM for these large sources and at late times would be of interest in relation to cluster cooling flows.

## 3 Implications for the Particle Content of Radio Sources

A long standing issue in understanding extragalactic radio sources is that of the nature of the energetic particle populations contained within these objects (e.g., De Young 2002). While relativistic electrons are clearly present, the nature of the neutralizing particles, whether protons or positrons, is still under debate. The relic radio bubbles appear to be in pressure equilibrium with the hot ICM, yet they do not emit x-rays at the same intensity or energy as does the ambient medium. The question then arises as to the nature of the material in the bubble that contributes the required energy density; it could be magnetic field, hot protons, or very hot positrons. Hot positrons could arise from the very low energy tail of a power-law distribution of electrons and positrons. If so, care has to be taken in choosing the low energy cutoff of such a distribution; it must supply the required pressure, but it must not be such that the electrons and positrons, themselves subject to Bremsstrahlung, will emit x-rays above the observed limits. For example, if a mildly relativistic electron-positron population is used with $`\gamma 1`$, as has been suggested for some compact electron-positron jets (e.g., Celotti & Fabian 1993), then pressure equilibrium with the surrounding ICM at $`p10^{10}`$ dyne cm-2 will give rise to a Bremsstrahlung emissivity per unit volume from both electrons and positrons that is roughly 10 percent of the radiation from the ambient ICM. This emission should lie within the sensitivity limits of the Chandra Observatory HRC but at the upper end of its energy bandpass.

An insight into this issue may arise from the possible need to reaccelerate the synchrotron emitting electron populations in the relic radio bubbles. The dynamic lifetimes are in excess of the radiative lifetimes (e.g., McNamara et al. 2001), and the requirement for reacceleration also emerges from the calculation of synchrotron radiation expected from the hydrodynamical simulations (Brüggen et al. 2002). Under the assumption that reacceleration is the relevant process, and if such reacceleration is to take place through internal shocks, Fermi acceleration, or turbulent acceleration via stochastic processes, then the fluid within the bubble must be able to sustain such motions and density inhomogeneities for a time long enough to provide significant reacceleration. If this fluid is an electron-positron gas, and if it provides enough internal energy to be in pressure balance with the hot external ICM, then this fluid is mildly relativistic. Such a hot gas has very high sound speeds, and as a result it will rapidly damp any discontinuities that produce shocks and will also not sustain large density or pressure inhomogeneities for significant times; i.e., for times in excess of many sound crossing times across the inhomogeneities. The exception to this is if significant cooling can be made to take place in density inhomogeneities, but in that case such areas would be extremely bright in x-rays. Hence it seems less likely that an electron-positron gas can both provide pressure equilibrium in the relic radio bubbles and also sustain the needed reacceleration of the radio emitting electron population. This is not the case for a electron-proton gas, which can be hot but not relativistic ($`T10^{10}`$ K) and rarefied ($`n10^4`$). This fluid will have low x-ray emissivity in the observed pass bands and will be in pressure equilibrium with the surrounding ICM. In this context it is interesting to note that thermalization of a mildly relativistic jet ($`v_j0.03c`$) through a termination shock at the ”working surface” will result in a temperature of $`10^{12}`$ K. Would such reacceleration processes produce observable inhomogeneities in the radio emission from the relic bubbles? As mentioned above, the emission seems homogeneous within the constraints of the low frequency observations. Inhomogeneities in emission would be expected if very large scale shock structures, comparable in size to the radio bubble, were present, but such structures would not be expected in an aged radio relic. A much more likely process is turbulent reacceleration, and this would occur more uniformly throughout the radio source volume on a ”fine grained” scale. Thus it would be unlikely to produce large scale intensity variations in the radio emission. (An alternative to reacceleration is that of ”storage” of energetic particles in regions of low magnetic field, with radiative losses taking place only in high field regions in an inhomogeneous magnetic field distribution. This process has been suggested as a mechanism for larger FR-II and FR-I radio sources (e.g., Eilek, Melrose & Walker 1997, Tribble 1994), and though it cannot be ruled out here, the radio emission from the relic radio bubbles does not display inhomogeneities or filamentary structures on the scale resolved by current radio data.)

If this reacceleration process is taking place as described, it implies that the electron-proton gas was present in the radio source at early times, since the very existence of the well defined and undisrupted relic radio bubbles implies that there has been no significant mixing with the ICM, and that magnetic fields have kept the internal and external fluids separated from the time of creation of the radio source cavity. This would imply that the electron-proton gas was either ejected from the nucleus of the parent galaxy or that the protons were entrained into the outflowing jet as it passed through the parent galaxy’s interstellar medium. If magnetic isolation also exists in the outwardly moving jet (which is more difficult to accomplish), then this would be evidence for an electron-proton gas being present in the jet very close to its point of origin.

## 4 Conclusions

Recent observations of relic radio ”bubbles” in the central regions of rich clusters imply that the lifetime of these objects is longer than that obtained from calculations of bubble disruption via the onset and development of hydrodynamic Rayleigh-Taylor and Kelvin-Helmholtz instabilities. Inclusion of the effects of the magnetic field in the ambient intracluster medium indicates that the creation of a tangential field on the surface of these bubbles as they expand can stabilize them against such instabilities for times comparable to their current lifetime. Calculation of the time of onset of instability under these conditions, together with calculation of the time required for the energy in the radio source to be dissipated as heat in the ICM, suggests that significant reheating of the ICM by radio sources may not commence until about $`5\times 10^8`$ years after the creation of the radio source. As this time is comparable to the cooling times in the centers of rich clusters, the role of radio sources as a means of reheating cooling flows may be more complex than earlier thought. Finally, the existence of relic radio bubbles, together with their likely pressure equilibrium with the surrounding intracluster medium, can be used to argue for the presence of a hot gas interior to the bubbles that is composed primarily of electrons and protons rather than electrons and positrons. The isolating effects of magnetic fields interior and exterior to the radio bubbles then suggests that the particle content of the original radio jet as it emerged from the parent galaxy was primarily electrons and protons.

I thank Brian McNamara for many useful discussions and a referee for helpful comments.


format: text

text:
# Theory of Melting and the Optical Properties of Gold/DNA Nanocomposites

## Abstract

We describe a simple model for the melting and optical properties of a DNA/gold nanoparticle aggregate. The optical properties at fixed wavelength change dramatically at the melting transition, which is found to be higher and narrower in temperature for larger particles, and much sharper than that of an isolated DNA link. All these features are in agreement with available experiments. The aggregate is modeled as a cluster of gold nanoparticles on a periodic lattice connected by DNA bonds, and the extinction coefficient is computed using the discrete dipole approximation. Melting takes place as an increasing number of these bonds break with increasing temperature. The melting temperature corresponds approximately to the bond percolation threshold.

The optical properties of metallic nanoparticles have been investigated intensively over the last decades . Recently, this investigation has expanded to include so-called functional metallic nanoparticles. Among these, there is a particular interest in DNA-modified gold nanoparticles. Such particles can form complex aggregates which may melt at high temperatures. The particles and their aggregates may have a variety of applications, e. g., in biological detection, which may be possible by using the optical and electrical sensitivity of the aggregates . Numerical model calculations of the optical properties of DNA modified gold nanoparticle aggregates show general agreement with experiments, including such features as the following: (i) for isolated gold nanoparticles in suspension, there is a strong ”surface plasmon” absorption in the visible; and (ii) this absorption maximum broadens and red-shifts when the cluster radius becomes comparable to the wavelength . However, some important physical details have not been explained. For example, in the gold/DNA nanoparticle system, the melting transition of a gold nanoparticle aggregate has a much narrower temperature width and occurs at a high temperature than that of a single DNA link . Also, the dependence of aggregate melting temperature on particle size has not yet been explained .

In this Letter, we model this novel “melting” transition of a gold/DNA nanoparticle aggregate. Our model accounts for most experimentally observed features, including (i) the small temperature width of the aggregate melting transition, compared to that of a single DNA link; (ii) the particle-size-dependence of the melting temperature $`T_m`$; and (iii) the temperature-dependence of the optical extinction coefficient $`C_{ext}(\lambda ,T)`$ at wavelength $`\lambda `$ and temperature $`T`$.

We first describe our model for the low-temperature morphology of the aggregate. In the simplest version of our model, the low-T cluster is taken simply as a collection of identical gold nanoparticles (each of radius $`a`$) which fill the sites of a simple cubic lattice of lattice constant $`d`$ ($`d>2a`$). The cluster is assumed to be a cube of edge $`L`$, containing $`N_{\mathrm{par}}=(L/d)^3`$ gold nanoparticles. We have also investigated the melting and the optical properties assuming that the low-T cluster is a fractal aggregate.

To describe the melting of this aggregate, we assume that each nanoparticle carries exactly $`N_s`$ single DNA strands. (This choice still retains the essential features of the melting.) The cluster exists at low $`T`$ because a chemical reaction converts two single strands $`S`$ into a double strand $`D`$. In this Letter, we adopt a simplified two-state model for this reaction, described by the relation

$$S+SD.$$
(1)
For short DNA strands (12-14 base pairs), this model describes melting well. Each gold nanoparticle has $`z`$ nearest neighbors ($`z=6`$ for a simple cubic lattice). Hence, there are $`N_{\mathrm{par}}z/2`$ bonds joining adjacent nanoparticles. We assume that a (temperature-dependent) fraction $`p(T)`$ of the single DNA strands form double strands by the reaction (1). In order to calculate the fraction $`p_{\mathrm{eff}}(T)`$ of bonds which contain at least one double strand, we adopt the following model. First, we assume that exactly $`N_s/z`$ of the single strands on a given nanoparticle are available to bond with any one of its z nearest neighbors. Thus, we assume that the maximum number of links that could be formed between any two particles is $`N_s/z`$. The probability that no link is formed is then taken to be

$$1p_{\mathrm{eff}}(T)=[1p(T)]^{N_s/z}.$$
(2)
The criterion for the melting temperature $`T_c`$ is that $`p_{\mathrm{eff}}(T_c)=p_c`$, $`p_c`$ being the bond percolation threshold for the lattice considered, at which an infinite connected path of double DNA strands first forms. (For example, $`p_c0.25`$ on a very large simple cubic lattice.) Thus, at percolation,

$$[1p(T_c)]^{N_s/z}=1p_c.$$
(3)
This equation implicitly determines $`T_c`$ in terms of $`z`$, $`p_c`$, and $`N_s`$.

$`p(T)`$ itself is determined by a chemical equilibrium between the single-strand and double-strand DNA molecules, which are attached to the gold nanoparticles. The chemical equilibrium condition corresponding to (1) is

$$\frac{[1p(T)]^2}{p(T)}=\frac{K(T)}{C_T}K^{}(T),$$
(4)
where $`K(T)`$ is a suitable chemical equilibrium constant, and $`C_T`$ is the molar concentration of single DNA strands in the sample (in our model all are attached to the nanoparticles). Since $`0<p(T)<1`$, the physical solution to eq. (4) is $`p(T)=1+\frac{1}{2}\left(K^{}\sqrt{K^{,2}+4K^{}}\right)`$.

Note that, according to eq. (3), $`p(T_c)`$ decreases with increasing $`N_s`$ and, hence, with increasing particle radius. Since $`p(T)`$ decreases monotonically with $`T`$, $`T_c`$ should thus be an increasing function of $`a`$, as reported in experiments.

In Fig. 1, we plot peff(T) for several $`a`$. We assume $`N_sa^2`$, set $`z=6`$ and use the experimental result that $`N_s=160`$ when $`a=8`$ nm . We have also assumed the simple van’t Hoff behavior $`K(T)=\mathrm{exp}[\mathrm{\Delta }G/k_BT]`$, with a Gibbs free energy of formation $`\mathrm{\Delta }G(T)=c_1(TT_M)+c_3(TT_M)^3`$, choosing the values of $`c_1`$, $`c_3`$, and $`T_M`$ to be consistent with experiments on these DNA molecules.

To calculate the $`T`$-dependent optical properties, we have considered the melting of two slightly different low-temperature aggregates. In the first, we assume that, at low $`T`$, the aggregate consists of a simple cubic collection of $`N_{\mathrm{par}}=(L/d)^3`$ gold nanoparticles on a simple cubic lattice, as described above, with all bonds occupied by DNA double-strand links. To generate a specific sample with a given peff(T), we randomly remove links with probability $`1p_{\mathrm{eff}}(T)`$, then identify the separate clusters, using a simple computer algorithm. We neglect gravitational forces, which may be important in some experimental circumstances. If the aggregate consists of two or more clusters, we simply place these aggregates in random positions and orientations within a larger bounding box (usually of edge 100$`d`$), taking care that the individual clusters do not overlap. The resulting geometry is shown schematically in Fig. 2. We have also carried out the same procedure to simulate the melting of a sample formed by reaction-limited cluster-cluster aggregation (RLCA) . A typical RLCA cluster is shown in Fig. 2 (d), and represents a possible fractal aggregate which might be produced by certain random growth processes at low $`T`$.

We calculate the optical properties of this sample using the Discrete Dipole Approximation (DDA). The sample is modeled as a collection of separate aggregates whose extinction coefficients are computed individually, then added. Each aggregate consists of many identical nanoparticles, which have complex frequency-dependent dielectric constant $`ϵ(\omega )`$, and polarizability $`\alpha (\omega )`$ related to $`ϵ(\omega )`$ by the Clausius-Mossotti equation $`ϵ(\omega )=1+4\pi n\alpha /[1(4\pi n/3)\alpha ]`$, where $`n=(d/L)^3`$. The resulting expressions for the induced dipole moment $`𝐩_i`$ of the $`i^{th}`$ sphere, and the corresponding expression for the extinction coefficient $`C_{ext}(k)`$ at wave number $`k`$, are given in Ref. .

In our case, each cluster consists of a number of DNA-linked individual gold nanoparticles. In our calculations, we do not include the optical properties of the DNA molecules, since these absorb primarily in the ultraviolet . We use tabulated values of the gold complex index of refraction, then calculate $`C_{ext}`$ for each cluster using the DDA. To improve the statistics, we average $`C_{ext}`$ for each cluster over possible orientations. We then sum the averaged extinction coefficients of all the individual clusters to get the total extinction coefficient of the suspension. This method is justified when the suspension is dilute.

In Fig. 3, we show $`C_{ext}(\lambda ,T)`$ versus $`\lambda `$ for the gold-DNA cluster at several values of peff(T), assuming $`a=20`$ nm. peff = 0 represents a dispersion of $`N_{\mathrm{par}}=(L/d)^3`$ individual nanoparticles, while peff = 1 represents a simple cubic lattice of connected nanoparticles. The calculated $`C_{ext}(\lambda ,T)`$ at each $`\lambda `$ changes strikingly for peff(T) $``$ pc. For small peff, there is a clear extinction peak near $`520`$ nm. This peak corresponds to the wavelength of the surface plasmon resonance (SPR) in individual Au nanoparticles. As peff increase, this peak first red-shifts, then greatly broadens, as the aggregate melts. The calculated peak is, however, shifted much more and is much broader than experiment. However, as shown in inset, if we assume that the low-$`T`$ aggregate is an RLCA fractal, the calculated peak shift is consistent with experiment .

In Fig. 4, we show $`C_{ext}(\lambda ,T)`$ at fixed $`\lambda =520`$ nm (close to the isolated-particle SPR), versus $`T`$, for several particle sizes, assuming that the low-T aggregate is a simple cubic cluster with $`N=1000`$ particles, as shown in Fig. 2. For each $`a`$, the extinction increases sharply at a characteristic $`T`$, corresponding to the melting of the aggregate for that $`a`$; at this $`T`$, the absorption due to the SPR increases sharply.

In the inset of Fig. 4, we compare $`C_{ext}(\lambda ,T)`$ for a regular and an RLCA cluster of particles of $`20`$ nm radius at $`\lambda `$ = 520 nm, both for $`N=1000`$ particles. Although the RLCA cluster has a slightly broader melting transition, as manifested in $`C_{ext}(\lambda ,T)`$, than does the regular lattice, both sets of data show a much sharper melting transition than that of a single DNA link. Also, although our normalized $`C_{ext}(\lambda ,T)`$ is calculated for the aggregates at 520 nm, we expect similar behavior at $`260`$ nm. (We have not carried out calculations at this $`\lambda `$ mainly because we have not included the DNA absorption properties.) In any case, the experimental melting curves at 260 and 520 nm are very similar.

Our calculated extinction coefficients are strikingly similar to recent experimental results. In particular, both experiment and calculation give a sharp increase in $`C_{ext}(\lambda ,T)`$ at fixed $`\lambda `$, as $`T`$ increases past a critical temperature, which we interpret as the melting temperature $`T_m`$. We also find, in agreement with experiment, that melting occurs over a much narrower range of $`T`$ in the aggregate than for a single bond, and that the melting occurs at higher $`T`$ for larger particles. The crucial point is that $`p_{\mathrm{eff}}(T)`$ is a much sharper function of $`T`$ than is $`p(T)`$, and this feature would not be affected by slight changes in the model (such as considering a body-centered-cubic rather than a simple cubic cluster). In the present model, the melting transition is in the universality class of bond percolation. This universality class might change because of the constant dissolving and reforming of the bonds of the aggregate (a possibility omitted from our model). This possibility should be studied further theoretically and experimentally. However, the optical properties may well be robust, in the sense that $`C_{ext}(T)`$ would be little affected even if this feature were included, since the key ingredient is the multiple links per bond as described above.

In summary, we have developed a simple model for the melting of aggregates of gold nanoparticles and DNA, and have calculated the T-dependent optical properties of these melting aggregates, using the DDA. We find that, at fixed $`\lambda `$, melting is accompanied by dramatic changes in the extinction coefficient $`C_{ext}(\lambda )`$. These calculated changes occur over a much narrower temperature range than that over which interparticle links themselves melt, the temperature width is smaller for particles with larger radius, and the melting occurs at higher temperatures for the larger nanoparticles than for the smaller ones. All these effects are in good agreement with experiment. It would be of great value if the predicted $`T`$-dependent structure of the clusters could be probed experimentally, e. g. by static light scattering, as has been done in other contexts. Finally, we mention that the present model can be applied to other similar nanocomposites.

This work has been supported by grant NSF DMR01-04987, by the U. S./Israel Binational Science Foundation, and by an Ohio State University Postdoctoral Fellowship awarded to S. Y. Park. We thank D. J. Bergman, C.-H. Kiang, A. A. Lazarides, and D. R. Nelson for valuable conversations.


text:
# Giant Pulses from PSR B0540-69 in the Large Magellanic Cloud

## 1 Introduction

PSR B0540$``$69 was discovered in the early 1980s by Seward, Harnden & Helfand (1984) using data from the Einstein X-ray Observatory. The pulsar is located inside the supernova remnant SNR 0540–693 in the Large Magellanic Cloud. It has a short rotation period ($``$50 ms) and a rapid spindown with a characteristic age of only 1500 yr. There are strong similarities between this system and the Crab. The pulsars have similar characteristics and the plerion around PSR B0540$``$69 is very similar to the Crab Nebula in size and energetics (Manchester, Staveley-Smith & Kesteven 1993). Following the discovery of the X-ray pulsations from PSR B0540$``$69, optical pulsations were soon detected (Middleditch & Pennypacker 1985), but the radio pulsations remained elusive. Finally, Manchester et al. (1993; hereafter M93) discovered a broad, weak radio pulse at 0.64 GHz. The flux density at that frequency is 0.4 mJy and the pulse duty cycle is more than 80% with the hint of a double profile. This is similar to the high energy profiles (de Plaa, Kuiper & Hermsen 2003) which also appear to show a broad, unresolved double structure. M93 failed to detect the pulsar at either 1.4 or 0.44 GHz.

In comparison, the radio profile of the Crab pulsar is complex. It consists of a so-called pre-cursor component and main and interpulse components. The pre-cursor has a very steep spectral index rendering it undetectable above about 1 GHz. The main and interpulse emission are in phase with the X-ray and $`\gamma `$-ray emission which both have narrow peaks. At frequencies above 5 GHz, further components have been detected (Moffett & Hankins 1999). The Crab pulsar also produces giant pulses when the emission can exceed the mean pulse energy by a factor of 1000 or more. This giant pulse emission occurs at the phases of the main pulse and interpulse but not at those of the pre-cursor or the ‘anomalous’ high frequency components (Lundgren et al. 1995; Hankins et al. 2003).

Giant pulses have now also been discovered from the millisecond pulsars PSRs B1937+21 (Cognard et al. 1996) and B1821–24 (Romani & Johnston 2001). Cognard et al. (1996) noticed that when pulsars are ordered by their magnetic field strength at the light cylinder, the pulsars with giant pulses are highly ranked. Romani & Johnston (2001) noted the similarity between the phases of the giant pulse emission and the X-ray emission in these pulsars and speculated that the giant pulse emission may also originate from near the light cylinder. McLaughlin & Cordes (2003) have very recently reported the non-detection of giant pulses from PSR B0540$``$69 through observations made at 0.66 and 1.52 GHz.

## 2 Observations and Data Reduction

Observations were made on 20 May 2001 (MJD 52049), with the Parkes 64-m radio telescope. We used the center beam of the 21-cm multi-beam system at an observing frequency of 1.38 GHz. The receiver had a system equivalent flux density of 27 Jy on cold sky. The back-end consisted of a filterbank system containing 512 channels per polarization, each of width 0.5 MHz for a total bandwidth of 256 MHz. The polarization pairs were summed to form the total power, each output was then sampled at 80 $`\mu `$s, one-bit digitized, and written to DLT for off-line analysis. A total of 8 hr of integration on PSR B0540$``$69 was obtained, split into 4 individual observations, each of duration 2 hr.

To search for giant pulses, it is not necessary to know the rotational frequency of the pulsar but it is necessary to know the dispersion measure (DM). As an initial trial, we used the DM of 146 cm-3pc given in M93. After applying the appropriate delay to each frequency channel, the 512 channels were summed. The mean and rms of groups of 8192 samples (0.65 s) were examined and those which showed obvious signs of interference were discarded. Our nominal 5$`\sigma `$ sensitivity in 80 $`\mu `$s is 0.9 Jy. In practice, even after the removal of high sigma points (clipping), we experienced substantially larger background fluctuations and the rms exceeded our expected rms by a factor of $``$2. The entire 8 hr observation (3.6$`\times 10^8`$ samples) was then searched for peaks. This technique was successfully used to discover giant pulses in PSR B1821–24 (Romani & Johnston 2001) and the giant micro-pulses in PSR B1706–44 (Johnston & Romani 2002). Two large amplitude, resolved pulses were discovered. We were able to determine the DM of these single pulses by subdividing the total bandwidth into 16 and fitting a straight line to their time-of-arrival across the 16 sub-bands. The DM of both pulses was determined to be 146.5$`\pm `$0.2 cm-3pc. Armed with this knowledge we then applied a DM of 146.5 cm-3pc to the entire data set.

We then used a more sophisticated search technique to detect weaker giant pulses. The DM introduces a broadening of 220 $`\mu `$s across each 0.5 MHz channel and the full width at half maximum of the strong giant pulses was $``$600 $`\mu `$s, or 8 samples. We therefore first smoothed the data using a top hat function with a width of 8 samples and then searched the data looking for 4 consecutive data points more than 3$`\sigma `$ from the mean.

## 3 Pulse Timing and Phase Comparisons

In order to obtain the correct pulsar phase for each time sample, an accurate timing solution for the pulsar is required. The latest published ephemeris for PSR B0540$``$69 from X-ray timing with the Rossi X-ray Timing Explorer (RXTE) does not cover the epoch of our observations (Zhang et al. 2001) and we therefore used an updated ephemeris provided to us by F. Marshall and W. Zhang from RXTE observations near MJD 52050, including overlap with our observing epochs. The peak of the X-ray profile at the solar system barycenter occurred 0.0366 s after MJD 52050 and at this epoch the spin frequency was 19.7886114394 s-1 and the frequency derivative was $``$1.875580$`\times 10^{10}`$ s-2. The position of the pulsar used in the ephemeris was that obtained from optical observations by Caraveo et al. (1992). Sample RXTE arrival times on the same dates as the radio observations agree with the adopted ephemeris, with uncertainties $`0.5`$ms. Using this timing solution we were therefore able to form an integrated profile of the entire 8 hr observation, and we assert that uncertainties in the X-ray ephemeris would be insufficient to broaden the pulse by more than $``$0.5 ms ($``$1% of the pulse period).

Direct comparison of the phases of the radio giant pulses with the X-ray profile is, however, subject to some uncertainty. We believe the absolute timing in the radio to be correct, as we have obtained arrival times for both the Crab pulsar and PSR B1937+21 which agree with independent radio ephemerides. However, the absolute X-ray timing is not secure and it is known that there are problems comparing RXTE timing of PSR B0540$``$69 and those made with the Chandra X-ray Observatory (A. Rots, private communication).

## 4 Integrated Profile

We failed to detect any significant integrated pulsed emission from the pulsar. The bottom panel of Fig. 1 shows the pulse profile with 128 phase bins after 8 hr of observations. The 1-$`\sigma `$ rms is 80 $`\mu `$Jy. As a check against the possibility that the X-ray ephemeris was incorrect we also performed a Fast Fourier Transform on each de-dispersed 2 hr dataset independently. No significant signals were detected near the expected pulsar frequency. Assuming the pulsar has a duty cycle of $``$10%, the upper limit on the continuum flux density at 1.38 GHz is therefore only 13 $`\mu `$Jy compared to 400 $`\mu `$Jy at 0.64 GHz (M93). This implies a spectral index steeper than –4.4 between the two frequencies, although this could be relaxed if the duty cycle of the profile is large and/or the 0.64 GHz flux density is incorrect. This is extremely high, even by pulsar standards, but it is comparable with the spectral index of the pre-cursor component in the Crab pulsar. The distance to the Large Magellanic Cloud is 49.4 kpc, hence the luminosity of the pulsar at 1.38 GHz is less than 32 mJy kpc2. Although this is low compared to the luminosity of other young pulsars, and a factor of 100 less than the luminosity of the Crab pulsar, it is significantly higher than that of the recently discovered PSR J0205+6449 in the SNR 3C58 (Camilo et al. 2002).

## 5 Giant Pulses

As described above, the initial search of the data detected two giant pulses. These are shown in the top panel of Fig. 1 where both the entire profile is shown and a more detailed view of the structure in the giant pulses. The widths of these pulses are entirely consistent with an intrinisic delta function convolved with the sampling rate (80 $`\mu `$s), the dispersion delay (220 $`\mu `$s) and an exponential tail with a 1/e time of 0.4 ms. This exponential broadening is likely the result of interstellar scattering. The two pulses are offset in phase by 0.03.

Using the analysis described in Section 2 (data smoothing plus searching for several large consecutive samples), we searched a range of DMs in order to determine the cut-off above which no spurious signals were detected. The best results were obtained using the criteria that four consecutive samples should exceed 3$`\sigma `$. We applied this analysis on the data at a DM of 146.5 cm-3pc and detected 15 giant pulses down to this sensitivity limit. Fig. 2 shows the flux density as a function of phase for the 15 detected giants superposed on the X-ray profile taken from De Plaa et al. (2003). There are two separate phase regions for the giant pulses, the first between phases 0.36 and 0.45 and the second between 0.61 and 0.66 phase. These phases should be considered preliminary in light of the discussion about absolute timing above. This is similar to the Crab giant pulses which also occur in two separate phase regions. However, the jitter in the phases of the giants in PSR B0540$``$69 is much larger than in the Crab; $``$5 ms compared to $``$0.5 ms. Fig. 3 shows the histogram of the flux densities obtained. Although we are dealing only with small number statistics, it seems clear that the slope is not as steep as the slope of –3.5 seen in the Crab pulsar (Lundgren et al. 1995) or that of –2.8 seen in PSR B1937+21 (Kinkhabwala & Thorsett 2000).

Two other points are worth noting. The first is that the brightest giant obtained has a flux density which exceeds the mean single pulse flux density by a factor $``$5000. This is a larger factor than any giant pulse seen in the Crab pulsar. All the giants detected here exceed the mean flux density by 700. Expressed in these terms, this rate of giants (1 in 35000) exceeds that of the Crab pulsar (1 in 350000; see Lundgren et al. 1995) by an order of magnitude. However, in terms of the emitted energy, the giants in PSR B0540$``$69 are significantly less powerful than in the Crab pulsar.

Finally, the non-detection of giant pulses in PSR B0540$``$69 at 1.5 GHz by McLaughlin & Cordes (2003) can probably be attributed to their observational setup. Their sampling rate was 0.6 ms, however, the dispersion smearing across each of their 5 MHz filters was $``$2 ms and this is sufficient to dilute the giant pulse enough for it to fall below their quoted detection threshold of 0.4 Jy. Their non-detection of giant pulses at 0.66 GHz would imply that either the spectral index of the giant pulses are shallower than –3, in contrast to the integrated profile, or that there is no giant pulse emission at this frequency. The latter seems unlikely given that both the Crab and PSR B1937+21 show giant pulse emission at frequencies well below 0.66 GHz. McLaughlin & Cordes (2003) do not comment on a detection or non-detection of the integrated pulse at this frequency.

## 6 Discussion

PSR B0540$``$69 is the 4th pulsar known to emit giant pulses and there are now two young pulsars and two millisecond pulsars which emit giants. We detected the giant pulses without detecting the underlying integrated pulsed emission. The four pulsars are in the top five pulsars when ranked by magnetic field at the light cylinder; only PSR B1957+20 from this list does not show giant emission. Clearly therefore this parameter is important to our understanding of the giant pulse phenomenon. Furthermore, all four objects have pulsed, non-thermal, X-ray emission. In the Crab, PSRs B1937+21 and B1821–24 the X-ray emission appears to be in phase with the giant radio emission. It is intriguing that for PSR B0540$``$69, the two giant pulse windows are $`0.5`$ out of phase with the X-ray pulse and that their spacing is similar to that of the high energy pulse components. Whether this indicates an error in the relative phasing between the radio and X-ray obserations or an association of the giant pulses with the magnetic pole opposite the high energy emission remains to be tested with independent phase comparisons.

Given the exponential broadening of the giant pulses at 1.38 GHz, it seems likely that the integrated profile at 0.64 GHz is dominated by scattering (M93). Scatter broadening scales with the inverse fourth power of the frequency and hence should be $``$20 times larger at 0.64 GHz than at 1.38 GHz. The 0.64 GHz profile is thus consistent with a double peaked profile separated by $``$0.3 phase and broadened by $``$10 ms of scattering. It is tempting to relate the location of this double peaked structure with the location of the giant pulses detected at the higher frequency. Further observations at 0.64 GHz are needed to confirm this. Furthermore, the scattering at high frequency would also explain the non-detection of the pulsar at 0.44 GHz reported by M93, where the scattering time would exceed the pulse period. These values also indicate that the effects of scintillation are not important in the giant pulse search toward PSR B0540$``$69, contrary to the suggestion of McLaughlin & Cordes (2003).

The DM is high compared to the other known pulsars in the Large Magellanic Cloud which have DMs in the range 67 – 107 cm-3pc (McConnell et al. 1991). A scattering time of $``$0.4 ms at 1.38 GHz is also large, especially when the scattering screen is likely to be located much closer to the pulsar (i.e. in the LMC) than to the Earth. The situation is analogous to that seen in the Vela pulsar, where both the DM and the scattering are much larger than expected for its distance. Both the excess DM and the scattering likely arise inside the Vela SNR itself where the line of sight intersects a dense clump of ionised material.

Expressed as a continuum flux density, the largest giant pulses have 1.4 GHz flux densities at a distance of 1 kpc of 4000, 200, 7 and 10 Jy for the Crab, PSRs PSR B0540$``$69, 1821–24 and 1937+21 respectively for roughly the same number of rotations of the pulsar ($``$106). The giant pulses from the two millisecond pulsars have very similar flux densities and these are significantly less than those of the young pulsars. From the point of view of searching for giant pulses, the peak flux is the important value. The giant pulses in the Crab pulsar and in PSR B1937+21 are known to have widths of order microseconds. Assuming all the giant pulses have durations of $``$$`\mu `$s, then the giant pulses have peak fluxes at 1.4 GHz at a distance of 1 kpc of 130, 11, 0.03 and 0.02 MJy for the Crab, PSRs PSR B0540$``$69, 1821–24 and 1937+21 which makes the differences between the young pulsars and the millisecond pulsars even more striking. In principle, therefore, the Parkes telescope is capable of detecting giant pulses from young pulsars at distances of $``$1 Mpc, and millisecond pulsars out to $``$50 kpc. In practise, however, lack of knowledge of the DM and scatter broadening both serve to weaken these limits. The prospects for future generation radio telescopes are excellent however. Even if, in young pulsars, the giant pulses only occur over the first $``$2000 yr, one might expect to see $``$20 per large galaxy and more in the nearby starburst galaxies. For millisecond pulsars, the faster objects which have high magnetic fields at the light cylinder, will produce giant pulses for a Hubble time. There therefore may be a large population of millisecond pulsars which produce giant pulses, enhancing the prospects for targetted searches of globular cluster in external galaxies. Detection of extra-galactic pulsars will yield the electron density content of the inter-galactic medium, important in understanding the nature of the dark matter (Maloney & Bland-Hawthorn 2001).

## 7 Summary

We have detected giant pulses from the young pulsar PSR B0540$``$69 in the LMC and yet have failed to detect the integrated flux from the pulsar after 8 hr observing. Four pulsars, including PSR B0540$``$69, are now known to have giant pulses and these four are in the top five when ranked by their magnetic field at the light cylinder. In the near future we plan to observe PSR B0540$``$69 at both lower and higher radio frequencies to determine the spectral index of the giant pulse emission and also to replicate the detection of the pulsar at 0.64 GHz by M93. We also plan a simultaneous radio and X-ray observation to both obtain accurate phasing between the radio and X-ray and to look for giant pulses in the X-ray.

The Australia Telescope is funded by the Commonwealth of Australia for operation as a National Facility managed by the CSIRO. We thank F. Marshall for providing sample RXTE pulse arrival times. We thank A. Lyne and D. Nice for providing timing ephemerides for the Crab and PSR B1937+21.


text:
# References

EFI 03-19

hep-ph/0305144

May 2003

A Comment on “Weak Phase $`\gamma `$ using Isospin Analysis and

Time-dependent Asymmetry in $`B^0K_S\pi ^+\pi ^{}`$111To be published as a Comment in Physical Review Letters

Michael Gronau222Permanent Address: Physics Department, Technion – Israel Institute of Technology, 32000 Haifa, Israel.

Enrico Fermi Institute and Department of Physics

University of Chicago, Chicago, Illinois 60637

> We comment on a recent suggestion for measuring $`\gamma `$ in $`B^0(t)K_S\pi ^+\pi ^{}`$. A difficulty is pointed out in relating electroweak penguin and tree amplitudes in $`B^+K^0\pi ^+\pi ^0`$, which is crucial for an implementation of this method. We clarify the necessary condition for such a relation.

PACS codes: 12.15.Hh, 12.15.Ji, 13.25.Hw, 14.40.Nd

In a recent interesting Letter  Nilendra Deshpande, Nita Sinha and Rahul Sinha propose to use $`BK\pi \pi `$ decays in order to determine the weak phase $`\gamma `$. Their method is based on the CP asymmetry in $`B^0(t)K_S\pi ^+\pi ^{}`$, and on an isospin triangle relation among the three amplitudes for $`B^+K^0(\pi ^+\pi ^0)_e,B^0K^0(\pi ^+\pi ^{})_e`$ and $`B^0K^0(\pi ^0\pi ^0)_e`$, in which the two pions are in an even angular momentum state. A crucial assumption of this method is that electroweak penguin and tree amplitudes contributing to $`B^+K^0(\pi ^+\pi ^0)_e`$ involve a common strong phase. Such a property was shown to hold in the SU(3) symmetry limit for the $`I=3/2`$ amplitude in $`BK\pi `$ , and in the isospin symmetry limit for the $`I=2B\pi \pi `$ amplitude .

The purpose of this short comment is to clarify the general condition under which tree and electroweak amplitudes in charmless $`B`$ decays can be related to each other. We will show that this condition is not fulfilled in the case studied in .

The effective Hamiltonian describing charmless $`\mathrm{\Delta }S=1`$ (or $`\mathrm{\Delta }S=0`$) decays  consists of current-current operators $`Q_1`$ and $`Q_2`$, QCD penguin operators $`Q_i,i=36`$, and electroweak penguin (EWP) operators $`Q_i,i=710`$. The operators $`Q_1`$ and $`Q_2`$, multiplying Wilson coefficients $`c_1`$ and $`c_2`$, respectively, and CKM coefficients $`V_{ub}^{}V_{us}`$ (or $`V_{ub}^{}V_{ud}`$), will be named tree operators. EWP operators involve CKM factors $`V_{tb}^{}V_{ts}`$ (or $`V_{tb}^{}V_{td}`$). The EWP operators $`Q_9`$ and $`Q_{10}`$ with the dominant Wilson coefficients, $`c_9`$ and $`c_{10}`$, have the same (V–A)(V–A) structure as the tree operators, and would have approximately the same matrix elements if they had also identical flavor SU(3) and isospin structure.

In order to find out when such a relation holds, one decomposes all four quark operators of the form $`(\overline{b}q_1)(\overline{q}_2q_3)`$ transforming as $`\overline{\mathrm{𝟑}}\mathrm{𝟑}\overline{\mathrm{𝟑}}`$ in terms of a sum of $`\overline{\mathrm{𝟏𝟓}}`$, $`\mathrm{𝟔}`$ and $`\overline{\mathrm{𝟑}}`$ . The representation $`\overline{\mathrm{𝟑}}`$ appears both symmetric ($`\overline{\mathrm{𝟑}}^{(s)}`$), and antisymmetric ($`\overline{\mathrm{𝟑}}^{(a)}`$) under the interchange of $`q_1`$ and $`q_3`$. The tree and electroweak parts of the $`\mathrm{\Delta }S=1`$ Hamiltonian are :

$$_T=\frac{G_F}{\sqrt{2}}V_{ub}^{}V_{us}\left[\frac{c_1c_2}{2}(\overline{\mathrm{𝟑}}_0^{(a)}+\mathrm{𝟔}_1)+\frac{c_1+c_2}{2}(\overline{\mathrm{𝟏𝟓}}_1+\frac{1}{\sqrt{2}}\overline{\mathrm{𝟏𝟓}}_0\frac{1}{\sqrt{2}}\overline{\mathrm{𝟑}}_0^{(s)})\right],$$
(1)
$$_{\mathrm{EWP}}=\frac{G_F}{\sqrt{2}}\frac{3V_{tb}^{}V_{ts}}{2}\left[\frac{c_9c_{10}}{2}(\frac{1}{3}\overline{\mathrm{𝟑}}_0^{(a)}+\mathrm{𝟔}_1)+\frac{c_9+c_{10}}{2}(\overline{\mathrm{𝟏𝟓}}_1\frac{1}{\sqrt{2}}\overline{\mathrm{𝟏𝟓}}_0\frac{1}{3\sqrt{2}}\overline{\mathrm{𝟑}}_0^{(s)})\right],$$
(2)
where subscripts denote the isospin of corresponding operators. Note that both the $`\mathrm{𝟔}`$ and $`\overline{\mathrm{𝟏𝟓}}`$ operators include a $`\mathrm{\Delta }I=1`$ component.

Eqs. (1) and (2) imply proportionality relations between EWP and tree operators transforming as $`\overline{\mathrm{𝟏𝟓}}`$ and $`\mathrm{𝟔}`$ :

$`_{\mathrm{EWP}}(\overline{\mathrm{𝟏𝟓}})`$ $`=`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{c_9+c_{10}}{c_1+c_2}}{\displaystyle \frac{V_{tb}^{}V_{ts}}{V_{ub}^{}V_{us}}}_T(\overline{\mathrm{𝟏𝟓}}),`$ (3)
$`_{\mathrm{EWP}}(\mathrm{𝟔})`$ $`=`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{c_9c_{10}}{c_1c_2}}{\displaystyle \frac{V_{tb}^{}V_{ts}}{V_{ub}^{}V_{us}}}_T(\mathrm{𝟔}).`$ (4)
The two proportionality constants are approximately renormalization scale independent, and are approximately equal in magnitudes but differ in sign ,

$$\frac{c_9+c_{10}}{c_1+c_2}\frac{c_9c_{10}}{c_1c_2}.$$
(5)
Therefore, EWP and tree amplitudes in $`B`$ decay processes which obtain contributions from either the $`\overline{\mathrm{𝟏𝟓}}`$ or the $`\mathrm{𝟔}`$ operator, but not from both, are proportional to each other and involve a common strong phase. This property does not hold when the two operators contribute because of the opposite signs in Eqs. (3) and (4). In this case EWP and tree amplitudes involve the sum and difference of two complex amplitudes and have different strong phases.

In the case of $`B(K\pi )_{I=3/2}`$ , the $`K`$ and $`\pi `$ are in an S-wave state, which is symmetric under an interchange of the two SU(3) octets. This state is a pure $`\mathrm{𝟐𝟕}`$. The only SU(3) operator which contributes to this transition is the $`\overline{\mathrm{𝟏𝟓}}`$. Consequently, the EWP and tree amplitudes are proportional to each other in the SU(3) approximation. The same holds true in the isospin symmetry limit for the EWP and tree amplitudes of $`B(\pi \pi )_{I=2}`$, since only the $`\overline{\mathrm{𝟏𝟓}}`$ contains a $`\mathrm{\Delta }I=3/2`$ component . On the other hand, in $`B^+K^0(\pi ^+\pi ^0)_e`$ studied in  the final state has $`I=3/2,S=1`$ and can be in a $`\mathrm{𝟐𝟕}`$ and in a $`\mathrm{𝟏𝟎}`$, to which the $`\mathrm{\Delta }I=1`$ components of both the $`\overline{\mathrm{𝟏𝟓}}`$ and the $`\mathrm{𝟔}`$ operators contribute. Hence, the condition for proportional EWP and tree amplitudes and for a common strong phase does not hold. Although this proportionality does not follow from symmetry considerations alone, it would be interesting to study possible dynamical assumptions which can lead to such a situation.

I thank Jonathan Rosner for useful discussions, and I am grateful to the Enrico Fermi Institute at the University of Chicago for its kind hospitality. This work was supported in part by the United States Department of Energy through Grant No. DE FG02 90ER40560.


text:
# Chapter 1 An Estimate Of The Primary Mass Of Cosmic Rays At 10¹⁸ eV As Inferred From Volcano Ranch Data.

## Chapter 1 An Estimate Of The Primary Mass Of Cosmic Rays At $`10^{18}`$ eV As Inferred From Volcano Ranch Data.

### Abstract

Accurate measurements of the lateral distribution of extensive air showers produced by cosmic rays of energy greater than $`10^{17}`$ eV were made in the 1970s by Linsley. At the time, the state of knowledge about the best hadronic interaction model to use to describe such data prevented conclusions about the cosmic ray mass from being drawn. We have used a modern model, qgsjet98, to infer the primary mass from these data, using the very careful records left by Linsley. We find that at a median energy of $`10^{18}`$ eV , the data are well described by an iron-dominated composition (88 $`\pm `$ 6(stat) $`\pm `$ 18 (syst))%. We discuss the systematic errors in this estimate that arise from model uncertainties and from the range of energies used in the work of Linsley. These data are used with the permission of the late John Linsley to whom this work is dedicated.

### 1. Introduction

Determining the mass composition above $`10^{17}`$ eV is an important and challenging measurement and has been tackled by several groups using different techniques. These results have not all been consistent in their conclusions. It is important to solve this problem because of its implications for cosmic ray models of origin, acceleration and propagation. We present a first interpretation of the precise Volcano Ranch measurements of the lateral distribution function using Monte Carlo tools that were not available when the data were recorded in the 1970s. The pioneering Volcano Ranch instrument was described by Linsley, in his various writings , with unusually detailed descriptions of his equipment, together with many shower details and a description of his data reduction methods.

### 2. Derivation of the primary mass composition

We have generated 1735 proton and iron showers with $`\mathrm{sec}\theta =1.01.5`$, and primary energy between $`10^{17}`$ eV and $`10^{19}`$ eV using the aires code (version 2.4.0), with the hadronic interaction generator qgsjet98. The results of the simulated showers were convolved with a simulation of the detector response made using geant4 to obtain the scintillator response . As a second step, each simulated shower was thrown up to 100 times on to a simulated Volcano Ranch array with random core positions in the range 0-150 m from the array center. A comparison between lateral distribution measurements  and average $`10^{19.1}`$ eV iron showers including the scintillator response of the detectors in the Volcano Ranch array configuration is presented in Fig 1 (left). The agreement of the lateral distributions is good and gives confidence in the procedures used.

A generalized version of the Nishimura-Kamata-Greisen (NKG) formula was used by Linsley to describe the lateral distribution :

$$S_{\mathrm{VR}}(r)=\frac{N}{r_m^2}\frac{\mathrm{\Gamma }(\eta \alpha )}{2\pi \mathrm{\Gamma }(2\alpha )\mathrm{\Gamma }(\eta 2)}\left(\frac{r}{r_m}\right)^\alpha \left(1+\frac{r}{r_m}\right)^{(\eta \alpha )}$$

From a subset of 366 showers detected with the array, the form of $`\eta `$ as a function of zenith angle $`\theta `$ and shower size $`N`$ was found to be  with $`\alpha =1`$:

$$\eta (\theta ,N)=a+b(\mathrm{sec}\theta 1)+clog_{10}(\frac{N}{10^8})$$
(1)
with $`a=3.88\pm 0.054`$, $`b=0.64\pm 0.07`$, and $`c=0.07\pm 0.03`$.
We estimate the primary mass by a maximum likelihood fit of the best linear combination of pure iron and proton generated samples to match the data sample assuming a bi-modal composition. The probability of observing a particular number of events $`d_i`$ in a particular bin is given by $`\mathrm{exp}^{f_i}f_i^{d_i}/d_i!`$ where $`f_i`$ is the predicted value for the number of events in this particular bin.
If we assume a bi-modal composition of proton and iron with proportions $`P_{Fe}`$ and $`P_p`$ then, $`f_i=C(P_{Fe}+P_p)`$ and $`C`$ is the overall normalization factor between numbers of data and Monte Carlo events. The estimates of the proportions $`P_j`$ are found by maximizing the total likelihood or its logarithm, $`\mathrm{ln}(\mathrm{\pounds })=\mathrm{\Sigma }d_i\mathrm{ln}(f_i)\mathrm{ln}(d_i!)f_i`$.

In Fig 2 we show the Monte Carlo pure iron sample (left) and pure proton sample (right) with the corresponding Volcano Ranch data points. As can be seen from Fig 2 (right), the tail at large $`\eta `$, in the comparison with iron, indicates that a light component must be included to fit the experimental data. The best fit gives a mixture with (88 $`\pm `$ 6)% of iron and a corresponding amount of protons. The resulting $`\eta `$ distribution is shown in Fig 2 (bottom) together with Volcano Ranch data. The systematic error arising from our lack of knowledge of the energy distribution of the events is estimated by repeating the fitting procedure with different energy spectra. From these analyses we estimate a systematic error of 12%. The systematic error introduced by the choice of hadronic model is estimated to be 14% .

### 3.  Discussion and Conclusions

In Fig 1 (right) we present the fraction of iron for Fly’s Eye, Agasa and Haverah Park using different hadronic interaction models. Also shown is the mean composition with the corresponding error we get for the Volcano Ranch energy range using qgsjet98 and an estimation of what we would get using qgsjet01 as it was done in Haverah Park . Cosmic rays at Volcano Ranch are found to be compatible with mean fraction (88 $`\pm `$ 6(stat) $`\pm `$ 18 (syst))% of iron in a bi-modal proton and iron mix, in the energy range $`10^{17.7}`$ eV to $`10^{19}`$ eV with mean energy $`10^{18}`$ eV.

The inconsistency between measurements of mass composition needs to be addressed in order to make a more solid conclusion on the origin, acceleration or propagation of cosmic rays. This work attempts to help answer these important questions by re-interpreting the pioneering Volcano Ranch experiment, but in doing so has in some ways deepened the mystery of the origin of these particles. Qualitatively there are significant differences between the estimates of composition at energies above $`10^{17}`$ eV. So the questions appear far from being resolved.

### 4. References

1. M. T. Dova et al, Proceedings of the XII ISVHECRI, CERN, 15-20 July 2002.

2. J. Linsley, Proc. 13th ICRC, Denver, (1973) 3212.

3. J. Linsley, Proc. 15th ICRC, Plovdiv,(1977) 56-62-89.

4. B. R. Dawson et al, Astropart. Phys. 9 (1999) 331.

5. M. Ave et al., Astropart. Phys. 19 (2003) 61.