text:
# Comparison of structural transformations and superconductivity in compressed Sulfur and Selenium

## I INTRODUCTION

The chalcogens exist in a rich variety of crystal structures at ambient pressures and upon compression. At low pressures, S, Se, and Te crystallize in different sequences of insulating two-fold-coordinated molecular or polymeric structures. With enough compression, however, all three elements transform into the same body-centered orthorhombic (bco) crystal structure. The bco structure is comprised of puckered layers with each site having a coordination number of four. With further compression, the more three-dimensional six-fold coordinated rhombohedral $`\beta `$-Po structure becomes stable in all three systems. At even higher pressures, Se and Te are observed to transform from the $`\beta `$-Po structure to the bcc structure, continuing the trend of increasing coordination and closer packing with pressure. Though is has been speculated that S should also undergo a $`\beta `$-Po $``$ bcc transition, recent first-principles calculations predict $`\beta `$-Po S to become unstable first to the relatively open simple-cubic (sc) structure. The calculations find simple-cubic S to be favored over a large range of pressures, with the bcc phase eventually becoming stable at very high pressures.

Electrical measurements show that S, Se, and Te are all metallic in the bco structure. Further, the metallic bco phases of S, Se, and Te are observed to be superconducting with transition temperatures of approximately 10, 5, and 3 K, respectively. Recent magnetic measurements show that the superconducting transition temperature in S jumps abruptly to 17 K upon transformation to the $`\beta `$-Po structural phase. This is among the highest transition temperatures observed in elemental solids. First-principles calculations indicate that this high transition temperature results from a moderate electron-phonon coupling combined with a large average phonon frequency, which sets the scale for the transition temperature. Earlier calculations had predicted the hypothetical bcc phase of S to be superconducting with a similar transition temperature of about 15 K near 550 GPa. As yet, no measurements of $`T_C`$ have been reported on the higher-pressure phases of the heavier chalcogens.

In this paper, we report density-functional calculations of the structural and superconducting transitions in S and Se in the $`\beta `$-Po and higher-pressure phases. Though the two materials are chemically similar, they display striking differences in their properties under pressure. In addition to the simple-cubic structure’s viability as a high-pressure phase for S but not Se, the stability of the fcc structure is another difference between the materials. An analysis of the contributions to the total energy shows that these structural differences can be understood in terms of differences between the S and Se cores. Although some of the trends in the calculated electron-phonon coupling and superconducting properties are similar in the two materials, the superconducting transition temperature is higher in S than in Se, even within the same structural phase. A comparison of the phonon spectra and electron-phonon coupling in these materials helps to explain the differences.

The remainder of this paper is organized as follows. In Sec. II, we compare and contrast the calculated high-pressure structural phase diagrams for S and Se. In Sec. III, the calculated vibrational spectra and electron-phonon coupling parameters for the high-pressure phases are presented and discussed. We also discuss ways in which our results can be tested and compared with diamond-anvil-cell-based experiments, including measurements of superconducting properties, transport coefficients, and optical conductivity. Concluding remarks are given in Sec. IV.

## II HIGH-PRESSURE STRUCTURAL PHASES

### A Method

Zero-temperature structural energetics are calculated using the plane-wave pseudopotential method within the local density approximation. The pseudopotentials are generated using the Troullier-Martins method. Kinetic-energy cutoffs of 70 Ry and 40 Ry are used for the plane-wave expansion of the Kohn-Sham orbitals in S and Se, respectively. The correlation potential is approximated using the Perdew-Zunger parameterization of the electron-gas results, and the nonlinearity of the exchange and correlation interaction between the core and valence charge densities is treated using a partial core. The Brillouin zone for each structure is sampled on Monkhorst-Pack meshes of at least $`20^3`$ points.

### B Results and discussion

The high-pressure crystal structures considered here for the chalcogens can all be described in terms of a single-atom rhombohedral unit cell with rhombohedral angles $`\alpha _r`$ of 60 (fcc), 90 (sc), approximately 104 ($`\beta `$-Po), and 109.47 (bcc). Figure 1 shows the total energy for S and Se calculated as a function of the rhombohedral angle at constant unit-cell volumes. The curves for the two materials are strikingly different, leading to different sequences of high-pressure phases. In particular, for the range of volumes shown in Fig. 1, the simple-cubic structure is never competitive in energy with the other structures considered for Se, while the fcc structure is never favorable for S.

For the largest unit-cell volume shown, corresponding to roughly 80 GPa in Se, the $`\beta `$-Po structure with $`\alpha _r104^{}`$ is lowest in energy. At a smaller volume of 76 a.u. (roughly 140 GPa), the local minimum for $`\beta `$-Po Se has disappeared, and the bcc structure is favored. At even smaller volumes, the minimum at $`\alpha _r=60^{}`$ drops below that at 109.5, indicating a bcc to fcc transition. Comparison of the enthalpies $`H=E_{\mathrm{total}}+pV`$ calculated for the different structures yields a transition pressure of 120 GPa for the $`\beta `$-Po$``$bcc transition, which is in reasonable agreement with the experimental value of 140 GPa. For the bcc$``$fcc transition, we predict a transition pressure of about 260 GPa.

For S, the $`\beta `$-Po structure, which is favored at a unit-cell volume of 56 a.u., is calculated to become unstable under compression with respect to the simple-cubic structure, and the simple-cubic structure subsequently becomes unstable to the bcc structure (not shown in Fig. 1). Using enthalpies calculated for the static lattices, we find transition pressures of 260 and 540 GPa for these structural transformations. Since we find that the structural phases of S differ significantly in their average phonon frequencies (see Sec. III), it is important to take into account zero-point contributions to the energies. Inclusion of zero-point energies shifts the estimated pressures in S to about 280 GPa for the $`\beta `$-Po$``$sc transition and about 500 GPa for the sc$``$bcc transition.

In terms of coordination number and packing, Se follows the expected trend of increasing coordination with pressure, going from six ($`\beta `$-Po) to eight (bcc) to twelve (fcc). With each transformation, Se adopts a more closely packed structure. Sulfur, on the other hand, transforms from the six-fold coordinated $`\beta `$-Po structure to the sc structure, which is also six-fold coordinated, but is more open. (The ratio of next-nearest to nearest-neighbor distances is 1.24 and 1.41 in the $`\beta `$-Po and sc structures, respectively.) While the more highly coordinated bcc structure is adopted at large compressions, the even more densely packed fcc structure is never favored. The key to these differences in the high-pressure phase diagrams of these isovalent elements lies in their cores.

Figure 2 compares the relative stability of the sc and bcc structures of S and Se as a function of volume. To facilitate the comparison, the volumes are measured with respect to the calculated equilibrium bcc volumes of $`V_0^{bcc}=98.6`$ and 125.9 a.u./atom for S and Se, respectively. The total energy can be divided into the one-electron contribution, which arises from the noninteracting kinetic energy and the electron–ion interaction, the electron–electron Coulomb energy, the exchange and correlation energies, and the ion–ion energy, which includes the Ewald energy and a term that accounts for the difference between the pseudopotential and the pure Coulomb potential of the ions. In both materials, the repulsive Hartree and ion–ion terms favor the more closely packed bcc structure, which has a more uniform distribution of both ionic and electronic charge. The attractive electron-ion Coulomb interaction favors the more open sc structure which tends to have very nonuniform charge distributions. The difference in total energy between the two structures is smaller than differences in the individual contributions to the total energy. Figure 2 shows that the one-electron contribution in S favors the sc structure significantly more strongly than it does in Se, thereby tipping the balance to stabilize the open sc structure. Further decomposition of the one-electron contribution shows that the kinetic energy favors the more uniform bcc structure. It is thus the electron-ion interaction that stabilizes the sc S phase.

The importance of the electron-ion interaction in producing different stable structures in S and Se can be understood in terms of the difference in their cores. In particular, since the Ne core of S contains only $`s`$ and $`p`$ electrons, the 3$`d`$ states in S have no orthogonality constraint with the core, resulting in a strongly attractive $`d`$ pseudopotential in the core region. On the other hand, in Se, which has $`d`$ states in the core, the repulsive Pauli component of the pseudopotential largely cancels the attractive Coulomb component in the core region, resulting in a relatively weak $`d`$ potential.

Under ambient conditions, the valence $`d`$ states are unoccupied in both S and Se, but with compression, the $`d`$ bands broaden and eventually cross the Fermi level. The electronic density of states for bcc and sc S are shown in Fig. 3. The volume of 50 a.u./atom corresponds to about 250 GPa, which is close to the pressure at which the sc phase is calculated to become stable. The thin lines show the $`d`$-projected density of states. Since the $`d`$ bandwidth varies roughly as $`C/d_{\mathrm{nn}}^5`$, where $`C`$ is the coordination number and $`d_{\mathrm{nn}}`$ is the nearest neighbor distance, the shorter bonds in the sc structure more than compensate for the lower coordination, leading to broader $`d`$ bands in the sc phase as compared to the bcc phase. With more $`s`$-$`d`$ transfer of electrons in the sc phase than in the bcc (or fcc) structures, the deep S $`d`$ pseudopotential becomes important in the energetics, stabilizing the sc phase. In Se, the sc phase similarly has larger $`d`$ occupation than the fcc and bcc phases, but because the $`d`$ potential is weak, this $`s`$-$`d`$ transfer does not lower the one-electron energy enough to stabilize the sc structure.

Similar reasoning helps explain the absence of the close-packed fcc structure in high-pressure S. Fig. 4 compares the relative stability of the fcc and bcc phases of S and Se as a function of volume. Again the different structural preferences of S and Se can be attributed to differences in the one-electron contributions to the total energy. With increasing compression, the one-electron energy in S increasingly favors the bcc structure, while in Se it increasingly favors the fcc structure. In both materials, the kinetic energy favors the more uniformly distributed fcc structure while the electron–ion term favors the bcc structure. With the deeper $`d`$ pseudopotential in S, the importance of the electron-ion term in S is enhanced, stabilizing the bcc structure. The difference in openness between the fcc and bcc structures is of course much smaller than that between the sc and bcc structures, so in comparisons of the fcc and bcc structures, the energy differences are smaller and the energy balance is more subtle.

This discussion of the differences between S and Se and the important role of the presence of $`d`$ states in the core suggests that the heavier chalcogen Te should be similar to Se in terms of high-pressure structural phases. Indeed, like Se, Te is observed to transform from the $`\beta `$-Po structure to the bcc structure with pressure. In our calculations, as well as earlier ones by Kirchhoff et al., this transition in Te is found close to the experimentally measured transition pressure of 27 GPa. In addition, we find that the close-packed fcc structure becomes stable above about 80 GPa. Finally, as in Se, the sc structure is calculated to be energetically unfavorable as a high-pressure phase of Te.

## III PHONON SPECTRA, ELECTRON-PHONON COUPLING, AND SUPERCONDUCTIVITY

### A Method

The density-functional linear response method is used to compute the phonon spectra for each structure at several different pressures. The dynamical matrices are computed using the scheme described in Ref. , where the self-consistent change in the Hamiltonian caused by ionic displacements is obtained by solving a Bethe-Saltpeter equation for the change in the charge density. The dynamical matrix is computed on a grid of $`6^3`$ phonon wavevectors $`𝐪`$ for $`\beta `$-Po S and $`8^3`$ wavevectors for sc and bcc S. For Se, $`8^3`$, $`10^3`$, and $`12^3`$ $`𝐪`$ points are sampled in the full $`\beta `$-Po, bcc, and fcc Brillouin zones, respectively.

The matrix element for scattering of an electron from state $`n𝐤`$ to state $`n^{}𝐤^{}`$ by a phonon with frequency $`\omega _{𝐪\nu }`$ and eigenvector $`\widehat{ϵ}_{𝐪\nu }`$ is given by

$$g(n𝐤,n^{}𝐤^{},𝐪\nu )=\sqrt{\frac{\mathrm{}}{2M\omega _{𝐪\nu }}}n𝐤|\widehat{ϵ}_{𝐪\nu }_𝐑V_{\mathrm{SCF}}|n^{}𝐤^{},$$
(1)
where $`_𝐑V_{\mathrm{SCF}}`$ is the gradient of the self-consistent potential with respect to atomic displacements. This scattering gives rise to a finite phonon linewidth,

$$\gamma _{𝐪\nu }=2\pi \omega _{𝐪\nu }[N(E_F)]^2|g_{𝐪\nu }|^2,$$
(2)
where $`N(E_F)`$ is the electronic density of states per spin at the Fermi level. The double brackets $`\mathrm{}`$ denote a doubly constrained Fermi surface average as defined in Ref. . This scattering process also contributes to the effective mass of the electrons via the mass enhancement parameter $`\lambda `$, obtained from the wavevector- and branch-dependent contributions

$$\lambda _{𝐪\nu }=\frac{\gamma _{𝐪\nu }}{\mathrm{}\pi N(E_F)\omega _{𝐪\nu }^2}$$
(3)
by summing over branches $`\nu `$ and averaging over wavevectors $`𝐪`$.

The electron-phonon spectral function, which measures the effectiveness of phonons of a given energy to scatter electrons on the Fermi surface, plays a central role in the Eliashberg strong-coupling theory of superconductivity. The spectral function is given by

$$\alpha ^2F(\omega )=\frac{1}{2\pi N(E_F)}\underset{𝐪\nu }{}\delta (\omega \omega _{𝐪\nu })\frac{\gamma _{𝐪\nu }}{\mathrm{}\omega _{𝐪\nu }}.$$
(4)
Within this framework, the mass enhancement parameter $`\lambda `$ is proportional to the inverse-frequency moment of the spectral function.

In the density-functional linear-response method, $`_𝐑V_{\mathrm{SCF}}`$ is computed in the process of determining the dynamical matrices. Therefore we calculate the electron-phonon matrix elements on the same grid of wavevectors used in the phonon calculations. For calculations of the phonon density of states and electron-phonon spectral function, which involve summing over q points throughout the Brillouin zone, a Fourier interpolation procedure is used to obtain the dynamical matrix and its dissipative part on a denser mesh of q points. The doubly constrained Fermi surface average of $`g`$ is computed using Bloch functions on dense meshes of at least $`30^3`$ k and k points in the full Brillouin zone, with delta functions at the Fermi level replaced by Gaussians of width of order 0.01 Ry chosen to reproduce the value of $`N(E_F)`$ obtained using the linear tetrahedron method.

### B Results and discussion

Tables I and II list parameters that characterize the phonon spectra and electron-phonon couplings in the high-pressure phases of S and Se. Included are the average phonon frequency $`\omega `$ and the logarithmic $`\lambda `$-weighted average phonon frequency $`\omega _{\mathrm{ln}}=\mathrm{exp}[_{𝐪\nu }\lambda _{𝐪\nu }\mathrm{ln}(\omega _{𝐪\nu })/_{𝐪\nu }\lambda _{𝐪\nu }]`$. Some of the values for $`N(E_F)`$ and $`\lambda `$ listed in Table I differ by a few percent from those published in our earlier paper on S because we use the more accurate tetrahedron method for calculating the electronic density of states in the present work.

Within each structural phase, compression stiffens the lattice and weakens the electron-phonon interaction, as indicated by the rise in $`\omega `$ and fall in $`\lambda `$. The pressure-dependence of $`\lambda `$ varies from structure to structure. In sc S, $`\lambda `$ varies approximately as $`1/\omega ^2`$, as might be expected from Eq. (3). In $`\beta `$-Po S, $`\beta `$-Po Se, and bcc Se, however, the pressure-dependence of $`\lambda `$ is much weaker. These differences can be understood by examining contributions to $`\lambda `$ from different phonon modes.

In some of the structural phases, it is possible to identify specific phonon modes that contribute strongly to $`\lambda `$. In the $`\beta `$-Po phase of Se, for example, there are strong anomalies in the phonon dispersion curves along the $`\mathrm{\Gamma }`$ to F direction (i.e., the $`(1\overline{1}0)`$ direction in terms of the reciprocal lattice vectors $`𝐛_i`$), as shown in Fig. 5. Similar but less pronounced anomalies exist in the $`\beta `$-Po phase of S. These anomalies are associated with nesting of the Fermi surface by the wavevector $`𝐪0.7(𝐛_1𝐛_2)`$, as confirmed by calculations of the geometric nesting factor $`\delta (𝐤𝐤^{}𝐪)`$. The strong Fermi surface nesting affects not only the phonon frequencies, but the linewidths as well. The combination of low frequencies and enhanced linewidths results in a large contribution to $`\lambda `$ from this part of the $`\beta `$-Po Brillouin zone. Anomalies are also present in the phonon spectra of bcc S and bcc Se. Figure 6 shows the dispersion curves for bcc Se at 120 GPa. While not correlated with nesting vectors, the low-frequency anomalies in the transverse branches along the $`\mathrm{\Gamma }`$ to N and $`\mathrm{\Gamma }`$ to H lines nevertheless give large contributions to $`\lambda `$.

The soft phonon anomalies in the bcc and $`\beta `$-Po phases are more weakly pressure-dependent than the overall phonon spectrum. The pressure dependence of $`\lambda `$ in these phases is strongly influenced by these persistent anomalies and is therefore weaker than the $`\omega ^2`$ dependence suggested by Eq. (3). In sc S, the electron-phonon coupling is not found to be dominated by any particular phonon mode, so the pressure dependence of $`\lambda `$ is to a good approximation determined by that of the average phonon frequency.

As can be seen in Tables I and II, transformations from structure to structure are usually accompanied by abrupt changes in the phonon and electron-phonon parameters. In comparing the various structural phases of S, the simple-cubic phase stands out in having the weakest electron-phonon coupling parameters. The spectral function for sc S, plotted in Fig.7, is strikingly different from those of bcc and $`\beta `$-Po S. Both the average and maximum phonon frequencies in the sc phase are much larger than those in the other phases. The higher frequencies in the sc phase are due to the shorter nearest-neighbor distances and the strong concentration of charge along the bonds in the open sc lattice. With these stiff bonds, the $`(1\overline{1}0)`$-polarized transverse branch along the (110) direction and both of the transverse branches along the (111) direction are nearly as high in energy as the longitudinal branches. This leads to a phonon density of states and a spectral function that are dominated by a high-frequency peak. The large phonon energy scale and the lack of low-energy modes that couple strongly to electrons, as evidenced by the relatively large value of $`\omega _{\mathrm{ln}}`$, combine to drive $`\lambda `$ down in the sc phase of S.

Among the high-pressure phases of Se, the bcc phase stands out in that $`\lambda `$ is largest in this phase. Although the average phonon frequencies are similar to those in the neighboring phases near the transition pressures, $`\omega _{\mathrm{ln}}`$ is significantly depressed in the bcc phase. This is a sign of the importance of the low-frequency anomalies already identified in the bcc Se phonon spectra. These anomalies enhance the spectral function at low energies, as shown in Fig.8, and they have a significant effect on $`\lambda `$ because they appear at low frequencies. Although $`\beta `$-Po Se also has soft phonon anomalies, there are fewer of them due to the lower symmetry of the lattice, and they have a smaller effect on $`\lambda `$. No comparable soft phonon anomalies are found in the fcc phase.

In comparing analogous structural phases of S and Se, the differences in electron-phonon coupling strength reflect not only differences in the phonon spectra, but also differences in the electron-phonon matrix elements. The stronger potential in S gives rise to a larger change in the bare potential $`_𝐑V_{\mathrm{bare}}`$ when atoms are displaced. Furthermore, because the $`d`$ electrons are more tightly bound in S, these electrons are less effective in screening $`_𝐑V_{\mathrm{bare}}`$. In combination, these two effects result in significantly larger matrix elements of $`_𝐑V_{\mathrm{SCF}}`$ in S than in Se. Combined with the smaller mass in S, this leads to larger linewidths. This is evident in both the $`\beta `$-Po and bcc phases, where the phonon frequencies in S are a factor of two or more higher than those in Se, yet the mass enhancement parameters in the two materials are close.

Our results can be compared to earlier density-functional-based studies of the electron-phonon interaction in $`\beta `$-Po and bcc Te and bcc S. Soft phonon anomalies very similar to those discussed here have been identified to be important in the electron-phonon coupling in $`\beta `$-Po and bcc Te. In addition, the electron-phonon mass enhancement parameter has been calculated to undergo a large jump at the $`\beta `$-Po to bcc transition in Te, just as we find for Se. For bcc S, our value of $`\lambda =0.70`$ at 550 GPa is somewhat larger than the value of 0.58 reported in Ref.  for a slightly higher pressure of 584 GPa. Since values of $`\lambda _{𝐪\nu }`$ are in reasonable agreement along high-symmetry directions, it is likely that the difference arises from a difference in sampling of wavevectors: In the present work, a uniform grid of points throughout the Brillouin zone is sampled, while in Ref. , $`\lambda _{𝐪\nu }`$ is calculated along a few high-symmetry directions and spherically averaged to estimate $`\lambda `$.

### C Making contact with experiments

The need to carry out experiments inside high-pressure cells limits the types of experimental probes available to investigate the electron-phonon interaction in the metallic phases of the chalcogens. In particular, the most common probe of $`\alpha ^2F`$ in superconductors, quasiparticle tunneling, is challenging because of the need to make well characterized tunnel junctions connected to leads that enter the high-pressure cell.

Measurement of superconducting properties such as the transition temperature $`T_C`$, the thermodynamic critical field $`H_C`$, the zero-temperature gap $`\mathrm{\Delta }_0`$, and the isotope effect on $`T_C`$, would be one avenue for probing the electron-phonon interaction. Since superconducting properties depend not only on the electron-phonon spectral function $`\alpha ^2F`$, but also on the Coulomb pseudopotential $`\mu ^{}`$, it is necessary to measure a combination of superconducting properties to get a handle on both the electron-electron and the electron-phonon parameters. In the weak coupling regime, $`T_C`$, $`\mathrm{\Delta }_0`$, and $`H_C`$ all essentially depend on the difference between $`\lambda `$ and $`\mu ^{}`$, rather than on $`\alpha ^2F`$ and $`\mu ^{}`$ separately. In this case, $`T_C`$ and the isotope effect on $`T_C`$, with their different $`\alpha ^2F`$ and $`\mu ^{}`$ dependences, would be the most useful combination of measurements.

Figure 9 shows our calculated superconducting transition temperatures for the high-pressure phases of S and Se for different values of $`\mu ^{}`$, obtained by using the calculated $`\alpha ^2F`$ functions as input to the Eliashberg equations. Available experimental data for bco Se, bco S, and $`\beta `$-Po S, are included. The theoretical data points in Fig. 9 correspond to $`\mu ^{}(\omega _{\mathrm{max}})=0.11`$, where $`\omega _{\mathrm{max}}`$ is the maximum phonon frequency. This representative value of $`\mu ^{}(\omega _{\mathrm{max}})`$ is chosen because it reproduces the experimentally measured $`T_C`$ of 17 K at 160 GPa in S. The error bars show the effect of varying $`\mu ^{}(\omega _{\mathrm{max}})`$ by $`\pm 0.02`$. While there is no reason to assume $`\mu ^{}`$ remains constant with changes in pressure or structure, or that $`\mu ^{}`$ in Se should be the same as in S, $`\mu ^{}`$ typically lies within the range of 0.1 to 0.14 for a wide range of materials. Furthermore, a rough upper bound for $`\mu ^{}`$ given by $`[\mathrm{ln}(E_F/\omega _{\mathrm{max}})]^1`$ suggests $`\mu ^{}`$ is less than about 0.16 in all the phases considered here.

The calculated transition temperatures follow the same trends as $`\lambda `$, with large increases in $`T_C`$ upon transition to both bcc Se and bcc S, and downward jumps in $`T_C`$ at the bcc$``$fcc transition in Se and the $`\beta `$-Po$``$sc transition in S. As reported earlier, our analysis attributes the large $`T_C`$ of 17 K in $`\beta `$-Po S not so much to strong electron-phonon coupling, but rather to a combination of moderate coupling and a large phonon energy scale. For analogous structural phases in S and Se, the primary reason for the higher transition temperatures in S is the stiffer lattice, which sets the energy scale for $`T_C`$.

The high phonon frequencies suggest that vertex corrections might be important in these materials. The other determining factor for the strength of vertex corrections is the magnitude of the electron-phonon interaction, which is moderate in both S and Se. Indeed, we estimate that the effect of the vertex corrections on $`T_C`$ is on the order of one percent, and can be neglected given the uncertainty in $`\mu ^{}`$.

Although we have used the Eliashberg formalism to calculate $`T_C`$ for all the phases, the coupling is sufficiently weak that the Allen-Dynes approximate formula for $`T_C`$ gives similar results. Furthermore, the ratios $`2\mathrm{\Delta }_0/k_BT_C`$ and $`\gamma T_C^2/H_C^2`$, where $`\gamma `$ is the linear coefficient of the specific heat, are calculated to be close to the BCS values of 3.53 and 1.68, respectively. For $`\beta `$-Po S at 160 GPa, for example, we find $`2\mathrm{\Delta }_0/k_BT_C`$ = 3.70 and $`\gamma T_C^2/H_C^2`$ =1.66, assuming $`\mu ^{}(\omega _{\mathrm{max}})`$ = 0.11.

As suggested above, measurement of the isotope effect on $`T_C`$ would be one way to obtain more information on the Coulomb parameter. For S at 160 GPa, we estimate that the isotope exponent $`\alpha =d\mathrm{ln}T_C/d\mathrm{ln}M`$ changes from 0.48 to 0.43 as $`\mu ^{}`$ is varied from 0.09 to 0.13, assuming isotopic masses $`M`$ of 32 and 36. To distinguish between $`\mu ^{}`$ = 0.09 and $`\mu ^{}`$ = 0.13 then requires being able to measure $`T_C`$ to an accuracy of better than 0.1 K.

Infrared spectroscopy can be used to probe $`\alpha ^2F`$ independent of $`\mu ^{}`$. In the normal state, the phonons contribute to the optical conductivity via their renormalization of the electronic quasiparticles and by phonon-assisted scattering processes. Panels (a) and (b) in Fig. 10 show the normal-state $`T=0`$ reflectance and optical conductivity calculated for different values of impurity scattering rates for $`\beta `$-Po S at 160 GPa. The expected Holstein structure, which occurs at frequencies near the phonon frequencies, is small because the coupling strength is relatively weak. Hence, it is likely that the extraction of $`\alpha ^2F`$ from phonon structure in the optical conductivity data would be challenging for this system. However, optical conductivity data could still be useful for estimating $`\lambda `$. The electron-phonon interaction narrows the Drude peak associated with impurity scattering by a factor of 1+$`\lambda `$ due to the renormalization of electronic quasiparticle energies and weights and transfers the low-frequency weight to the mid-infrared region, as illustrated in Fig. 10(c). With an independent determination of the plasma frequency $`\omega _p`$, the width and height of the Drude peak could be used to estimate $`\lambda `$.

Perhaps a simpler way to determine $`\lambda `$ from experiments is via the dc conductivity. At temperatures on the order of the Debye temperature, the phonon-limited electrical resistivity is linear in temperature with a slope proportional to $`\lambda _{\mathrm{tr}}=2d\omega \alpha _{\mathrm{tr}}^2F(\omega )/\omega `$. Like $`\alpha ^2F`$, the transport spectral function $`\alpha _{\mathrm{tr}}^2F`$ measures the effectiveness of phonons to scatter electrons on the Fermi surface, but it is weighted to take into account the change in direction of the electron velocity. Experimental measurements of the electrical resistivity have been reported for the bco and lower pressure phases of S and Se. Extensions of these measurements to higher pressures would test our results for the electron-phonon coupling strength. Of particular interest is the linear coefficient of the resistivity in $`\beta `$-Po S. We calculate $`\lambda _{\mathrm{tr}}=0.78`$ at 160 GPa, close to our result for $`\lambda `$. We believe this is the most direct way to verify our theoretical prediction that, despite the large $`T_C`$, the electron-phonon coupling is not particularly strong in this phase.

## IV CONCLUSIONS

We have presented density-functional based calculations on the pressure-induced structural phase transitions in S and Se. Although there are some similarities in the high-pressure phase diagrams of these isovalent elements, there are also striking differences. With increasing compression, Se adopts a sequence of ever more closely packed structures ($`\beta `$-Po$``$bcc$``$fcc), while S favors more open structures ($`\beta `$-Po$``$sc$``$bcc). These differences can be understood in terms of the deeper $`d`$ pseudopotential in S arising from the lack of $`d`$ states in the S core.

All the high-pressure phases of S and Se are calculated to have moderate electron-phonon coupling strengths and superconducting transition temperatures in the range of 0.5-20 K. Structure and pressure dependences of the coupling strength and $`T_C`$ can be understood in terms of changes in the phonon spectra. In particular, the large observed $`T_C`$ in $`\beta `$-Po S is calculated to arise from a combination of a moderate mass enhancement parameter $`\lambda `$ and a large overall phonon energy scale. The large increase in $`\lambda `$ and $`T_C`$ predicted upon transformation to the bcc phase in both S and Se is due to the presence of soft modes in the bcc phonon spectra that strongly couple to electrons.

We suggest that measurements of the electrical resistivity as a function of temperature, the transition temperature $`T_C`$, and the isotope effect on $`T_C`$ would be the most promising avenues for testing our results on the electron-phonon interaction in these compressed phases. Further, $`T_C`$ measurements could be used to search for the predicted sc phase of S and the predicted fcc phase of Se since our results indicate that $`T_C`$ should change abruptly when these structural transformations take place. The calculated transition pressures for both sc S and fcc Se are within the range of current diamond-anvil-cell experiments.

## V ACKNOWLEDGMENTS

We thank R. J. Hemley and V. V. Struzhkin for valuable discussions. This work was supported by the National Science Foundation under Grant DMR-9973225 and by the NPACI.


format: text

text:
# Isotropization of Ultra-High Energy Cosmic Ray Arrival Directions by Radio Ghosts

## 1 Introduction

### 1.1 The Problem of UHECR Isotropy

The upper end of the cosmic ray spectrum, at total energies above the (as yet unobserved) Greisen-Zatsepin-Kutzmin (GZK) cut-off ($`4\times 10^{19}`$ eV - \[Greisen, 1966\], \[Zatsepin and Kuzmin, 1966\]) represents a challenge for particle-physicists and astrophysicists alike. The nature of the sources of these ultra-high energy cosmic rays and their distance scale are still unknown. Only our own galactic disk can be ruled out at present as a major source site, as a compatible anisotropy has not been observed by any of the experiments sensitive to the UHECR energy range \[Takeda, 1999, Bird et al., 1999, Medina Tanco and Watson, 1999\].

Although photons, neutrinos, or some unknown particle cannot be disregarded, the muon to electron ratio measured for extensive air showers \[Halzen et al., 1995, Hayashida, 1996\] points to hadrons as the primaries hitting the upper atmosphere and triggering the cascades. Neutrons with relativistic factor $`\gamma 10^{11}`$ decay into protons after a path of $`1`$ Mpc. Heavy nuclei, on the other hand, may lose as much as $`24`$ nucleons per Mpc \[Puget et al., 1976, Cronin, 1992\] due to photodisintegration in interactions with the cosmic microwave and infrared backgrounds. Therefore, although new estimates of the infrared background \[Malkan and Stecker, 1998\] seem to indicate that the previous photodisintegration rates may be overestimated by a factor of $`10`$ \[Stecker, 1998\], UHECR, unless galactic, are very likely light nuclei, probably mainly protons. We will assume the latter in the remaining of this work.

Depending on the large scale configuration of the intergalactic magnetic field (IGMF), field values as low as $`B_{IGMF}10^9`$ G \[Kronberg, 1994a\] may be expected. Consequently, proton gyroradius at $`E10^{20}`$ eV can be of the order of $`10^2`$ Mpc. This represents a potentially interesting perspective for astronomy since UHECR should, under favorable conditions, point to their sources opening a new astronomical window. This would allow not only the direct study of their sources, but also the measurement of the intervening intergalactic and galactic (particularly halo) magnetic fields.

Relevant to this potential is the nature of the production mechanism involved. So far, both bottom-up and top-down UHECR production mechanisms have been devised and constitute perfectly viable scenarios. The first set of models invoke some kind of acceleration agent that is always spatially associated with baryonic matter, closely traced by luminous matter. Top-down mechanisms, on the other hand, resort to the decay of relics from GUT phase transitions in the early Universe \[Bhattacharjee and Sigl, 2000\]. Under the assumption that relics aggregate like dark matter (DM) and that their decay products are susceptible of inelastic interactions with the CMB, then any observed anisotropy should trace the nearby distribution of DM. Therefore, the UHECR flux should have two components, one originated inside the galactic halo and another extragalactic still roughly traced by luminous matter. \[Dubovsky and Tinyakov, 1998\] showed that, under general conditions, the halo component would dominate the extragalactic flux by at least two orders of magnitude. This is only true, however, in the unrealistic case of dark matter uniformly distributed in intergalactic space. Nevertheless, dark matter aggregates strongly and tends to be overabundant, by factors of $`10^2`$, in the center of galaxy clusters when compared to its abundance in the halos of isolated galaxies. It can therefore be shown that \[Medina Tanco, 2000a\], in a sample of 47 events (as AGASA’s above $`4\times 10^{19}`$ eV), and assuming Virgo as the only source of extragalactic events, 3-7 events should originate in Virgo and arrive inside a solid angle of approximately the size of the cluster. This could give rise to a slight anisotropy that correlates with the SGP when combined with the almost isotropic flux originated in a large galactic halo.

Photo-pion production interactions with the cosmic microwave background (CMB) restrict the source region to nominally $`D\stackrel{}{<}100`$ Mpc \[Waxman et al., 1997, Hillas, 1998\]. Therefore, the signatures of the characteristic inhomogeneity of the local universe should be clearly observable in the angular arrival distribution of UHECR regardless of the production mechanism. In spite of the low statistics available at present, the latter does not seem to be the case \[Takeda, 1999, Medina Tanco, 2000b\].

If UHECR are charged particles originated inside the large scale structures in the nearby Universe, then any isotropization must be due to intervening magnetic fields. This could happen either inside the galactic halo, like in the presence of a magnetized galactic wind \[Ahn et al., 2000\], by the interaction with highly structured IGMF inside walls and filaments \[Medina Tanco, 1998b\] or by the scattering off magnetic irregularities permeating the intergalactic medium.

In the present work we analyze the latter possibility considering radio ghosts \[Enßlin, 1999\] as scattering centers of UHECR protons.

### 1.2 The Existence of Radio Ghosts

Active galaxies eject large amounts of radio emitting plasma – short: radio plasma – into their environment. There it forms the typical cocoons of radio galaxies. The radio emission results from synchrotron emission of a population of relativistic electrons in the radio plasma’s magnetic fields. It is possible to estimate the minimal energy density of the electron and magnetic components of the radio plasma required in oder to produce the observed emissivity. This minimum is given by rough energy equipartition between electron population and the magnetic fields. The resulting minimal pressure is typically of the order of the environmental thermal pressure, indicating that relatively strong magnetic fields are present.

After a cosmological short time of $`10^710^8`$ years the radio luminosity of the cocoon decays strongly due to radiation and expansion energy losses of the electron population, or since the central engine of the radio galaxy stopped its activity. Although undetectable by our instruments, the radio plasma is still present in the IGM inside a fossil radio cocoon, a so called radio ghost. The subsequent evolution of the radio plasma is unclear, since observationally poorly constrained. It can be expected that the strong magnetic fields of the ghosts allow it to resist erosion by subsonic turbulence. The fossil radio cocoon should be kinematically decoupled from the ballistic motion of the parent galaxy and follow mostly the flow pattern of the embedding material. From this one would expect the ghosts to have a cosmological distribution comparable to that of the galaxies. And one would further expect that the oldest ghosts are swept into clusters of galaxies by the flow of structure formation. But buoyancy of the probably very light radio plasma can produce some relative motion between ghosts and the IGM gas. That buoyancy is operating on radio plasma is impressively demonstrated by a recent simultion of a buoyant rising radio blob in a cluster center, which reproduces radio and X-ray observations of M87 strikingly well \[Churazov et al., 2000\].

This might allow the ghosts in clusters or filaments of galaxies to ascend to larger radii, until they get stopped by freely infalling matter at the accretion shock. It depends crucially on the topology of this accretion shock surface if the ghosts are able to escape from gravitationally bound structures as clusters and filaments of galaxies, or not. It is therefore difficult to predict the spatial distribution of ghosts. In this work we consider cases where the ghosts are distributed as the galaxies, where they are more, and where they are less clustered.

There are several observations supporting their existence and future experiments might be able to see further signatures of ghosts.

* The progenitor of radio ghosts, radio cocoons of active radio galaxies, are observed and well studied. Since the sharp boundaries of radio cocoons indicate that these objects are able to resists their very turbulent birth, it is reasonable to assume that they also survive the much weaker ambient turbulence.
* A radio ghost can be revived to radio emission by the compression in a shock wave. This should lead to regions of diffuse synchrotron emission within clusters of galaxies, where accretion or merger shock waves are expected. Such objects are indeed observed, and are named cluster radio relics. Models for such processes can be found in Enßlin et al. \[Enßlin et al., 1998a\] and Enßlin & Krishna (2000).
* Another class of diffuse radio emitter in galaxy cluster might also be connected to ghosts: cluster radio halos. There is the speculative possibility that the short-lived radiating electrons are replenished by the decay of charged pions produced in hadronic collisions of long-lived cosmic ray protons, which were released from their confinement in radio ghosts \[Enßlin, 1999\].
* Although at low energies, the electron population within ghosts should remain mostly relativistic. Comptonization of the CMB photons by this electron population should produce CMB distortions, which have a spectral signature comparable to that of a grey absorber. This signature is weak with an optical depth of the order of $`10^7`$ \[Enßlin and Kaiser, 2000\], but principally observable. The energetically up-scattered CMB photons might also be detectable.

For this work, to estimate the deflection of UHECR by ghosts, most of their above discussed properties are not of any importance. Only the magnetic fields are required, and their existence is nearly guaranteed by the existence of radio galaxies, and the extremely low magnetic diffusivity in extragalactic plasmas. Therefore, in the context of this work, the term ‘radio ghost’ can be simply translated into ‘remnant magnetic fields of a former radio cocoon’. Even if turbulence is able to disrupt ghosts, the shear motions should increase the magnetic energy content of the fields, leading to an enhanced efficiency in UHECR deflection. We summaries, that the existence of regions in the IGM filled with remnant magnetic fields of former radio galaxies – as it is assumed in this work – is difficult to be doubted.

## 2 Properties of Radio Ghosts

### 2.1 The Distribution Function of Radio Ghosts

The properties of a radio ghost depend on three factors: the state of the radio plasma, as released by a radio galaxy (RG), the turbulent history of the environment, and the present environment of the ghost. We are mainly interested in ghosts located in filaments of galaxies, where turbulence is less strong compared to galaxy clusters. In our simplistic approach, we neglect therefore any environmentally induced changes in the ghost properties. We assume that a ghost was originally produced with a total energy $`E_{\mathrm{gh}}`$, and this energy did not change significantly over cosmological times. The ghost consists of magnetic fields and relativistic particles (electrons and protons), which should be roughly in energy equipartition. The ghost’s internal energy density

$$\epsilon _{\mathrm{gh}}=\epsilon _B+\epsilon _{\mathrm{CR}}2\epsilon _B$$
(1)
is related to it’s internal pressure via the relativistic equation of state, which the CR population and the tangled magnetic fields do follow:

$$P_{\mathrm{gh}}=\frac{1}{3}\epsilon _{\mathrm{gh}}\frac{2}{3}\epsilon _B.$$
(2)
Since the ghost is in pressure equilibrium with its environment $`P_{\mathrm{gh}}=P`$, the energy released by the RG can be written as the sum of the internal energy of the ghost and the volume work required to produce a cavity in the IGM for the ghost:

$$E_{\mathrm{gh}}=\epsilon _{\mathrm{gh}}V_{\mathrm{gh}}+PV=4PV_{\mathrm{gh}}.$$
(3)
The relevant properties of a ghost, as volume $`V_{\mathrm{gh}}`$ and magnetic field strength $`B_{\mathrm{gh}}`$, are therefore determined by $`P`$ and $`E_{\mathrm{gh}}`$:

$$V_{\mathrm{gh}}(E_{\mathrm{gh}},P)=\frac{E_{\mathrm{gh}}}{4P}$$
(4)
$$B_{\mathrm{gh}}(P)=\sqrt{12\pi P}$$
(5)
Thus we need a model for the IGM pressure. The enthalpy of the gas is increased during structure formation by compression of the flow, which is driven by the evolving gravitational potential111During structure formation the gas follows the energy equation
$$\frac{d}{dt}\left(\frac{1}{2}v_{\mathrm{gas}}^2+\frac{\gamma }{\gamma 1}\frac{P_{\mathrm{gas}}}{\varrho _{\mathrm{gas}}}+\mathrm{\Phi }_{\mathrm{grav}}\right)=\frac{\mathrm{\Phi }_{\mathrm{grav}}}{t},$$
(6)  if heating and cooling processes can be neglected. This equation can be integrated, if we approximate the potential to be time independent, yielding
$$\frac{1}{2}v_{\mathrm{gas}}^2+\frac{\gamma }{\gamma 1}\frac{P_{\mathrm{gas}}}{\varrho _{\mathrm{gas}}}=(\mathrm{\Phi }_{\mathrm{grav}}\mathrm{\Phi }_{\mathrm{grav},\mathrm{o}}).$$
(7)  Whenever the kinetic energy is small a unique relation between the potential and the specific enthalpy is given.  . This picture neglects non-gravitational heating and cooling processes. Cooling is only important in the center of dense clusters of galaxies. Whenever the gas is subsonic, the pressure is related to the gas density via

$$P(n_{\mathrm{gas}})=P_\mathrm{o}(n_{\mathrm{gas}}/n_{\mathrm{gas},\mathrm{o}})^\gamma ,$$
(8)
where $`\gamma =5/3`$ is the adiabatic index of the gas, otherwise the pressure is lower. We neglect deviations from this equations assuming that only a small fraction of gas has a sonic or higher velocity. This is justified, since supersonic gas becomes soon shocked when it hits the next accretion shock around a cluster or filament of galaxies.

Unfortunately, the distribution of gas in the local Universe is not observed with current instruments. But the number density $`n_{\mathrm{gal}}`$ of the galaxy distribution can be used to give a rough idea of the gas and pressure distribution. We assume that the baryon fraction of the matter density is approximately the same in all structures ($`\varrho _{\mathrm{gas}}/\varrho _{\mathrm{matter}}=\mathrm{const}`$), and that the bias between optical galaxies and matter can be approximated by

$$n_{\mathrm{gal}}=n_{\mathrm{gal},\mathrm{o}}(n_{\mathrm{gas}}/n_{\mathrm{gas},\mathrm{o}})^{b_{\mathrm{gal}}}.$$
(9)
The galaxy bias can be found to be $`b_{\mathrm{gal}}=1.6`$ by linearizing this expression and comparison with the observed linear bias for optical galaxies \[Peacock, 1999\]. Thus we get a rough proportionality between galaxy density and pressure:

$$P(n_{\mathrm{gal}})=P_\mathrm{o}(n_{\mathrm{gal}}/n_{\mathrm{gal},\mathrm{o}})^{\gamma /b_{\mathrm{gal}}}n_{\mathrm{gal}}^{1.042}.$$
(10)
This relation can be normalized by using the observed properties of clusters of galaxies. A typical density and temperature of a cluster is $`n_e=10^3\mathrm{cm}^3`$ and $`kT=5`$ keV, leading to $`P_o=P_{\mathrm{cluster}}=10\mathrm{keV}\mathrm{cm}^3`$. We therefore get a ghost magnetic field strength of $`25\mu `$G inside of galaxy clusters and $`0.54\mu `$G in galaxy filaments (assuming $`n_e=10^5\mathrm{cm}^3`$).

The distribution of ghost energies is assumed to follow a Schechter distribution

$$n_{\mathrm{gh}}(E_{\mathrm{gh}})dE_{\mathrm{gh}}=N_{\mathrm{gh}}\left(\frac{E_{\mathrm{gh}}}{E_{}}\right)^{\alpha _{\mathrm{gh}}}\mathrm{exp}\left(\frac{E_{\mathrm{gh}}}{E_{}}\right)\frac{dE_{\mathrm{gh}}}{E_{}},$$
(11)
for $`E_{\mathrm{gh}}>E_{\mathrm{min}}`$, $`n_{\mathrm{gh}}(E_{\mathrm{gh}})=0`$ otherwise. Both the normalization factor $`N_{\mathrm{gh}}`$ and the upper cutoff $`E_{}`$ are allowed to depend on the position, parameterized by the galaxy density ($`N_{\mathrm{gh}}=N_{\mathrm{gh}}(n_{\mathrm{gal}})`$ and $`E_{}=E_{}(n_{\mathrm{gal}})`$).

It is convenient to introduce the ghost-volume-filling factor

$$\mathrm{\Phi }_{\mathrm{gh}}=_{E_{\mathrm{min}}}^{\mathrm{}}𝑑E_{\mathrm{gh}}V_{\mathrm{gh}}(E_{\mathrm{gh}},P)n_{\mathrm{gh}}(E_{\mathrm{gh}},n_{\mathrm{gal}}),$$
(12)
and the average energy density in ghosts

$$ϵ_{\mathrm{gh}}=\frac{4}{3}\epsilon _{\mathrm{gh}}\mathrm{\Phi }_{\mathrm{gh}}.$$
(13)
The factor $`\frac{4}{3}`$ appears here, since our definition of the energy of ghosts included the volume work done by the expanding radio lobe. This part of the ghosts energy, which is $`\frac{1}{4}`$ of its total energy budget, is located in the IGM in the form of heat, kinetic, and gravitational energy.

The normalization factor can therefore be expressed as

$$N_{\mathrm{gh}}=\frac{ϵ_{\mathrm{gh}}}{E_{}\mathrm{\Gamma }(2\alpha _{\mathrm{gh}},E_{\mathrm{min}}/E_{})},$$
(14)
where $`\mathrm{\Gamma }`$ is the incomplete Gamma-function:

$$\mathrm{\Gamma }(a,x):=_x^{\mathrm{}}𝑑tt^{a1}e^t.$$
(15)
The total number density in ghosts reads

$$n_{\mathrm{gh}}=_{E_{\mathrm{min}}}^{\mathrm{}}E_{\mathrm{gh}}n_{\mathrm{gh}}(E_{\mathrm{gh}})=\frac{ϵ_{\mathrm{gh}}}{E_{}}\frac{\mathrm{\Gamma }(1\alpha _{\mathrm{gh}},E_{\mathrm{min}}/E_{})}{\mathrm{\Gamma }(2\alpha _{\mathrm{gh}},E_{\mathrm{min}}/E_{})}.$$
(16)
Finally we assume that the energy density in ghosts scales with the galaxy density to some power $`b_{\mathrm{gh}}`$

$$ϵ_{\mathrm{gh}}(n_{\mathrm{gal}})=ϵ_{\mathrm{gh},\mathrm{o}}(n_{\mathrm{gal}}/n_{\mathrm{gal},\mathrm{o}})^{b_{\mathrm{gh}}}.$$
(17)
$`b_{\mathrm{gh}}=1`$ would correspond to constant energy output per galaxy. This might be unrealistic, since the older, and bigger galaxies in the denser regions are expected to have produced more radio plasma during their lifetime. We therefore also adopt $`b_{\mathrm{gh}}=2`$ in order to have a ghost distribution more clustered than the galaxies.

On the other hand, buoyant motion of the ghost might have allowed a significant fraction to escape from clusters and filaments of galaxies. This could have let to a distribution of ghosts which is less clustered than the galaxies. We mimic this case by also using $`b_{\mathrm{gh}}=1/2`$.

### 2.2 Guessing the Parameters

The parameters, which are still undetermined, are $`ϵ_{\mathrm{gh},\mathrm{o}}`$, $`\alpha _{\mathrm{gh}}`$, $`E_{\mathrm{min}}`$, and $`E_{}(n_{\mathrm{gal}})`$.

The progenitors of ghosts are RGs. We therefore try to use the observed radio luminosity function (RLF) of radio galaxies to get an idea about the shape and normalization of the ghost size distribution. The present day spatially averaged number density of ghosts is simply given by the time integral of their birth rate

$$\overline{n}_{\mathrm{gh}}(E_{\mathrm{gh}})=𝑑t\dot{\overline{n}}_{\mathrm{gh}}(E_{\mathrm{gh}},t).$$
(18)
The birth rate is given by the number density of visible radio galaxies divided by their typical lifetime $`t_{\mathrm{visible}}`$. We assume that the time the galaxy is active $`t_{\mathrm{active}}`$ and the time the galaxy is visible $`t_{\mathrm{visible}}`$ are the same for all galaxies ($`t_{\mathrm{active}}=t_{\mathrm{visible}}=t_{\mathrm{rg}}=310^7`$ yr). This is a typical value out of the observed range of ages of large expanded radio lobes, which are of the order of 100 Myr \[Alexander and Leahy, 1987, Komissarov and Gubanov, 1994, Venturi et al., 1998, Slee and Roy, 1998\].Although these fairly robust age estimates are based on interpreting the observed radio spectral steepening in terms of the inverse Compton losses of the relativisitic electrons against the ubiquitous cosmic microwave background photons uncertainties and scatter are likely up to one order of magnitude. We note, that the assumed age only affects the distribution of the total ghost energy budget onto small and large ghosts sizes, but not the budget itself.

We further assume that the jetpower $`q_{\mathrm{jet}}`$ and radio luminosity $`L_\nu `$ of a RG does not change during its active phase. These are very coarse assumptions. But we hope they are justified since we use them only to get a feeling for the ghosts size distribution. Radio galaxies with jetpower $`q_{\mathrm{jet}}`$ produce therefore ghosts with the energy

$$E_{\mathrm{gh}}=\frac{1}{2}q_{\mathrm{jet}}t_{\mathrm{rg}}$$
(19)
with a birth rate of

$$\dot{\overline{n}}_{\mathrm{gh}}(E_{\mathrm{gh}},t)=2\frac{n_{\mathrm{jet}}(q_{\mathrm{jet}})}{t_{\mathrm{rg}}}\frac{dq_{\mathrm{jet}}}{dE_{\mathrm{gh}}}=4\frac{n_{\mathrm{jet}}(q_{\mathrm{jet}},t)}{t_{\mathrm{rg}}^2},$$
(20)
where $`n_{\mathrm{jet}}(q_{\mathrm{jet}},t)dq_{\mathrm{jet}}`$ is the jetpower distribution function of RGs. The number two appears since RG have typically two radio lobes.

The jetpower distribution function can be derived from the radio luminosity function of RG with the additional assumption that there is a unique relation between radio luminosity and jetpower. Enßlin et al. (1997) have fit a power-law relation to the jetpower derived by Rawlings and Saunders (1991) of a sample of radio galaxies. Their values are based on minimum energy arguments and age estimates. The real energy of the radio plasma can easily be much higher than the minimal energy estimate by some factor $`f_{\mathrm{power}}>1`$ due to the presence of relativistic protons, low energy electrons or deviations from equipartition between particle and field energy densities. A rough estimate of $`f_{\mathrm{power}}`$ can be derived from observations of radio lobes embedded in the intra-cluster medium (ICM) of clusters of galaxies. They show a discrepancy of the thermal ICM-pressure to the pressure in the radio plasma following minimal energy arguments by a factor of $`510`$, even if projection effects are taken into account \[Feretti et al., 1992\]. Since also a filling factor smaller than unity of the radio plasma in the radio lobes can mimic a higher energy density we chose $`f_{\mathrm{power}}=3`$ in order to be conservative.

The jetpower-radio luminosity correlation at $`\nu =2.7`$ GHz then reads

$$q_{\mathrm{jet}}=a_\nu (L_\nu /(\mathrm{Watt}\mathrm{Hz}^1h_{50}^2))^{b_\nu }f_{\mathrm{power}}$$
(21)
for which $`b_\nu =0.82\pm 0.07`$, $`a_\nu =10^{45.2826.22b_\nu \pm 0.18}\mathrm{erg}\mathrm{s}^1h_{50}^2`$ \[Enßlin et al., 1997\]. This relation allows us to translate the observed radio luminosity function into a jetpower distribution function. We adopt different radio-luminosity functions parameterized by Dunlop & Peacock \[Dunlop and Peacock, 1990\] and integrate Eq. 18 using Eqs. 20 and 21 in an Einstein de Sitter Cosmology up to redshift $`z=10`$ (the result depends only weakly on the upper cutoff due to the strong decline of the used RLFs beyond $`z=3`$, the pure luminosity evolution (PLE) RLF is only integrated up to $`z=4`$). We further calculate the ghost distribution function for $`b_\nu =0.7,0.82,1`$ in order to show the dependence on the uncertainties. The results can be seen in Fig. 1.

Fig. 1 shows that $`\alpha _{\mathrm{gh}}2`$, which we adopt in the following. A reasonable lower cutoff is then $`E_{\mathrm{min}}=610^{57}`$ erg (corresponding to $`L_\nu =10^{23}\mathrm{Watt}\mathrm{Hz}^1`$). The ghost distribution does not need to have a real cutoff there, but the part which is accessible to radio observations of the parent radio galaxies ends there. In order to conservative we do not extrapolate to lower energies.

The radio luminosity function itself becomes dominated by starburst galaxies at smaller radio luminosities. This gives a second contribution to the magnetization of the IGM from magnetized galactic winds. An extrapolation of the strongly redshift dependent, recent starburst history of galaxies to high redshifts ($`z=10`$) indicates that a substantial cosmological volume might be magnetized by such winds \[Kronberg et al., 1999\]. However, the fields resulting from this process are expected to be tangled on small, galactic scales, and therefore less effective scattering centers for UHECR as radio ghosts.

An integration of the ghost energy distribution functions shown in Fig. 1 gives the total energy output of radio galaxies per volume. The results do not depend on the assumed lifetime $`t_{\mathrm{rg}}`$ of radio galaxies, they scale linearly with $`f_{\mathrm{power}}`$, and depend only moderately on $`b_\nu `$. Detailed results are given in Tab. 1. The energy input in form of radio plasma is roughly $`\overline{ϵ}_{\mathrm{gh}}=10^{66\mathrm{}67}\mathrm{erg}\mathrm{Gpc}^3`$, it is distributed over $`\overline{n}_{\mathrm{gh}}=10^{7\mathrm{}8}`$ ghosts per Gpc3, which have an average energy of $`\overline{\epsilon }_{\mathrm{gh}}=10^{58\mathrm{}59}`$ erg. Note that the X-ray background, which is believed to be dominated by AGN emission, corresponds to an injection energy of $`310^{67}\mathrm{erg}\mathrm{Gpc}^3`$ (Chokshi & Turner 1992; Soltan 1982). Either the X-ray energy losses of AGNs exceed the radio plasma release, or the jetpower is strongly underestimated here and $`f_{\mathrm{power}}=10\mathrm{}30`$ would be more realistic. Even such a strong energy input into the IGM seems to be consistent with present the day limit to the Comptonizarion of the IGM \[Enßlin and Kaiser, 2000, Enßlin et al., 1998b\].

The next parameter we have to fix for our model is the maximal ghost energy. $`E_{}(n_{\mathrm{gal}})`$ can be roughly estimated from observations as shown in the following. The progenitor of large ghosts in rare regions should be giant radio galaxies (GRGs). Their linear size distribution seems to have a sharp cutoff at 3 Mpc $`h_{50}^1`$ (Ishwara-Chandra & Saikia 1999). There is only one GRG significantly exceeding this: 3C236 has a linear size of 5.7 Mpc $`h_{50}^1`$. This could be – for example – due to a very low pressure in the environment of 3C236. For a typical galaxy filament environment, we assume the linear size of GRG to have a cutoff at 3 Mpc $`h_{50}^1`$. The diameter of one of the two radio lobes is roughly $`\frac{1}{3}`$ of the linear size of the whole GRG. We adopt therefore a maximal radius of our (spherically approximated) ghosts in the dilute environment of galaxy filaments of $`r_{}=0.5`$ Mpc $`h_{50}^1`$. The distribution function of ghost sizes

$$f_\mathrm{o}(r)dr=\frac{3(r/r_{})^{23\alpha _{\mathrm{gh}}}}{\mathrm{\Gamma }(1\alpha _{\mathrm{gh}},E_{\mathrm{min}}/E_{})}\mathrm{exp}\left(\frac{r^3}{r_{}^3}\right)\frac{dr}{r_{}}$$
(22)
has indeed a very sharp cutoff at $`r_{}`$, as required by the observation of GRG. This $`r_{}`$ is related to $`E_{}`$ via Eqs. 3 and 2:

$$E_{}=\frac{16}{3}\pi r_{}^3P\frac{32}{9}\pi r_{}^3\epsilon _B=\frac{4}{9}r_{}^3B_{\mathrm{gh}}^2$$
(23)
Equipartition (or minimum energy) magnetic field strength of GRG are of the order of a few $`\mu `$G (Ishwara-Chandra & Saikia 1999), so that we adopt a conservative value of $`B_{\mathrm{gh}}=1\mu `$G for ghosts located in galaxy filaments. From this we find that $`E_{}=1.610^{60}\mathrm{erg}`$, and that a corresponding RG had ejected twice as much energy ($`3.210^{60}\mathrm{erg}`$), which seems to be moderate. Comparison with Fig. 1 shows that ghosts with $`10100`$ times higher energy are expected from the radio counts (if $`t_{\mathrm{rg}}`$ is as big as $`310^7`$ yr). Such ghosts are also required since a $`10^{60}`$ erg ghosts in a dense cluster environment ($`n_\mathrm{e}=10^3\mathrm{cm}^3B_{\mathrm{gh}}=25\mu `$G) has a diameter of $`0.1`$ Mpc, which is clearly below the observed maximal size of radio lobes inside clusters. In order to allow a diameter of $`0.5`$ Mpc inside a cluster one needs $`E_{}10^{62}`$ erg. We therefore assume that the maximal energy of a ghost scales linearly with the gas density:

$$E_{}=10^{62}\mathrm{erg}\frac{n_\mathrm{e}}{10^3\mathrm{cm}^3}=10^{62}\mathrm{erg}\left(\frac{n_{\mathrm{gal}}}{n_{\mathrm{gal},\mathrm{o}}}\right)^{1/b_{\mathrm{gal}}},$$
(24)
where $`n_{\mathrm{gal},\mathrm{o}}`$ is the galaxy density in the center of a typical cluster of galaxies. Such a scaling might be justified due to the fact that the most massive galaxies, which probably host the most powerful AGNs, are the central elliptical galaxies within clusters.

For the purpose of our simulation it also makes sense to define the average total cross section for the assemble of ghosts:

$$\overline{\sigma }_{\mathrm{gh}}=𝑑rf(r)\pi r^2=\pi r_{}^2\frac{\mathrm{\Gamma }(\frac{5}{3}\alpha _{\mathrm{gh}},E_{\mathrm{min}}/E_{})}{\mathrm{\Gamma }(1\alpha _{\mathrm{gh}},E_{\mathrm{min}}/E_{})}$$
(25)
### 2.3 Constructing the Ghost Distribution

The scaling relations in Eqs. 10, 17, and 24 are normalized for a reference density equal to the central galaxy density of a galaxy cluster. We estimate $`n_{\mathrm{gal},\mathrm{o}}`$ by identifying the peaks in our galaxy distribution. The energy density of ghosts at our reference gas density $`n_{\mathrm{gal},\mathrm{o}}`$ (inside clusters) is given by

$$ϵ_{\mathrm{gh},\mathrm{o}}=\frac{\overline{ϵ}_{\mathrm{gh}}V}{𝑑V(n_{\mathrm{gal}}/n_{\mathrm{gal},\mathrm{o}})^{b_{\mathrm{gh}}}},$$
(26)
where $`\overline{ϵ}_{\mathrm{gh}}=10^{66\mathrm{}67}\mathrm{erg}\mathrm{Gpc}^3`$ and V is the volume of our galaxy catalog. Whenever we use the CfA redshift catalog we exclude volumes close to the galactic plane in our estimate of V and the above integral due to the under-sampling of galaxies there.

If the ghosts would get disrupted and mix with the ambient gaseous medium a significant magnetization of the IGM would result. Assuming complete mixing, we estimated the resulting magnetic field strengh as a function of the local cosmic matter density in a few scenarios ($`\overline{ϵ}_{\mathrm{gh}}=10^{66\mathrm{}67}\mathrm{erg}\mathrm{Gpc}^3`$, $`b_{\mathrm{gh}}=0.5`$ and $`b_{\mathrm{gh}}=1`$). The field strength as a function of density is plotted in Fig. 12, and the volume fraction of the different density regimes in Fig. 13. Such magnetic fields are not excluded by Faraday rotation measurements, which give limits of $`B_{\mathrm{IGMF}}<nG`$ for a homogeneous background field with Mpc coherence-scales \[Kronberg, 1994b\], or $`B_{\mathrm{IGMF}}<\mu `$G if cosmological fields are concentrated in the large scale structures like sheets and filaments of galaxies \[Ryu et al., 1998a, Blasi et al., 1999\]. The implied values for cluster of galaxies of $`B_{\mathrm{ICM}}0.3\mathrm{}3\mu `$G are also consistent with Faraday rotation measurements of magnetic fields. Note, that shear flows due to cosmic structure formation can significantly increase the strength of IGM magnetic fields \[Dolag et al., 1999, Roettiger et al., 1999\].

### 2.4 The Radio Ghosts Produced by Our Galaxy

The formation of the $`M_{\mathrm{bh}}=2.610^6M_{}`$ black hole in the center of our galaxy (Eckart et al. 1997) was probably accompanied by dissipative processes as radiation and ejection of relativistic plasma in a temporary AGN. We assume that a fraction $`ϵ_{\mathrm{diss}}=0.1`$ of the rest mass energy of the accreted matter was dissipated and that $`ϵ_{\mathrm{plasma}}=0.5`$ of this energy went into the relativistic plasma. The radio ghosts formed by the remnant of this radio plasma got therefore an energy budget of

$$E_{\mathrm{gh}}=\frac{ϵ_{\mathrm{plasma}}ϵ_{\mathrm{diss}}}{1ϵ_{\mathrm{diss}}}M_{\mathrm{bh}}c^2=2.610^{59}\mathrm{erg}.$$
(27)
The term $`1ϵ_{\mathrm{diss}}`$ in the denominator corrects for the mass discrepancy between accreted mass and final black hole mass due to dissipation of energy. Eqs. 1, 2, and 3 imply that the filled volume is given by

$$V_{\mathrm{gh}}=\frac{3E_{\mathrm{gh}}}{8\epsilon _B}=0.083\mathrm{Mpc}^3\left(\frac{B}{\mu \mathrm{G}}\right)^2.$$
(28)
Comparing this to the gyroradius $`r_\mathrm{g}`$ of a $`10^{20}\mathrm{eV}`$ cosmic ray proton shows the possible importance of the ghosts of our own galaxy: $`V_{\mathrm{gh}}66r_\mathrm{g}^3(B/\mu \mathrm{G})`$.

But it cannot be assumed that the ghosts produced by our own galaxy are located in the direct vicinity of our galaxy since the observed peculiar velocity of 162 km/s with respect to our own galaxy group (Rauzy & Gurzadyan 1998) indicates differential motion between the Milky Way and the local IGM. Radio ghosts produced by other members of the local group may be nearby today, but this is already included into our statistical approach described in Sec. 2, where the density of ghosts scales with the number density of galaxies.

## 3 UHECR Flux at Earth

### 3.1 Scattering by Radio Ghosts

Equation (5) can be used to calculate the intensity of the magnetic field inside ghosts. The estimation of the topology of the internal field is, however, a much more difficult matter.

The generation of any ghost was a turbulent event fed by the parent radio galaxy. This process certainly produced vortices at all scales. Ghosts, however, have been dynamically decoupled afterwards and, therefore, smaller scale wiggles should have been dissipated over time. Hence, most of the remaining power in the field probably lays at large scales at present.

In Fig. 2 we exemplify one possible internal magnetic field topology. The shown magnetic field has null perpendicular component at the surface and is therefore decoupled from the environment field. The amplitude of the fluctuations obeys a power law power spectrum of spectral index $`\xi =4`$.

In order to analyze the scattering capability of ghosts and to estimate how critical is the uncertainty in the actual topology of the field to the present study, we performed numerical simulations using spherical radio ghosts, similar to the one in Fig. 2. Different power spectral indexes $`\xi `$ and radius were considered. We illuminate the ghosts with a beam of test particles, protons, with a power law energy spectrum, $`dN/dEE^3`$, and energy $`E>4\times 10^{19}`$ eV. The angular distribution of the scattered particles is shown in Fig. 3.a for ghosts of different sizes and $`\xi =5/3`$, $`3`$ and $`4`$. The radii of the ghosts are given in terms of the average gyroradius of the injected UHECR spectrum. As a comparison, an isotropic scattering particle distribution is also shown in the same figure. It can be seen that ghosts are very efficient at scattering UHECR regardless of power spectral index. Forward scattering becomes dominant only at the low end of the ghost size distribution, when the ghost’s radius becomes comparable to the average gyroradius of the particles. Figure 3.b is the same as figure 3.a but for a monoenergetic proton injection at $`E=10^{20}`$ eV. The same ghosts radii as in fig. 3.a (physical units) were used, and so they are down by a factor of two in units of UHECR gyroradius in fig. 3.b. The stronger dominance of forward scattering for the smallest ghosts can be appreciated, as well as for the flattest fluctuation spectrum ($`\xi =5/3`$). Nevertheless, as previously argued, small scale magnetic turbulence inside a typical ghost should have been dissipated over time, and a steep spectrum of magnetic fluctuations is expected.

Therefore, given the current uncertainties in radio ghost shape, size, magnetic field topology and spatial distribution, isotropic scattering is, very likely, an acceptable assumption and it will be used in what follows.

### 3.2 Numerical Model

We use Monte Carlo numerical simulations to track UHECR propagation through the intergalactic medium and to evaluate their arrival distribution at Earth.

We start from the basic assumption that most UHECR are protons of extragalactic origin whose sources aggregate spatially as either luminous or cold DM.

Charged particles, even at the extreme energies considered, are coupled to the intervening magnetic fields. From this point of view, the propagation region can be divided in several components: (a) the sources and their immediate neighborhood, (b) the intergalactic medium, (c) the galactic halo, (d) the galactic disk, (e) the heliosphere and (f) the magnetosphere.

The close environment of the sources can be included inside the sources just by re-defining them as large as needed. This should be acceptable in general since, the newly defined source, would still be below the resolution capability of any current experiment.

The heliospheric and magnetospheric fields, despite their relatively large values, are confined to small volumes and, from the point of view of deflection, are negligible.

The galactic disk magnetic field is one to two orders of magnitude smaller than the latter fields, but fills regions at least 8 orders of magnitude larger. However, exception made of lines of sight crossing directly the galactic bulge or a relatively narrow strip surrounding the galactic plane, the disk field has probably little impact on UHECR deflection.

Consequently, the two main factors responsible for blurring pointing information, are the intergalactic (IGMF) and the galactic halo magnetic field.

Few is known regarding the large scale structure of the IGMF and only sparce rotation measurements are available to constrain its intensity \[Kronberg, 1994b, Vallee, 1997\]. Basically, two extreme scenarios can be imagined that satisfy the observational constraints. In one of them, the field is compressed inside walls and filaments, where it attains high intensity ($`0.11\mu `$G) and high degree of correlation (over scales of up to tens of Mpc), leaving unmagnetized voids behind \[Ryu et al., 1998b\]. We will call this the laminar IGMF model. A second possibility is that of a more evenly distributed field, whose intensity scales somehow with matter and is coherent on scales $`L_c1\mathrm{Mpc}(B/\mathrm{nG})^2`$. A numerical model for the latter case is an ensemble of cells of characteristic size $`L_c`$, with uniform field inside but randomly oriented with respect to adjacent cells. The galaxy density is used to scale $`B`$ and $`L_c`$. We call this model the cellular IGMF.

The first IGMF scenario can affect propagation in a radical way, with UHECR mainly constrained to the interior of large scale structures and drifting along walls and filaments. The observable flux depends critically on the location of the observer and on the details of the actual topology of the field inside the GZK-sphere. Isotropy can be attained in a rather natural way, but the results are strongly model dependent. The propagation of UHECR in this IGMF has been treated in an earlier work \[Medina Tanco, 1998b\], and will not be considered here.

A similar uncertainty exists around the topology of the magnetic field inside the galactic halo. Nevertheless, unless there is a large scale magnetized galactic wind that structures this field and maintains a large azimuthal component up to large galactocentric distances, the halo should be unimportant regarding isotropization. If a magnetized galactic wind does exist on the scale of several hundreds of kpc, then it would be the dominant factor in the determination of the UHECR flux at Earth \[Ahn et al., 2000, Biermann et al., 2000\]; therefore, extragalactic propagation details would be mostly irrelevant.

Therefore, we restrict the present analysis to a cellular IGMF and neglect the galactic halo.

The same procedure as in \[Medina Tanco, 1997, Medina Tanco, 1998a\] is used in the description of the cell-like spatial structure of the IGMF. The cell size is given by the correlation length, $`L_cB_{IGMF}^2(r)`$. The intensity of the IGMF, in turn, scales with gas density as $`B_{IGMF}n_{gas}^\eta (r)`$ and the proposed IGMF value at the Virgo cluster ($`10^7`$ G, \[Arp, 1988\]) is used as the normalization condition. We use $`\eta =0.5`$, as a compromise between a frozen-in field ($`\eta =2/3`$) and Valleé’s (1997) estimate ($`\eta 0.35`$), which may be too flat due to the assumed values for the magnetic field in superclusters and larger scales. Nevertheless, tests have been conducted for different values of $`\eta `$ covering the previous interval, and the scaling is not critical to our conclusions.

The formulation of section 2 is used to define the size and spatial distribution of radio ghosts. This requires the knowledge of the gas density distribution inside the simulation volume. Actually, as galaxies are easier to survey over large volumes, equation(9) is used to transform between galaxy and gas density distributions. There are, however, serious an unavoidable bias and sampling problems inherent to galaxy surveys. Therefore, we relay on both galaxy surveys and cold DM large scale structure simulations to perform independent evaluations.

The sources of UHECR are distributed according to either the galaxy or cold DM distributions respectively.

The 1999 version of the CfA catalog \[Huchra et al., 1992\] is used to characterize the galaxy distribution.

Cold DM simulation data are from Springel et al. (in preparation). They carried out simulations that mimic the Local Universe. The initial conditions of these simulations have been constrained by the redshift survey of IRAS galaxies. As a result, the simulations develop the same local large-scale structure (e.g., the Great Attractor and Cetus Wall; clusters like Virgo and Coma are also found at the right place).

The scenario is completed by the introduction of radio ghosts according to the formulation of section 2.

Relativistic test particles (UHECR protons) are injected at the sources with a spectrum $`dN/dEE^2`$ and propagated through the intergalactic magnetic field up to the detector on Earth. Adiabatic energy losses due to redshift, pair production and photo-pion production due to interactions with the cosmic microwave background radiation (CMBR) are also included \[Berezinsky and Grigor’eva, 1988, Achterberg et al., 1999\].

The flux at Earth can be divided into two components: (a) a direct  radiation field, constituted by particles that fly from source to detector without encountering ghosts and (b) a diffuse  radiation field, comprising particles which underwent at least one encounter with a radio ghost.

### 3.3 Different Scenarios

Fig. 4 is meant as a control. It shows the Aitoff projection of the two-dimensional arrival probability density (galactic coordinates with the anti-galactic center at the center of the figure), for sources distributed according to nearby luminous matter (CfA catalog, \[Huchra et al., 1992\]) inside $`100`$ Mpc and no radio ghosts. The same procedure as in \[Medina Tanco, 1997, Medina Tanco, 1998a\] is used in the description of the intergalactic magnetic field (IGMF): a cell-like spatial structure, with cell size given by the correlation length, $`L_cB_{IGMF}^2(r)`$. The intensity of the IGMF, in turn, scales with luminous matter density, $`\rho _{gal}`$ as $`B_{IGMF}\rho _{gal}^{0.3}(r)`$ \[Vallee, 1997\] and the observed IGMF value at the Virgo cluster ($`10^7`$ G, \[Arp, 1988\]) is used as the normalization condition. The mask covers the plane of the galaxy, where the actual distribution of galaxies is not well known due to obscuration by dust. The curved, thick line is the celestial equator. Northern hemisphere is the sky patch to the right, enclosed by that line.

Superimposed on the figure are the available events with $`E>4\times 10^{19}`$ eV observed by AGASA (47 events \[Takeda, 1999\]), Haverah Park (27 \[Reid and Watson, 1980\]), Yakutsk (24 \[Afanasiev, 1995\]) and Volcano Ranch(6 \[Linsley, 1980\]).

The arrival probability contours trace roughly the local large scale structure. Distinguishable observational signatures should be expected towards the region of the Southern branch of the supergalactic plane (to be observed in the near future by the Auger experiment) at $`l45^o`$, the lines of sight to the more distant Pisces-Perseus wall and Perseus cluster and, very prominently, towards a large area surrounding the Virgo and Ursa Major clusters.

The actually observed distribution of UHECR is clearly much more isotropic than what one would expect under the implicit assumptions in Fig. 4.

In figures 5, 6 and 7, we show how the UHECR arrival distribution function would be modified by the presence of radio ghosts in the intergalactic medium for different values of $`b_{gh}`$ and $`\alpha _{gh}=2.0`$. The results are insensitive to $`\alpha _{gh}`$ in the range of interest. The most important parameter is $`b_{gh}`$, defined in equation (17), which tells how are ghosts distributed with respect to galaxies. $`b_{gh}=1`$ implies the same spatial distribution as galaxies and larger values a more clustered distribution.

It can be seen from Fig. 5 and 6, that clustered distributions of ghosts produce no noticeable effects in the observed UHECR flux at Earth.

In fact, if these angular distributions were observed with the same exposure in declination as AGASA’s after 7.5 years of integration, the observed (ensemble averaged) amplitude of the first harmonic \[Linsley, 1975\] would vary only by $`5\%`$ when going from Fig. 4 to 6. This is basically because a distribution of scatterers that is more clustered than the sources is, in practical terms, equivalent to increase the physical size of the sources, which is irrelevant from the point of view of all sky anisotropies.

The result is different when the scatterers are distributed in a larger volume than the sources. This is the case in figure 7, where $`b_{gh}=0.5`$. In this scenario it is assumed that e.g. radio ghosts buoy out of cosmological structures, creating thick halos around walls and filaments, permeating voids to some extent. This diminishes considerably the direct component (Fig. 7a) and accounts for a considerable increase in the diffuse component (Fig. 7b) which becomes, by far, dominant. The composite flux (Fig. 7c) still shows a smooth, large scale gradient towards the region of the Virgo cluster but is much more isotropic that any of the previous scenarios.

A more quantitative picture can be obtained from Fig. 8. An ensemble of 1000 independent samples, of 50 UHECR each, was built by observing the arrival distribution function in Fig. 7c with an exposure equivalent to AGASA’s. The resultant distribution function of the amplitude of the first harmonic is shown, in comparison with the actually observed values by Yakutsk, AGASA, Haverah Park and Volcano Ranch (see, \[Medina Tanco and Watson, 1999\]). It can be seen that $`b_{gh}=0.5`$ is the upper limit for the clustering of ghosts that is compatible with the observed anisotropy. However, it is also clear that that an even smoother distribution of ghosts in intergalactic space would be much more acceptable. Increasing the volume average ghost energy density by a tolerable factor of 3 would help to move the results in the right direction, but it wouldn’t be a determinant factor by itself. Rephrased in other words, and looking at Fig. 9, the mean-free path of UHECR for interactions with ghosts in the intergalactic medium, should be reduced to less than a few Mpc at any location in the propagation region, i.e., walls and filaments as well as voids, to obtain a level of isotropy consistent with observation.

Whether such a smooth spatial distribution of radio ghosts is physically reasonable is a very complicated matter which we leave open. The possibility that radio ghosts are able to escape by buoyancy from the denser environments in which they were produced exists. However, specific large scale structure hydrodynamic simulations should be performed to decide whether buoyancy is actually operational and, furthermore, a deeper knowledge of the composition and mechanical properties of the radio plasma is necessary in order to decide if ghosts are stable enough to survive as an entity under these conditions.

As was mentioned previously, there are uncertainties associated with unavoidable biases and sampling incompleteness associated with galaxy surveys, as is the case with the CfA catalog used up to here. This problem can be specially critical in this analysis, since the results depend on an absolute normalization of the density as a function of depth into the local universe. To check the extension of the distortions occurring in our previous analysis, we repeated our calculations using the distribution of cold DM, calculated by large scale structure hydrodynamic simulations (Springel, in preparation). They carried out simulations that mimic the Local Universe, developing the observed local large-scale structure.

Fig. 10 and 11 show the arrival flux for $`\alpha _{gh}=2.0`$ and $`b_{gh}=2.0`$ and $`0.5`$ respectively for the simulation using the CDM model. The resultant UHECR flux is plotted in an Aitoff projection in supergalactic coordinates. The large spot near the center of Fig. 10 is the cosmic ray image from the Virgo cluster. The central band is the supergalactic plane. Note that the supergalactic plane is not such a neatly defined structure in UHECR with average energy $`E\stackrel{}{>}4\times 10^{19}`$ eV. It is clear that, regarding isotropization of the UHECR flux, the same kind of effect is present in these new scenarios: the dominant signature from Virgo at $`b_{gh}>1`$ is strongly diluted at $`b_{gh}0.5`$.

## 4 Discussion

The, to some extent unexpected, isotropy in the UHECR flux observed particularly by Yakutsk and AGASA, is a very strong constraint to production and propagation models alike. Most of the scenarios proposed in the literature should produce a sizable anisotropy as either extragalactic luminous or dark matter is normally associated with the invoqued particle sources. If that is really the case, then it is our understanding of the topology and intensities of the intervening magnetic fields that is critically incomplete. The problem could reside inside the network of walls and filaments, in the interior of the large surrounding voids, or even in our nearby environment, namely the Galaxy halo.

In the present work we explore the possibility that radio ghosts, blobs of magnetized radio plasma remnant from past periods of activity in radio galaxies, being able to scatter UHECR in the intergalactic medium. Such a process could, in principle, degrade the direct incoming flux from the sources and build up a diffuse UHECR component large enough to be responsible for the observed degree of isotropy.

Our results show that, over the most conservative region of the radio ghost parameter space, such isotropization is not possible. This stands not from an inability of ghosts to scatter UHECR, but mainly from the fact that, under general conditions, ghosts should tend to cluster more strongly than the sources of the particles.

If, however, radio ghosts are able to buoy out into the surroundings of the dense large scale structures and into voids, while surviving the process, UHECR isotropy could be obtained in those cases in which the mean free path for interactions with ghosts is reduced below some few Mpc all over the propagation region ($`b_{gh}<0.5`$).

Finally, one can question the basic assumption of this work, that radio ghosts are able to survive for cosmological times as intact objects. If radio ghosts mix with the ambient medium a significant magnetization of the IGM should result, which can be sufficiently strong in order to explain observed magnetic fields in clusters of galaxies. These magnetic fields would have short coherence length, and therefore do not violate Faraday rotation measurements. The ability of the fields to deflect UHECR particles would be increased compared to the case of intact, non-mixing ghosts for the following reasons: the average deflection of particles in a ghosts vicinity is proportional to $`V_{\mathrm{gh}}B_{\mathrm{gh}}`$ (cross section $`V_{\mathrm{gh}}^{2/3}`$ times scattering angle, the latter is proportional to $`B_{\mathrm{gh}}`$ times $`V_{\mathrm{gh}}^{1/3}`$). For a fixed magnetic energy budget $`E_{B_{\mathrm{gh}}}=V_{\mathrm{gh}}B_{\mathrm{gh}}^2/(8\pi )`$ of the ghosts, this implies that the scattering efficiency scales with $`V_{\mathrm{gh}}^{1/2}`$ and therefore increases with increasing size. For fixed magnetic flux ($`B_{\mathrm{gh}}V^{2/3}`$) the scattering efficiency inreases with the size as $`V_{\mathrm{gh}}^{1/3}`$. This demonstrates that relaxing the hypothesis of compact ghosts makes this scenario even more powerful.

Acknowledgments
The authors wish to thank Dr. Simon White, Dr. Volker Springel, Dr. Avishai Dekel, Dr. Gerard Lemson and others for kindly making available the DM distribution from their large-scale-structure numerical simulation of the local Universe using a constrained realization of the initial perturbation field. It is a pleasure to acknowledge the constructive comments of the anonymous referee. GMT is partially supported by the Brazilian agencies FAPESP and CNPq. He thanks the MPA for the hospitality during his visit. TAE was temporary supported by the National Science and Engineering Research Council of Canada (NSERC).

Figure Captions:

Fig. 1. Energy distribution of radio ghosts multiplied by $`E_{\mathrm{gh}}^2`$ for better display. The lines are labeled with $`b_\nu `$, the exponent of the radio luminosity-jetpower correlation. The solid lines are derived from the pure luminosity evolution (PLE) radio luminosity function of Dunlop & Peacock (1990). The dashed lines from their radio luminosity function No. 2 from the MEAN-z data (RLF2). This latter luminosity function has an exponential cutoff at luminosities of $`10^{28}`$ Watt Hz-1. The curves are plotted for the range of radio luminosities between $`10^{23}`$ to $`10^{29}`$ Watt Hz-1. A horizontal line in this diagram corresponds to $`\alpha _{\mathrm{gh}}=2`$. E.g. the PLE with $`b_\nu =1.0`$ has a low energy index of $`\alpha _{\mathrm{gh}}=1.75`$.

Fig. 2. Example of a possible internal radio ghost magnetic field topology. The vectorial representation of the component of the magnetic field parallel to the plane of a slice, cut through the middle of the ghost, is shown. The length of the arrows scales linearly with magnetic field intensity. The magnetic field is tangent to the surface of the ghost. The amplitude of the fluctuations obeys a power law power spectrum of spectral index $`\xi =4`$. See text for further details.

Fig. 3. (a) Angular distribution of scattered UHECR for ghosts of different radius (in units of proton gyroradius) and values of $`\xi `$. The incident particles have a power law spectrum, $`dN/dEE^3`$ for E $`>4\times 10^{19}`$ eV. For comparison, the thick line shows an isotropic angular distribution. Forward scattering becomes dominant only as the ghost radius becomes comparable to the UHECR gyroradius. (b) Same as (a) but for protons at a fixed energy, E$`=10^{20}`$ eV. The same ghosts are used as in (a) and so their radii are down by a factor of approximately two in units of particle gyroradius.

Fig. 4. UHECR propagation without radio ghosts. Aitoff projection of the arrival probability density. Galactic coordinates are used with the antigalactic center at the center of the figure. UHECR sources are distributed according to nearby luminous matter (CfA catalog). Data points correspond to cosmic rays observed by AGASA, Haverah Park, Yakutsk and Volcano Ranch with $`E>4\times 10^{19}`$ eV. Clearly, the observations are more isotropic than what should be expected from the model.

Fig. 5. UHECR received at Earth for $`\alpha _{gh}=2.0`$, $`b_{gh}=2.0`$ and sources distributed as luminous matter.

Fig. 6 UHECR received at Earth for $`\alpha _{gh}=2.0`$, $`b_{gh}=1.0`$ and sources distributed as luminous matter.

Fig. 7 UHECR received at Earth for $`\alpha _{gh}=2.0`$, $`b_{gh}=0.5`$ and sources distributed as luminous matter. (a) UHECR that encounter no ghost from source to Earth; (b) UHECR that underwent at least one scattering off a RG. Comparing this figure with Fig. 4, 5 and 6, it can be seen that ghosts must have a more extended distribution than the sources ($`b_{gh}<0.5`$) in order to isotropize the UHECR flux at Earth. A tolerable (3-fold) increase in average ghosts energy density goes in the right direction but is not a determinant factor in itself.

Fig. 8 Amplitude of the first harmonic for UHECR with $`E>4\times 10^{19}`$ eV for $`\alpha =2.0`$ and $`b_{gh}=0.5`$. 1000 independent samples of 50 UHECR each were selected from the distribution in Fig. 7c with the same sensitivity in declination as AGASA’s. The amplitudes of the first harmonic obtained for these samples are shown in the histogram. As a comparison, the amplitudes obtained by AGASA, Haverah Park, Yakutsk, and Volcano Ranch \[Medina Tanco and Watson, 1999\] are also shown. It can be seen that, in order to isotropize the UHECR flux, the condition $`b_{gh}<0.5`$ must be fulfilled. That is, ghosts should be much more smoothly distributed in space than the galaxies that originated them.

Fig. 9 UHECR mean free path \[Mpc\] as a function of $`b_{\mathrm{gh}}`$, the bias between galaxies and ghost distributions, and gas density relative to the reference level. Ghosts are able to isotropize the UHECR flux only for $`b_{gh}\stackrel{}{<}0.5`$, for which the mean free path is smaller than a few Mpc both at voids and structures.

Fig. 10 UHECR received at Earth for $`\alpha _{gh}=2.0`$, $`b_{gh}=2.0`$ and sources distributed as $`\mathrm{\Lambda }`$CDM. Supergalactic coordinates are used in this figure. The strong excess in the arrival direction distribution corresponds to the Virgo cluster.

Fig. 11 UHECR received at Earth for $`\alpha _{gh}=2.0`$, $`b_{gh}=0.5`$ and sources distributed as $`\mathrm{\Lambda }`$CDM. Again, isotropization is obtained only for a ghost population that is much more extended than the UHECR sources. Supergalactic coordinates are used in this figure as in Fig. 10.

Fig. 12 Distribution of injected magnetic field strength as a function of background density, in the case that all radio ghosts completely mix with the ambient medium. The volume fractions of the different density regimes can be seen in Fig. 13.

Fig. 13 The volume fractions of the different density regimes ($`dV/d\mathrm{log}(n_{\mathrm{gas}}/n_{\mathrm{gas},0})/𝑑V`$) computed from the simulated dark matter distribution. The median density is marked. A count in cell box with size of $`1.5`$ Mpc was used.