text:
# 1 Introduction

## 1 Introduction

Magnetism in strongly correlated electron systems has been a subject of extensive studies for a long time. In many examples, the energy gain from electron kinetic energy term is a crucial driving force for stabilizing magnetically symmetry broken states. Even in the Mott insulating state, this is equally true. Although the itinerancy of electrons is lost by the Coulomb interaction in the Mott insulator, electrons still find an optimized way to gain the energy from the kinetic energy term through a virtual process, where the strong coupling expansion becomes an adequate description under totally suppressed coherent motion of electrons. Anderson’s mechanism of the superexchange interaction thus proposed  has long been one of the most fruitful concepts in the condensed matter physics. The double exchange mechanism proposed by de Gennes , Zenner  and Anderson and Hasegawa  is another example, where a ferromagnetic (FM) metal is stabilized by the kinetic energy gain under a strong Hund’s rule coupling between degenerate orbitals.

Superconductivity, on the other hand, has been analyzed mainly from a different view since the dramatic success of BCS theory . Various types of pairing interactions mediated by bosonic excitations have been considered to be responsible for the mechanisms. To say nothing of the electron-phonon interaction, roles of AFM and FM paramagnons and excitons have been studied in various contexts. In these proposals, the superconducting ground state is the consequence of pair formation by the mediated attraction in the presence of coherent and metallic quasiparticles.

In this report, we discuss in more general perspective the importance of the kinetic origin for various symmetry breakings. We particularly discuss the kinetic mechanism of the superconductivity under the suppressed coherence of quasiparticles near the Mott insulator. The mechanism of superconductivity we discuss in this paper is conceptually different from the conventional magnetic mechanism discussed in the spin fluctuation theories and in the $`t`$-$`J`$ model, because the origin is not in the mediated bosons but in the direct kinetic process as in the superexchange for the mechanism of antiferromagnetism. Although the superexchange interaction is the only process for the energy lowering by the kinetic energy term in the Mott insulating phase, this constraint is lost in the presence of doped carriers and various other processes may equally or even more importantly contribute to the kinetic energy gain. The existence of various other processes is in fact easily confirmed by the strong coupling expansion of the Hubbard-type models, where in the second order, namely in the order of $`t^2/U`$, it not only generates the superexchange interaction but also other two-particle terms including the pair hopping processes. Here in the strong coupling expansion, the onsite Coulomb repulsion $`U`$ is assumed to be larger than the kinetic transfer term $`t`$. These second-order terms usually have secondary importance in the presence of the kinetic energy term in the first order of $`t`$ in the normal metal. However, as similarly to the superexchange interaction in the Mott insulating phase, they become the primary origin of the kinetic energy gain even in metals when the single-particle processes are suppressed by some reason but the amplitude of the second-order terms are retained. In the study on the strong-coupling expansion of the d-p model, it was shown that the d-wave pair-hopping process appears in the third order, although the superexchange interaction appears only in the fourth order . In the next section we first discuss the origin of the unusual suppression of coherent single-particle motion in the proximity of the Mott insulator before considering the mechanism of superconductivity.

High-temperature cuprate superconductors show a variety of unusual properties in the normal state . Among others, we concentrate on two remarkable properties widely observed in the cuprates for the purpose of examining the relevance of the kinetic mechanism as the driving mechanism of the high-Tc superconductivity. One of the remarkable properties in the cuprates is the flat dispersion in single-particle excitations. Angle resolved photoemission spectra (ARPES) in Y and Bi based high-Tc cuprates show an unusual dispersion which is far from weak correlation picture . The dispersion around $`(\pi ,0)`$ and $`(0,\pi )`$ is extremely flat beyond the expectations from usual van Hove singularities. The flat dispersion also shows rather strong damping.

The other remarkable property we discuss is the pseudogap phenomenon observed in the underdoped region . It is observed both in spin and charge excitations in which gap structure emerges from a temperature $`T_{\mathrm{PG}}`$ well above the superconducting transition point $`T_\mathrm{c}`$. The gap structure is observed in various different probes such as NMR relaxation time, the Knight shift, neutron scattering, tunneling, ARPES, specific heat, optical conductivity, and DC resistivity. The ARPES  data have revealed that the pseudogap starts growing first in the region around $`(\pi ,0)`$ and $`(0,\pi )`$ from $`T=T_{\mathrm{PG}}`$ much higher than $`T_c`$. Therefore, the pseudogap appears from the momentum region of the flattened dispersion and it is likely that the mechanism of the pseudogap formation is deeply influenced from the underlying flatness. The superconducting state itself also shows a dominant gap structure in this flat spots, $`(\pi ,0)`$ and $`(0,\pi )`$, due to the $`d_{x^2y^2}`$ symmetry. In fact, the pseudogap structure above $`T_c`$ appears continuously to merge into the $`d_{x^2y^2}`$ gap below $`T_\mathrm{c}`$. To understand the superconducting mechanism and the origin of the high transition temperatures, a detailed understanding of the physics taking place in the flat dispersion region is required.

The emergence of the flat dispersion around $`(\pi ,0)`$ and $`(0,\pi )`$ has also been reported in numerical simulation results rather universally in models for strongly correlated electrons such as the Hubbard and $`t`$-$`J`$ models on square lattices . As we see in the next section, this criticality is interpreted from a strong proximity of the Mott insulator where strong electron correlation generates suppressed dynamics and coherence.

## 2 Kinetic Pairing Derived from Electron Differentiation

In a simple picture, the correlation effects emerge as the isotropic mass renormalization, where the Coulomb repulsion from other electrons makes the effective mass heavier. This effect was first demonstrated by Brinkman and Rice  in the Gutzwiller approximation and refined in the dynamical mean field theory .

In the numerical results on a square lattice, as discussed above, the correlation effects appear in more subtle way where the electrons at different momenta show different renormalizations. When the Mott insulator is approached and the doping concentration becomes small, the mass renormalization generally becomes stronger. However, once the renormalization effect gets relatively stronger in a part of the Fermi surface, it is further enhanced at that part in a selfconsistent fashion because the slower electrons become more and more sensitive to the correlation effect. This generates critical differentiation of the carriers depending on the portion of the Fermi surface.

On square lattices, the stronger renormalization happens around ($`\pi ,0`$) and ($`0,\pi `$). A part of this anisotropic correlation effect concentrating near ($`\pi ,0`$) and ($`0,\pi `$) is intuitively understood from the carrier motion under the background of AFM correlations. As we see a real space picture in Fig.1a, the carrier motion in the diagonal directions does not disturb the correlations due to the parallel spin alignment, while the motion in horizontal and vertical directions strongly disturbs the AFM backgrounds as expressed as the wavy bonds and the motion itself is also disturbed as a feedback. Such strong coupling of charge dynamics to spin correlations causes flattening and damping of electrons around ($`\pi ,0`$) and ($`0,\pi `$), but not around the diagonal direction ($`\pm \pi /2,\pm \pi /2`$). The anisotropic renormalization effect eventually may generate a singularly flat dispersion on particular region of the Fermi surface, which accepts more and more doped holes in that region due to the enhanced density of states.

The transition to the Mott insulator is then governed by that flattened part, since the carriers reside predominantly in the flat region. The criticality of the metal-insulator transition on the square lattice is thus determined from the doped carriers around the flat spots, $`(\pi ,0)`$ and $`(0,\pi )`$. The hyperscaling relation becomes naturally satisfied because singular points on the momentum space govern the transition. In fact, the hyperscaling relations are numerically supported in various quantities and shows agreements with experimental indications. For example, the electronic compressibility critically diverges as $`\kappa 1/\delta `$ with decreasing doping concentration $`\delta `$ while the Drude weight is unusually suppressed as $`D\delta ^2`$ . The coherence temperature (the effective Fermi temperature) is also scaled as $`T_F\delta ^2`$ and indicates unusual suppression. In more comprehensive understanding, all the numerical data are consistent with the hyperscaling relations with a large dynamical exponent $`z=4`$ for the metal-insulator transition . Such large exponent opposed to the usual value $`z=2`$ for the transition to the band insulator is derived from the slower electron dynamics even at $`T=0`$ generated by the flat dispersion.

We, however, should keep in mind that the relaxation time of quasiparticles and the damping constant of magnetic excitations do not have criticality at the transition point to the Mott insulator. A general remark is that the relaxation time is critical only in the case of the Anderson localization transition and not in the case of the transition to the Mott insulator. The DC transport properties and magnetic relaxation phenomena are contaminated by such noncritical relaxation times $`\tau `$ and are influenced by the carriers in the other portion than the flat part because the flat part has stronger damping and less contributes to the DC properties. Large anisotropy of $`\tau `$ masks the real criticality and makes it difficult to see the real critical exponents in the $`\tau `$-dependent properties. Relevant quantities to easily estimate the criticality is the $`\tau `$ independent quantities such as the Drude weight and the compressibility.

Near the metal-insulator transition, the critical electron differentiation and selective renormalization may lead to experimental observations as if internal degrees of freedom of the carriers such as spin and charge were separated because each degrees of freedom can predominantly be conveyed by carriers in different part of the Fermi surface. Another possible effect of the electron differentiation is the appearance of several different relaxation times which are all originally given by a single quasiparticle relaxation in the isotropic Fermi liquids, but now depend on momenta of the quasiparticles.

Another aspect one might ask in connection to the relevance of the flat part to the metal-insulator transition is the observed level difference between ($`\pi ,0`$) and ($`\pm \pi /2,\pm \pi /2`$) in the undoped and underdoped cuprates . The level at ($`\pi ,0`$) is substantially far from the Fermi level, and at a first glance, it does not have a chance to contribute to the metal-insulator transition. However, it has been clarified  that this level difference is likely to be absent before the d-wave interaction channel starts growing where the flat part indeed governs the criticality, while it is developed by renormalization of the single-particle level accompanied by the d-wave interaction in the lower energy scale. It is remarkable that this renormalization exists even in the insulator.

If the mass renormalization would happen in an isotrpic way as in the picture of Brinkman and Rice, the renormalization can become stronger without disturbance when the insulator is approached. However, if the singularly renormalized flat dispersion emerges critically only in a part near the Fermi surface but the whole band width is ratained, that flattened part has stronger instability due to the coupling to larger energy scale retained in other part of the momentum space. The instability can be mediated by local and incoherent carrier motion generated from two-particle processes derived in the strong coupling expansion . The local two-particle motion is given in the order of $`t^2/U`$ with the bare $`t`$ while the single-particle term is renormalized to $`t^{}`$ which can be smaller than $`t^2/U`$ at the flattened portion. The instability of the flat dispersion was studied by taking account such local and incoherent terms in the Hubbard and $`t`$-$`J`$ models . The inclusion of the two-particle terms drives the instability of the flat part to the superconducting pairing and the formation of the d-wave gap structure. In fact, even at half filling, the two-particle process stabilizes the d-wave superconducting state and reproduces the basic feature of the pseudogap formation  observed in the BEDT-TTF compounds .

The paired bound particles formed from two quasiparticles at the flat spots have different dynamics from the original quasiparticle. In fact, when the paired singlet becomes the dominant carrier, the criticality changes from $`z=4`$ to $`z=2`$, resulting in the recovery of coherence and kinetic energy gain . It generates a strong pairing interaction from the kinetic origin. This pairing mechanism is a consequence of suppressed single-particle coherence and electron differentiation due to strong correlations.

The instability of the flat dispersion coexisting with relatively large incoherent process was further studied . It has turned out that promotion of the above scaling behavior and the flat dispersion offers a way to control potential instabilities. Even when a flat band dispersion is designed near the Fermi level by controlling lattice geometry and parameters, it enlarges the critical region under the suppression of single-particle coherence in the proximity of the Mott insulator mentioned above. In designed lattices and lattices with tuned lattice parameters, it was reported that the superconducting instability and the formation of the spin gap have been dramatically enhanced .

## 3 Pseudogap Phenomena in Cuprates as Superconducting Fluctuations

As is mentioned in §1, the pseudogap in the high-Tc cuprates starts growing from the region of the flat dispersion. When the single-particle coherence is suppressed, the system is subject to two particle instabilities. As clarified in §2, the superconducting instability in fact grows. However, the AFM and charge order correlations are in principle also expected to grow from other two-particle (particle-hole) processes and may compete each other. In particular, the AFM long-range order is realized in the Mott insulator and its short-range correlation is well retained in the underdoped region. Therefore, to understand how the superconducting phase appears in the underdoped region, at least competition of dSC and AFM correlations has to be treated with underlying suppressed coherence in the region of ($`\pi ,0`$) and ($`0,\pi `$). The authors have developed a framework to treat the competition by employing the mode-mode coupling theory of dSC and AFM fluctuations where these two fluctuations are treated on an equal footing .

It should be noted that the strong dSC pairing interaction is resulted from a highly correlated effect with electron differentiation while the critical differentiation has not been successfully reproduced from the diagrammatic approach so far. Then, within the framework of the mode-mode coupling theory, at the starting point, we have assumed the existence of correlation effects leading to the flattened dispersion and the $`d`$-wave pair hopping process. The AFM and dSC fluctuations are predominantly generated by the contributions from the quasiparticle excitations in the flattened regions ($`\pi ,0`$) and ($`0,\pi `$). These fluctuations are treated in a set of selfconsistent equations with mode couplings of dSC and AFM. From the selfconsistent solution, the pseudogap formation is well reproduced in a region of the parameter space. The pseudogap emerges when the mode coupling between dSC and AFM is repulsive with a severe competition and dSC eventually dominates at low temperatures. Such competition suppresses $`T_c`$, while above $`T_c`$ it produces a region where pairing fluctuations are large. This region at $`T_{PG}>T>T_c`$ shows suppression of $`1/T_1T`$ and the pseudogap formation around ($`\pi ,0`$) and ($`0,\pi `$) in $`A(k,\omega )`$. These reproduce the basic feature of the pseudogap phenomena experimentally observed in the underdoped cuprates. The pseudogap formation is identified as coming from the superconducting fluctuations. The momentum dependence shows that the pseudogap formation starts around ($`\pi ,0`$) from higher temperatures and the formation temperature becomes lower with increasing distance from ($`\pi ,0`$) . All of the above reproduce the experimental observations.

We, however, note a richer structure of the gap formation observed in the transversal NMR relaxation time $`T_{2G}`$ and the neutron resonance peak. One puzzling experimental observation is that the pseudogap structure appears in $`1/T_1T`$ , while in many cases $`1/T_{2\mathrm{G}}`$, which measures $`\mathrm{Re}\chi (Q,\omega =0)`$ at $`Q=(\pi ,\pi )`$, continuously increases with the decrease in temperature with no indication of the pseudogap. In addition, the so called resonance peak appears in the neutron scattering experiments . A resonance peak sharply grows at a finite frequency below $`T_\mathrm{c}`$ with some indications even at $`T_c<T<T_{\mathrm{PG}}`$. This peak frequency $`\omega ^{}`$ decreases with lowering doping concentration implying a direct and continuous evolution into the AFM Bragg peak in the undoped compounds. The neutron and $`T_{2\mathrm{G}}`$ data support the idea that the AFM fluctuations are suppressed around $`\omega =0`$ but transferred to a nonzero frequency below $`T_{\mathrm{PG}}`$.

To understand these features, a detailed consideration on damping of the magnetic excitations is required. With the increase in the pairing correlation length $`\xi _d`$, the pseudogap in $`A(k,\omega )`$ is developed. Since the damping is mainly from the overdamped Stoner excitations, the gap formation in $`A(k,\omega )`$ contributes not only to suppress growth of AFM correlation length $`\xi _\sigma `$ but also to reduce the magnetic damping because, inside the domain of the $`d`$-wave order, the AFM excitations are less scattered due to the absence of low-energy quasiparticle around ($`\pi ,0`$). If the quasiparticle damping is originally large around ($`\pi ,0`$), the damping $`\gamma `$ can be reduced dramatically upon the pseudogap formation. Under this circumstance, our calculated result reproduces the resonance peak and the increase in $`1/T_{2G}`$ with lowering temperature at $`T>T_c`$ in agreement with the experimental observations in YBa2Cu3O6.63, YBa2Cu4O8 and some other underdoped compounds .

A subtlety arises when the damping around ($`\pi /2,\pi /2`$) starts contributing. This is particularly true under the pseudogap formation. If contributions from the $`(\pi /2,\pi /2)`$ region would be absent, the damping of the magnetic excitation would be strongly reduced when the pseudogap is formed around $`(\pi ,0)`$ as we mentioned above. However, under the pseudogap formation, the damping can be determined by the Stoner continuum generated from the $`(\pi /2,\pi /2)`$ region and can remain overdamped. This process is in fact important if the quasiparticle damping around the $`(\pi /2,\pi /2)`$ region is large as in the case of La 214 compounds . The formation of the pseudogap itself is a rather universal consequence of the strong coupling superconductors. However, the actual behavior may depend on this damping. If the damping generated by the $`(\pi /2,\pi /2)`$ region is large, it sensitively destroy the resonance peak structure observed in the neutron scattering experimental results.

## 4 Conclusion and Discussion

Electron critical differentiation is a typical property of the proximity of the Mott insulator. The flattening of the quasiparticle dispersion appears around momenta $`(\pi ,0)`$ and $`(0,\pi )`$ on square lattices and determines the criticality of the metal-insulator transition with the suppressed coherence in that momentum region of quasiparticles. Such coherence suppression subsequently causes an instability to the superconducting state when a proper incoherent kinetic process is retained. The d-wave superconducting state is stabilized from such retained microscopic process derived from the strong correlation expansion. The origin of the superconductivity is ascribed to the kinetic energy gain.

By assuming the d-wave channel and the presence of strongly renormalized flat quasi-particle dispersion around the $`(\pi ,0)`$ region, we have constructed the mode-mode coupling theory for the AFM and $`d`$SC fluctuations. The pseudogap in the high-$`T_\mathrm{c}`$ cuprates is reproduced as the region with enhanced $`d`$SC correlations and is consistently explained from precursor effects for the superconductivity. The existence of the flat region plays a role to suppress the effective Fermi temperature $`E_F`$. This suppressed $`E_F`$ and relatively large local pair hopping process both drive the system to the strong coupling superconductor thereby leading to the pseudogap formation. The pseudogap formation is also enhanced by the AFM fluctuations repulsively coupled with dSC fluctuations.

Several similar attempts have also been made to reproduce the pseudogap phenomena observed in the cuprates. A common conclusion inferred from these calculations including those by the present authors is that the pseudogap is reproduced when the d-wave channel is explicitly assumed  while crucial features such as the pseudogap formation around ($`\pi ,0`$) cannot be well reproduced if one tries to derive the superconducting channel itself from the AFM spin fluctuations . This difficulty is summarized: if one desires to stabilize a strongly fluctuating pseudogap region well above $`T_c`$ generated originally by the magnetic interaction, one has to treat the strong-coupling superconductivity and hence even more strong-coupling magnetic interaction. In such a circumstance, it is dificult to escape from the magnetic instability before the pseudogap formation. The difficulty is naturally interpreted from our analysis: The kinetic d-wave channel we derived in the strong coupling expansion is not ascribed to the magnetic origin and is not contained neither in the spin-fluctuation mechanism nor in the $`t`$-$`J`$ model. This superconducting channel has a comparable amplitude to the magnetic one, namely $`t^2/U`$. Furthermore the channel is enhanced by the flattened dispersion and the electron critical differentiation. Such effects are far beyond the one-loop level of the weak coupling approach. The change in the quasiparticle dispersion in the portion of the Fermi surface is not well reproduced for the moment and the momentum dependent selfenergy appears to be significantly underestimated in the existing diagrammatic evaluations.

The strong coupling Hamiltonian was considered before  and also in terms of the pairing arising from the attractive superexchange interaction . The superconductivity from direct kinetic origin was also discussed before in a specific model  as well as in the interlayer tunneling mechanism , where interlayer kinetic energy gain is required to stabilize the superconductivity. If we identify the intralayer charge incoherence observed in the normal state as the crossover phenomena above the unusually suppressed Fermi temperature in the proximity of the Mott insulator, and notice that no quantum coherence is reached even within the layer, the kinetic origin of the superconductivity is found solely in a two-dimensional plane.

Under the suppression of coherence in the proximity of the Mott insulator, the second-order process in the strong-coupling expansion may also drive other type of ordered state in addition to the superconductivity. In particular, because the charge compressibility is enhanced due to the flattened dispersion, the charge fluctuation may be strongly enhancced and the charge ordering and stripes can be triggered both by the kinetic energy gain and the intersite Coulomb interaction. When the competition between the superconductivity and the charge ordering becomes serious, presumable repulsive mode coupling of them may become another origin of the pseudogap formation, although such evidence is not clear in the present experimental results of the cuprates.

The presend theory of the kinetic superconductor predicts a specific form of the kinetic energy gain to be seen in the optical conductivity  and the single particle spectra . Qualitatively, it is expected that the in-plane kinetic energy starts gained even in the pseudogap region. The energy gain dominantly coming from the $`(\pi ,0)`$ and $`(0,\pi )`$ region of the single particles must have a significant doping dependence associated to the change in the dynamical exponent $`z`$ from 4 to 2.


format: text

text:
# Evolution of Buoyant Bubbles in M87

## 1 Introduction

The giant elliptical galaxy M87 continues to be the subject of numerous experimental and theoretical studies. It is located at, or near, the center of the very nearby ($``$ 18 Mpc, 1 arcminute corresponds to $``$5 kpc) irregular Virgo cluster and is at the center of a weak cooling flow with an estimated mass deposition rate of some 10-40 $`M_{}`$ per year (e.g. Peres et al. 1998). The radio source Virgo A (3C274), associated with the galaxy M87, is well known because of the spectacular jet which is observed both in radio and optical bands (see e.g. Biretta 1999 for a review). In addition to the jet, which is contained within the central 2 kpc or $`30^{\prime \prime }`$ region, there is a lower surface brightness radio halo extending up to a distance of $``$40–50 kpc from the nucleus (e.g. Kassim et al. 1993, Böhringer et al. 1995, Owen, Eilek and Kassim 1999, 2000).

The X–ray emission is strongly peaked at the position of the nucleus of M87 and is largely symmetric around it. However, it was found in the Einstein HRI observations (Feigelson et al. 1987) and then in the ROSAT data (Böhringer at al. 1995) that there are departures from a symmetric model and these departures resemble the morphology of some prominent features in the outer radio lobes. These findings initiated numerous discussions on the interaction of the radio and X–ray emitting plasmas (e.g. Feigelson et al. 1987, Böhringer et al. 1995, Böhringer et al. 1999, Owen, Eilek and Kassim 1999,2000, Harris et al. 1999). In a number of other cooling flow clusters with strong radio halos surrounding the dominant galaxy (e.g. Perseus – Böhringer at al. 1993, McNamara et al. 1996, Churazov et al. 2000, Fabian et al. 2000a, A4059 – Huang and Sarazin 1998, Hydra A – McNamara et al. 2000) there is a clear anticorrelation of the radio and X–ray emitting plasmas. This makes the case of M87 especially interesting.

We discuss below one particular model which is able to reproduce qualitatively both the radio and X–ray morphologies observed in the central region surrounding the M87 galaxy.

### 1.1 Radio and X–ray morphology

There is a long history of M87 radio halo observations, starting with Mills (1952) shortly after the initial discovery of the radio source Virgo A (Bolton, Stanley and Slee 1949). Because of the very bright compact source, the dynamic range of the images was always a problem. Recently, Owen et al. (2000) presented a new high resolution, high dynamic range map of the halo of M87 at 327 MHz (see Figure 1). For the sake of further discussion this map has been rotated by 90 degrees to have East up and North to the right. In this map several distinct structures are discernible. The high brightness centre represents the well-studied inner lobe structure (oriented approximately east-west) with the famous jet pointing west north-west (approximately bottom–right for the orientation of images adopted in Fig.1). The highly structured outer halo is much fainter. It consists of the ear-shaped (torus-like) eastern bubble, the much less well-defined western bubble, both of which are connected to the central emission by a ‘trunk’, and the two very faint, almost circular, emission regions northeast and southwest of the centre.

Owen et al. (2000) find a peak surface brightness in the region of the eastern bubble of about 53 mJy/beam while the circular regions range from 14 mJy/beam to about 25 mJy/beam. Rottmann et al. (1996) show that all of the described structures are also present at 10.6 GHz. They also show that the outer halos have a radio spectral index in the range from $`2.8`$ to $`2.5`$ between 4.8 GHz and 10.6 GHz. The spectral index is flatter by a factor of about 1.5 between 330 MHz and 1.5 GHz. This suggests that radiative energy losses of the synchrotron emitting particles lead to a significant steepening of the radio spectrum at about 3 GHz. Unfortunately, the maps used by Rottmann et al. (1996) have low resolution since they are obtained with single dish observations and so the spectral index distribution of the halo and its connection to the detailed structures, e.g., bubbles and attached “trunks”, is not known.

In Fig.1 (lower left) we show the X–ray surface brightness distribution of the same region (and similar orientation) as the radio map. The data are the combined 200 ksec ROSAT/HRI image, adaptively smoothed to suppress high frequency noise. As pointed out by a number of authors (Feigelson et al. 1987, Böhringer at al. 1995, Owen, Eilek and Kassim 1999, Harris et al. 1999) there is evidence for a correlation between X–ray and radio emitting features. The simplest explanation of this correlation is that the excess X–ray emission is due to inverse Compton scattering of the cosmic microwave background photons by the same relativistic electrons which produce the synchrotron radio emission (Feigelson et al. 1987). However, ROSAT/PSPC observations have shown that the excess emission has a thermal spectrum (Böhringer at al. 1995) and the X–ray emitting gas in these region has a lower temperature than that in the ambient regions.

Thus, we summarize three observational facts which we are trying to explain below:

* There are prominent “torus–like” features in the radio image.
* There is a correlation (but not one to one) of the X–ray and radio bright regions.
* Excess X–ray emission, associated with the radio bubbles and attached “trunks”, arises from thermal gas whose temperature is less than that of the ambient X-ray emitting gas.

### 1.2 Analogy with atmospheric explosions

Concentrating on the “torus–like” radio features one can note its striking similarity with some evolutionary stages of hot buoyant bubbles formed by a powerful (e.g. nuclear) explosion in the Earth’s atmosphere. The transformation of the initially spherical bubble into a torus is a common and well known property of buoyant bubbles lacking strong surface tension (e.g. Walters and Davison 1963, Onufriev 1967, Zhidov et al. 1977). Morphologically similar structures resembling ”mushrooms” appear in Rayleigh-Taylor unstable configurations: as the fluid rises through the ambient medium, Kelvin-Helmholtz instabilities create the torus-like head of the ”mushroom”. Images of various laboratory experiments where such structures were observed and discussion of the relevant physical mechanisms are widespread (see e.g. Batchelor 1967, Turner 1973, Inogamov 1999). The idea that buoyancy plays an important role in the evolution of the radio lobes in galaxies was first proposed by Gull and Northover (1973) and has been used to estimate the life time of the radio lobes in M87 by Böhringer et al. (1995). The role of buoyancy was also discussed by Baum and O’Dea (1991) in application to PKS 0745-191. The similarity of the M87 “torus–like” features to the Rayleigh–Taylor mushrooms was pointed out by Churazov et al. (2000) and to a subsonic vortex ring by Owen, Eilek and Kassim (2000).

Further pursuing the analogy with powerful explosions, we note that during the transformation of a bubble into a torus some ambient gas is captured and uplifted by the rising bubble/torus. In the cooling flow surrounding M87, the entropy of the X–ray emitting gas rises quickly with distance from the centre, and therefore the emissivity of the gas increases very steeply if it remains in pressure balance with the ambient gas (Böhringer et al. 1995, Nulsen 1997). This may qualitatively explain the correlation of the radio and X–ray emitting plasmas and naturally accounts for the thermal nature of the excess thermal emission.

In the case of an atmospheric explosion, the final stage of the hot bubble occurs when the bubble reaches a height at which the density of the ambient gas is equal to the density of the bubble. Then (neglecting oscillations near the equilibrium position) the bubble expands laterally (forming a “pancake”) to occupy a thin layer in the atmosphere having the same density. The bubbles/tori in M87 may share the same fate. In a spherically symmetric potential, the bubble will try to fill a segment (equipotential surface) of a sphere. The largest distinct features in the radio map (see Fig.1) could be just those late-stage bubbles.

The sketch of the possible overall source structure based on the analogy with a powerful atmospheric explosion is shown in Fig.1. Here the black region in the center denotes the inner radio lobes, the gray mushrooms correspond to the buoyant bubbles already transformed into tori and the gray lens-shaped structures are the “pancakes” formed by the older bubbles.

Buoyant bubbles and their role in the complex X–ray and radio morphology of M87 have been discussed previously e.g. in Böhringer et al. (1995), Nulsen (1997), Böhringer (1999), Owen, Eilek and Kassim (2000). Below we follow the line of arguments suggested in Churazov et al. (2000) and present hydrodynamic calculations of the radio lobe evolution to illustrate the qualitative picture described above.

## 2 Method

In the next sections we describe the method, initial conditions, and assumptions used in the numerical simulations of bubbles in a hot gaseous atmosphere. We first detail the method and the use of “tracer” particles to track the motion of hot gas which is entrained in the rising bubble. We next describe our approach to calculation of the synchrotron radio emission from the evolving bubble.

### 2.1 Hydrodynamic simulations

The simulations were obtained using the ZEUS-3D code which was developed especially for problems in astrophysical hydrodynamics (Stone & Norman 1992a). The code uses finite differencing on an Eulerian grid and is fully explicit in time. It is based on an operator-split scheme with piecewise linear functions for the fundamental variables. The fluid is advected through a mesh using the upwind, monotonic interpolation scheme of van Leer (1977). For a detailed description of the algorithms and their numerical implementation see Stone & Norman (1992a, b). To study the spectral aging of the relativistic gas, the ZEUS code was modified to follow the motion of ‘tracer’ particles which are advected with the fluid.

In our simulations we employed an ideal gas equation of state and we ignored the effects of magnetic fields and rotation. The cooling time is sufficiently long compared to the time scales considered here since, for an electron density of 0.02 cm<sup>-3</sup> and temperature of the order of 2 keV (see Fig.2) the cooling time is $`5\times 10^8`$ years. This is an order of magnitude longer than the typical duration of the run of $`5\times 10^7`$ years. We therefore completely neglected cooling. The simulations were computed on a spherical grid in two dimensions and were performed on an IBM RS/6000 cluster. The computational domain spans 51 kpc in radius and $`\pi /2`$ rad in angle and was covered by 200 $`\times `$ 200 grid points (in the $`r`$\- and $`\theta `$-direction).

A model for the mass and initial temperature was adopted from Nulsen & Böhringer (1995) which is shown in Fig.2. The gas density distribution was then found by assuming hydrostatic equilibrium to maintain an initially static model. A spherical bubble was created with a radius of $`r_\mathrm{b}=5`$ kpc at a distance of $`d=9`$ kpc from the gravitational center. It was made buoyant by reducing its density with respect to the background by a factor of 100 and simultaneously raising the temperature by the same factor. During its subsequent evolution, the gas is treated as a single fluid which is assumed to obey a polytropic equation of state with $`\mathrm{\Gamma }=5/3`$. The bubble was filled with 2000 tracer particles and was placed on the $`\theta =0`$ axis along which it rose in the gravitational potential. Subsequently, we rotated the two-dimensional output about this axis (azimuthally). To test the assumption of rotational symmetry, we performed some 3D simulations with reduced resolution and found that this assumption is fairly robust as far as the global motions and morphologies are concerned. As we discuss in detail below, buoyancy deforms the bubble and drives it through the ambient medium as shown in Fig.3. Shear instabilities cause the formation of tori which separate from the main bubble. This process, also called “vortex shedding”, has been observed and studied extensively (see Norman et al. 1982).

Finally, we should address some issues related to the accuracy of these kinds of finite-difference hydrodynamic simulations. While the code can simulate large-scale mixing due to Rayleigh-Taylor and Kelvin-Helmholtz instabilities, it does not include real particle diffusion. Any observed diffusion is therefore entirely numerical. The boundary between the bubble and the ambient medium becomes less sharp as the simulation proceeds due to discretization errors in the advection scheme. For a test of the advection algorithm in the ZEUS code see Stone & Norman (1992a,b). In simple advection tests, it was found that, during the advection of a sharp discontinuity over a grid of 200 zones, the discontinuity is spread over 3-4 grid cells. Therefore, the small features in our simulation are likely to be affected by these advection errors whereas the larger features are not.

Second, numerical viscosity is also responsible for suppressing small-scale instabilities at the interface between the bubble and the cooler, surrounding, X-ray emitting gas. To assess the effects of numerical viscosity, we have repeated our simulations on grids with 100 $`\times `$ 100 grid points. From our experiments, we can conclude that “global parameters” such as the position and size of the bubble as well as the presence of “toroidal” structure are relatively insensitive to the resolution. The detailed small scale morphology does depend on the resolution and the initial conditions.

Finally, we stress that the simulations presented here are not intended to provide a detailed simulation of the radio/X-ray structure surrounding M87. Rather they are intended as a guide, a toy model, to assist our interpretation of the complex emission (see discussion in Section 3.1).

### 2.2 Simulation of the radio emission

The observations of radio synchrotron emission from the buoyant bubbles in M87 imply the presence of a relativistic, magnetized plasma. However, the simulations presented in previous sections are based on a purely hydrodynamic scheme with a single, non-relativistic fluid. To estimate the radio flux expected from the simulation, we make the following assumptions.

The relativistic plasma and the associated magnetic field, the ‘relativistic fluid’, are confined within small volumina intermixed with the non-relativistic, thermal plasma governing the dynamics. This implies that the relativistic fluid may mix with the thermal fluid on macroscopic scales but not on microscopic ones. This is very similar to the behavior of air bubbles in water which are confined by surface tension (e.g. Zhidov et al. 1977).

We assume the magnetic field to be tangled on scales small compared to the size of the volumina of relativistic fluid. This allows us to treat the magnetic field as part of the relativistic fluid. The total energy density inside the small volumina is therefore the sum of the energy density of the magnetic field, $`u_\mathrm{B}`$, and that of the relativistic particles, $`u_\mathrm{e}`$. We assume that the small volumina of relativistic fluid are in pressure equilibrium with the thermal fluid throughout the simulation, i.e. $`3p_{\mathrm{th}}=u_\mathrm{B}+u_\mathrm{e}`$. We also assume that there is no significant re-acceleration of relativistic particles in the bubble during its evolution.

At the start of the simulation, we assume some initial value of the magnetic field (the same value for every small volumina of relativistic particles within the bubble) and a power law distribution for the relativistic electrons (in terms of the particle Lorentz factor), $`\gamma `$, $`n_\mathrm{e}(t_\mathrm{o})d\gamma =n_\mathrm{o}\gamma ^pd\gamma `$. At the start of the simulation $`u_\mathrm{e}`$ is therefore given by

$$u_\mathrm{e}=m_\mathrm{e}c^2n_\mathrm{o}_{\gamma _{\mathrm{min}}}^{\gamma _{\mathrm{max}}}\gamma ^p\left(\gamma 1\right)𝑑\gamma ,$$
(1)
where we assume $`\gamma _{\mathrm{min}}=1`$, $`\gamma _{\mathrm{max}}=10^6`$. For a given initial magnetic field the normalization $`n_\mathrm{o}`$ of the particle distribution is then adjusted such that the entire relativistic fluid is in pressure equilibrium with the thermal fluid.

For values of $`p2`$, the energy density of relativistic particles is dominated by low energy particles. The life time of such particles, due to radiative losses, is much longer than the duration of the simulation (and no Coulomb losses are present since we assume that thermal gas is not mixed with relativistic particles on microscopic scales). The motion of the bubble (see Section 3) is found to be subsonic even with respect to the cluster gas. We therefore assume that the energy density of the magnetic filed and relativistic particles changes adiabatically according to the changes of the pressure of the surrounding thermal gas: $`u_\mathrm{B}u_\mathrm{e}p_{\mathrm{th}}`$. Note, however, that our hydrodynamic simulations assume an adiabatic index of 5/3 for the nonrelativistic gas, while the relativistic fluid expands and contracts adiabatically with an adiabatic index $`\mathrm{\Gamma }_\mathrm{r}=4/3`$. This implies that, during the expansion of a buoyant bubble, the fraction of the volume occupied by the relativistic fluid increases relative to the volume of the thermal fluid as $`V_\mathrm{r}V_{\mathrm{th}}^{5/4}`$. Strictly speaking this means that accounting for relativistic fluid with a different adiabatic index will not affect the dynamics of the bubble evolution as long as the fraction of the volume occupied by the relativistic fluid is small. We will see in the following that this assumption is probably invalid, but we believe that it does not very severely affect the results.

The energy losses of the relativistic particles are due to adiabatic expansion, synchrotron radiation and inverse Compton scattering of the CMB. These losses depend on the pressure history of each volumina of relativistic particles. The pressure history can be recovered using the tracer particles. For a given region within the bubble the change of the Lorentz factor $`\gamma `$ of a given electron due to adiabatic expansion or contraction of the relativistic fluid and radiative losses can be written as

$`{\displaystyle \frac{d\gamma }{dt}}`$ $`=`$ $`a{\displaystyle \frac{\gamma }{p_{\mathrm{th}}(t)}}{\displaystyle \frac{dp_{\mathrm{th}}(t)}{dt}}`$ (2)
$``$ $`{\displaystyle \frac{4\sigma _\mathrm{T}}{3m_\mathrm{e}c}}\gamma ^2\left[u_\mathrm{B}(t_\mathrm{o}){\displaystyle \frac{p_{\mathrm{th}}(t)}{p_{\mathrm{th}}(t_\mathrm{o})}}+u_{\mathrm{IC}}\right],`$
where $`a=(\mathrm{\Gamma }1)/\mathrm{\Gamma }=1/4`$.

This can be integrated (e.g. Kaiser, Dennett-Thorpe & Alexander 1997) to give the relation between the Lorentz factor of the electron at time $`t`$ and the initial Lorentz factor of the same electron

$$\gamma _0=\frac{\gamma p_{\mathrm{th}}(t)^a}{p_{\mathrm{th}}(t_0)^a+\gamma b},$$
(3)
with

$`b={\displaystyle \frac{4\sigma _\mathrm{T}}{3m_\mathrm{e}c}}{\displaystyle _{t_\mathrm{o}}^t}\left[u_\mathrm{B}(t_\mathrm{o}){\displaystyle \frac{p_{\mathrm{th}}(t)}{p_{\mathrm{th}}(t_\mathrm{o})}}+u_{\mathrm{IC}}\right]p_{\mathrm{th}}(t)^a𝑑t.`$ (4)
The value of $`b`$ is calculated from the pressure history of the tracer particles. The integral in this expression must be evaluated numerically using the discrete time steps of the simulation.

Finally, the adiabatic change of the volume in which the relativistic particles are contained implies $`n_\mathrm{o}(t)=n_\mathrm{o}(t_\mathrm{o})[p_{\mathrm{th}}(t)/p_{\mathrm{th}}(t_\mathrm{o})]^{1/\mathrm{\Gamma }}`$. The energy distribution of relativistic electrons at time $`t`$ is then

$$n_\mathrm{e}d\gamma =n_\mathrm{o}(t)\gamma _0^p\frac{d\gamma _0}{d\gamma }d\gamma =n_\mathrm{o}(t_\mathrm{o})\frac{p_{\mathrm{th}}(t)\gamma _\mathrm{o}^{2p}}{p_{\mathrm{th}}(t_\mathrm{o})\gamma ^2}d\gamma .$$
(5)
The energy distribution of the relativistic particles combined with the strength of the magnetic field yields the monochromatic synchrotron emissivity (e.g. Longair 1994)

$$j_\nu =\frac{4}{3}\sigma _\mathrm{T}cu_\mathrm{B}_{\gamma _{\mathrm{min}}}^{\gamma _{\mathrm{max}}}\gamma ^2\varphi (\gamma ,\nu )n_\mathrm{e}𝑑\gamma ,$$
(6)
where $`\varphi (\gamma ,\nu )`$ is the dimensionless spectrum emitted by a single electron.

To obtain the synchrotron surface brightness from the simulation, we assume rotational symmetry about the polar-axis and integrate the emissivity along lines of sight through the resulting 3-dimensional model.

## 3 Results and discussion

With the foundation and assumptions described in the previous section, we describe the results of two numerical simulations. We first present the results from the two numerical simulations with different initial conditions for the bubbles. We compute the velocity and position versus time of the rising bubble. For the first simulation, we compute the synchrotron radio emission from the bubble for two orientations, one aligned with the symmetry axis in the plane of the sky and the second withthe symmetry axis aligned 45 with respect to the plane of the sky. The two orientations differ in the initial conditions for the magnetic field strength in order to reproduce the observed radio spectrum. We derive the X-ray appearance of the bubble from the first simulation, also for the two different orientations. Finally, we discuss the “pancake” phase and calculate the largest distance to which bubbles can rise as a function of their initial location and volume fraction of trapped, hot gas.

### 3.1 Bubble configuration

Fig.3 shows snapshots of the temperature distribution during the evolution of the bubble. One can identify two stages of the bubble evolution. The initially spherical bubble flattens and develops a “cone” at the rear, which is filled with entrained gas. The bubble then transforms into a torus which at later stages may fragment into smaller structures, but which continue to show a global toroidal structure. Note that ambient gas (in particular the gas captured during the bubble $``$ torus transformation) occupies the central part of the rising structure. Note also that in the last stages (shown in Fig.3) the coldest temperatures are found not at the center of the cluster but within the rising structure. We stress again that no radiative cooling was included in the simulations and these low temperatures regions are associated with gas which has been uplifted from the central regions and has expanded adiabatically to match the ambient gas pressure at the current location of the bubble.

Of course the particular shape of the torus/mushroom structure (e.g. two tori at different distances from the center in Fig.3) is the result of our choice of initial and boundary conditions and assumed axial symmetry. For example, Fig.4 shows snapshots of the temperature distribution for slightly different initial conditions (with an initial temperature contrast, between the ambient gas and the bubble, which is a factor of two larger than in the previous run). Note the changes in the detailed morphology of the rising bubble compared to Fig.3. In general one would expect a large variety of structures to be formed, especially in 3D. But the torus-like geometry is a generic feature of rising bubbles in the absence of strong surface tension.

### 3.2 The velocity of the rising bubble

The simplest velocity estimate of the rising bubble can be obtained by equating the ram pressure and buoyancy forces acting upon a bubble (e.g. Gull and Northover 1973). The buoyancy force is obviously:

$`F_\mathrm{b}=Vg(\rho _\mathrm{a}\rho _\mathrm{b}),`$ (7)
where $`V`$ is the bubble volume, $`g`$ is the gravitational acceleration (we assume the ambient gas is in hydrostatic equilibrium), $`\rho _\mathrm{a}`$ and $`\rho _\mathrm{b}`$ are the mass densities of the ambient and the bubble gas respectively. The ram pressure (drag) is:

$`F_dC{\displaystyle \frac{1}{2}}Sv^2\rho _\mathrm{a},`$ (8)
where $`S`$ is the cross section of the bubble. The numerical coefficient $`C`$ (drag coefficient) depends on the geometry of the bubble and the Reynolds number. Thus the terminal velocity of the bubble is

$`v\sqrt{g{\displaystyle \frac{V}{S}}{\displaystyle \frac{2}{C}}{\displaystyle \frac{\rho _\mathrm{a}\rho _\mathrm{b}}{\rho _\mathrm{a}}}}\sqrt{g{\displaystyle \frac{V}{S}}{\displaystyle \frac{2}{C}}}.`$ (9)
Here the factor of $`(\rho _\mathrm{a}\rho _\mathrm{b})/\rho _\mathrm{a}`$ can be dropped if the bubble density is low compared to the ambient gas density. The expression for the terminal velocity can be further rewritten using the Keplerian velocity at a given distance from the cluster center: $`v\sqrt{(r/R)(8/3C)}v_\mathrm{K}`$, where $`r`$ is the bubble radius, $`R`$ is the distance from the center and $`v_\mathrm{K}=\sqrt{gR}`$ is the Keplerian velocity. For a solid sphere moving through an incompressible fluid the drag coefficient $`C`$ is of the order 0.4–0.5, for Reynolds numbers in the range $`10^3`$$`10^5`$ (e.g. Landau and Lifshitz 1963). The drag coefficient for a rising bubble should of course be different. First of all, it is not a solid sphere, but a rather complicated structure, resembling a smoke ring in air. Such rings may travel large distances through the air with low drag. Secondly the Mach number of the rising bubble in our simulations is $``$0.6–0.7 and the compressibility of the gas is important. This would tend to increase the drag coefficient, e.g., for a sphere moving at $`M0.7`$ the drag coefficient is $`C0.6`$. Thirdly, in a stratified medium an additional contribution to the drag may come from the excitation of internal gravity waves. In our simulations, the typical Keplerian velocity was $``$ 400 km s<sup>-1</sup>. Fig.5 shows the position of the bubble front as a function of time. The solid line shows the expected position of the feature moving with the constant velocity of $``$ 390 km s<sup>-1</sup>. The effective drag coefficient which can be estimated from the simulations, assuming the above formula for the rise velocity, is $`C0.75`$. We note here that this value may in turn be affected by the numerical resolution adopted here. However, since the velocity of the bubble depends on the square root of the drag coefficient, formula (9) still yields crude order of magnitude estimates of the bubble velocity. Thus, a large and strongly underdense bubble will rise with a velocity comparable to the Keplerian velocity.

In the case of M87, the bubble’s environment may not be very underdense compared to the ambient gas (since no clear X–ray holes have been seen) and the velocity of the bubble may be somewhat smaller. However, we believe that the above estimate of the velocity is approximately correct since the velocity depends only weakly (as the 0.5 power) on the bubble size and gas density gradients. This velocity defines a time scale of $``$ few $`10^7`$ years for the evolution of the bubble to the stage seen as a bright torus in Fig.1.

### 3.3 The rising bubble in synchrotron emission

In Figure 6 we show the simulation results at 327 MHz. The radio emission was calculated for the hydrodynamic simulations shown in Fig.3. Upper and lower rows correspond to different viewing angles. The symmetry axis is perpendicular to the line of sight in the upper row (case 1) in this figure and at an angle of $`45^{}`$ to the line of sight (case 2) in the lower row. From the radio observations of M87 we note that the well-defined eastern bubble is observed at a distance of roughly 25 kpc from the source centre. This corresponds to the third panel in the upper row of Figure 6, i.e. bubble age of 42 Myrs, while for the simulation projected at an angle of 45 (lower row), the fourth panel, at an age of 67 Myrs, gives the correct projected distance. Note here that these distances are for the primary bubble/vortex ring. Figure 6 clearly shows the presence of a secondary vortex ring which is closer to the centre. The appearance of this secondary ring is fairly dependent on the particular choice of initial conditions (see Section 2 and Figures 3,4). In the case studied here, the synchrotron surface brightness of the secondary vortex ring exceeds that of the primary bubble. No such structure is observed in M87 and so we concentrate in the following discussion on the primary bubble.

For an initial power law distribution of relativistic electrons in a constant magnetic field, a significant steepening of the synchrotron spectrum will be observed for frequencies higher than the break frequency, $`\nu _{\mathrm{break}}`$, which is given by (e.g. Leahy 1991)

$$\nu _{\mathrm{break}}=2.52\times 10^6\frac{B}{\left(B^2+B_{\mathrm{CMB}}^2\right)^2t^2}\mathrm{GHz},$$
(10)
where $`B`$ is the strength of the magnetic field in $`\mu `$G and $`t`$ is the time (in Myrs) for which the plasma was exposed to this field. $`B_{\mathrm{CMB}}`$ is the strength of the magnetic field in $`\mu `$G equivalent to the energy density of the CMB, $`u_{\mathrm{IC}}`$. For the redshift of M87, $`u_{\mathrm{IC}}4.2\times 10^{13}`$ erg cm<sup>-3</sup> and so $`B_{\mathrm{CMB}}3.2\mu `$G. From equation 10 the time $`t`$ can be expressed as a function of the break frequency and the strength of the magnetic field. For any given break frequency, the time $`t`$ is maximal if $`B(t_{\mathrm{max}})=B_{\mathrm{CMB}}/\sqrt{3}2\mu `$G. We have seen that the observed radio spectral indices (Section 1.1) require $`\nu _{\mathrm{break}}3`$ GHz, for which a maximum age is $`t_{\mathrm{max}}10^8`$ years. This is comparable with the time span of our simulations. For example, if the symmetry axis is at angle of $`45^{}`$ with the line of sight (case 2 above), then the front of the bubble will reach a projected distance of $``$25 kpc after 67 Myrs. This simple estimate shows that we require a rather low value of the initial magnetic field strength (comparable to $`B(t_{\mathrm{max}})2\mu `$G) in order to avoid having the $`\nu _{\mathrm{break}}`$ lie at a frequency lower than $`3`$ GHz.

For case 1, we used $`p=2.3`$ and an initial magnetic field strength of 6.5 $`\mu `$G. We then calculated the synchrotron emission of the bubble as described in Section 2.2. After $``$ 42 Myrs, i.e. when the bubble is $``$25 kpc away from the center, the spectral index of the primary bubble is $`2.7`$ from 4.7 – 10.6 GHz and $`1.0`$ from 330 – 1500 MHz which compares well with the observed values (Rottmann et al. 1996, Section 1.1).

For case 2 the magnetic field inside the bubble must be lower since here the bubble is older when its upper edge reaches the correct projected height from the cluster centre. We set $`B=2.9\mu `$G and $`p=2.3`$. For this case the high frequency spectral index is found to be $`3.0`$ while at lower frequencies it is $`1.0`$.

The assumption of pressure equilibrium determines the initial energy density of the relativistic particles. The energy density of the magnetic field is $`1.7\times 10^{12}`$ ergs cm<sup>-3</sup> for case 1 and $`3.4\times 10^{13}`$ ergs cm<sup>-3</sup> for case 2. This is much lower than the energy density of the thermal cluster gas and the particles energy density is therefore $`u_\mathrm{e}3p_{\mathrm{th}}4.2\times 10^{10}`$ ergs cm<sup>-3</sup> for both cases discussed. This implies a severe departure from the equipartition condition of $`u_\mathrm{B}u_\mathrm{e}`$ often assumed for synchrotron emitting plasmas. The initial volume filling factor of the radio emitting plasma in the initial bubble can then be adjusted to provide the required surface brightness of the final bubble. The peak surface brightness predicted in the region of the primary vortex ring is 55 mJy/pixel in case 1 (see Owen et al. 2000, Section 1.1), however, the required volume filling factor of the relativistic plasma is very close to unity. For case 2, even for a filling factor of unity, the peak surface brightness of the bubble is only 2.4 mJy/pixel. Furthermore, we neglected any possible contribution from the relativistic ions (the so-called ”k” factor) and thermal plasma to the pressure inside the bubble, which would further reduce its synchrotron emissivity. We also note here that volume filling factors $``$1 pose a serious problem for our assumption that the relativistic fluid does not influence the overall dynamics of the buoyant bubble.

The problems discussed above arise from the life time of the synchrotron emitting particles being short compared to the dynamical time scale of the structure. Gull & Northover (1973) first noted this lifetime problem for buoyant models. Below, we discuss various solutions.

From Figure 6, we find that the diameter of our primary bubble is smaller than the observed diameter of the eastern bubble in M87. As we showed in Section 3.2, the rise time of the bubble is correlated with its size. It is therefore quite possible that the eastern bubble in M87 is somewhat younger than our simulations suggest. Such a younger age would allow a stronger magnetic field in the bubble and could explain the observed spectrum, which in turn would allow a smaller volume filling factor for the relativistic fluid.

The buoyant bubble may also be filled mainly by a weak magnetic field interspersed with regions of a much higher field strength but small volume filling factor. The relativistic electrons may then survive for a much longer time in the mainly weak magnetic field and only diffuse slowly into the regions of greater field strengths which are responsible for most of the observed synchrotron emission (Eilek, Melrose & Walker 1997). Subsonic turbulence can also enhance the energy density of the magnetic field in regions with initially weak field (e.g. Eilek, Owen and Zhou 1999).

If the magnetic field in the buoyant bubbles is stronger than that estimated above and rather homogeneous, then in situ acceleration of relativistic particles is required. In this case, the age of the buoyant bubble cannot be constrained from radio observations. The two latter mechanisms were also discussed by Owen et al. (2000) to explain the radio emission of the outer M87 halo.

The fundamental conclusion of these simulations is that one can reproduce the observed radio brightness and the break in the spectrum of the torus–like feature in M87 assuming that (i) the bubble, filled with the relativistic plasma, is in pressure equilibrium with the thermal gas, (ii) the initial magnetic field strength is relatively low so as to prolong the lifetime of the electrons and (iii) no in situ reacceleration of particles or generation of magnetic field occurs during the bubble evolution. However, the equipartition condition is sacrificed and the parameter space which satisfies these assumptions is rather tight. It is likely that the real situation is much more complex than was assumed in the simple model discussed above.

### 3.4 The rising bubble in X–ray emission

As is clear from Fig.3, the bubble captures ambient gas during its transformation to a torus and then uplifts this gas to larger radii. Since all motions in our simulations are subsonic, this uplifted gas remains in approximate pressure equilibrium with the ambient gas at larger radii. The volume emissivity of this gas can be compared with the volume emissivity of the ambient gas:

$`ϵ_u=n_u^2\mathrm{\Lambda }(T_u)`$
$`ϵ_a=n_a^2\mathrm{\Lambda }(T_a),`$ (11)
where $`n`$ and $`T`$ are the density and temperature of the uplifted and ambient gas respectively, $`\mathrm{\Lambda }(T)`$ is the emissivity in a given energy band. If we assume adiabatic evolution of the entrained and uplifted gas then:

$`n_u=n_{a,0}\left({\displaystyle \frac{P_a}{P_{a,0}}}\right)^{1/\mathrm{\Gamma }}=n_a{\displaystyle \frac{T_a}{T_{a,0}}}\left({\displaystyle \frac{P_a}{P_{a,0}}}\right)^{1/\mathrm{\Gamma }1}`$ (12)
$`T_u=T_{a,0}\left({\displaystyle \frac{P_a}{P_{a,0}}}\right)^{11/\mathrm{\Gamma }},`$ (13)
where the subscript $`0`$ stands for the values of the uplifted gas at its initial location where it was captured by the rising bubble. Thus the uplifted gas will have an emissivity which is higher than the ambient gas by a factor:

$`F={\displaystyle \frac{ϵ_u}{ϵ_a}}=\left({\displaystyle \frac{T_a}{T_{a,0}}}\right)^2\left({\displaystyle \frac{P_a}{P_{a,0}}}\right)^{2/\mathrm{\Gamma }2}\times {\displaystyle \frac{\mathrm{\Lambda }(T_u)}{\mathrm{\Lambda }(T_a)}}.`$ (14)
The enhancement factor $`F`$ (excluding the cooling functions) is shown in Fig.7. One can see that in a cooling flow (density decreases with distance while temperature increases), the emissivity of the uplifted gas greatly exceeds the emissivity of the ambient gas if the displacement is large. The ratio $`\mathrm{\Lambda }(T_u)/\mathrm{\Lambda }(T_a)`$ (when the temperatures are near 1 keV and the soft X–ray band is considered) further increases the brightness of the uplifted gas with respect to the ambient gas.

To illustrate how the X–ray surface brightness is affected by the a radio emitting plasma bubble and lumps of uplifted cold gas, we calculated the projected surface brightness distribution (for the ROSAT energy band from 0.5 to 2 keV) for several stages of the bubble evolution. When calculating these maps, we assumed axial symmetry and assumed that the axis of the bubble is either perpendicular to the line of sight (upper row of images in Fig.8) or that the angle between the line of sight and bubble axis is 45 (lower row of images in Fig.8). Fig.8 shows the relative deviation of the X–ray surface brightness from the unperturbed value of the surface brightness at the same distance from the center. As expected the “trunk” of the mushroom is brighter than the surrounding regions.

Fig.9 illustrates the effect of entrainment of cold gas as the bubble rises. Shown in this figure are the distributions of tracer particles at the initial and final (after $``$67 Myrs) stages of the bubble evolution. The tracer particles fill the entire volume outside the bubble. The color assigned to each particle characterizes its initial distance from the cluster center. For example, “green” and “blue” particles seen in the upper torus have been uplifted from the central region of the galaxy or cluster atmosphere.

Abundance gradients are frequently observed in cluster atmospheres around central bright galaxies (e.g. Matsumoto et al. 1996). If our interpretation of the X–ray features as uplifted lumps of gas is correct, then one would expect the abundance in these regions to be higher than in the surrounding regions. This can be tested with Chandra and XMM.

### 3.5 The pancake stage

If the lumps of radio and X–ray emitting plasmas in the bubble are bound together (e.g. by magnetic stresses as in the multiphase cooling flow picture of Nulsen (1986)) then one can calculate the maximum radius to which the whole bubble can rise for given volume fractions of “entrained” gas and (weightless) radio emitting plasma (see eq. (6) in Churazov et al. 2000). Fig.10 shows the dependence of the final position of the bubble on the initial position and fraction of the bubble volume occupied by entrained gas which was captured near the initial location of the bubble. At this final location, the average mass density of the bubble (consisting of separate lumps of radio and X–ray emitting plasma) is equal to the density of the ambient gas. At the final position the bubble will spread along the equipotential (isodensity) surface and may form a pancake–like structure. The X–ray appearance of the bubble at this stage crucially depends on whether mixing of radio and X–ray emitting plasmas is microscopic or macroscopic (Böhringer et al. 1995, Nulsen 1997)

A simple simulation demonstrating pancake formation is presented in Fig.11. In this simulation, we assumed that the bubble is initially filled with gas having the same entropy as the cluster gas at a distance of $``$25 kpc. This is done by setting the temperature of the gas inside the bubble a factor of $``$2 higher than the temperature of the surrounding cluster gas. The density in the bubble was set a factor of 2 lower in order to maintain pressure equilibrium. The resulting modest contrast in entropy between the bubble and the surrounding gas approximately models the situation when the filling factors of the radio emitting plasma and thermal plasma in the bubble are comparable as considered in the previous section. The upper row of Fig.11 shows the entropy distribution at the beginning of the simulation and after $`2.6\times 10^8`$ years. In the lower panel, we show the distribution of tracer particles, which initially are located within the bubble, at the same times. The two black curves show, respectively, the boundary of the bubble at the beginning of the simulation and the region in the cluster where the cluster gas has the same entropy as the gas initially inside the bubble. It is clear that at the end of the simulation, tracer particles are concentrated along a surface where cluster gas has a similar entropy to that of the gas in the original bubble. Pancake–like structures are already visible after $`10^8`$ years from the start of the simulations. Oscillations near the neutral buoyancy point were also observed.

Strongly underdense bubbles may rise to much larger distances from the cluster center before reaching a neutral buoyancy point (see Fig.10). We did not simulate the evolution of strongly underdense bubbles to such late stages, but qualitatively the picture should be the same.

The largest structures visible in Figure 1 are the two almost circular low surface brightness emission regions to the northeast and southwest of the centre. Owen et al (2000) suggest that these are spherical bubbles of radio plasma seen in projection. In this case, they may receive substantial energy input from the inner lobes which contain the active jet. Alternatively one can identify the circular emission regions seen in the radio images with earlier bubbles that have reached their isodensity distance and are now in the process of spreading and forming pancake-like structures. We stress again that this requires that either relativistic and thermal particles are mixed on microscopic scales or the lumps of relativistic plasma and the ambient medium are bound together.

Owen et al. (2000) have noted the sharp boundary of the M87 halo which is traced by all observations at different frequencies and they argue that this excludes a dominant diffusive transport process for the relativistic plasma. In the proposed model the radius and sharpness of the boundary is naturally explained by the height to which buoyant bubbles can rise. Also the complicated structure of the outer halo may in this model reflect several layers of partly fragmented bubbles.

The time required for pancake formation at large distances from the cluster center is at least several times longer than that required for a strongly buoyant bubble to travel such distances. Therefore, it is extremely difficult for relativistic particles to survive (see 3.3). Either reacceleration or exposure of relativistic particles to stronger magnetic fields seems to be inevitable to explain the radio emission from the outer radio halo if it is indeed the result of lateral spreading of a bubble which has reached the neutral buoyancy point.

If, on the contrary, lumps of radio emitting plasma are not bound to the captured lumps of thermal gas, then radio bubbles will continue to rise along the pressure gradient and will eventually leave the central cooling flow region. Overdense lumps of captured gas will separate from the bubbles and fall back towards the central region. One possible example of overdense material having fallen back toward the central galaxy may be found in A1795. Cowie et al. (1983) reported a remarkable $`45^{\prime \prime }`$ (41 kpc long for the Hubble constant of 100 km s<sup>-1</sup> Mpc<sup>-1</sup>) filamentary emission line (H$`\alpha `$ \+ \[NII\]) chain originating near the central cD galaxy.

### 3.6 Total energy budget

The X–ray luminosity of the whole cooling flow region around M87 is about $`10^{43}`$ erg s<sup>-1</sup>. On the other hand, estimates of the total energy supplied by the jet to the ICM give a much larger value of $`10^{44}`$ erg s<sup>-1</sup> (e.g. Owen, Eilek and Kassim 2000). If all this energy can be dissipated into the thermal gas in the cooling flow region then the cooling flow should instead be a “heating” outflow (Owen, Eilek and Kassim 2000), at least at the present moment. The situation may be transient so that inflow and outflow periods alternate depending on the strength of the AGN activity (e.g. Tucker and David 1997, Binney 1999, Owen, Eilek and Kassim 2000). We speculate below on the consequences of a large input of power in the form of relativistic plasma into the cooling flow region.

Initially the radio lobe, inflated by the jet, expands supersonically (e.g. Heinz, Reynolds, Begelman 1998). Later, as the size of the radio lobe grows, the expansion becomes subsonic. As we discussed above (Section 3.2), a large and strongly underdense bubble would rise with approximately the Keplerian velocity of $`400`$ km s<sup>-1</sup>. Assuming that this velocity limits the growth of the original bubble due to the jet input power $`L`$, one can write: $`\frac{\mathrm{\Gamma }}{\mathrm{\Gamma }1}P_04\pi r^2v=L`$, where $`P_0`$ is the ambient gas pressure, $`\mathrm{\Gamma }=4/3`$ is the adiabatic index of the relativistic gas in the bubble, $`r`$ is the size of the bubble. Substituting $`L10^{44}`$ erg s<sup>-1</sup>, $`P_02\times 10^{10}`$ erg cm<sup>-3</sup> and $`v400`$ km s<sup>-1</sup> we find $`r5`$ kpc, which is approximately the size of the inner lobes surrounding the jets. If the subsonic expansion phase is longer than the duration of the initial strongly supersonic expansion, then most of the energy deposited by the AGN will go into the enthalpy of the expanded lobe and not into shock heating the cluster gas.

The inner lobes break loose (during a time $`tr/v10^7`$ years) and form bubbles which then subsonically rise through the cluster atmosphere (Gull and Northover 1973). At any given moment, thermal gas tries to settle in the potential of the cluster in order to form a nondecreasing entropy profile. Rising bubbles produce convective motions in the gas and mix lumps of thermal gas with different entropies. Using tracer particles added to the thermal gas (see Fig.9) we estimated that some 3–4$`\times 10^7M_{}`$ of cool, thermal gas has been uplifted from below 15 kpc to above 35 kpc by the bubble after $`8\times 10^7`$ years. Much larger masses of gas undergo smaller displacements by that time. For example, $`1.3\times 10^8M_{}`$ has been lifted from below 14 kpc to above 25 kpc. Although these are very crude estimates, which are sensitive to the particular initial and boundary conditions adopted in our simulations, one can expect that if the bubbles are launched approximately every $`10^7`$ years then they are able to transport some 1–10 $`M_{}`$ of cold gas per year from the center to the periphery of the cooling flow. Interestingly this value is comparable to the mass deposition rate in M87 estimated from the X–ray data. One can further list several important consequences of the convection induced by the rising bubbles:

* Due to uplifted cold gas, the X–ray surface brightness distribution should be shallower than in the absence of bubbles. The entropy profiles also flatten.
* If the bubbles are formed and are rising systematically over some preferred direction then the cold gas will be continuously uplifted over this direction while it will be flowing inward in other directions.
* Thermal instabilities can effectively take place within an environment of extensive convective motions. As a result, distributed mass deposition from cooling gas can naturally occur.
* Overdense lumps of thermal gas, uplifted by the rising bubbles, may fall back (when separated from the “carrier” buoyant bubble) to the central region thereby producing filamentary structures, stretching radially from the galaxy.

Even if the bubbles are not mixed on microscopic scales with the thermal plasma and are able to eventually leave the cooling flow region, most of the energy has to be deposited in the thermal gas (if there is a strong contrast in pressure between the inner and outer regions of the cooling flow). For a bubble moving at a constant terminal velocity the dissipation of energy is characterized by the change of the bubble enthalpy $`\frac{\mathrm{\Gamma }}{\mathrm{\Gamma }1}PV(P/P_0)^{\frac{\mathrm{\Gamma }1}{\mathrm{\Gamma }}}`$, where $`P`$ is the pressure of the cluster thermal gas. The bubble energy goes into sound waves, internal waves, turbulent motion in the wake, potential energy of the uplifted cold gas etc. Sound waves can carry energy away from the central region, while most of the other channels would eventually result in heating the cooling flow region. If the jet power is indeed $`10^{44}`$ erg s<sup>-1</sup>, then even a modest 10% efficiency of local dissipation into heat should be enough to exceed the radiative cooling of the gas. The gas then heats up, but this heating is distributed through a large volume. Since no shock waves are envisaged in this picture one would not expect to find localized, hot, compressed regions.

If strong energy input from AGN is maintained for a very long period (and dissipation into heat is efficient), then a core with a flat entropy profile will be formed and the size of this core will grow with time. We note here that the rate of energy extraction from the bubble depends on the pressure gradient (see above expression for enthalpy). Therefore, if the gas is heated to a temperature comparable to the depth of the local potential well, the dissipation efficiency will decrease in this region. This effect may provide a self-regulation for heating in the central part of the cooling flow region.

Eventually, strong heating is capable of completely eliminating the cooling flow, unless the power of an AGN is in turn regulated by the conditions in the cooling flow region as discussed by Tucker and David (1997), Binney (1999), Owen, Eilek and Kassim (2000). In M87, heating at a rate of $`5\times 10^{43}`$ ergs s<sup>-1</sup> would eliminate the cooling flow for $`10^8`$ years. This is comparable to the age of the largest radio structures. On the other hand, the complicated radio halo morphology suggests that periods of high AGN activity have alternated with the periods of low AGN power. We note here that dissipation of the subsonic motions into heat may take a longer time than the characteristic bubble rise time and therefore the actual gas heating rate may be related to the AGN activity averaged over a long period of time.

## 4 Conclusions

Buoyant bubbles of cosmic rays, slowly rising through the cooling gas, can qualitatively explain the complicated radio and X–ray morphology of the central 40 kpc region around M87. Torus–like features, seen in the radio, may be similar to the “mushrooms” which appear in Rayleigh–Taylor unstable configurations: as the fluid rises through the ambient medium, Kelvin-Helmholtz instabilities create torus-like “mushroom” heads. The excess X–ray emission trailing some prominent radio features could be due to cold gas captured by the rising bubbles and uplifted to large distances from the central source.

Rising bubbles produce strong convection in the cooling flow region. Convection flattens X–ray surface brightness and entropy profiles in the cooling flow. Convection may promote distributed mass deposition in the cooling flow and the formation of the filaments. If more than 10% of the energy supplied by the jet is dissipated into heat in the cooling flow region, then large cores with flat entropy profiles may be formed.

New Chandra and XMM observations of cooling flow clusters (e.g. McNamara et al. 2000, Fabian et al. 2000a) demonstrate that the interaction of the radio sources and thermal gas is widespread in these objects. Furthermore, the mass of cooling gas at temperatures less than $``$1 keV, derived from the spectra, seems to be significantly smaller than expected from the simplest cooling flow models (Böhringer et al. 2001, Tamura et al. 2001, Peterson et al. 2001). Possible explanations for this discrepancy are summarized by Fabian et al. (2000b; see also references therein). As we discussed above, heating by the buoyant bubbles may play an important role in the thermal balance in cooling flows. One of the predictions of the buoyant bubble picture, namely a correlated morphology of the colder uplifted gas and radio emission, seems to be in a broad agreement with recent XMM results (Belsole et al. 2001).

Our model is of course highly oversimplified and is not intended to closely reproduce the M87 environment. Clearly more detailed simulations are needed, perhaps in 3D, which would also follow the formation of the initial bubble (possibly with continuous energy input) and include radiative cooling. However, despite the simple treatment of bubbles in cluster and galaxy atmospheres presented above, we believe that this approach indicates areas of particular interest.

We are grateful to the editor Steven N. Shore and the referee for important comments and suggestions. We thank Alexei Kritsuk, Nail Inogamov, Christine Jones, Nail Sibgatullin and Henk Spruit for useful discussions. This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center.


text:
# The Size Distribution of the Kuiper Belt Objects

## 1 Introduction

After the first KBO was discovered in 1992 (Jewitt & Luu 1993), it is known that the outer solar system beyond Neptune is populated by small bodies. The total number of KBOs is increasing quickly and the general picture of the outer solar system is changed. There are three classes of KBOs according to the current observations. One of the population are the Resonant KBOs. They are trapped in 3:2 mean-motion resonance with Neptune at 39.4 AU. About one-third of the total population are in this class. The second class mainly occupies the region between 41 AU and 46 AU with small eccentricities. They are called the Classical KBOs and about two-third of KBOs belong to this class. Finally, several KBOs which behave very differently from the above two were discovered since 1996 (Luu et al. 1997) and form a new class. They are Scattered KBOs and they are moving on the large semi-major axises, highly eccentric, and inclined orbits. These KBOs of new dynamical class could originate from the scattering by Neptune.

One of the most important properties of KBOs is the size distribution because it tells the possible total number of KBOs and also gives the hint of density profile of Solar nebula. To extract the size distribution from the luminosity function, Jewitt et al. (1998) used Monte Carlo models to simulate the survey under the necessary assumptions about the geometric albedos and distance distributions of the KBOs. In principle, they assume a power-law size distribution as $`r^q`$ and determine the best value of $`q`$ by reproducing the luminosity function from the models. They found $`q=4.0\pm 0.5`$ and claimed this value of $`q`$ is consistent with $`q3.5`$ from what was expected from a collisional production function in Dohnanyi (1969). Kenyon & Luu (1999) did some planetesimal accretion calculations in a single annulus at 35 AU. They found that all models produce two power-law size distributions, $`q=2.5`$ for radii $`0.33`$ km and $`q=3`$ for radii $`13`$ km and the results are nearly independent of the initial mass in the annulus.

Different from the approach in Jewitt et al. (1998), we first get the number of KBOs against absolute magnitude histogram directly from the data of current known KBOs. Then, we convert this histogram into the size distribution. In addition to the size distribution for all current known KBOs, the Plutinos and Classical KBOs’ size distributions are studied separately. These results would be in Section 2 and there are concluding remarks in Section 3.

## 2 The Results

We get the data of current known KBOs from the daily updated List of Transneptunian Objects at the website: http://cfa-www.harvard.edu/iau/lists/TNOs.html and Figure 1 is the plot of eccentricity against semi-major axis for all KBOs from these data. There are 282 KBOs in the list on 27th July 2000. We regard the Plutinos to be those objects with semi-major axises from 38.1 AU to 40.5 AU. The Scattering KBOs are those objects with eccentricities $`e>0.3`$ and semi-major axis $`a>40.5`$ AU. All the rest are the Classical KBOs. Therefore, 63 of these current known KBOs are Plutinos, 212 objects belong to Classical KBOs and 7 objects are Scattering KBOs.

We use the table of converting absolute magnitudes to diameters at the same website to get the diameters. However, we transform the values of diameters in the table to be the values for albedo = 0.04. The histogram for number of KBOs against the absolute magnitude is then ploted. This plot can be transformed to the histogram for number of KBOs against the diameters, i.e. the size distribution easily.

### 2.1 All Known KBOs

Figure 2 is the result for all current known KBOs and it shows that there is a peak around 230 km. For those objects larger than this peak value, the size distribution follows

$$N(s)=\frac{150}{(s/200)^4}$$
(1)
approximately, where $`s`$ is the diameters (km) of the objects. For those objects smaller than the peak value, the size distribution becomes

$$N(s)=90(s/200)^3.$$
(2)
It is encouraging that those objects larger than 230 km follow the same size distribution as in Jewitt et al. (1998), i.e. $`q=4`$ power-law. However, the existence of this peak around 230 km and the different size distribution for objects smaller than the peak value has never been mentioned and discussed before. One might argue that this current observational result is not complete enough to confirm this peak. The size distribution would be completely different after all 100 km size KBOs are discovered. However, considering that all KBOs around or larger than 100 km are all within the detection limit and thus the detection probability is proportional to the number of objects. Equation 1 gives us $`N(100)/N(230)28`$ and thus the current known 100 km size KBOs should be about 28 times more than 230 km size KBOs if all sizes of KBOs follow the size distribution of Equation 1. It would be a puzzle that the current number of 230 km size KBOs is larger than the number of 100 km size KBOs if the size distribution is described by single power-law only.

Therefore, it is obviously that the peak we discover from the data of current known KBOs should be real and Equation 2 might be a good approximation for objects between 70 km to 230 km, though probably not for objects smaller than 70 km due to current detection limit.

### 2.2 Plutinos

As suggested by Malhotra (1995), the resonance sweeping due to the migration of Neptune might explain the trapping of KBOs into the 3:2 resonance. During this resonance sweeping process, it is possible that not all KBOs of different sizes can be trapped and become the Plutinos. Therefore, the size distribution for Plutinos might give us information about this.

According to the histogram in Figure 3, we found that for the objects larger than 150 km, the size distribution can be fitted by a $`q=4`$ power-law

$$N_\mathrm{p}(s)=\frac{17}{(s/200)^4},$$
(3)
for the objects smaller than 150 km, we have:

$$N_\mathrm{p}(s)=60(s/200)^4.$$
(4)
Comparing these results with Figure 2, it seems that the distribution is narrower and sharper around the peak and the peak is about 150 km. This might imply that there is a particular size, i.e. about 150 km, which is preferable to the resonance trapping process. This result of Plutinos’ size distribution might be important for the dynamical calculations of resonance sweeping.

### 2.3 Classical KBOs

Different from the Resonant and Scattering KBOs, the Classical KBOs might be the better tracer for the primordial Solar nebula because there were not trapped or scattered by Neptune and thus closer to where they were formed. Therefore, it would be interesting to know the size distribution for the Classical KBOs only.

According to the histogram in Figure 4, we found that for the objects larger than 230 km, the size distribution can be fitted by

$$N_\mathrm{c}(s)=\frac{140}{(s/200)^4},$$
(5)
for the objects smaller than 230 km, we have:

$$N_\mathrm{c}(s)=65(s/200)^3.$$
(6)
These results are pretty much the same as the results for all KBOs because at least two-third of all KBOs are in fact the Classical KBOs.

## 3 Concluding Remarks

We have used the current known data of KBOs to study their size distribution. We found that there is a peak around 230 km and therefore the size distribution is unlikely to be described by one single power-law. What does this observational result imply ?

The size distribution is related to the density of the Solar nebula because during the formation stage of KBOs, the amount of available material determined the final size and therefore the size distribution of KBOs. In fact, the initial density of the Solar nebula is one of the most important factors to lead the evolution of the Solar system. One possibility is to use the current mass distribution of planets to set up the ‘minimum mass Solar nebula’ model, then rescale to a higher density value but keep the same density profile because there must be certain amount of mass loss during the planet formation. From these assumptions, we can calculate that there was more than 20 Earth mass between 30 to 50 AU. However, the total mass of discovered KBOs is only about 0.07 Earth mass from the results in Section 2 or up to 0.26 Earth mass by Jewitt et al. (1998). Because of this contradiction, one might argue that the ‘minimum mass Solar nebula’ model could be wrong, especially for the outer Solar system.

One way to approach this problem is to study the formation of KBOs. Stern (1995) claimed that the timescale to produce the QB1s is longer than the disc lifetime and suggested that the present-day disc is not representative of the ancient structure. Kenyon & Luu (1999) did some accretion calculations and conclude that the ‘minimum mass solar nebula’ would be enough to form large KBOs in a reasonable timescale. Thus, from these calculations of Stern (1995), Kenyon & Luu (1999), we do need the outer Solar system to be at least as massive as the ‘minimum mass Solar nebula’ model to form the current observed KBOs through the accretion. Therefore, the outer Solar system between 30 to 50 AU was at least two order of magnitude more massive and the depletion in this region definitely happened during the formation of the Solar system.

On the other hand, the accretion calculations in Kenyon & Luu (1999) told us the size distribution for QB1s can approximately described by single power-law when these QB1s were formed. We therefore come to the conclusion that the size distribution should have changed during the depletion.

Jewitt (1999) suggested that for collisions to have removed 99 percent of the initial mass, one would have to postulate a steep initial size distribution ($`q5`$ in which objects large enough to escape commotion carry only 1 percent of the mass. This suggestion implies that the size distribution can be different in the past and also that the larger objects can survive easier than the smaller objects during the depletion. It is therefore possible that the peak we discovered was formed during this depletion process. A complete simulation to combine the accretion process and the orbital dynamics as in Jiang, Duncan & Lin (2000) might help the understanding of this depletion and thus the evolution of size distribution.

## Acknowledgments

Ing-Guey Jiang wishes to take this opportunity to acknowledge the suggestions on Solar System Dynamics from Scott Tremaine and the discussion with David Jewitt, Eugene Chiang, Li-Chin Yeh, Doug Lin during their visit at Academia Sinica, Institute of Astronomy and Astrophysics. He also thanks the hospitality and discussion of Martin Duncan and Doug Lin during his trips to Queen’s University and UC Santa Cruz.

Cheng-Pin Chen acknowledges the financial support of Undergraduate Summer Student Program from Academia Sinica, Institute of Astronomy and Astrophysics.

## REFERENCES

Dohnanyi, J. 1969, J. Geophys. Res., 74, 2531

Jewitt, D. 1999, Annual Review of Earth and Planetary Science, 27, 287

Jewitt, D., Luu, J. X. 1993, Nature, 362, 730

Jewitt, D., Luu, J. X. & Trujillo, C. 1998, AJ, 115, 2125

Jiang, I.-G., Duncan, M. J. & Lin, D. N. C. 2000, in progress

Kenyon, S. J. & Luu, J. X. 1999, AJ, 118, 1101

Luu, J. X., Marsden, B. G., Jewitt, D., Trujillo, C. A., Hergenrother, C. W., Chen, J., Offutt, W. B. 1997, Nature, 387, 573

Malhotra, R. 1995, AJ, 110, 420

Stern, S. A. 1995, AJ, 110, 856