text:
# Structure Shocks as a Source of Cosmic Rays in Clusters

## 1. Introduction

At least two lines of reasoning lead us to examine the properties of energetic particle acceleration at structure formation shocks. First, virtually all the gas in filaments, groups and clusters111To simplify the discussion below we do not distinguish groups from clusters. has at some time passed through one, or more probably, several structure shocks (e.g., Quilis et al 1998; Cen & Ostriker 1999; Miniati et al. 2000). Structure shocks involve very diffuse plasmas, so will be “collisionless”. Thus, they are likely to accelerate high energy particles (call them “cosmic rays” or CRs) through the so-called “diffusive shock acceleration”, provided there is a weak magnetic field present (e.g., Blandford & Eichler 1987). As discussed by Blasi at this meeting and outlined briefly below, clusters should be good reservoirs of CR protons, so that over time they will accumulate, possibly even to a level where they can contribute significantly to the cluster pressure (e.g., Berezinsky et al 1997). CRs have a softer equation of state than nonrelativistic thermal plasma, but CR protons are effectively immune to nonadiabatic cooling. Consequently, they can change the thermodynamic properties of the ICM. In that event it becomes important to include them in consideration of cluster dynamics, especially in cooling flows (see, e.g., the contribution by Ryu et al in these proceedings).

The second rationale for understanding particle acceleration at structure shocks comes from diffuse nonthermal emissions seen in at least some clusters. The most compelling such evidence, known for some time, is the existence of diffuse cluster radio halos and so-called “radio relic” sources (e.g., Feretti & Giovannini 1996; Feretti 1999). The detailed properties of these two classes of radio source are different in some respects, such as polarization and location inside the clusters, but both involve substantial volumes in their host clusters. The radio halos tend to be centered on the cluster cores and are unpolarized, while the relics are most likely to be found on the perimeters of clusters and can be highly polarized. Both types of radio sources result from synchrotron radiation by substantial populations of $`>`$$``$GeV electrons.

As discussed at this meeting by Fusco-Femiano, X-ray observations now show convincing evidence for diffuse, nonthermal hard X-ray emission in at least Coma and A2256 (e.g., Rephaeli et al 1999; Fusco-Femiano et al 2000). Again, this implies nonthermal electron populations. Here, however, as discussed by several other speakers at this workshop, the origin of that emission and the energy of the electrons are less certain. If the emission is nonthermal bremsstrahlung the electrons are only a little more energetic than the thermal electrons responsible for the soft X-ray continuum. On the other hand, one of the prime candidates for the hard X-ray excesses is inverse-Compton scattered CMB photons, again involving roughly GeV electrons. So far, the evidence requires only nonthermal electrons, since $`\pi ^0`$ decays coming from inelastic p-p collisions have not yet been detected. However, for particle acceleration in normal plasmas we should expect energetic hadrons as well. If the accelerators behave similarly to galactic CR accelerators the energy carried by hadronic CRs is likely to be one to two orders of magnitude greater than for electrons, in fact.

To facilitate what follows it may be useful here to review very quickly the key issue of CR longevity in clusters, since that largely controls the needs for extended accelerators of CRs. First, it has long been noted that electrons responsible for the observed radio halos have such short lifetimes to radiative losses that they cannot possibly fill a cluster from a single point source (e.g., Jaffe 1977). Supposing for example, that the cluster magnetic field is $`1\mu `$G, then electrons radiating at 1 GHz have Lorentz factors, $`\gamma 2\times 10^4`$. For Lorentz factors above a few hundred the dominant energy losses will come from inverse-Compton scattering in this case (e.g., Sarazin 1999), leading to lifetimes $`t_r4\times 10^{12}\gamma ^1`$ yrs, or about $`2\times 10^8`$ yrs in this case. It is simple to show that if we fix the observed radio frequency at 1 GHz, this is about the maximum lifetime of the relevant radiating electrons against combined inverse-Compton and synchrotron radiation.

Although in free flight relativistic electrons could cross a cluster in a few million years, diffuse radio emissions from the clusters and Faraday rotation through some clusters reveal the clear presence of weak magnetic fields (see the contributions by Feretti, Clarke, and Kronberg at this meeting, for example). The observations point to a tangled, perhaps turbulent field. MHD fluctuations will severely restrict their propagation. As an example consider Bohm diffusion, corresponding to scattering on saturated field fluctuations, so that a particle mean free path approximates the particle gyro radius; that is, $`D_B=\frac{1}{3}cr_g`$, for relativistic particles. This would allow GeV particles to diffuse less than a kiloparsec in a Hubble time. Bohm diffusion is a limiting case, and advective motions will surely carry CR electrons farther than that by a large factor even in the much smaller radiative lifetime of the electrons. But, the very small range of these electrons remains valid under any reasonable set of circumstances. The main consequence of this result is that CR electrons responsible for observed diffuse cluster radio and probably X-ray and predicted $`\gamma `$-ray emissions must be continuously replenished somehow. The observed hard X-ray luminosity of Coma exceeds $`10^{43}`$erg s-1, placing a lower bound on the replenishment rate in that cluster.

Note as well that hadronic CRs, and protons, in particular, are confined similarly by MHD wave scattering. In that case, on the other hand, radiative losses are negligible up to extremely high energies, as are losses due to inelastic collisions with the CMB and the cluster thermal plasma (e.g., Berezinsky et al 1997). So, the key consequence is that CR protons below about $`10^{15}`$eV = 1 PeV will be essentially locked to the ICM forever, once they are introduced.

The importance of shock heating to the ICM makes diffusive shock acceleration (DSA) an immediate candidate for production of cluster CRs, since DSA can transfer some tens of percents of the energy dissipated at the shock to the CR population (e.g., Blandford & Eichler 1987). For a typical cluster shock that would correspond to a CR energy input rate $`0.1\rho u^3R^210^{45}`$ erg/sec (e.g., Jones et al 2001). Over a Hubble time this could amount to as much as $`10^{63}`$ergs of CR energy. We emphasize from the start that our aim is not to argue for structure shocks as the only source of diffuse CRs in clusters, but that it is one very likely to be there and to be important. Others have estimated that radio galaxies can contribute a similar CR energy flux (e.g., En$`ß`$lin et al 1997). Intense starburst activity during galaxy formation has also been proposed (e.g., Völk et al 1996). Turbulent acceleration in the ICM may also play a significant role, especially as a means to reaccelerate CRs introduced by some other means (e.g., Jaffe 1977; Eilek & Wetherall 1999; Brunetti et al 2001).

In the remainder of this paper we will outline the role of structure shocks in determining the conditions in the ICM, review briefly the relevant properties of shock acceleration physics and discuss the results of some numerical simulations of structure formation that include treatment of CR acceleration via DSA, followed by advective transport and relevant energy loss mechanisms. Finally, since the ultimate test of such calculations is their match to observed cluster properties, we have computed “synthetic observations” of nonthermal emissions from simulated clusters, so describe some of those results here, as well. More extensive discussions of most of these issues are contained in our cited works, as well.

## 2. The Role of Shocks in Cluster Formation

Figure 1 illustrates the rich distribution of shock surfaces that are likely to be associated with large scale structure formation. For reference a volume rendered image of thermal X-ray emissivity is included to help locate clusters. These images are taken from one of our numerical simulations using a grid-based N-body/hydro scheme (Ryu et al 1993). The hydro part of the code uses a “TVD”, Riemann-solver scheme that cleanly captures even relatively weak shocks inside only 2-3 zones, which each span about 100 kpc in the simulation shown. The volume displayed is a little less than 20 $`h^1`$Mpc on a side at z = 0, and was extracted from a full simulation box 75 $`h^1`$ Mpc on a side. For this SCDM model $`\mathrm{\Omega }_M=1`$, $`\mathrm{\Omega }_B=0.13`$, $`h=0.5`$ and $`\sigma _8=0.6`$, which produces a cluster population at the current epoch consistent with observations.

Structure shocks are most commonly characterized as either “accretion shocks”, if they result from infall of diffuse, intergalactic material onto the perimeter of a cluster, or “merger shocks” if they result from collision between two clusters. A quick glance at Figure 1 shows that this is an overly simplified picture. For example, collisions between flows in filaments can lead to shocks, and accretion shocks are often hard to distinguish from merger shocks, given the complexity that accompanies the accumulation of mass in regions where clusters are forming. Dissipation at these shocks provides the basic heating of the ICM, although other processes, including feedback from star formation may also be important contributors.

Structure shocks come in a wide range of strengths, most generally indicated by Mach number. That depends especially on the temperature of the inflowing gas, since relative flow velocities tend generally to be of order $`10^3`$ km $`s^1`$. To the extent that the interacting material has been virialized the Mach numbers should be small, of course, while cooler material entering a cluster from a filament, for example, may pass through much higher Mach number shocks. Miniati et al (2000) and Miniati (2002a) have examined the histories of shock heating as revealed through cosmological simulations. Integrating over cosmic time they found that moderate strength shocks with Mach numbers roughly in the range $`2<M<7`$ contribute most to heating of the ICM. Such shocks are capable of accelerating CRs efficiently, so deserve close scrutiny in that context.

It is important to remember that, since CRs are effectively tied to the ICM up to pretty high energies, the integrated shock history of the ICM determines the character of the CR proton population (as well as their secondary products). So, one would not expect CRs produced in this way to be only associated with recent merger events, for example. On the other hand, also keep in mind that electrons lose energy so quickly that we can expect to see them only very close to where they have been accelerated or, in the case of secondary electrons, produced via decay of $`\pi ^\pm `$. In the context of structure shocks, we then should find electronic emissions either in association with current shocks or as secondaries associated with the accumulated CR proton population and the thermal baryons. Synchrotron radiation depends sensitively on magnetic field strength, so those emissions will also be heavily weighted towards regions of the strongest fields, wherever they form.

## 3. Diffusive Shock Acceleration

The DSA paradigm depends on the ability of energetic charged particles to pass through the dissipative layer in a shock mostly unimpeded, but also on the existence of sufficiently strong scattering to cause the particles to propagate diffusively upstream and downstream of the shock with respect to the bulk flow. That requires a weak magnetic field, which may not exist in primordial matter. On the other hand the first stars, galaxies and perhaps shocks should have seeded the ICM with magnetic flux. Fields $`>0.1\mu `$G should be adequate to accelerate particles in structure shocks to very high energies on reasonable time scales (e.g., Kang et al 1996). If a shock is planar on the scale of the particle scattering lengths and involves a simple, sharp transition, then the classic “test particle” solution for the momentum distribution of energetic particles is a power law, $`f(p)p^q`$, independent of any other details with a slope given by

$$q=\frac{3r}{r1}\frac{4M^2}{M^21}4(1+\frac{1}{M^2}),$$
(1)
where $`r`$ is the compression ratio of the shock and the arrow corresponds to a $`\gamma =\frac{5}{3}`$ gas. This leads to the well-known behavior, $`q4`$, in the strong shock limit.

For that spectral form the energy and pressure in this population diverge logarithmically as the momentum, $`p\mathrm{}`$, so it was recognized a long time ago that diffusive shock acceleration could in principle become quite efficient (e.g., Axford 1982). Already for $`M=3`$, equation 1 gives $`q=4.5`$, while for $`M=5`$ the slope is $`q=4.17`$. So, even moderate shocks can transfer substantial fractions of the kinetic energy flux into CRs if the seed population is adequate and the maximum energy is relativistic. The particle acceleration time scales with the diffusion time across the shock, $`t_d=D/u_s^2`$. Again employing a Bohm diffusion model we can estimate that the time to accelerate CRs to a PeV in a structure formation shock with speed $`u_s10^3`$km/sec can be as little as $`10^7`$yr when the magnetic field $`1\mu `$G. The seed CR population can come from previous sources or particles injected at the shock from the thermal plasma. The latter process is not well understood, but is certainly nonlinear and probably depends on the orientation of the magnetic field. Within our current understanding injection is sometimes modeled as a “thermal leakage” in which a small fraction of the downstream thermalized population is successful in crossing the postshock turbulence and the shock itself back into the approaching plasma. From there it can be returned downstream with an energy boost factor $`u_s/c`$ to begin the acceleration process. Theoretical models (e.g., Malkov 1998) and numerical simulations (Kang et al 2002) suggest that a simple way to model the injection physics is to match the thermal and isotropic nonthermal particle distributions at a particle speed that is several times the characteristic postshock thermal particle speed (Kang et al 2002), defining an “injection momentum”. In the simulations described below we will employ that model with injection momentum set to 2.6 times the postshock thermal peak. Further details can be found in Miniati 2001. This will generally lead to an injected fraction of CRs $`10^4`$, consistent with more sophisticated numerical simulations of the diffusive shock acceleration process (Kang et al 2002).

We note as well with this level of CR injection that the acceleration efficiency of shocks with $`M`$ 3-5 should be on the order of 10-30% (e.g., Kang et al 2002). That estimate is obtained from nonlinear simulations of CR modified shocks that include feedback between the CR population, bulk flow upstream and downstream of the shock as well as the injection process itself. Still, this is sufficiently small that shock modifications due to CR pressure gradients are modest and the CR momentum distribution is not greatly modified from the test particle form. This is convenient, since it allows us to make reasonable estimates of CR properties in clusters using test particle models for acceleration.

## 4. Cosmological Simulations Including CR Acceleration at Shocks

### 4.1. Methods

A number of authors have explored analytically the acceleration of CRs by one or more processes (e.g., Kang et al 1996; Eilek & Wetherall 1999; Blasi 2000; Brunetti et al 2001) in clusters. Yet, given the complex histories of clusters, the absence of stationary structures and likely wide ranges of properties among individual clusters, numerical simulations offer a powerful alternative tool to gain a more complete picture. With this in mind we have carried out several structure formation simulations that incorporate an efficient numerical scheme for transport of CRs, including diffusive acceleration at shocks, as well as adiabatic losses for protons and radiative and Coulomb losses for electronic components. The scheme is Eulerian, uses a finite volume approach to advection in coordinate and momentum space, and models the momentum distribution as a piecewise power law function to allow a coarse momentum grid. To conserve computational effort the simulations described consider CRs only with energies below 1 PeV, so that we can reasonably assume acceleration at shocks is instantaneous compared to dynamical time steps and that spatial diffusion of CRs in smooth flows can be ignored. Dynamical feedback of CRs is neglected; that is, the CR population is treated in the test particle limit. Thus, CRs emerging from shocks take a spectral form governed by equation 1 for injected CRs or whenever that spectrum is flatter than the incident CR spectrum.

Electrons are split into a “primary” population derived from thermal leakage at shocks and a “secondary” population derived from $`\pi ^\pm `$ decays following inelastic p-p collisions. Those same collisions produce $`\gamma `$-rays via $`\pi ^0`$ decay, and that distribution is also computed. Since injection of primary electrons at shocks is not well understood, we simply scale it as a fraction, $`R_{e/p}`$ of the proton injection described above. Using galactic CRs for a reference, we may expect $`R_{e/p}0.01`$. The proton, primary and secondary electron populations are all evolved separately. Our full numerical method is an extension of that described in Jones et al (1999) and has been explained in detail by Miniati (2001). The Eulerian cosmology code is the same as that employed in the simulation shown in Figure 1, and outlined in §2.

### 4.2. Some Results: CR Properties

Detailed discussions of results from these simulations can be found in Miniati et al (2001a,b) and Miniati (2002a,b). Here we mention briefly a few of the salient results, especially as they relate to themes of this meeting.

A single structure formation simulation leads to many clusters that have formed naturally from random initial fluctuations. While we pay a price in reduced resolution when we simulate a large volume, we gain an unbiased sample of objects whose properties can be explored from multiple perspectives. The simplest question we can ask from these simulations is an estimate of the total CR population we can expect. Figure 2 shows estimates of the fraction of individual cluster core pressures contributed by CR protons resulting from diffusive acceleration at structure shocks. These data were obtained from a “concordance model” $`\mathrm{\Lambda }`$CDM simulation with $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_B=0.04`$, $`h=0.67`$ and $`\sigma _8=0.9`$ (Miniati 2002a,b). (Similar results were reported for an analogous SCDM simulation, Miniati et al 2001a.) The simulation shown used a computational box 50 $`h^1`$ Mpc on a side in comoving coordinates with $`512^3`$ grid zones and $`256^3`$ dark matter particles.

The mean CR pressure contribution in these clusters is about 20% in the cores, but there is a very wide scatter. These are first estimates, since the model for diffusive CR acceleration was applied in the test particle limit and the injection rates are uncertain. Nonetheless, the numbers are consistent with expectations from more sophisticated nonlinear simulations of individual shocks in the range of Mach numbers mentioned earlier that are most important for heating of cluster ICMs. There is no clear trend of $`P_{CR}/P_{tot}`$ with cluster temperature. The scatter seen in Figure 2 is real and reflects the influence of varied shock histories among individual clusters. We emphasize again that $`P_{CR}`$ depends on the time integrated shock history of the plasma in a cluster rather than a single event, such as a recent merger.

Reflecting that fact, inspection of the spatial distribution of CR protons in the simulations shows them to be distributed roughly similar to the thermal plasma (Miniati et al 2001a). However, since the CRs are continuously injected at accretion shocks, which lie mostly outside cluster cores, the CR distributions are broader than the X-ray emitting thermal gas. From these behaviors we can immediately anticipate a key property of secondary electron CRs as well as associated $`\gamma `$-rays from $`\pi ^0`$ decays. Both will be concentrated in cluster cores, since they will result from a convolution of CR and thermal proton distributions. That is distinctly different from the distribution of primary CR electrons and their emissions, which are concentrated close to contemporary shocks. This last point was mentioned earlier and comes from the very short energy loss time scales of electrons. While shocks can and do penetrate into cluster cores, they are more common, stronger and perhaps more obvious outside the cores.

### 4.3. Some Results: Nonthermal Emissions

The spatial and momentum distributions of the CRs can be used to compute expected emissions from such processes as inverse-Compton scattering of the CMB, bremsstrahlung and, since our simulations also include passive magnetic fields, synchrotron radio emissions. From inelastic p-p scattering products we include emissions from secondary $`e^\pm `$ and $`\gamma `$s from $`\pi ^0`$ decay, as well.

Since there has been much recent discussion about nonthermal X-rays and $`\gamma `$-rays from clusters, we focus our attentions there. If the primary electron to proton ratio in cluster CRs is similar to that in galactic CRs, so that primary electrons constitute something like a percent of the total primary CRs, then the average $`\gamma `$-ray luminosities in our simulated clusters are produced in comparable measure by inverse-Compton emission from primary electrons and $`\pi ^0`$ decay. The primary electron inverse-Compton contribution is especially important in extended regions outside cluster cores, where strong shocks are most likely, while $`\pi ^0`$ decays can dominate at high energies in cluster cores.

Figure 3 illustrates the X-ray to $`\gamma `$-ray spectral energy distribution computed from one $`T_x4`$ keV cluster in the simulation of Miniati (2002a,b). IC emission of primary and secondary electrons (actually $`e^\pm `$), nonthermal bremsstrahlung and $`\gamma `$-rays from $`\pi ^0`$ decay are included. The emissions were integrated inside a spherical volume of radius 5 $`h^1`$ Mpc.

Primary electron inverse-Compton and $`\pi ^0`$ emissions dominate as mentioned, while nonthermal bremsstrahlung is largely unimportant. Note at high energies that the spectra of inverse-Compton emissions from both electronic components and $`\pi ^0`$-decay $`\gamma `$s are similar over a substantial range. All reflect the spectra of primary protons, with adjustments for radiative losses in the electrons and for pion generation and decay cross sections in the other case. Below a GeV, however, the $`\pi ^0`$ contribution drops due to production threshold and phase space limits, while above a a few TeV the primary electron spectrum cuts off since electron energy loss rates dominate DSA acceleration rates.

Examination of the spatial distributions of emissions of the different contributions uncovers some important details of how structure shocks may reveal themselves, and suggests some strategies for testing models for particle acceleration. Figure 4 separates the emissions from the clusters shown in Figure 3 into contributions from the cluster core (2 Mpc diameter) and from the “outskirts” defined as an annulus between 1 Mpc and 3 Mpc radius. For this plot the cluster was assumed to be at the distance of the Coma cluster, taken as d = 105 Mpc. The figure illustrates how emissions associated with secondary CRs products concentrate towards cluster cores, while emissions from from the outer regions are dominated by primary electrons. Spectral distributions are distinctly different in the two cases, so it should be possible in principle to identify emissions from the hadronic CR component in the cluster cores. That observational challenge is also currently underway in efforts to establish the hadronic CR population in galactic supernova remnants (e.g., Aharonian et al 1994). The right panel in the figure shows the integrated flux from the different parts of the selected cluster in comparison to various detection limits.

The above discussion also is relevant to our ability to understand the origins of diffuse radio halos and relic sources. Recall that the former are found primarily in cluster cores, while the latter tend to lie in cluster outskirts. In our simulations core regions do emit synchrotron emissions that are produced mostly by secondary electrons. As discussed in Miniati et al 2001b the radio luminosities and polarization properties are similar to observed halos, and the radio luminosity is a very steep function of cluster temperature, consistent with the comparative rarity of observed halos. The main handicap in trying to explain radio halos entirely in terms of secondary decay products is that this explanation does not appear to explain naturally observations of radio spectral steepening (e.g., Deiss et al 1997). However, more sophisticated models including other contributions to the electron population or its acceleration might solve this problem. The other issue sometimes raised for secondary electron radio halo models is whether the required electron population is consistent with limits to the inverse-Compton X-ray and $`\gamma `$-ray fluxes in Coma, for example. For very weak fields the inverse-Compton flux from the electron population required to account for radio halo in Coma would become excessive (Blasi & Colafrancesco 1999; Miniati et al 2001a). On the other hand if the magnetic field is in excess of a few $`\mu `$G the required electron population to explain the synchrotron source is reduced to a level leading to inverse-Compton luminosities consistent with current limits on the $`\gamma `$-ray flux, and below the expected $`\pi ^0`$ flux at high energies, in fact.

The simulations also produce diffuse radio sources resembling radio relics in some of the same outer regions where the inverse-Compton emissions from primary electrons are prominent (Miniati et al 2001b). These are immediate postshock volumes where primary electrons still remain energetic and where magnetic fields are relatively strong. The simulations predict that these regions should be highly polarized, since the local magnetic fields tend to become aligned with the shock faces. These properties all seem consistent with those seen in radio relic sources.

Finally we illustrate in Figure 5 the individual inverse-Compton $`\gamma `$-ray luminosities from primary accelerated electrons computed in a $`\mathrm{\Lambda }`$CDM simulation of Miniati (2002a,b). The solid line is a least squares power law fit to the simulation data, which has a slope 2.6, in agreement with what is expected from cluster scaling relations (Miniati 2002a). The simulated cluster luminosities were integrated inside a 5 $`h^1`$ Mpc radius spherical volume. For the illustration a relative electron shock injection rate of $`R_{e/p}=0.01`$ was used. Computed fluxes scale directly with that parameter. Several EGRET upper limits for nearby clusters are shown in the figure for comparison. The estimated luminosities from the simulations are consistent with the current $`\gamma `$-ray non-detections of clusters, but they also suggest that $`\gamma `$-ray fluxes may be large enough to be seen by the next generation of experiments.

## 5. Summary

Shocks are a ubiquitous consequence of cosmic structure formation, and they play an essential role in heating of cluster media . Virtually all of the gas in cluster media has been processed by at one or more shocks of at least moderate strength. Since these shocks involve highly tenuous ionized media, they are collisionless in nature, so will not fully thermalize the plasmas passing through them. One likely consequence is the acceleration of relativistic particles, or cosmic rays. We have begun to explore through numerical simulations the roles that particle acceleration in structure shocks may play. Our current conclusions are:

$``$The shocks that are primarily responsible for heating cluster ICMs can be efficient particle accelerators, possibly generating nonthermal proton pressures on the order of 10% or more of the total virial pressure in cluster cores.

$``$Cluster ICMs are very good reservoirs for energetic protons, so the cosmic ray populations there reflect the full history of the cluster medium more than any single event, such as a recent merger.

$``$Relativistic electrons at most energies have loss lifetimes so short that they must either be accelerated locally or be secondary products from energetic hadronic cosmic rays in order to have populations great enough to account for detectable nonthermal emissions.

$``$There are two main regions for production of nonthermal radiation in clusters; the X-ray bright core and the outskirts where strong shocks are most likely.

$``$Inverse-Compton emission and $`\pi ^0`$ decays dominate the production of $`\gamma `$-rays in typical clusters. Inverse-Compton from primary electrons dominates in the outskirts, provided electrons are injected in proportion to ions comparably to the galactic cosmic rays; that is, so that a fraction of a percent of the energy flux through shocks is transferred to electrons.

$``$In cluster cores $`\gamma `$-ray emission above a few hundred MeV should be dominated by $`\pi ^0`$ decays, whereas inverse-Compton emission from secondary electrons dominates in these regions at lower $`\gamma `$-ray energies.

$``$Primary and secondary electrons may also contribute substantially to nonthermal radio synchrotron emissions in clusters. Primary electron emissions are confined to volumes close to contemporary shocks, so should be seen mostly in cluster outskirts, contributing to radio relic sources. Secondary electronic emissions should again be concentrated in cluster cores, contributing to radio halos.

#### Acknowledgments.

TWJ and EJH have been supported by NASA through grant NAG5-10774, by the NSF through grant Ast00-71167 and by the University of Minnesota Supercomputing Institute. FM acknowledges support from the Max-Planck-Gesselschaft Rechenzentrum in Garching. DR and HK were supported by grant No. R01-1999-00023 from the Korea Science & Engineering Foundation. We thank the organizers of this meeting for their hard work and for providing a provocative and illuminating program.

## References

Aharonian, F., A., Drury, L. O’C., & Völk, H. J. 1994, A&A, 285, 645

Axford, W. I., Leer, E., & McKenzie, J. F. 1982, A&A, 111, 317

Berezinsky, V. S., Blasi, P. & Ptuskin, V. S. 1997, ApJ, 487, 529

Blandford, R. D. & Eichler, D. 1987, Phys. Rept., 154, 1

Blasi, P. 2000, ApJ, 532, L9

Blasi, P. & Colafrancesco, S. 1999, Astropart. Phys., 12, 169

Brunetti, G., Feretti, L. & Giovannini, G. 2001, MNRAS, 320, 365

Cen, R. & Ostriker, J. P. 1999, ApJ, 514, 1

Deiss, B. M., Reich, W., Lesch, H., & Wielebinski 1997, A&A, 321, 55

Dennison, B. 1980, ApJ, 239, L93

Eilek, J. & Wetherall, J. C. 1999, in Diffuse Thermal and Relativistic Plasma in Galaxy Clusters, ed: H. Böhringer, L. Feretti & P. Schuecker, Garching: MPE, p 249

En$`ß`$lin, T. A., Biermann, P. L., Kronberg, P. P. & Wu, X.-P. 1997, ApJ, 477, 560

Feretti, L. 1999, in Diffuse Thermal and Relativistic Plasma in Galaxy Clusters, ed: H. Böhringer, L. Feretti & P. Schuecker, Garching: MPE, p 3

Feretti, L. & Giovannini, G. 1996, in IAU Symposium 175, Extragalactic Radio Sources, ed: R. D. Ekers, C. Fanti, & L. Pardrielli, Dordrecht: Kluwer, p 333

Fusco-Femiano, R. et al, 2000, ApJ, 534, L7

Jaffe, W. J. 1977, ApJ, 212, 1

Jones, T. W., Ryu, D. & Engel, A. 1999, ApJ, 512, 105

Jones, T. W., Miniati, F., Ryu, D., & Kang, H. 2001, in High Energy Gamma-Ray Astronomy, ed: F. A. Aharonian & H. J. Völk, New York: AIP, p 436

Kang, H., Ryu, D. & Jones, T. W. 1996, ApJ, 456, 422

Kang, H., Jones, T. W. & Gieseler, U. D. J. 2002, ApJ(in press)

Malkov, M. A. 1998, Phys. Rev. E, 58, 4911

Miniati, F. 2001, Comp. Phys. Comm, 141, 17

Miniati, F. 2002a, astro-ph/0203014, MNRAS(submitted)

Miniati, F. 2002b MNRAS(submitted)

Miniati, F., Jones, T. W., Kang, H. & Ryu, D. 2001b, ApJ, 562, 233

Miniati, F., Ryu, Kang, H. & Jones, T. W. 2001a, ApJ, 559, 59

Miniati, F., Ryu, D., Kang, H., Jones, T. W., Cen, R. & Ostriker, J. P. 2000, ApJ, 542, 608

Quilis, V., Ibáñez, & Sáez, D. 1998, ApJ, 502, 518

Rephaeli, Y., Gruber, D. E. & Blanco, P. 1999, ApJ, 511, L21

Ryu, D., Ostriker, J. P., Kang, H. & Cen, R. 1993, ApJ, 414, 1

Sarazin, C. L. 1999, ApJ, 520, 529

Völk, H. J., Aharonian, F. A., & Breitschwerdt, D. 1996, Space Science Reviews, 75, 279


format: text

text:
# Behaviour of 𝑆-Wave poles Near Threshold and the Scalar Meson Nonet Below 1 GeV

## Introduction

In a theory with quarks and mesons one can study strong interactions through meson-meson scattering. Observed spectra may be interpreted in terms of quark-antiquark or more complicated systems . Necessarily, spectra must be confirmed in several different experiments, and compiled in spectroscopic tables. Then, for a lowest-order approximation of strong interactions, confinement models may be constructed. Their usefulness can be measured by the model’s achievements when adjusting its parameters to experiment.

For further refinements of strong-interaction models, there are several possible directions. One way is to reanalyse experiments in order to adjust the data to new insights , in particular when popular models  are not capable of reproducing specific experimental results . Another way is to adapt the model by including more sophisticated interactions , though often at the cost of more model parameters. One may also just ignore the spectroscopy which is distilled from experiment, and, instead of comparing resonance positions and widths to information from the spectroscopic tables, compare the model’s predictions directly, if possible, to experimental scattering cross sections and phase shifts .

Here we will discuss the eternally disputed  low-lying nonet of $`S`$-wave poles in meson-meson scattering cross sections, which are predicted by a four-parameter model . It should thereby be noted that the model’s parameters have already been fitted to the $`J^P=1^{}`$ $`c\overline{c}`$ and $`b\overline{b}`$ spectra, as well as to $`P`$-wave meson-meson scattering data .

First, let us outline the motivation for our work. For the interaction in the vicinity of a resonance in meson-meson scattering, one may consider quark-exchange or quark-pair-creation processes, as depicted in Fig. (1). When the intermediate $`q\overline{q}`$ system is close enough to a genuine bound state of confinement, then the system will resonate, resulting in a resonance in meson-meson scattering. Another picture for the same phenomenon, shown in Fig. (2), is that a mesonic quark-antiquark system (M) gets a self-energy correction from a virtual meson loop.

Either picture describes the same physical situation, namely a mesonic resonance or bound state , but in a rather different way. Our aim is to merge both pictures in one model.

## 1 Quark exchange

In the quark-exchange picture we obtain a resonance in the particular partial-wave meson-meson-scattering cross section which agrees with the quantum numbers of the intermediate $`q\overline{q}`$ system.

Such a process may be described by scattering phase shifts of the form

$$\text{cotg}\left(\delta _{\mathrm{}}(s)\right)\frac{E_R\sqrt{s}}{\mathrm{\Gamma }_R/2},$$
(1)
where $`E_R`$ and $`\mathrm{\Gamma }_R`$ represent the central invariant meson-meson mass and the resonance width, respectively.

However, formula (1) is a good approximation for the scattering cross section only when the resonance shape is not very much distorted and the width of the resonance is small. Moreover, the intermediate state in such a process is essentially a constituent $`q\overline{q}`$ configuration that is part of a confinement spectrum (also referred to as bare or intrinsic states), and hence may resonate in one of the eigenstates. This implies that the colliding mesons scatter off the whole $`q\overline{q}`$ confinement spectrum of radial, and possibly also angular excitations, not just off one single state . Consequently, a full expression for the phase shifts of formula (1) should contain all possible eigenstates of such a spectrum as long as quantum numbers are respected. Let us denote the eigenvalues of the relevant part of the spectrum by $`E_n`$ ($`n=0`$, $`1`$, $`2`$, $`\mathrm{}`$) and the corresponding eigenstates by $`F_n`$. Then, following the procedure outlined in Ref. , we may write for the partial-wave phase shifts the more general expression

$$\text{cotg}\left(\delta (s)\right)=\left[I(s)\underset{n=0}{\overset{\mathrm{}}{}}\frac{\left|_n\right|^2}{\sqrt{s}E_n}\right]^1\left[R(s)\underset{n=0}{\overset{\mathrm{}}{}}\frac{\left|_n\right|^2}{\sqrt{s}E_n}\mathrm{\hspace{0.33em}1}\right].$$
(2)
In $`R(s)`$ and $`I(s)`$ we have absorbed the kinematical factors and the details of two-meson scattering, and moreover the three-meson vertices. The details of formula (2) can be found in Ref. .

For an approximate description of a specific resonance and in the rather hypothetical case that the three-meson vertices have small coupling constants, one may single out, from the sum over all confinement states, one particular state (say number $`N`$), the eigenvalue of which is nearest to the invariant meson-meson mass close to the resonance. Then, for total invariant meson-meson masses $`\sqrt{s}`$ in the vicinity of $`E_N`$, one finds the approximation

$$\text{cotg}\left(\delta (s)\right)\frac{\left[E_N+R(s)\left|_N\right|^2\right]\sqrt{s}}{I(s)\left|_N\right|^2}.$$
(3)
Formula (3) is indeed of the form (1), with the central resonance position and width given by

$$E_RE_N+R(s)\left|_N\right|^2\text{and}\mathrm{\Gamma }_R\mathrm{\hspace{0.33em}2}I(s)\left|_N\right|^2.$$
(4)
In experiment one observes the influence of the nearest bound state of the confinement spectrum, as in classical resonating systems. Nevertheless, formula (3) is only a good approximation when the three-meson couplings are small. Since the coupling of the meson-meson system to quark exchange is strong, the influence of the higher- and lower-lying excitations is not negligible.

In the other hypothetical limit, namely of very large couplings, we obtain for the phase shift the expression

$$\text{cotg}\left(\delta (s)\right)\frac{R(s)}{I(s)},$$
(5)
which describes scattering off an infinitely hard cavity.

The physical values of the couplings come out somewhere in between the two hypothetical cases. Most resonances and bound states can be classified as stemming from a specific confinement state . However, some structures in the scattering cross section stem from the cavity which is formed by quark exchange or pair creation . The most notable of such states are the low-lying resonances observed in $`S`$-wave pseudoscalar-pseudoscalar scattering .

## 2 Meson loops

¿From the discussion in Sec. (1) one may conclude that to lowest order the mass of a meson follows from the quark-antiquark confinement spectrum. It is however well-known that higher-order contributions to the meson propagator, in particular those from meson loops, cannot be neglected.

Virtual meson loops give a correction to the meson mass, whereas decay channels also contribute to the strong width of the meson. One obtains for the propagator of a meson the form

$$\mathrm{\Pi }(s)=\frac{1}{s\left(M_{\text{confinement}}+\mathrm{\Delta }M_{\text{meson loops}}\right)^2},$$
(6)
where $`\mathrm{\Delta }M`$ develops complex values when open decay channels are involved.

For the full mass of a meson, all possible meson-meson loops have to be considered. A model for meson-meson scattering must therefore include all possible inelastic channels as well. Although in principle this could be done, in practice it is not manageable, unless a scheme exists dealing with all vertices and their relative intensities. In Ref.  relative couplings have been determined in the harmonic-oscillator approximation assuming $`{}_{}{}^{3}P_{0}^{}`$ quark exchange. However, further kinematical factors must be worked out and included.

## 3 The spectrum

The complete model consists of an expression for the $`K`$ matrix, similar to formula (2), but extended to many meson-meson scattering channels, several constituent quark-antiquark channels, and more complicated transition potentials , which at the same time and with the same set of four parameters reproduces bound states, partial-wave scattering quantities, and the electromagnetic transitions of $`c\overline{c}`$ and $`b\overline{b}`$ systems .

The $`K`$ matrix can be analytically continued below the various thresholds, even the lowest one, with no need of redefining any of the functions involved, in order to study the singularities of the corresponding scattering matrix. Below the lowest threshold, these poles show up on the real $`\sqrt{s}`$ axis, and can be interpreted as the bound states of the coupled system, to be identified with the stable mesons. For the light flavours one finds this way a nonet of light pseudoscalars, i.e., the pion, Kaon, eta, and eta. For the heavy flavours, the lowest-lying model poles can be identified with the $`D(1870)`$, $`D_s(1970)`$, $`\eta _c(1S)`$, $`J/\psi (1S)`$, $`\psi (3686)`$, $`B(5280)`$, $`B_s(5380)`$, $`\mathrm{{\rm Y}}(1S)`$, $`\mathrm{{\rm Y}}(2S)`$, and $`\mathrm{{\rm Y}}(3S)`$.

Above the lowest threshold, the model’s partial-wave cross sections and phase shifts for all included meson-meson channels can be calculated and compared to experiment, as well as the inelastic transitions. Singularities of the scattering matrix come out with negative imaginary part in the $`\sqrt{s}`$ plane. To say it more precisely: out of the many singularities in a rather complex set of Riemann sheets, some come close enough to the physical real axis to be noticed in the partial-wave phase shifts and cross sections. In fact, each meson-meson channel doubles the number of Riemann sheets, hence the number of poles. Consequently, with ten scattering channels one has for each eigenvalue of the confinement spectrum 1024 poles in 1024 Riemann sheets, out of which usually only one contains relevant poles in each $`\sqrt{s}`$ interval in between the thresholds. Those can be identified with the known resonances, like the $`\rho `$ pole in $`\pi \pi `$ scattering, or the $`K^{}`$ pole in $`K\pi `$ scattering. However, there might always be a pole in a nearby Riemann sheet just around the corner of one of the thresholds, which can be noticed in the partial-wave cross section. The study of poles is an interesting subject by itself .

Once the four model parameters have been adjusted to the experimental phase shifts and cross sections, the pole positions can be determined and compared to tables for meson spectroscopy. But usually no new information is gained from such a comparison. Here we will report on the singularity structure of the lowest poles in $`S`$-wave meson-meson scattering.

## 4 Scattering-matrix poles

In the hypothetical case of very small couplings for the three-meson vertices, we obtain poles in the scattering matrix that are close to the eigenvalues of the confinement spectrum. Let us denote by $`M_1`$ and $`M_2`$ the meson masses, and by $`\mathrm{\Delta }E`$ the difference between the complex-energy pole of the scattering matrix and the energy eigenvalue, $`E_N`$, of the nearby state of the confinement spectrum. Using formula (4), we obtain

$$\mathrm{\Delta }E\left\{R(s)iI(s)\right\}\left|_N\right|^2.$$
(7)
We may distinguish two different cases:

(1) $`E_N>M_1+M_2`$ (above threshold),

(2) $`E_N<M_1+M_2`$ (below threshold).

When the nearby state of the confinement spectrum is in the scattering continuum, then $`\mathrm{\Delta }E`$ has a negative imaginary part and a real part, since both $`R(s)`$ and $`I(s)`$ of formula (7) are real, and $`I(s)`$ is moreover positive. The resonance singularity of the scattering matrix which corresponds to this situation is depicted in Fig. (3).

Notice that the resonance singularity is in the lower half of the complex-energy plane (second Riemann sheet), as it should be.

When the nearby state of the confinement spectrum is below the scattering threshold, then $`\mathrm{\Delta }E`$ has only a real part, since $`I(s)`$ turns purely imaginary below threshold, whereas $`R(s)`$ remains real. The bound-state singularity of the scattering matrix corresponding to this situation is depicted in Fig. (4).

Note that the bound-state pole is on the real axis of the complex-energy plane, as it should be.

## 5 Threshold behaviour

Near the lowest threshold, as a function of the overall coupling constant, $`S`$-wave poles behave very differently from $`P`$\- and higher-wave poles. This can easily be understood from the effective-range expansion  at the pole position. There, the cotangent of the phase shift equals $`i`$. Hence, for $`S`$ waves the next-to-lowest-order term in the expansion equals $`ik`$ ($`k`$ represents the linear momentum related to $`s`$ and to the lowest threshold). For higher waves, on the other hand, the next-to-lowest-order term in the effective-range expansion is proportional to $`k^2`$.

Poles for $`P`$ and higher waves behave in the complex $`k`$ plane as indicated in Fig. (5$`b`$). The two $`k`$-plane poles meet at threshold ($`k=0`$). When the coupling constant of the model is increased, the poles move along the imaginary $`k`$ axis. One pole moves towards negative imaginary $`k`$, corresponding to a virtual bound state below threshold on the real $`\sqrt{s}`$ axis, but in the wrong Riemann sheet. The other pole moves towards positive imaginary $`k`$, corresponding to a real bound state.

For $`S`$-wave poles, the behaviour is shown in Fig. (5$`a`$). The two $`k`$-plane poles meet on the negative imaginary $`k`$ axis. When the coupling constant of the model is slightly increased, both poles continue on the negative imaginary $`k`$ axis, corresponding to two virtual bound states below threshold on the real $`\sqrt{s}`$ axis. Upon further increasing the coupling constant of the model, one pole moves towards increasing negative imaginary $`k`$, thereby remaining a virtual bound state for all values of the coupling constant. The other pole moves towards positive imaginary $`k`$, eventually passing threshold ($`k=0`$), thereby turning into a real bound state of the system of coupled meson-meson scattering channels. Hence, for a small range of hypothetical values of the coupling constant, there are two virtual bound states, one of which is very close to threshold. Such a pole certainly has noticeable influence on the scattering cross section.

Although we are not aware of any experimental data that could confirm the above-described threshold behaviour of poles, we suspect this to be possible for atomic transitions in cavities. Unfortunately, it does not seem likely that in the near future similar processes can be studied for strong coupling.

## 6 The low-lying nonet of $`𝑺`$-wave poles

The nonet of low-lying $`S`$-wave poles behave as described in Sec. (4), with respect to variations of the model’s overall coupling constant. However, they do not stem from the confinement spectrum, but rather from the cavity. For small values of the coupling, such poles disappear into the continuum, i.e., they move towards negative imaginary infinity , and not towards an eigenstate of the confinement spectrum as in Fig. (3).

In Fig. (6) we study the hypothetical pole positions of the $`K_0^{}(730)`$ pole in $`K\pi `$ $`S`$-wave scattering. The physical value of the coupling constant equals 0.75, which is not shown in Fig. (6). A figure for smaller values of the coupling constants can be found in Ref. . The physical pole in $`K\pi `$ isodoublet $`S`$-wave scattering, related to experiment , comes out at $`727i263`$ MeV in Ref. . Here we concentrate on the threshold behaviour of the hypothetical pole movements in the complex $`k`$ and $`\sqrt{s}`$ planes. Until they meet on the axis, which is for a value of the coupling constant slightly larger than 1.24, we have only depicted the right-hand branch.

In the left figure we observe how the poles arrive on the imaginary $`k`$ axis, and then continue to move along that axis. One of the poles moves upwards, initially describing a virtual bound state, and crossing the real $`k`$ axis for a value of the coupling constant slightly larger than 1.30. The other pole moves downwards, remaining a virtual bound state for further increasing values of the coupling constant.

In the right figure the same pole has been depicted in the $`\sqrt{s}`$ plane. Here, the situation is more confusing, since the pole positions are in the same interval of energies. The pole which corresponds to the one moving downwards along the imaginary $`k`$ axis moves to the left on the real $`\sqrt{s}`$ axis. Its positions as a function of the coupling constant are indicated by solid circles. The pole which moves upwards along the imaginary $`k`$ axis initially moves towards threshold and then turns back, following the former pole, but in a different Riemann sheet. The positions of the latter pole are indicated by open circles. In the inset we try to better clarify its motion. Notice that, since we took 0.14 GeV and 0.50 GeV for the pion and the Kaon mass, respectively, we end up with a threshold at 0.64 GeV.

It is interesting to notice that in a recent work of M. Boglione and M.R. Pennington  also a zero-width state is found below the $`K\pi `$ threshold in $`S`$-wave scattering. Here, we obtain such a state for unphysical values of the coupling.

In Fig. (7) we have depicted the movement of the $`a_0(980)`$ pole in $`S`$ wave $`I=1`$ $`KK`$ scattering (threshold at 1.0 GeV) on the upwards-going branch. One observes a very similar behaviour as in the case of $`K\pi `$ scattering, but with two important differences, to be described next.

The $`K_0^{}(730)`$ poles meet on the real $`\sqrt{s}`$ axis only 16 MeV below threshold (see Fig. 6), and for a value of the coupling constant which is well above the physical value of 0.75, whereas the $`a_0(980)`$ poles meet 238 MeV below threshold, when the coupling constant only equals 0.51. At the physical value of the coupling constant, the $`a_0(980)`$ pole is a real bound state some 9 MeV below threshold.

But there is yet another difference. Whereas the $`K\pi `$ channel represents the lowest possible scattering threshold for the $`K_0^{}(730)`$ system, $`KK`$ is not the lowest channel for the $`a_0(980)`$. In a more complete description, at least all pseudoscalar meson-meson loops should be taken into account. One of these is the $`\eta \pi `$ channel, which has a threshold well below $`KK`$. Consequently, upon including the $`\eta \pi `$ channel in the model, the pole cannot remain on the real $`\sqrt{s}`$ axis, but has to acquire an imaginary part in a similar way as shown in Fig. (3). In Ref.  we obtained a resonance-like structure in the $`\eta \pi `$ cross section, representing the physical $`a_0(980)`$. The corresponding pole came out at $`962i28`$ MeV.

For the $`f_0(980)`$ system the situation is very similar to that of the $`a_0(980)`$. Assuming a pure $`s\overline{s}`$ quark content , we obtain for the variation of the corresponding pole in $`KK`$ $`I=0`$ $`S`$-wave scattering a picture almost equal to the one shown in Fig. (7). However, only in lowest order the $`KK`$ channel could be considered the lowest threshold for the $`f_0(980)`$ system. In reality $`s\overline{s}`$ also couples to the nonstrange quark-antiquark isosinglet through $`KK`$, and hence to $`\pi \pi `$ . This coupling is nevertheless very weak, which implies that the resulting pole does not move far away from the $`KK`$ bound state. In Ref.  we obtained a resonance-like structure in the $`\pi \pi `$ cross section representing the physical $`f_0(980)`$. The corresponding pole came out at $`994i20`$ MeV.

At lower energies, we found for the same cross section a pole which is the equivalent of the $`K_0^{}(730)`$ system, but now in $`\pi \pi `$ isoscalar $`S`$-wave scattering. This pole at $`470i208`$ MeV may be associated with the $`\sigma `$ meson, since it has the same quantum numbers and lies in the ballpark of predicted pole positions in $`\sigma `$ models (for a complete overview of $`\sigma `$ poles, see Ref. ).

We do not find any other relevant poles in the energy region up to 1.0 GeV.

## 7 Conclusions

We have shown that the poles of the $`a_0(980)`$ and $`f_0(980)`$ belong to a nonet of scattering-matrix poles. The lower-lying isoscalar pole and the isodoublet poles in the complex-energy plane have real parts of 0.47 GeV and 0.73 GeV, respectively, and imaginary parts of 0.21 GeV resp. 0.26 GeV. Whether these poles represent real physical resonances  is not so relevant here. Important is that the $`a_0(980)`$ and $`f_0(980)`$ are well classified within a nonet of scattering-matrix poles with very specific characteristics, different from those of the poles stemming from confinement, like the confinement-ground-state nonet of scalar mesons $`f_0(1370)`$, $`a_0(1450)`$, $`K_0^{}(1430)`$, and $`f_0(1500)`$. The latter poles vary as a function of the coupling constant exactly the way indicated in figure (3). For vanishing coupling constant they end up on the real $`\sqrt{s}`$ axis at the positions of the various ground-state eigenvalues of the confinement spectrum, which are the light-flavour $`{}_{}{}^{3}P_{0}^{}`$ states at 1.3 to 1.5 GeV .

The low-lying $`S`$-wave poles related to the cross sections in the $`f_0(470)`$, $`K_0^{}(730)`$, $`f_0(980)`$, and $`a_0(980)`$ regions move to negative imaginary infinity in the $`\sqrt{s}`$ plane for decreasing values of the coupling. The pole positions for the physical value of the coupling are well explained by their threshold behaviour. Whether or not these poles have large imaginary parts, leading to large widths and strong resonance distortion, depends in a subtle way on the thresholds and couplings of the various relevant scattering channels .

As to the nature of the light scalar mesons, which has recently been discussed in Refs. , we can only remark that in a many-coupled-channel model each of the channels contributes to the states under the resonance, not just one specific channel.

Comparing to experiment cross sections and phase shifts that follow from modelling strong interactions is a particularly useful strategy, since then the model gains independence with respect to the criteria for omitting experimental results from the spectroscopic tables . Moreover, this strategy has another advantage. Phase shifts do not easily follow from confinement models, nor does Nature provide a bound-state spectrum of mesons. Confinement models and spectroscopic tables go well together, but do not shed sufficient light on strong interactions. Both can be avoided when following the strategy outlined above. Nonetheless, it would be useful to have some coordination in the availability of experimental results on meson-meson scattering.

Acknowledgement: We wish to thank Frieder Kleefeld for useful discussions on the distribution of $`S`$-matrix poles in the complex $`E`$ and $`k`$ planes.

This work has been partly supported by the Fundação para a Ciência e a Tecnologia of the Ministério da Ciência e da Tecnologia of Portugal, under contract numbers POCTI/35304/FIS/2000, and CERN/P/FIS/40119/2000.


text:
# Rearrangement collisions between gold clustersSupported by the Fondo Nacional de Investigaciones Científicas y Tecnológicas (FONDECYT, Chile) under grant #8990005 and #1010988 and DIPUC-Chile.

## 1 Introduction

The study of nanostructures has recently attracted widespread interest among theoretical and experimental physicists and chemists, and because of its many applications has also come to the forefront of technology rev-heer . On the theoretical side ab initio procedures are now capable of providing incisive insights into the properties of these systems. In addition, novel and sophisticated nanostructure fabrication, manipulation and measurement techniques have given impetus to experiment, and reliability to a large amount of experimental data clev97 ; scha97 ; tayl92 ; koga98 ; spas98 ; palp98 . Moreover, the technological applications on a variety of devices has strongly stimulated activity in the nanostructure field since they can be considered as building blocks of novel nanostructured materials and devices whet96 ; andr96 ; mirk96 ; aliv96 .

In particular, metallic clusters provide an interesting subject of study for at least two reasons: i) clusters constitute intermediate systems between isolated atoms and molecules, on the one extreme, and bulk solids on the other (i.e. they constitute genuine mesoscopic systems); and ii) often they exhibit an interesting phenomenology of their own. Gold clusters have received widespread attention during the last two decades, both experimentally clev97 ; scha97 ; tayl92 ; koga98 ; spas98 ; palp98  and theoretically zhao94 ; hand94 ; garz96 ; habe97 ; wang01 ; garz98 ; belb98 ; barn99 ; sole00 ; wils00 ; liyi00 ; hakk00 ; gron00 ; rome02 .

In principle one does expect ab initio procedures to be the definitive tool to handle this type of systems; however, it is not always feasible (or at least practical) to implement ab initio calculations. In fact, since we are interested in structures with a fairly large number of atoms, or arranged in several different interacting nanostructures, ab initio computations first become very time consuming and, in the end, impractical. Much the same happens when attempting to obtain a detailed description of the long time evolution of small systems, or to describe their properties when subject to a large variety of external conditions.

Here we intend to develop an adequate description of colliding gold clusters, a phenomenon which falls into one or more of the categories described in the preceding paragraph. In this paper we employ mainly classical molecular dynamics, an alternative computation scheme that has proven to be quite reliable sole00 . Classical molecular dynamics (MD) is a valid option in this case, since it allows to considerably reduce the computation time (relative to MD ab initio calculations), and/or to increase substantially the number of particles that can be handled. Obviously, there is a price to pay; as examples of this cost we mention that the method is at best semi-classical, that the detailed electronic structure is ignored, and that phenomenological potentials, adjusted to fit bulk properties, are used. In spite of these shortcomings the MD technique can reliably be used to compute the properties dominated by the ionic contribution, which is the case for the phenomena treated in this paper. This is especially true when MD is checked, in some physical limit, against first principle calculations. An alternative, but equally convenient procedure, is to close in on a solution, for example for geometrical optimization via MD, to be followed by first principle calculations. The synergy between ab initio and MD thus allows to significantly reduce the resources that are required and to expand the set of problems amenable to treatment.

More than a decade ago the bombardment with gold clusters of metallic surfaces was investigated experimentally matt86  and theoretically eche90 ; shap90 ; shap91 . Moreover, the dramatic energy accommodation that occurs in cluster-cluster collisions, which is crucial to understand the growth mechanism during the early stages of particle formation, was investigated around the same time by Blastein et al. blas92 . We focus our interest on the dynamics of gold cluster collisions. In particular we study the structure and symmetry of small Au clusters, and the dynamics of atom–cluster and cluster–cluster collisions, for different values of the impact parameter $`b`$ and as a function of center of mass energy $`E`$.

This paper is organized as follows: after this Introduction we describe, in Sec. 2, both the ab initio and EA methods employed in the computations. In Sec. 3 we provide results, for a large variety of cases, of the implementation of the codes, for the cluster structures, their symmetries and the dynamics of the collision process. Finally, in Sec. 4 we draw conclusions and close the paper.

## 2 Simulation method

### 2.1 Ab-initio method

In order to asses the quality of the semi-empirical embedded atom (EA) procedure used in the context of our classical molecular dynamics simulations, we start contrasting EA against first principles calculations. Thus, ab initio geometrical optimization was carried out within the Car-Parinello approach carpar85 , in the framework of the density functional theory, using gradient corrections in the PBE implementation PBE . Gradient corrected functionals have been adopted in recent theoretical studies of geometrical optimization of metallic clusters, mainly because they are more accurate than the local density functional, even though there still is some controversy in the literature on the ground state geometry of small Au clusters habe97 ; part90 ; semi97 . The calculations were performed only at the $`\mathrm{\Gamma }`$-point of the Brillouin zone, using norm conserving pseudopotentials trou91 . The wavefunctions were expanded in plane waves, with an energy cutoff of 60 Ry. We have explicitly checked that, with this energy cutoff, the structural properties of our system are well converged. The box used in the calculations was always at least three times larger than the cluster diameter. The results obtained, as well as their EA counterparts, are given in Fig. 2 and in Table Rearrangement collisions between gold clustersthanks: Supported by the Fondo Nacional de Investigaciones Científicas y Tecnológicas (FONDECYT, Chile) under grant #8990005 and #1010988 and DIPUC-Chile.. The ab initio optimized geometries are illustrated in Fig. 1. The agreement between ab initio, EA and Wilson and Johnston wils00  (WJ) procedures is quite satisfactory for large clusters (N$``$10) and provides a reasonable basis to trust the EA calculations that constitute the core of the present paper.

### 2.2 Embedded atom method

The interatomic interaction between gold atoms is modeled using a semi-empirical (EA) potentials daw84 ; foil86 . On the basis of these EA potentials we obtain the average binding energy per atom, $`E_b`$, which is later minimized to yield the optimal cluster geometry. The latter is achieved using a Monte Carlo procedure, for which we adopted as starting configurations, for the different cluster sizes, the geometries found by Wilson and Johnston wils00  (WJ). Once the optimal geometry is established several static properties, like nearest neighbor distances and angles, and the average coordination number, are readily evaluated.

After the various different clusters are properly characterized we use classical molecular dynamics (MD) to simulate the cluster–cluster scattering process. Many body EA semi-empirical potentials are used throughout. To integrate the equations of motion we implement the Verlet velocity algorithm, with a 1 femtosecond time step. Since the collision fragments heat up as a consequence of the scattering process it is necessary to cool them down; to do so we simply rescale the temperature by 1 % every 1000 steps, during a total of 100 000 time steps, to reach a final energy of 63 % of the initial one. Finally, the collision fragments are carefully scrutinized to extract the physical information we are looking for.

## 3 Results and Discussion

### 3.1 Structures

The first issue we address is to check the reliability of our MD procedure as compared with alternative methods. In Table Rearrangement collisions between gold clustersthanks: Supported by the Fondo Nacional de Investigaciones Científicas y Tecnológicas (FONDECYT, Chile) under grant #8990005 and #1010988 and DIPUC-Chile. we compare the average binding energy $`E_b`$ of the lowest energy configurations, and the average nearest neighbor distance $`R`$, of clusters built with different numbers of Au atoms. On the one hand we have ab initio results obtained within the framework of density functional theory, and on the other the empirical potential results of WJ wils00  (who used the Murrell–Mottram murr90 ; cox99  potential) as well as the EA values that we obtained.

It is quite apparent that the EA estimates for $`E_b`$ differ considerably from the ones found by WJ, and also with the more trustworthy ab initio results. However for small clusters the geometrical parameters are not so satisfactory. This small N (where N is the number of atoms in the cluster) error margin is not unexpected, since the EA potential has been adjusted to fit bulk properties and cannot be expected to fully succeed in systems where N is tiny. However, and also as expected, the situation improves as N increases, which is precisely what is borne out by Table Rearrangement collisions between gold clustersthanks: Supported by the Fondo Nacional de Investigaciones Científicas y Tecnológicas (FONDECYT, Chile) under grant #8990005 and #1010988 and DIPUC-Chile., where we display results for the range of clusters sizes from N=3 to 12. The same trend is observed in Fig. 2, where we plot $`E_b`$ versus N. We notice that ab initio calculations were performed for N = 4, 5 and 6 planar and three dimensional structures. It is apparent that as the cluster size increases, specially when the number of atoms is larger than 10, that both the ab initio and EA values tend towards the bulk value.

The average nearest neighbor distances $`R`$ that we compute also exhibit larger errors than those of WJ. In spite of the fact that we employ the WJ cluster configurations as the starting point for our calculations, but in which we use a binding energy obtained from a different potential, we derive geometrical structures which differ from those of WJ. However, once again, increasing N yields compatible results. For example, for a 3-atom cluster the difference in nearest neighbor distances amounts to 15%, but already for a 6-atom cluster it reduces to only 8%. These average distances are illustrated in Fig. 3, where we plot $`R`$ as a function of the number of atoms in the cluster, N. In contrast with Fig. 2 the plot of Fig. 3 is not smooth, but shows abrupt variations between two successive values of N. Despite this roughness the tendency of the WJ and our EA plots is to approach each other as N$`1`$. Moreover, our results are in good agreement with the ab initio ones obtained by Wang wang01 .

### 3.2 Symmetries

In addition to the binding energies and interatomic distances the cluster symmetry is a relevant characteristic and, in the context of gold cluster topologies, the Jahn-Teller effect is also an important element. Ab initio calculations predict a C2v symmetry for Au3 and Au4, while EA yields C2v and D2d, and the Murrell–Mottram potential used by WJ yields D3h and Td symmetries, respectively. These differences are quite apparent in the pair correlation function $`g(r)`$ plotted in Fig. 4 .

The second difference in the binding energy is defined by

$$\mathrm{\Delta }_2E_b(N)=2E_b(N)E_b(N1)E_b(N+1),$$
(1)
and gives an indication of the stability of a cluster with respect to disproportionation rev-heer , as well as its ionic hardness. A plot of $`\mathrm{\Delta }_2E_b(N)`$ as a function of N is given in Fig. 5, where we observe a good agreement of our EA values with those reported by Wilson et al. wils00 , both for the position and magnitude of the hardness peaks. The maxima (minima) of $`\mathrm{\Delta }_2E_b(N)`$ imply that there are values for which it is more difficult (easier) to add an atom to the cluster. Moreover, as N $`1`$ the plot becomes quite smooth, which constitutes an indication that the cluster can incorporate an additional atom without major hindrance.

### 3.3 Collisions between a single Au atom and a Au cluster

Next we report the results of our simulations of the collisions between a single Au atom and a variable size Au cluster, for several values both of the initial per atom energy $`E`$ and of the impact parameter $`b`$. The precise details that describe the collision process are as follows: at time $`t=0`$ the centers of mass of the atom and the cluster are placed along the $`x`$-axis, $`\pm 10`$ Å away from the origin, respectively. The atom is located a distance $`b`$ away from the $`x`$-axis, along the $`y`$-axis on the $`xy`$ plane. The principal symmetry axes of the various clusters are aligned perpendicular to the direction of motion, that is parallel to $`y`$. The impact parameter $`b`$ is varied between 0 (head-on collision) and 7 Å; the latter corresponds to the distance where the interaction potential effectively vanishes, since the average radius of a cluster with $`12N14`$ is less than 3 Å and the interaction is cut off at 5.5 Å. The projectile energy is varied between 0.1 and 1.5 eV per atom, in steps of 0.2 eV (8 values in all). The maximum energy $`E=1.5`$ eV per atom corresponds to approximately one half of the cluster binding energy. After the collision takes place the resulting fragments are stabilized, by gradual cooling through the rescaling of the internal velocities. Finally, we analyze the data characterizing the collision fragments for several special cases.

We consider two categories: i) the scattering of a single gold atom against clusters with N$`=12`$, $`13`$ and $`14`$ atoms, which is dealt with here; and, ii) the scattering of a variable size projectile (N$`=12`$, $`13`$ y $`14`$) on a variable size target with a similar number of Au atoms, which is presented in 3.4. Throughout we use the concepts of low and high energies, and small and large impact parameters. Low energies are defined to be in the range 0.1$`E`$0.7 eV and large within 0.7$`E`$1.5 eV. Similarly, small impact parameters cover the range $`0b3`$ Å, and large is defined as $`3b7`$ Å. The upper bounds on $`E`$ and $`b`$, as mentioned above, are related to the binding energy of the cluster and the distance at which the projectile does not interact with the target, respectively.

One atom on 12 Because of the rich variety of results, and to facilitate their understanding by the reader, we have chosen to illustrate them by means of figures. In particular Fig. 6 describes the one gold atom collision with a 12 atom cluster. In this and ensuing figures the columns correspond to different values of the energy per atom (in eV), while the rows correspond to several different impact parameters in Å units. Each entry characterizes a collision on the basis of the following symbols: single atoms, dimers and trimers are represented by dots, two dots joined by a line, and three dots that form a triangle, respectively. When several of these collision fragments are generated we denote their number by a factor in front of the corresponding symbol. If the fragment contains four or more atoms the symbolic representation is a circle. Finally, if the number of atoms in the target is not altered after the collision, we stress this fact representing it by a square with the original number of target atoms in its interior, as long as N$`4`$.

In Fig. 6 it is noticed that for low energies and small impact parameter ($`0b3`$ Å), projectile and target fuse into a 13 atom cluster. When the energy is increased to $`0.9E1.5`$ eV, keeping the impact parameter fixed, coalescence is observed for a few cases, while more often 11 or 9 atom clusters, and one or two dimers, respectively, are generated. For larger impact parameters, $`4b7`$, fusion is present in the low energy region, but for a few cases target and projectile remain unaltered. Finally, for large $`b`$ and $`E`$ values we observe a few cases of coalescence, some 11 atom clusters plus a dimer and many instances (denoted by squares) in which projectile and target size do not change.

A particularly interesting scattering process, illustrated in Figs. 7, occurs for an energy of 4 eV and impact parameter 1.8 $`<b<`$ 2.2 Å. When the projectile approaches the target (Fig. 7a) it attracts the nearest (lowest in Fig. 7b) atom, but without removing it from the cluster. This generates a large energy transfer, which in turn induces large amplitude vibrations in the “lowest” atom, as the projectile leaves the scene (Fig. 7c). As a consequence of these large amplitude vibrations this particular atom does overcome the energy barrier and ends up at the center of the cluster (Fig. 7c). We have estimated a lower bound for this energy barrier of 0.03 eV. Quite remarkably this lower symmetry cluster has an energy slightly smaller than the fully symmetric 12 atom cluster we accepted above as the stable configuration.

In fact, the original 12($`I_h`$) cluster has a binding energy $`E_b=3.03`$ eV, while the energy of the less symmetric configuration 12(C5v) equals $`E_b=3.09`$ eV, which is 2% lower. Thus, we are faced with the question of why this lower energy configuration was not obtained in the minimization process we reported above. The explanation of this apparent contradiction is related to the fact that the energy barrier the atom has to overcome, to shift to the center of the cluster, is rather large and cannot be achieved in a minimization process that starts with an icosahedron and allows only small displacements from the original equilibrium positions. It is worth mentioning that this asymmetric structure was also asigned minimum energy by Wilson and Johnston wils00 . However, they argue that the bond compression of the 12(C5v) structure generates a repulsion strong enough to destabilize it in favor of the icosahedron.

One atom on 13 When a Au atom collides with the “magic number” N=13 cluster the results are not much different from the previous 1 on 12 case. However, there is a larger number of coalescence cases, many instances of dimer formation, and also instances are observed where there is no change in the number of atoms of projectile and target. For example, when $`E=0.9`$ eV, there is fusion for $`b2`$, dimer plus a 12 atom cluster generation for $`b=3`$, fusion for $`b=4`$, no change $`b=5`$, and fusion for $`b=6`$ and $`b=7`$ Å.

One atom on 14 In this case the results fall into only two categories: either target and projectile fuse or they remain unaltered. Fusion prevails for values of $`b5`$ Å. For $`b>5`$ Å there is coalescence only for very small energies $`E0.3`$, and no change in the number of atoms in projectile and target for larger energies. This small $`b`$ large $`E`$ behavior can be understood as follows: for low impact parameter the collision gives rise to violent cluster vibrations and deformations, which precludes the trapping of the projectile. Instead, larger $`b`$ values induce less drastic cluster deformations and sufficient attraction to fuse projectile and target.

One atom on 12: cluster rotation In all preceding cases the cluster symmetry axis was taken to be parallel to the $`y`$ axis, that is, perpendicular to the initial projectile velocity. Now we align the cluster symmetry axis parallel to the $`y`$ axis. A qualitative change is observed: single atoms as a result of the collision. In fact, for low energies and small $`b`$ there is fusion. However, for $`E>0.7`$ eV several alternatives are observed: either coalescence, or 8, 9 and 11 atom cluster formation accompanied by the creation of dimers and single atoms. When $`b>3`$ Å we obtain 11 atom clusters plus a dimer or a pair of isolated atoms. Finally, for the largest $`b=7`$ Å value no atomic reordering is observed. However, we notice that in general the overall structure of the results is equivalent to the perpendicularly oriented cluster impact discussed above, which allows us to concentrate on cluster collisions without paying much attention to their relative spatial orientation.

### 3.4 Collisions between two Au clusters

Now we turn to the problem of the collision of two clusters of a few, but in general different, number of atoms.

One dimer on a 12 atom cluster

Fusion is obtained for low energies ($`E0.7`$ eV) and practically all values of the impact parameter $`b`$. In many instances the end result is a dimer and a 12 atom cluster, which are the outcome after a complex dynamic interaction that, finally, yields a reconstruction into the two original clusters. For $`E>0.7`$ eV and small $`b`$ we notice a diversity of results: single atoms plus clusters of 5, 6, 7, 8, 10 and 11 atoms are obtained, which reflects the fact that the more complex the projectile the richer the variety of collision fragments.

One 12 atom cluster on another 12 atom cluster  Fusion is observed only for low energies ($`E=0.1`$ eV per atom) over the whole range of impact parameters $`0b7`$ Å. For large energies and small impact parameter collisions, either large fragments plus a couple of dimers or trimers, or dimers and single atoms are generated. For large $`E`$ and large $`b`$, collisions without cluster size rearrangement are predominant.

One 13 atom cluster on a 12 atom cluster  Again there is coalescence for $`E`$ 0.1 eV and $`0b7`$ Å, while for large $`E`$ and large $`b`$ collisions without rearrangement predominate. In the $`0b`$ 5 Å and $`E0.5`$ eV region large and medium size fragments plus trimers, dimers and single atoms are produced. Total break up of the cluster is seen almost exclusively for small $`b`$ and large $`E`$ collisions.

One 14 atom cluster on a 12 atom cluster  Again here the results are quite similar to the previous ones. For small $`b`$ there is fusion, while for large $`b`$ and large $`E`$ the collision does not modify the size of the colliding clusters. Moreover, for a large region of parameter space ($`0b5`$ Å and $`0.5E1.5`$ eV) a whole variety of fragments does result: large fragments plus dimers, medium size fragments and finally, in for the largest energies, total cluster breakup into small pieces.

One 13 atom cluster on a 13 atom cluster

The collision of two “magic number” clusters yields the rich variety of results illustrated in Fig 8. It is readily noticed that there are many notable exceptions to the general trends observed in the preceding cases, and which are only present for this particular case. This is specially noticeable in the upper right of the figure, which illustrates the parameter values for which dimers, trimers and clusters of 6, 7, 8, 9, 11, 12, 14, 15, 16, 17 and 24 atoms are generated. However, again fusion is observed for low energy collisions ($`E0.1`$ eV) for all $`b`$ values, as well as no rearrangement collisions for large impact parameters, which are the final outcome of a complex dynamic interaction that in the end yields a reconstruction into two clusters equal to the original ones.

One 14 atom cluster on a 13 atom cluster

A low energy coalescence region, as well as a large impact parameter zone where after the collision projectile and target rearrange into their original structures, is again obtained. Also, and just as in previous cases, for large $`E`$ and small $`b`$ a diversity of fragments (dimers, trimers, and 5, 6, 7, 8 and 9 atom clusters) do result. In addition there are several cases with a single Au atom exchange between projectile and cluster, such that the original (14, 13) pair is converted into a (15, 12) one. However, in this same $`b`$ and $`E`$ region when three, instead of two, collision fragments are generated, always a dimer (and not a single atom) is created. This dimer originates in the 14 atom projectile, yielding in the end a 12 atom and a 13 atom cluster plus a dimer, as collision fragments.

## 4 Conclusions

The dynamics of the collision process of gold clusters has been investigated by means of classical molecular dynamics in combination with embedded atom (EA) potentials. First, the reliability of the EA potentials was confirmed by comparison with ab initio values, finding that EA is in good agreement for the cluster sizes we considered, an agreement which improves as N increases. Next, structural characteristics and the symmetry of the various Au clusters were obtained and contrasted with published results wils00 .

Several type collisions were investigated, finding regions of coalescence, fragmentation and scattering. Which of these outcomes actually occurs depends mainly on the values of the projectile energy $`E`$ and the impact parameter $`b`$. Coalescence is dominant for low energies ($`E<0.7`$ eV) and small impact parameters ($`b<3`$ Å). Simple scattering, with no change in the size and structure of the colliding clusters, prevails for large $`E`$ and $`b`$ values. For large energies and small impact parameters fragmentation and scattering are generally the case. For large $`E`$ and large $`b`$ scattering is the most probable outcome. When the cluster does breaks up the main collision products, apart from large fragments, are dimers. On the other hand cluster coalescence provides a viable mechanism to generate larger cluster sizes. It is also of interest that the collisions themselves turn out to be rather insensitive to the relative orientation of the projectile and target main symmetry axes, and that cluster collisions can generate metastable structures, which are usually not accesible due to the existence of a potential barrier.

| Cluster | Method | Symmetry | $`E_b`$ \[eV\] | $`R`$ \[Å\] | $`v`$ \[Å3\] |
| --- | --- | --- | --- | --- | --- |
| 3 atoms | EA | C2v | 2.405 | 1.937–2.438 | 3.317 |
| | CP | C2v | 2.182 | 2.607–2.753 | 4.225 |
| | WJ | D3h | 1.759 | 2.780 | 4.133 |
| 4 atoms | EA | D2d | 2.313 | 1.904–2.874 | 3.996 |
| | CP (a) | D2d | 2.234 | 2.570–3.120 | 5.857 |
| | CP (b) | D2h (planar) | 2.517 | 2.700 | 2.346 |
| | CP (c) | C2v (planar) | 2.504 | 2.648 | 29.759 |
| | WJ | Td | 2.178 | 2.793 | 5.004 |
| 5 atoms | EA | Oh | 2.685 | 2.507 | 7.287 |
| | CP (a) | Oh | 2.464 | 2.801 | 4.672 |
| | CP (b) | C2v (planar) | 2.652 | 2.715 | 7.706 |
| | WJ | Oh | 2.382 | 2.801 | 4.290 |
| 6 atoms | EA | Oh | 2.825 | 2.536 | 5.767 |
| | CP (a) | Oh | 2.523 | 2.834 | 8.663 |
| | CP (b) | D3h (planar) | 2.870 | 2.712 | 4.365 |
| | WJ | Oh | 2.574 | 2.789 | 7.667 |
| 8 atoms | EA | D2d | 2.937 | 2.577 | 5.512 |
| | CP | D2d | 2.857 | 2.837 | 7.385 |
| | WJ | D2d | 2.736 | 2.790 | 6.792 |
| 10 atoms | EA | D4d | 2.999 | 2.594 | 9.794 |
| | CP | D4d | 2.914 | 2.829 | 36.225 |
| | WJ | D4d | 2.837 | 2.778 | 30.854 |
| 12 atoms | EA (a) | Ih | 3.027 | 2.604 | 15.190 |
| | EA (b) | C5v | 3.089 | 2.693 | 21.142 |
| | CP (a) | C5v | 2.879 | 2.896 | 33.032 |
| | CP (b) | D4h | 2.962 | 2.814 | 8.780 |
| | WJ | Ih | 2.886 | 2.776 | 18.395 |