text:
# Nuclear masses: evidence of order–chaos coexistence

## Abstract

Shell corrections are important in the determination of nuclear ground–state masses and shapes. Although general arguments favor a regular single–particle dynamics, symmetry–breaking and the presence of chaotic layers cannot be excluded. The latter provide a natural framework that explains the observed differences between experimental and computed masses.

PACS numbers: 21.10.Dr, 24.60.Lz, 05.45.Mt

Different approaches have been developed to reproduce the systematics of the observed nuclear masses. Some of them are of microscopic origin, while others, more phenomenological, are inspired from liquid drop models or Thomas–Fermi approximations. Often, the total energy is expressed as the sum of two terms

$$U(Z,N,x)=\overline{U}(Z,N,x)+\stackrel{~}{U}(Z,N,x),$$
(1)
where $`Z`$ and $`N`$ are the proton and neutron numbers, respectively, and $`x`$ represents a set of parameters that define the shape of the atomic nucleus. The first term $`\overline{U}`$ describes the bulk (macroscopic) properties of the nucleus, and contains all the contributions that vary smoothly with proton and neutron numbers. Its typical value is close to 8 MeV/$`A`$ throughout the periodic table, where $`A=N+Z`$ is the mass number. The second term $`\stackrel{~}{U}`$ describes shell effects, related to the shape–dependent microscopic fluctuations of the nuclear wavefunctions . A method to incorporate these effects in a consistent manner was originally proposed by Strutinsky .

Global nuclear mass calculations have been pursued over the years with increasing precision . Despite the numerous parameters contained in the models, the accuracy of the results obtained and the quality of the predictions are impressive. Different models and potentials yield similar results, and give an accuracy of $`5\times 10^4`$ for a medium–heavy nucleus whose total (binding) energy is of the order of 1000 MeV.

One may ask whether the difference between measured and computed masses has a particular significance, since it seems to be quite model independent. Among other possibilities, a natural explanation would come from many-body effects not included in the mean–field scheme. Our purpose here is to show that there is a very natural and appealing dynamical explanation. We will argue that $`\stackrel{~}{U}`$ may be splitted into two parts, $`\stackrel{~}{U}=\stackrel{~}{U}_{reg}+\stackrel{~}{U}_{ch}`$, and that the present calculated masses are able to correctly reproduce only $`\stackrel{~}{U}_{reg}`$. The two contributions $`\stackrel{~}{U}_{reg}`$ and $`\stackrel{~}{U}_{ch}`$ originate in regular and chaotic components of the motion of the nucleons, respectively. Although our calculations are done within a mean–field approximation, the final result for the fluctuations produced by the chaotic part of the motion are in fact of a much more general validity, and may be interpreted as arising from the residual interactions.

With essentially no free parameters, this dynamical symmetry–breaking mechanism provides a quantitative description of the basic observations on nuclear masses. First, the amplitude of $`\stackrel{~}{U}_{ch}`$ is much smaller than that of $`\stackrel{~}{U}_{reg}`$. Second, the typical size of the fluctuations $`\stackrel{~}{U}_{reg}`$ have a constant amplitude as a function of the number of nucleons. The RMS of $`\stackrel{~}{U}_{reg}`$ is found to be of about 3 MeV, which is in good agreement with the typical size of the difference $`\delta U=U_{exp}\overline{U}`$ between the experimental and the computed bulk properties of the nuclei . Third, the RMS of $`\stackrel{~}{U}_{ch}`$ is of order $`0.5`$ MeV, in good agreement with the typical size of the difference $`\delta U=U_{exp}U_{calc}`$ between the experimental and the calculated masses including shell effects . The mass-number dependence of this difference is also well reproduced.

From a semiclassical point of view, shell effects are interpreted as modulations in the single–particle spectrum produced by the periodic orbits of the corresponding classical dynamics . Given the single–particle energy levels $`E_j(x)`$ computed from some Hamiltonian $`H(x)`$, the level density $`\rho (E,x)=_j\delta \left[EE_j(x)\right]`$ is approximated by

$$\rho (E,x)=\overline{\rho }(E,x)+\stackrel{~}{\rho }(E,x).$$
(2)
The quantity $`\overline{\rho }`$ is the average density of states, whereas the oscillating part is expressed as

$$\stackrel{~}{\rho }=2\underset{p}{}\underset{r=1}{\overset{\mathrm{}}{}}A_{p,r}(E,x)\mathrm{cos}\left[rS_p(E,x)/\mathrm{}+\nu _{p,r}\right].$$
(3)
The sum is over all the primitive periodic orbits $`p`$ (and their repetitions $`r`$) of the single–particle Hamiltonian. Each orbit is characterized by its action $`S_p`$, stability amplitude $`A_{p,r}`$, and Maslov index $`\nu _{p,r}`$.

The shell correction to the nuclear mass is computed by inserting the oscillatory part of the density of states into the expression of the energy, $`\stackrel{~}{U}(x,A,T)=𝑑EE\stackrel{~}{\rho }(E,x)f(E,\mu ,T)`$, with $`f`$ the Fermi function. In a semiclassical expansion, the leading order of the integral comes from the energy dependence of the action, and the dependence of the prefactors can be ignored. Setting moreover the temperature $`T`$ to zero, one obtains

$$\stackrel{~}{U}(x,A)=2\mathrm{}^2\underset{p}{}\underset{r=1}{\overset{\mathrm{}}{}}\frac{A_{p,r}}{r^2\tau _p^2}\mathrm{cos}(rS_p/\mathrm{}+\nu _{p,r}),$$
(4)
where $`\tau _p=S_p/E`$ are the periods of the periodic orbits. The classical functions entering this expression are evaluated at the Fermi energy $`E_F`$, related to the mass number and shape parameters through the condition $`_0^{E_F}\overline{\rho }(E,x)𝑑E=A`$. According to Eq. (4), each periodic orbit produces a modulation or bunching of the single–particle states on an energy scale $`h/\tau _p`$. The presence of fluctuations in the total energy is therefore a very general phenomenon that occurs for an arbitrary Hamiltonian, irrespective of the nature of the corresponding classical dynamics. However, their importance (i.e., their amplitude) depends strongly on the properties of the dynamics, and in particular on its chaotic or regular character . In order to see this explicitly, we compute the variance of the fluctuations of the energy. The average is performed over a mass number window around a given nucleid, as when analyzing experimental data. The size of the window is taken to be small with respect to macroscopic quantities like the Fermi energy, but large compared to the typical oscillatory scale of $`\stackrel{~}{U}`$. This average is denoted by brackets. From Eq. (4) we obtain for the variance a double sum over periodic orbits labeled by the indices $`p`$ and $`p^{}`$. By considering the dependence of the density of periodic orbits with their period, it has been shown  that the dominant contribution to the double sum is given by the diagonal terms $`p=p^{}`$, $`r=r^{}`$. The resulting sum is convergent. Using moreover the semiclassical definition of the form factor,

$$K(\tau )=2h_0^{\mathrm{}}𝑑ϵ\mathrm{cos}(ϵ\tau /\mathrm{})\stackrel{~}{\rho }(Eϵ/2)\stackrel{~}{\rho }(E+ϵ/2)$$

whose diagonal part is

$$K_D(\tau )=h^2\underset{p,r}{}A_{p,r}^2\delta \left(\tau r\tau _p\right),$$
(5)
the energy variance may be written as

$$\stackrel{~}{U}^2\frac{\mathrm{}^2}{2\pi ^2}_0^{\mathrm{}}\frac{d\tau }{\tau ^4}K_D(\tau ).$$
(6)
This is the basic equation that provides the starting point for our analysis. An accurate evaluation of $`\stackrel{~}{U}^2`$ requires the knowledge of the periodic orbits, and therefore of the Hamiltonian. However, we do not want to make a calculation for a specific system. Instead we are interested in making general statements based on the statistical behavior of the single–particle orbits. That behavior depends on the regular or chaotic nature of the dynamics.

General considerations suggest that the single–particle motion in the nucleus should be dominated by regular orbits. This is based on the following arguments. There are several factors that are important when considering the amplitude of the shell effects as given by Eq. (4). The first one is the prefactor $`A_{p,r}/\tau _p^2`$. It implies that the more stable and short the periodic orbit is, the larger its contribution. The second is associated with the interference properties. Shells are wave effects: they are determined by a superposition of oscillatory contributions from different periodic orbits. The third factor, closely related to the previous one, is the phase space structure of the periodic orbits, i.e., their degeneracy. It depends on the existence of conserved quantities. In fully chaotic dynamics (where only the energy is conserved), the periodic orbits are isolated and unstable. In contrast, if there are as many conserved quantities as degrees of freedom, the motion is regular (or integrable) and periodic orbits come in degenerate families, all orbits having the same properties (action, stability, etc) within a family. Since many orbits contribute coherently, this degeneracy leads to an enhancement of the size of the shell corrections. In its equilibrium state, and for a given number of nucleons, the nucleus will adapt its shape in order to minimize its energy. This minimization procedure, that takes also into account the contributions from $`\overline{U}`$, favors regular single–particle dynamics that produce strong effects.

Equation (6) allows to estimate the typical fluctuations arising from a regular dynamics. In the range $`\tau _{min}\tau \tau _H`$, where $`\tau _{min}`$ is the period of the shortest periodic orbit in the semiclassical sum (4) and $`\tau _H=h\overline{\rho }`$ is the Heisenberg time, the form factor of the regular single–particle levels is given by $`K_D(\tau )=\tau _H`$ . This approximation is not valid for short times, of the order of $`\tau _{min}`$. The reason is that for $`\tau \tau _{min}`$, Eq. (5) produces system specific delta function peaks. If for simplicity we ignore this and extrapolate the statistical behavior down to $`\tau \tau _{min}`$, and moreover use $`K_D(\tau )=0`$ for $`\tau <\tau _{min}`$ (cf. Eq. (5)), then Eq. (6) leads to

$$\stackrel{~}{U}_{reg}^2=\frac{1}{24\pi ^4}gE_c^2.$$
(7)
The quantity $`E_c`$ is the energy associated with the shortest periodic orbit,

$$E_c=h/\tau _{min},$$
(8)
whereas the parameter $`g`$ measures this energy in units of the single–particle mean level spacing $`\delta =1/\overline{\rho }`$

$$g=E_c/\delta =\tau _H/\tau _{min};$$
(9)
$`g`$ counts the number of single–particle states on the scale $`E_c`$. To compute $`E_c`$ we need the period of the shortest orbit. Its length will typically be two or three characteristic nuclear dimensions. For simplicity we assume it to be two. Then, for a flat mean–field, like the Woods–Saxon potential, we have $`E_c=\pi E_F/k_Fr`$, where $`r`$ is the nuclear radius and $`k_F`$ is the Fermi wave vector. Since $`r1.1A^{1/3}`$ fm, we arrive at

$$E_c=\frac{77.5}{A^{1/3}}\mathrm{MeV}.$$
(10)
$`E_c`$ is the largest oscillatory energy scale where coherent bunching effects occur in any thermodynamic quantity. This result should be contrasted to the traditional shell–effect estimate based on the harmonic oscillator potential, $`E_c=\mathrm{}w_040/A^{1/3}`$ MeV .

Since $`\delta 2E_F/3A=25/A`$ MeV, Eq. (9) leads to

$$g=\pi A^{2/3}.$$
(11)
From these results we obtain for the typical fluctuations $`\sigma _{reg}=\sqrt{\stackrel{~}{U}_{reg}^2}`$ of a nucleus with a regular single–particle dynamics the value

$$\sigma _{reg}=2.84\mathrm{MeV},$$
(12)
an expression independent of the number of nucleons. This estimate is in good agreement with the $`3`$ MeV observed when $`\overline{U}`$ is subtracted out from the experimental or calculated values. The use of the real nuclear mean level spacing $`\delta `$ to compute (11) (i.e., including spin and isospin degrees of freedom) amounts to treating the different contributions as uncorrelated, and therefore provides a lower bound.

In deriving Eq. (12) a regular single–particle motion of the nucleons has been assumed. Even though we have stressed that a regular motion favors the minimization of the total energy, in the full many–body problem nothing guarantees the perfect integrability of the semiclassical dynamics. If some dynamical symmetries are broken, chaotic components will coexist with regular single–particle motion. The consequences of their presence on the behavior of the nuclear masses will now be discussed.

From a semiclassical point of view, the most simple approximation that can be made in the generic case of a mixed dynamics, when regular orbits coexist with chaotic layers, is to split the sum over periodic orbits in Eq. (4) into two terms, one from the regular part, another from the chaotic part. The shell energy is now written as

$$\stackrel{~}{U}=\stackrel{~}{U}_{reg}+\stackrel{~}{U}_{ch}.$$
(13)
The two terms are, from a statistical point of view, independent, $`\stackrel{~}{U}_{reg}\stackrel{~}{U}_{ch}=0`$. This happens because, as already pointed out, the dominant contribution to the energy comes from the short orbits, with $`\tau _p\tau _H`$. Since the orbits contributing to each term are different, the cross products vanish by the averaging procedure (assuming that the actions of the orbits are incommensurable). The decomposition (13) identifies the difference between measured and computed masses, that have a typical size of $`0.5`$ MeV, with the oscillatory contribution of the periodic orbits lying in the chaotic components.

We are therefore led to evaluate the variance of the shell corrections that originate in the chaotic layers. This can be done from the general Eq. (6). Again, our purpose is to make a statistical estimate valid for a generic chaotic component, with no reference to a particular system. In the range $`\tau _{min}\tau \tau _H`$ the form factor of chaotic single–particle levels is given by $`K_D(\tau )=2\tau `$ . It coincides with the random matrix prediction . This result is not valid for times of the order of $`\tau _{min}`$. Extending nevertheless this behavior down to $`\tau \tau _{min}`$, and imposing again $`K_D(\tau )=0`$ for $`\tau <\tau _{min}`$, we now arrive at

$$\stackrel{~}{U}_{ch}^2=\frac{1}{8\pi ^4}E_c^2.$$
(14)
Using Eq. (10), the typical size of the fluctuations $`\sigma _{ch}=\sqrt{\stackrel{~}{U}_{ch}^2}`$ of a chaotic component in the single–particle dynamics is given by

$$\sigma _{ch}=\frac{2.78}{A^{1/3}}\mathrm{MeV}.$$
(15)
This expression has to be compared with the RMS of the difference $`\delta U`$ between the experimental and the computed masses, given in Ref. . The comparison is shown in Fig. 1.

The agreement is extremely good. The amplitude in Eq. (15) is uncertain up to an overall factor of say, $`2`$. It can be varied by increasing slightly the period of the shortest orbit (we have chosen the shortest possible one; any modification will diminish $`\sigma _{ch}`$) and by the inclusion of spin and isospin (this increases $`\sigma _{ch}`$ by a factor $`2`$ if these components are treated as uncorrelated). The $`A`$ dependence is very well fitted in the region $`A\stackrel{>}{}75`$, with deviations observed for lower mass numbers. This is in agreement with the limited accuracy of the Strutinsky corrections and in general of semiclassical theories for light nuclei.

There are several features of the present theory that make its predictions reliable. First of all, Eq.(14) contains only one physical parameter, the period of the shortest chaotic periodic orbit, which is a function of $`A`$ because the size of the nucleus increases with the mass number. But it has no dependence on the relative size of the chaotic region, i.e., on the fraction of phase space occupied by chaotic motion. Without this quite remarkable and important fact we would have been forced to estimate that fraction, something that is hardly possible with present knowledge despite the efforts in this direction since the pioneering work of Ref.. This is in contrast to what happens for the regular regions. The enhancement factor $`g`$ present in Eq. (7) – but not in Eq. (14) – depends on the corresponding regular average level density $`\overline{\rho }_{reg}`$. The evaluation of $`g`$ in Eq. (11) assumes that the nuclear single–particle levels are regular, and the good qualitative agreement with the experimental results of $`\stackrel{~}{U}_{reg}^2`$ suggests that indeed most of the phase space is occupied by regular trajectories. Mixed systems present therefore a peculiar structure of shell effects. The amplitude of the fluctuations of the regular phase space regions is proportional, through the factor $`g`$, to the square root of the regular average level density, $`\sigma _{reg}\sqrt{\overline{\rho }_{reg}}`$. In contrast, the amplitude of the fluctuations associated with the chaotic motion is constant, independent of their phase space volume (aside from the transient nearly–integrable regime, not considered here, where the chaotic layers are formed). In the extreme case of a fully chaotic dynamics, $`\overline{\rho }_{reg}=0`$ and $`\stackrel{~}{U}=\stackrel{~}{U}_{ch}`$.

Equation (14) is in fact quite robust. It is not only independent of the chaotic phase–space volume, but is also valid for arbitrary dimensions. This fact suggests a very natural origin of these fluctuations, that we haven’t discussed yet. Although our analysis is based on a single–particle picture, it can be extended to the full many–body phase space. In that space, and in a semiclassical picture, to first approximation the point representing the system follows a very simple trajectory driven by the regular mean–field that dominates the motion. On top of that, it is likely that the residual interactions, not taken into account in that approximation, induce chaotic motion. The presence of chaotic orbits would then introduce additional long–range modulations in the regular single–particle density of states. Eq. (14) would then still be valid to evaluate the amplitude of those modulations, with $`\tau _{min}`$ the period of the shortest chaotic orbit in the multidimensional space. Rough estimates indicate that $`\tau _{min}`$ is comparable to the three–dimensional non–interacting period. Therefore Eq. (15) presumably gives an estimate of the mass fluctuations arising from neglected many–body effects.

To summarize, we have investigated shell effects of nuclear ground states within a unified framework, namely periodic orbit theory. We show that a dynamical symmetry–breaking mechanism that introduces chaotic layers in the motion of the nucleons produces additional shell corrections to the total energy. The intensity of the effect is, to first approximation, independent of the size or the topology of the layers. It is governed by a single parameter, $`E_c`$, which is proportional to the inverse time of flight of a nucleon across the nucleus at Fermi energy. The comparison of the typical size of the chaotic fluctuations is in very good agreement with the deviations between computed and experimental values, for the amplitude as well as for their dependence with the mass number. In contrast, the fluctuations produced by the regular components are independent of the mass number, and are proportional to the regular phase–space volume. The picture we are suggesting is coexistence of order and chaos as produced, for instance, by residual interactions. It manifests in nuclear low–energy properties. Further evidence should be given. In particular, our theory predicts autocorrelations in the total energies , as well as the effect of the presence of chaotic layers on the total level density. This deserves further investigation.

We are indebted to W. J. Swiatecki whose questions and persistence led to the present investigation, and to G. Bertsch and N. Pavloff for enlightening discussions and comments. The Laboratoire de Physique Théorique et Modèles Statistiques is an Unité de recherche de l’Université Paris XI associée au CNRS.


format: text

text:
# Instabilities in Molecular Dynamics Integrators used in Hybrid Monte Carlo Simulations

## 1 INTRODUCTION

The subject of instabilities in molecular dynamics (MD) integrators is not new. Studies of the phenomenon for the case of lattice QCD have been reported in . In  the lattices used were small ($`4^4`$ sites). The onset of instabilities was investigated as a function of the integration step size ($`\delta \tau `$) for heavy and light quark masses. For both masses studied the instabilities set in at some critical step size. The critical value of $`\delta \tau `$ was smaller for the lighter mass than for the heavier one. However, in both cases $`\delta \tau `$ was still so large that the energy change ($`\delta H`$) from the step size induced errors along an MD trajectory was sufficiently great to make simulations impractical in the unstable region of parameter space anyway.

Our study  using larger lattices and quark masses that are relatively light (by current standards) encountered instabilities at step sizes that are small enough to be relevant for large scale simulations. Our motivation to summarise these results is to highlight the impact these instabilities may have on future simulations using MD integrators as components, in particular HMC.

This article is organised as follows. In section 2 we demonstrate the instability in the simple harmonic oscillator (SHO) system. At this point we will restate the hypothesis of  to explain the onset of the instability. We present our main results in section 3 and show that they are completely consistent with this hypothesis. Finally in section 4 we will draw our conclusions.

## 2 SHO RESULT AND HYPOTHESIS

In the case of the SHO system a leapfrog update can be written as:

$$\left[\begin{array}{c}q\left(t+\delta \tau \right)\\ p\left(t+\delta \tau \right)\end{array}\right]=\left[\begin{array}{cc}1\frac{1}{2}\left(\omega \delta \tau \right)^2& \omega \delta \tau \\ \omega \delta \tau +\frac{1}{4}\left(\omega \delta \tau \right)^3& 1\frac{1}{2}\left(\omega \delta \tau \right)^2\end{array}\right]\left[\begin{array}{c}q\left(t\right)\\ p\left(t\right)\end{array}\right]$$
(1)
where $`q`$ and $`p`$ are the coordinates and conjugate momenta, $`t`$ is MD time and $`\omega `$ is the angular frequency.

The eigenvalues of the matrix in (1) are

$$\lambda _\pm =e^{\pm i\mathrm{cos}^1\left(1\frac{1}{2}\left(\omega \delta \tau \right)^2\right)}$$
(2)
When $`\omega \delta \tau <2`$ the inverse cosine in (2) is real, the eigenvalues are complex, and the phase space trajectories describe ellipses. However, when $`\omega \delta \tau >2`$, the inverse cosine becomes imaginary, and one eigenvalue is greater than one. At this point the phase space trajectories become hyperbolic, and the energy change along a trajectory grows exponentially.

The hypothesis of  is that the short distance modes of an asymptotically free quantum field theory behave like a collection of loosely coupled SHO modes. It is then anticipated that the MD integration of lattice QCD will go unstable when the highest frequency mode goes unstable.

We note that in the SHO example $`\omega `$ plays the role of the MD force $`F`$. Hence we surmise that for QCD the instability will occur when the combination $`F\delta \tau `$ reaches a critical value. Furthermore, we expect that the fermionic contribution to the force will increase in some manner with the inverse quark mass. This suggests that as the quark masses are decreased, the $`F\delta \tau `$ term will reach its critical value for smaller values of $`\delta \tau `$. This indeed is exactly what is reported in  and .

## 3 NUMERICAL RESULTS

Our numerical computations were carried out using 10 lattice configurations from a UKQCD lattice computation . The lattices were picked evenly spaced from a larger ensemble. The lattice volume was $`16^3\times 32`$ sites, and the physical parameters for the ensemble were $`\beta =5.2`$ for the gauge coupling, $`\kappa =0.1355`$ for the hopping parameter and $`c_{\mathrm{SW}}=2.0171`$ for the $`O(a)`$–improvement coefficient . These parameters correspond to a regime where the ratio of the pseudoscalar to vector masses is $`\frac{m_\pi }{m_\rho }0.6`$ .

We evolved the configurations for unit length MD trajectories using a variety of step sizes at a variety of $`\kappa `$ values. Along each trajectory we computed $`\delta H`$ as well as the 2-norms and $`\mathrm{}`$ norms of the gauge and fermionic contributions to the MD force ($`F_g_2`$,$`F_f_2`$, $`F_g_{\mathrm{}}`$, and $`F_f_{\mathrm{}}`$) averaged over all the time steps along the trajectory.

In figure 1 we show the results of fitting the fermionic force norms (averaged over configurations) to the fit ansatz:

$$F=C(am)^\alpha $$
(3)
in order to check that the forces increase in some inverse manner with the quark mass. To calculate $`am`$ we used the formula

$$am=\frac{1}{2}\left(\frac{1}{\kappa }\frac{1}{\kappa _{\mathrm{crit}}}\right)$$
(4)
where $`\kappa _{\mathrm{crit}}`$ is the critical value of $`\kappa `$ corresponding to the massless limit of the pion. In our fits we have left $`\kappa _{\mathrm{crit}}`$ as a free parameter, however to determine $`am`$ for the axes of figure 1 we used the value $`\kappa _{\mathrm{crit}}=0.13633`$, determined from spectroscopy on the ensemble .

One can see from figure 1 that the values of $`\alpha `$ are negative for both the 2–norm and the $`\mathrm{}`$–norm indicating that the forces do indeed grow inversely with $`am`$. The fits for $`\kappa _{\mathrm{crit}}`$ are consistent with the spectroscopic determinations indicating that the forces will diverge as $`am0`$.

In figure 2, we plot the variation of the force norms and $`\delta H`$ against the step size $`\delta \tau `$ at a fixed value of $`\kappa =0.1355`$ for a single configuration. It can be seen that at a value of about $`\delta \tau =0.0105`$, the fermionic forces start increasing and that $`\delta H`$ increases exponentially indicating the onset of the instability. The gauge contributions to the force norms seem to show no change in behaviour.

The $`\mathrm{}`$–norms grow at a larger rate than the 2–norm. As the $`\mathrm{}`$–norm is the largest absolute value contribution to the fermionic force, it plays a similar role to the highest frequency mode of a corresponding set of loosely coupled SHOs whereas the 2–norm of the force can be likened to the behaviour of the corresponding bulk modes. Hence this result is consistent with the earlier hypothesis, that the highest frequency modes drive the system unstable.

In figure 3 we show the behaviour of $`F_f_{\mathrm{}}`$, $`F_g_{\mathrm{}}`$ and $`\delta H`$ against the hopping parameter $`\kappa `$ for $`\delta \tau =0.01`$ and $`\delta \tau =0.012`$. When $`\delta \tau =0.01`$ the system is completely stable for all values of $`\kappa `$ (for this single configuration); however, for $`\delta \tau =0.012`$ the fermionic force norm rises sharply for $`\kappa =0.1355`$ accompanied by a large increase in $`\delta H`$ signalling the onset of the instability. This again is consistent with our earlier hypothesis that instabilities set in when the $`F\delta \tau `$ term reaches a critical value.

## 4 CONCLUSIONS AND DISCUSSION

We have exhibited the onset of an instability in the leapfrog MD integration scheme, used in HMC simulations and inexact simulation algorithms. The instability is manifest in the leapfrog scheme for a simple SHO system. We have also presented numerical evidence demonstrating that the instability is present also in lattice QCD systems. We have shown that our data is consistent with the hypothesis that the instability occurs when the $`F\delta \tau `$ term in the leapfrog equations reaches some critical value, and that the fermionic contribution to the force increases with decreasing quark mass. We have shown that as quark masses become lighter the instability sets in at smaller values of $`\delta \tau `$. Furthermore, for simulations with light dynamical quarks the instability sets in at values of $`\delta \tau `$ that are small enough to be relevant to large scale numerical simulations.

We anticipate that the major stumbling block for simulations with lighter quark masses will not come from “exceptional” configurations, where the value of $`\kappa `$ may become supercritical for some configurations, but rather from the onset of the instability which we expect to set in before $`\kappa `$ can become supercritical.

The instability problem is exponentially bad, but can be controlled easily by reducing $`\delta \tau `$. We are concerned for simulations with inexact algorithms, where the instability may be hard to detect as there is no $`\delta H`$ to monitor, and full finite step size extrapolations are seldom made.

Finally we wish to mention that neither the use of higher order integration schemes nor the use of double precision arithmetic are expected to alleviate the instability problem .

## 5 ACKNOWLEDGEMENTS

We acknowledge support from PPARC under grant number GR/L22744.


text:
# The dark side of star formation in galaxy clusters: spectroscopic follow-up of clusters observed with ISOCAM

## 1. Introduction

Since the classical study of Butcher & Oemler (1984) who put in evidence an increasing fraction of blue, presumably star-forming, galaxies as a function of cluster redshift, several authors have tried to estimate the importance of the cluster environment on the evolution of cluster galaxies. Recent studies with sample of distant clusters (Dressler et al. 1999, Balogh et al. 2000) pointed out that the star formation rate (SFR) per cluster galaxy appears to be lower than that in similar types of galaxies in the surrounding field. The main drawback of these studies is that the SFR rely on measurements of optical line fluxes (usually \[OII\] or in a few cases H$`\alpha `$) which suffer from strong dust extinction. Since dust obscuration depends on object types and on environmental conditions, it is impossible to correct it on average. On the other hand, although galaxies affected by extinction show typical spectral features (see Dressler et al. 1999), it is quite impossible to quantify the total SFR on the basis of the optical data alone. Mid-IR surveys, along with radio centimetric fluxes, provide the best estimates of the dust obscured SFR (see e.g. Chary & Elbaz 2001) when far-IR are not available. For this reason a sample of 10 clusters at different redshift ($`0.1<z<0.9`$) has been observed with the ISOCAM camera on-board ISO. We present here the optical/mid-IR properties of the galaxies detected at 15 $`\mu `$m in the A1689 cluster (see Fadda et al. 2000 for ISOCAM data) and preliminar results of the VLT follow-up program of J1888, a cluster at $`z=0.56`$ deeply observed with ISOCAM.

## 2. Hidden star formation in A1689

Using available images (HST and NTT) and about 100 spectra measured with the NTT (Duc et al. 2001), we studied membership, morphologies and spectral types of the objects detected by ISOCAM in A1689. Our analysis confirms the high fraction of blue galaxies in the cluster ( already detected by Butcher & Oemler 1984) relying on a secure cluster membership of the galaxies (spectroscopic and photometric redshifts are known for all the galaxies observed). Most of the 15 $`\mu `$m sources are luminous, blue, emission-line or morphologically disturbed galaxies, i.e. a population of galaxies usually associated with the “Butcher–Oemler” effect. However, $``$30% of the 15 $`\mu `$m sources do not show any sign of star-formation activity in their optical spectrum. More of 70% of the emission-line galaxies in our spectroscopic sample are detected at 15 $`\mu `$m and all the galaxies classified as dusty starbursts ( “e(a)” type in Dressler et al. 1999) are 15 $`\mu `$m sources. On the contrary, none of the galaxies with a post-starburst optical spectrum have been detected at 15 $`\mu `$m.

Since the AGN activity is very low in A1689, the 15 $`\mu `$m flux is a reliable tracer of the dust-obscured star formation activity (see Chary & Elbaz 2001). Comparing the SFR based on the mid-IR (15 $`\mu `$m flux) and optical ( \[OII\] line flux), we found that for galaxies with mid-IR emission the ratio SFR(IR)/SFR(Opt) is very high and ranges between 10 and 100, being the highest among “e(a)” galaxies (see Figure 1). The median SFR(IR) is 2 $`M_{}yr^1`$, while the median SFR(Opt) of the \[OII\] detected galaxies is 0.2 $`M_{}yr^1`$. Although we are not observing luminous infrared galaxies (the highest total IR luminosity $`6.2\times 10^{10}L_{}`$ is measured in a “e(a)” galaxy), we conclude that a significant portion of the star formation activity is visible only using mid-IR data. At least 90% of the star formation is missed when estimated from the \[OII\] line. A comparison with the field, essential to investigate the effects of the environment on the star formation of the cluster galaxies, will be possible only when a statistical sample of mid-IR galaxies for the coeval field of A1689 will be available.

## 3. Preliminary results in J1888

J1888, a cluster at $`z=0.56`$, has been surveyed by ISOCAM at the same depth of the Hubble Deep Field covering a strip of $`14.5^{}\times 3.5^{}`$ crossing the cluster center in order to explore possible variation of the mid-IR properties with the cluster-centric distance (see Figure 2). We have recently obtained about 120 spectra with VLT/FORS1 and a Ks image with NTT/SOFI, which complements B and R images obtained with the 2.2m ESO telescope. So far we have redshifts for 55 out of 70 sources detected at 15 $`\mu `$m and 22 out of 47 sources detected at 6.7 $`\mu `$m.

The redshift distribution of the cluster, which appears very loose in the optical images, peaks at $`z=0.56`$ but the field is highly polluted by background galaxies. Considering the range of magnitudes of the cluster members (21.4 $`<B<`$25.3), 55% of the galaxies do not belong to the cluster. Considering only the 15 $`\mu `$m sources, only 30% of the ISOCAM sources are cluster members.

The color-magnitude relation is not obvious to determine even for cluster members. Quite surprisingly there is a large number of “red outliers” ($`BR>2`$) which belong to the cluster and most interestingly several of them are detected at 15 $`\mu `$m. The presence of “red outliers” is not uncommon in clusters at $`z0.50.6`$ (see e.g. Margoniner & de Carvalho 2000) but in this case any contamination by field galaxies is ruled out. Contrary to the case of A1689, we have almost no mid-IR emitters among blue cluster members.

On the basis of VLT spectra we did a classification of galaxies according to the scheme of Dressler et al. (1999). As in the case of A1689, we detect at 15 $`\mu `$m all the dusty starbursts (“e(a)” type) and no post-starburst galaxies. Compared to A1689, we detect a larger fraction of spiral-like spectra which is not unexpected for clusters showing the “Butcher-Oemler” effect.

A detailed reanalysis of the J1888 ISOCAM data is now in progress using the new method developed by Lari et al. (2001) in parallel with the analysis of the deep ISOCAM surveys in the region of the Lockman Hole. The comparison of the data from this survey, whose redshift distribution peaks at redshift of 0.6 close to that of J1888, will be fundamental to study the influence of cluster environment on the evolution of galaxies.

## References

Balogh, M. L., Morris, S. L., Yee, H. K. C. et al. 1999, ApJ, 527, 54

Butcher, H. & Oemler, A. Jr. 1984, ApJ, 285, 426

Chary, R. & Elbaz, D. 2001, ApJ, 556, 562

Dressler, A., Smail, I., Poggianti, B. M. et al., 2000, ApJS, 122, 51

Duc, P.-A., Poggianti, B.-M., Fadda, D. et al., 2001, A&A, submitted

Fadda, D., Elbaz, D., Duc, P.-A. et al. 2000, A&A, 361, 827

Lari, C., Pozzi, F., Gruppioni, C. et al., 2001, MNRAS, 325, 1173

Margoniner, V. E. & de Carvalho, R. R. 2000, AJ, 119, 1562


text:
# Adaptive Random Walks on the Class of Web Graphs

## I Introduction

It has been recognized recently that the key role in the dynamics of complex systems is played by the evolving networks of various structures . Evolution of complex networks representing catalytic reactions, cell metabolic processes, relaxation in disordered systems, protein folding, and realistic ecological, social, and information networks  are guided by variety of principles inherent to each particular network. Technically, the study of the evolving networks includes theory of graphs, as opposed to regular lattices, which are appropriate mathematical objects in the study of driven cellular automata with diffusive dynamics . Both, the dynamics of nodes and arcs (links) contribute to the evolution of the network. Thus, network growth implies that a complex architecture of links emerges, which is peculiar to the network, as an effect of the individual microscopic rules at work.

Many complex networks self-organize into a scale-free structures of links (for a recent review see ). The appearance of a scale-free structure is not accidental, having that networks grow from stochastic rules, that assumes certain degree of disorder on the microscopic scale. This rules result in the occurrence of dynamic fractal structures on the growth time-scale , which is then manifested in the emergent hierarchical organization of node degrees at large evolution times. In growing networks the preference attachment such as in socially motivated linking, e.g., in science citation network , was shown to be one of dominant mechanisms  leading to the scale-free structure of links . However, in many other growing networks the preference attachment alone has no capacity to take into account the intimate relationship between structural complexity of the networks and their evolution. This implies dynamically changing wiring diagram, in sharp contrast to the rigid link structure in the science citation network, for example. Here we concentrate on this type of networks with the world-wide Web as a prominent example, in which constantly updated links yield rapidly changing wiring diagram .

On the other hand, the occurrence of various universality classes—characterized by the same scaling exponents—suggests that only certain details of the microscopic dynamics are relevant for the universal scaling behavior of complex networks, in analogy to critical behavior in physical systems in equilibrium . Although the theoretical background of the universality in complex dynamic systems and networks is still missing, numerical simulations on specific models  and master equation calculations  suggest that the following local properties (alone or in the appropriate combinations) of the dynamics are relevant for the emergent universal behavior: (a) Directed vs symmetric links (looked from a node in the network); (b) Fixed vs variable wiring diagram and time scale of its variations compared to the growth time scale; (c) Type of preference linking; and (d) Presence of constraints in the microscopic growth rules .

Recently we have introduced a model of growing directed graph with bias updates and biased attachment of links, motivated by the conduct of agents in the real world-wide Web . It was shown that the model captures minimum relevant properties for the dynamics of the Web, leading to a satisfactory quantitative comparison of the universal scaling quantities with the ones measured in the real Web . The relevant control parameter of the model $`\beta `$, defined as the average ratio of updated and added links per time step (thus measuring the degree of rewiring in the graph) can vary in the range $`(0,\mathrm{})`$.

The dynamic processes such as random walk specifically designed to the Web graphs reveal important information on network’s relaxation upon triggering and suggest how these networks can be efficiently explored. In Ref.  we introduced a random walk with adaptive walk strategy, that utilizes the locally available information on the underlying graph structure. In previous works  we have discussed the structure and dynamics of the world-wide Web, representing a directed graph in this class with rapidly changing wiring diagram: Estimated  value of $`\beta 34`$ from the scaling agreement with the observed behavior .

In this work we extend the study to the entire class of graphs generated by varying the control parameter $`\beta `$ in the physical range $`(0,\mathrm{})`$. In the limit of extreme rewiring ($`\beta \mathrm{}`$) some of the scaling features of the class disappear. On the other end, for $`\beta \beta _c(N)0.1`$ the graph undergoes a rigidity transition, where majority of links in the graph remains fixed in time. Among the networks between these two limits are the world-wide Web and networks representing a potential range of naturally occurring metabolic networks and catalytic reactions in open environment . Specifically, we study (1) how the global graph structure varies with the control parameter $`\beta `$ in the entire range before the rigidity transition, and (2) how these structural properties influence the relaxation of the network and the dynamics of the random walk with variable walk strategies. In addition, study of the random-walk dynamics reveal how the local clusters participate in the global behavior of the graph. We briefly discuss the potential application of the random-walk path in communication processes on these class of graphs.

The organization of the paper is as follows: In Sec. II we present the growth rules of the class of Web graphs. We grow large networks of $`5\times 10^6`$ nodes to determine the distributions of out- and in-degree for the range of values of the control parameter $`\beta `$. In Sec. III we define the adaptive random walk on this class of graphs and study the structure of connected subgraphs scanned by an ensemble of such walkers, which we then compare with the statistics of a standard random walk on the same graphs. In Sec. IV we determine the survival probability distribution of these two types of random walk and their access time to a given hierarchy level. Sec. V is devoted to a short summary and conclusions.

## II Growth and Structure of the Web Graph

A directed Web graph is grown from the dynamic rules proposed in Ref. , which are based on the criteria (a)–(c) listed above. (Potential effects of aging and other constraints on the microscopic dynamic rules are left out of this work). A growth step is defined by adding a node and subsequently creating and/or removing links in the entire network, that results in total increment $`M(t)`$ links at time step $`t=i`$. Naturally, we take $`M\overline{M(t)}>0`$, suggesting that, in the average, the number of links increase at each time step. The links are updated in the following way. An outgoing link originates from a node $`n<i`$ with the probability

$$p_1(n,i)=\frac{\alpha M+q_{out}(n,i)}{(1+\alpha )Mi},$$
(1)
which is determined by the number of outgoing links $`q_{out}(n,i)`$ previously created at that node. The link points towards the node $`k<i`$ with the probability

$$p_2(k,i)=\frac{\alpha M+q_{in}(k,i)}{(1+\alpha )Mi},$$
(2)
where $`q_{in}(k,i)`$ is the current number of incoming links collected at target node up to the step $`i`$. It is assumed that at the time of addition of a node $`i`$ to the network $`q_{out}(i,i)=q_{in}(i,i)=0`$. In addition a fraction $`f_0(t)\alpha M(t)`$ of new links are outgoing links from the new added node $`i=t`$. Whereas the remaining $`f_1(t)(1\alpha )M(t)`$ links are the updated links at preexisting nodes in the network. Hence, the relevant parameter in the model is the ratio of updated and added links at each time step, i.e., $`\beta f_1(t)/f_0(t)=(1\alpha )/\alpha `$, (independent of the actual increment $`M(t)`$), which measures the degree of rewiring (or flexibility) of the graph. The Eqs. (1)-(2) are motivated by the conduct of the agents who create links in the real Web: A new link is is likely to be created by the most active agents in the network, and, similarly, the most probable target node will be the one that already attracted many links. For consistency, we use $`M=1`$ throughout this work .

The numerical simulations in Ref.  demonstrated that the rules in Eqs. (1)-(2) are compatible with the observed scaling behavior of the real Web, when the single control parameter $`\beta `$ is set to a given value $`\beta 3`$, corresponding to $`\alpha 0.2`$ within estimated error bars. The dynamically emerging out- and in-degree cumulative distributions are given by the power-law functions

$$P(q_{out})q_{out}^{(\tau _{out}1)};P(q_{in})q_{in}^{(\tau _{in}1)},$$
(3)
in the asymptotic region for large number of links (i.e., large networks). The distributions are given in Fig. 1 for various values of the control parameter, showing the $`\alpha `$-dependent scaling exponents $`\tau _{out}`$ and $`\tau _{in}`$.

The asymptotic scaling behavior of the in-degree distribution that is generated by the shifted-linear preference rule in Eq. (2) with $`\tau _{in}=2+\alpha `$ was proved exactly in Ref. , that agrees completely with our results in Fig. 1. Recently, a similar analysis using rate equation approach was done in Ref.  for a directed graph in which in addition the probability of creating a link $`C(n,k)`$ from the preexisting node $`n`$ to node $`k`$ depends on the out-degree of the node $`n`$ and in-degree of the target node $`k`$. It was demonstrated that when $`C(n,k)`$ is given by a product of two shifted-linear functions of the corresponding degrees, the emergent distributions of in- and out-degree are power-law type and statistically correlated . According to the above rules in Eqs. (1)-(2), the probability for link creation $`C(n,k)`$ in our model is $`C(n,k)=p_1(n,i)\times p_2(k,i)`$, which is exactly the type of function compatible with a double-power law distributions for out- and in-degrees, in agreement with numerical results in Fig. 1. Emergence of statistical correlations between in- and out-degree distributions is discussed in Ref. . The model of Ref.  is slightly different from ours, involving three parameters that can not be reduced to a single control parameter, making it difficult to extract an exact value of the exponent $`\tau _{out}`$ in terms of our control parameter $`\beta `$ (or $`\alpha `$). The observed behavior shown in Fig. 1 suggests that $`\tau _{out}`$ increases approximately linearly with $`\alpha `$ for $`0.1<\alpha 0.5`$ (i.e., $`9>\beta 1`$). However, $`\tau _{out}`$ increases rapidly in the region below $`\beta 1/2`$, suggesting that the distribution of out-degree becomes exponential when the rigidity of the graph increases towards a finite critical value. Note that in the ideally rigid graph $`\beta 0`$ the out-degree distribution is trivial: each node has exactly $`M`$ links that remain fixed in time.

Study of the distribution of connected components done in Ref.  suggests that for $`\beta `$ larger than some critical value $`\beta >\beta _c(N)`$ a giant connected component occurs, similar to the real Web . By decreasing the control parameter $`\beta `$ the network structure changes gradually, eventually undergoing the phase transition into a structure without a giant component . The critical value of $`\beta `$ depends on the network size $`N`$ and in the cases studied here it is close to $`\beta _c1/12`$. Loss of scaling in the out-degree distribution at $`1/3>\beta >\beta _c(N)`$, while the in-degree still exhibits scaling behavior, demonstrates the power of cooperation in the network: Although the rules (1) and (2) formally look alike, the dynamical variations in $`q_{out}(n,i)`$ are due to local force, whereas the increase of $`q_{in}(k,i)`$ with time is made by a collective effect. Consequently, certain dynamical processes which are related to the out-degree structure might change their character before the network becomes rigid. Next we study random walk processes on this class of graphs.

## III Random Walks with Variable Strategy

We define two types of the random walks: a standard random walk, and a walk that learns its strategies from the information stored at a visited node on the graph, as defined below. First we grow the graph of $`N`$ nodes using the rules in Eqs. (1)-(2), and then start a walk from a randomly selected node, say node $`n`$. In the simple case, here called a naive random walk (NRW) for the reason to be clear soon, the walker moves along one of the $`q_{out}(n,N)`$ outgoing links of that node, selecting the link with equal probability

$$w(n)=1/q_{out}(n),$$
(4)
where we write $`q_{out}(n)q_{out}(n,N)`$, assuming for simplicity that the network does not grow during the walk time. In this way the walker selects a target node, say $`k`$, and moves there making one time step of the walk. Note that this standard random walk on the graph of the structure given in Eq. (3) is more complex compared to the case of hierarchical graphs with a constant branching ratio, in that it moves in an environment with variable structure of the relevant out-degree. An adaptive random walk (ARW), on the other hand, makes move selections with the statistical weight, which is correlated with the linking strategies of the visited node. In this case we assume that a link $`\mathrm{}`$ from node $`n`$ to node $`k`$ has a weight defined by

$$w_{\mathrm{}}(n)p_2(k,N),\mathrm{\Sigma }_{\mathrm{}=1}^{q_{out}(n)}w_{\mathrm{}}(n)=1,$$
(5)
where $`p_2(k,i=N)`$ is given in Eq. (2) with the normalization indicated in Eq. (5). Hence, the adaptive random walker uses the same type of strategy that was used earlier by the visited node, thus moving preferably along “hot” links of the visited node. Note that the probability $`w_{\mathrm{}}(n)`$ is not identical to $`p_2(k,i)`$ both because it is evaluated at a later time and it is normalized. In this way, the adaptive walker exploits full information about local graph structure, which is done by the out-degree of the visited node and in-degrees of all with that node connected target nodes . The naive random walker, on the contrary, uses only the first part of information. As a quantitative measure of different move strategies it is useful to consider the first-step statistics . We show in Fig. 2 the comparative statistics of the first jump of an ensemble of $`N_a`$ walkers with the naive and adaptive strategies, respectively, simulated on the same graph and using the same set of random numbers. It is shown in Fig. 2 that an adaptive search of the direction by the first jump is guided by the graph structure: the slope of the distribution of the first-jump distances (measured in node ranks) resembles the in-degree distribution. For the ensemble of this size the out-degree distances made in the first jump are correlated with the in-degree distances in view of Eq. (5). In the random search the distribution is still a power-law, however, the graphs structure is much less recognized (slope of the distribution of the first jumps is by one lower than the in-degree distribution of the graphs). Also much weaker, if any, correlations between in- and out-degree distances are observed. The walker continues as long as it is visiting nodes with nonzero number of out-links $`q_{out}>0`$.

### A Structure of connected subgraphs

A non-interrupted walk scans a set of connected nodes on the graph. The set of visited nodes is thus a connected subgraph smaller than a connected component that is usually searched by the Web crawl . An ensemble of random walkers thus can be used to scan the local structure of connected subgraphs, which in principle differs from the global graph structure.

The “distance” traversed by the walker on the graph is best measured in the node-rank differences, i.e., distance in the levels of hierarchy, in which the graph structure is nontrivial. Spatial distribution of nodes is, in principle, random, whereas the temporal fractal structure  occurs on the growing time scale. In Fig. 3 we show the time-integrated probability distributions of the node-rank distances for adaptive random walk on the graph with $`N=10^4`$ nodes and varying the control parameter $`\beta `$. An ensemble of $`N=20\times 10^4`$ walkers was employed. In Fig. 4 the corresponding distributions obtained by the ensemble of naive random walkers are shown that are simulated in the same conditions as the adaptive random walkers in Fig. 3.

The distributions of connected subgraphs found both by the ensemble of adaptive, $`W(\mathrm{\Delta }q)`$, and the ensemble of naive, $`R(\mathrm{\Delta }q)`$, random walkers are power-law type ($`\kappa =out,in`$):

$$W(\mathrm{\Delta }q_\kappa )(\mathrm{\Delta }q_\kappa )^{\delta _\kappa };R(\mathrm{\Delta }q_\kappa )(\mathrm{\Delta }q_\kappa )^{\rho _\kappa }.$$
(6)
Distinction between the respective distributions for in- and out-degree distances are resembling the underlying graphs structure. On the other hand, the difference in the scaling exponents, for instance in $`\delta _{out}`$ and $`\rho _{out}`$ and similarly with $`\delta _{in}`$ and $`\rho _{in}`$ suggest that ARW and NRW are two qualitatively distinct types of processes. Namely, both $`\delta _{out}>2`$ and $`\delta _{in}>2`$, whereas $`\rho _{out}<2`$ and $`\rho _{in}<2`$ indicating that the average distance in the case of the naive random walk diverges  with $`\mathrm{\Delta }q\mathrm{}`$ at $`N\mathrm{}`$, whereas the adaptive random walk remains in the confined areas-distances, $`<\mathrm{\Delta }q><\mathrm{}`$ when $`N\mathrm{}`$. These properties have important consequences, for instance for the disease spreading  along the random walk path on the graph.

In the case of ARW the similarity between the $`W(\mathrm{\Delta }q_{out})`$ and $`W(\mathrm{\Delta }q_{in})`$ and the corresponding distributions in Fig. 1 in the scaling region is striking. By increasing the probability $`\alpha `$ of a link outcoming from the new-added node, fast decrease of the distribution $`W(\mathrm{\Delta }q_{out})`$ in Fig. 3 makes the exponent $`\delta _{out}`$ difficult to measure for $`\alpha >0.6`$ . At a crossover value close to $`\beta 1/3`$ the distribution approaches the exponential form. At $`\alpha =1`$ all nodes have exactly one out-link, rendering this distribution trivial.

## IV Survival and Access Time of the Walk

Next we study the survival probability of a walk and the average access time of a walker to a given level of hierarchy. These properties are relevant for the relaxation of the graph and message passing to a given distance in hierarchy on the graph, respectively.

The access time is defined as a number of steps necessary for a walker to make a given distance here measured in the difference of the hierarchy levels of the departing and accessed node. In Fig. 5 the average access time normalized to the size of the ensemble $`N_a`$ is plotted against the difference in the hierarchy levels between departing and visited nodes for two values of the parameter $`\beta `$ and distances measured both by in- and out-degree. A remarkable feature is the power-law behavior of the access time . In general, the scaling form

$$<t_{acccess}>(\mathrm{\Delta }q_\kappa )^{\theta _\kappa }f(t/\mathrm{\Delta }q_\kappa )$$
(7)
applies both for adaptive and naive random walks when $`\beta `$ is large. However, the exponents $`\theta _\kappa `$ ($`\kappa `$ stands for “in,A” and “out,A” for the adaptive and similarly “in,R” and “out,R” for random walk) are both larger than 2 in the case of the adaptive search of the destination, whereas $`1<\theta _\kappa <2`$ for the naive random walk. Hence, apart from the small distances, the average access time of the adaptive random walk is much shorter compared to the standard random walk. This makes the adaptive random walk a good basis for its potential application as an efficient algorithm  for message-passing to a given level of hierarchy. In the case of catalytic reactions, this property of the adaptive random walk suggests that in the average fewer number of reactions occur along the path before a targeted reaction at a given hierarchy level is reached. Accordingly, the survival probability, $`P_s(t)`$ shown in Fig. 6, represents two different types of relaxation processes for large graph flexibility (large $`\beta `$): a nearly exponential relaxation in the case of ARW, and stretch-exponential for the NRW, respectively.

The difference between adaptive and random search for the target gradually diminishes with increasing rigidity of the graph (i.e., decreasing the degree of rewiring $`\beta `$), making the wandering on the graphs similar in both cases when the critical point is approached $`\beta \beta _c(N)`$. In Fig. 6 we show the probability $`P_s(t)`$ that a walk survives for $`t`$ steps on the graph, when the control parameter $`\beta `$ is varied. At the critical value ($`\beta _c0.081`$ for the graph sizes used in this simulations), the survival probabilities of the two types of walks almost coincide statistically. In the limit $`\beta =0`$, i.e., when $`\alpha =1`$ in Eqs. (1)-(2), each node has exactly one out-link, rendering the both walk strategies redundant. As we already pointed out, this limit corresponds to the graph with fixed links, belonging to another universality class compared to the Web graph that we discuss here.

## V Discussion and Conclusions

We have studied a class of directed graphs with variable wiring diagram caused, for instance, by frequent updates of the out-links in the world-wide Web while the graph evolves. The relevant control parameter that defines the class is the average ratio of updated vs. fixed links, $`\beta (1\alpha )/\alpha `$. In the case of the Web $`\beta `$ is estimated  in the range 3–4, i.e., to each added link in the Web comes in the average 3 to 4 links that connect preexisting nodes at each time step of the evolution. Other members of this class of graphs can be sought in the catalytic biochemical reactions, in which nodes are substrates (such as ATP, $`\mathrm{CO}_2,\mathrm{H}_2\mathrm{O}`$, etc.) and links are reactions among these substrates . Since the backward reactions, although present, are usually driven by much smaller rates than forward reactions , the biochemical reactions are represented by directed hierarchical graphs, with ATP as a top connected node. The presence of catalysts or enzymes selects appropriate reaction with a high accuracy. Investigations and potential applications of catalytic reactions among complex molecules is a great challenge of sciences in the future . The control parameter in this case can be roughly identified as the concentration of the catalyst, that can vary depending on the type of the reactants and/or catalyst.

Here we investigated within the growth model of Ref.  how the structure of the out- and in-connections of the graph vary with the control parameter and how that structure influences the relaxation processes on the graph. We summarize the properties that appear as direct consequences of the flexibility of the graph:

(1) In addition to in-degree, the out-degree distribution appears to be a power-law distribution with a new scaling exponent; Correlations between out- and in-degree spontaneously develop;

(2) These structural properties of the graphs determine the properties of the random walks on these graphs. Richness of the doubly-hierarchical connections allows us to define various types of the random walk strategies. Specifically, the adaptive random walker fully utilizes local information on both in- and out-connectivity of a visited node and adapts its moves to follow the node linking preferences. On the other hand, a standard random walker makes use of the out-degree of a visited node searching only the possible ways out. For large graph flexibility the adaptive strategy proves effective in the search of connected subgraphs and in reducing the access time to a given hierarchy level. Quick relaxation of the graph represented by the ARW path in this range of parameter indicates a short list of reactions before a targeted reaction occur.

The dynamic processes on this class of graphs reflect the underlying graph structure, leading to a variety of scale-free properties with distinct functional dependences on the in- and out- node degrees. Summary of all scaling exponents for varying flexibility (or rigidity) of the graph is shown in Fig. 7. To understand the theoretical basis of these dependences an analytical study in terms of rate equations of the random-walk dynamics is necessary.

(3) By increasing the graphs rigidity (reducing degree of rewiring $`\beta `$) the advantage of the strategy “reading full information on the way” is lost: Survival of the adaptive random walk becomes statistically closer to the one of the naive walk strategy; Access time gradually approaches the Poisson character; At a critical flexibility $`\beta _c0.1`$ the graph structure undergoes the rigidity percolation transition to another universality class characterized with predominantly fixed wiring diagram, in which the random walk dynamics becomes trivial.

(4)The other limit $`\beta \mathrm{}`$, representing the extreme flexibility, is formally interesting as a possible course of the world-wide Web evolution. In this limit the evolving network at each growth stage shares some similarity with the static random graphs! We detect certain numerical instability in the distributions that can be related to the changes in the linking probabilities Eqs. (1)-(2) when $`\alpha 0`$. After a large number $`N_1`$ nodes are added, these probabilities are $`N_1`$-dependent as $`p_\kappa \alpha /(1+\alpha )N_1+1/(1+\alpha )N_1^{(11/\tau _\kappa )}`$, where we assumed that maximum degree varies with the network size approximately as $`q_\kappa N_1^{1/\tau _\kappa }`$, (as before, $`\kappa `$ stands for “in” or “out”). In the random-graph theory the leading inverse linear term of linking probability is compatible with the occurrence of cycles (i.e., triangles, squares, pentagons, etc.) of all sizes, whereas the additional inverse sublinear term introduces admixture of full graphs . By letting $`\alpha 0`$ the leading linear term disappears, inducing sudden change in the structure of subgraphs—now being only full graphs. Now, having $`\tau _{in}2`$ for small $`\alpha `$ leads to $`p_{in}1/N_1^{1/2}`$, that is compatible with the occurrence of fully connected pentagons. Whereas $`\tau _{out}`$3–2.5, leading to $`p_{out}1/N_1^{2/3}`$$`1/N^{3/5}`$ allowing appearance of at most fully connected squares. This suggests why the exponent of the out-degree should stay larger than the one of in-degree distribution when $`\alpha 0`$. Numerical simulations suggest (see Fig. 7) that $`\tau _{out}2.5`$, whereas $`\tau _{in}2`$ in this limit.

Random walk path may serve as a communication channel in modeling of packet transport on this class of graphs. In addition to channels, the communication processes depend on the properties of packets to be transported and of driving rate resulting in load at individual nodes. By varying these parameters on a hierarchical tree it was shown in Ref.  that a continuous transition to a congested regime occurs, that exhibits several universal features. Looked in this context, our study concentrates on the properties of the communication channels only, and we do not specify any details on the character of packets. The simulation conditions thus correspond to idealized low driving rate below the transition point. Apart from the occurrence of cycles (closed loops), the Web graph differs from the hierarchical tree in that linked nodes are preferably at large distances of hierarchy, in contrast to the tree structure where only next hierarchy levels are coupled. We have demonstrated that due to the complex double-hierarchical structure of the Web graph the channel selection can be done using various strategies, each of them resulting in a nontrivial statistical properties of the channels. How these channel properties would influence the communication processes when more realistic packets and driving conditions are considered remains for the future study.

In conclusion, we have demonstrated that flexible wiring diagrams in the class of directed graphs induces a number of structural and relaxation properties that are crucial both for evolution of these graphs and for design of dynamic processes for exploring their structure.

###### Acknowledgements.

This work was supported by the Ministry of Education, Science and Sports of the Republic of Slovenia. I thank to Vyatcheslav Priezzhev for helpful suggestions.


text:
# Galaxy Cores as Relics of Black Hole Mergers

## 1 Introduction

Ferrarese et al. (1994) and Lauer et al. (1995) divide elliptical galaxies into two classes based on their nuclear properties, which Lauer et al. call “core” and “power-law” galaxies. Core galaxies exhibit a definite break in the surface brightness profile at some radius $`R_b`$; inward of this break, the logarithmic slope gently decreases in a manner that mimics a constant-density core. Power-law galaxies show essentially a single power-law profile throughout their inner regions, $`\mathrm{\Sigma }(R)R^\mathrm{\Gamma }`$, $`\mathrm{\Gamma }0.8\pm 0.2`$. The brightest galaxies, $`M_V<21`$, are exclusively core galaxies while faint galaxies, $`M_V>16`$ always exhibit power laws; galaxies of intermediate luminosity can exhibit either type of profile (Gebhardt et al., 1996). While the two categories were initially seen as distinct, nonparametric deprojection revealed that even the “core” galaxies exhibit power laws in their central space densities, $`\rho r^\gamma `$, with $`\gamma <\mathrm{\hspace{0.33em}1}`$ (Merritt & Fridman, 1995). Power-law galaxies have $`1<\gamma <\mathrm{\hspace{0.33em}2.5}`$ (Gebhardt et al., 1996). Furthermore the distribution of de-projected slopes is essentially continuous as a function of galaxy luminosity in the larger samples now available (Ravindranath et al., 2001; Rest et al., 2001).

Here we assume that the steep central density cusps of faint ellipticals and bulges, $`\rho r^2`$, are characteristic of the earliest generation of galaxies, and ask: How do the low-density cores associated with bright galaxies form? An appealing hypothesis links cores to nuclear black holes (BHs): in a galactic merger, the BHs will fall to the center of the merger remnant and form a bound pair, releasing their binding energy to the surrounding stars (Begelman, Blandford & Rees, 1980; Ebisuzaki, Makino & Okumura, 1991). High-resolution $`N`$-body simulations verify that this process can convert a steep power-law cusp, $`\rho r^2`$, into a shallow power-law cusp, $`\rho r^1`$, within the radius of gravitational influence of the BHs (Milosavljević & Merritt, 2001). Successive mergers would presumably lower the density of the core still more. In this model, power-law galaxies are those which have not experienced a major merger since the era of peak BH growth, or which have re-generated their cusps via star formation (Milosavljević & Merritt, 2001).

Preliminary tests of the cusp-disruption model were presented by Faber et al. (1997) and Milosavljević & Merritt (2001). The former authors plotted core properties (break radius, core luminosity) versus global properties in a sample of 19 early-type galaxies and noted a rough proportionality. In the paradigm investigated here, core properties should correlate more fundamentally with BH mass, since the mass of stars ejected by a decaying binary BH is expected to be of order the BHs’ mass. Milosavljević & Merritt (2001) used the new empirical relation between galaxy velocity dispersion and BH mass, the $`M_{}\sigma `$ relation (Ferrarese & Merritt, 2000; Gebhardt et al., 2000), to estimate BH masses in the Faber et al. sample. They found that a rough dynamical estimate of the “mass deficit” – the mass that would need to be removed from an initially $`r^2`$ density cusp in order to produce the observed profile – correlated well with $`M_{}`$.

In this paper we present the most careful test to date of the BH merger hypothesis for the formation of galaxy cores. We use nonparametric deprojection to compute the mass deficit in a sample of galaxies with high-resolution imaging data from HST (§2). We find a strong correlation between this mass and the mass of the nuclear BH; typical ejected masses are $``$ several $`M_{}`$ (§3). We argue (§4) that this result is consistent with the formation of cores via hierarchical mergers of galaxies containing pre-existing BHs. Cusps of non-interacting dark matter particles should behave in the same way as cusps of stars in response to heating by binary BHs, and we argue that the damage done to stellar cusps by this mechanism is a reasonable guide to the damage that would be done to dark matter cusps.

## 2 Data and Method

Our data set is drawn from a sample of 67 surface brightness profiles of early-type galaxies observed with HST/WFPC2 by Rest et al. (2001), and three additional galaxies: NGC 4472 and 4473, observed with WFPC1 by Ferrarese et al. (1994); and a WFPC2 F547M image of M87 (Jordán et al., 2002). The Rest et al. sample was selected from the set of all early-type galaxies with radial velocities less than 3400 km s<sup>-1</sup>, absolute $`V`$-band magnitudes less than $`18.5`$, and absolute galactic latitude exceeding 20 degrees. From this sample we excluded 13 galaxies for which central velocity dispersions were not available in literature; as discussed below, velocity dispersions were needed to compute BH masses in most of the galaxies. For specifics of the image-data reduction we refer the reader to the sources cited above. Surface brightness profiles used in this study were major-axis profiles. We applied a crude correction for the apparent ellipticity of the galaxies by multiplying the volume-integrated quantities, defined below, by $`(1ϵ_b)`$, where $`ϵ_b`$ is the ellipticity of the isophote at the break radius.

The intrinsic luminosity profiles $`\nu (r)`$ were obtained by deprojecting the PSF-deconvolved surface brightness profiles $`\mathrm{\Sigma }(R)`$ using the non-parametric MPL technique (Merritt & Tremblay, 1994). We opted for one-step deprojection via maximization of the penalized likelihood functional

$$_\lambda [\nu ]=\underset{i}{}\frac{\left(\mathrm{\Sigma }_iP_i[\nu ]\right)^2}{\left(\mathrm{\Sigma }_i^{err}\right)^2}\lambda _0^{\mathrm{}}\left[\frac{d^2\mathrm{log}\nu }{(d\mathrm{log}r)^2}\right]^2d\mathrm{log}r.$$
(1)
The first term compares the observed surface brightness data $`\mathrm{\Sigma }_i\mathrm{\Sigma }(R_i)`$ to the projections $`P_i[\nu ]`$ of the intrinsic luminosity estimate $`\nu `$ given by the operator

$$P_i[\nu ]=2_{R_i}^{\mathrm{}}\frac{\nu (r)rdr}{\sqrt{r^2R_{i}^{}{}_{}{}^{2}}}.$$
(2)
The second term in equation 1 is the penalty function which assigns zero penalty to any power-law $`\nu (r)`$; hence our estimate of $`\nu `$ is unbiassed if $`\nu `$ is an unbroken power law and should be minimally biassed if $`\nu `$ is approximately a power law. The relative strength of the penalty term is regulated through the parameter $`\lambda `$ which was chosen by eye and equal to $`0.01`$ for all galaxies; integrated quantities like the mass deficit defined below are only weakly dependent on $`\lambda `$.

Distances for 26 of the galaxies were drawn from the SBF survey (Tonry et al., 2001). For the remaining 34 galaxies we adopted distances computed assuming a pure Hubble expansion with $`H_0=80\text{ km s}^1\text{ Mpc}^1`$ corrected for Virgo-centric infall (Rest et al., 2001). Luminosity densities were converted to mass densities $`\rho (r)=\mathrm{{\rm Y}}\nu (r)`$ using the individual mass-to-light ratios $`\mathrm{{\rm Y}}_V`$ quoted in Magorrian et al. (1998) or their best-fit relation $`\mathrm{log}(\mathrm{{\rm Y}}_V/\mathrm{{\rm Y}}_{,V})=1.11\pm 0.33+(0.18\pm 0.03)\mathrm{log}(L_V/L_{,V})`$ for galaxies not included in that study.

We define $`\gamma d\mathrm{log}\rho (r)/d\mathrm{log}(r)`$ as the local, negative logarithmic slope of the deprojected density profile. Power-law galaxies are defined as those in which $`\gamma 2`$ at all radii; typically the profiles of such galaxies show no clear feature that can be identified as a “break radius” and are unlikely candidates for cusp destruction by binary BHs. In the remaining 35, “core” galaxies, the slope varies from $`\gamma >2`$ at large radii to $`\gamma <2`$ at small radii; the radius at which the slope crosses $`\gamma =2`$ in the positive sense ($`d\gamma /dr>0`$) is called here the “break radius” $`r_b`$ (Figure 1). This definition has little in common with the more standard definition based on fitting of the surface brightness profile to an ad hoc parametric function. In four galaxies the slope crosses $`\gamma =2`$ in the positive sense at more than one radius and thus the definition of $`r_b`$ is ambiguous. In such cases, we select the crossing toward larger radius from the largest dip of the slope below $`\gamma =2`$.

We define the mass deficit as the difference in integrated mass between the deprojected density profile $`\rho (r)`$ and a $`\gamma =\gamma _0=2`$ profile extrapolated inward from the break radius:

$$M_{\mathrm{def}}4\pi (1ϵ_b)_0^{r_b}\left[\rho (r_b)\left(\frac{r}{r_b}\right)^{\gamma _0}\rho (r)\right]r^2𝑑r.$$
(3)
Our choice of an $`r^2`$ density profile to characterize the “undisrupted” core is to a certain extent arbitrary; adiabatic growth of BHs can produce cusps with $`1.5<\gamma <\mathrm{\hspace{0.33em}2.5}`$ depending on initial conditions, and the faintest ellipticals with measured cusp slopes exhibit a similar range of $`\gamma `$’s (e.g., Gebhardt et al. (1996)). To test the sensitivity of our results to the assumed initial profile, we repeated the analysis using fiducial slopes of $`\gamma _0=1.75`$ and $`1.5`$. Values of $`\gamma _0>2`$ were found to exclude all but a few galaxies.

BH masses for a few of the galaxies in our sample are available from spatially-resolved kinematical studies (Merritt & Ferrarese, 2001). For all other galaxies we estimated $`M_{}`$ via the $`M_{}\sigma `$ relation,

$$M_{}1.4\times 10^8M_{}\left(\frac{\sigma _c}{200\text{ km s}^1}\right)^{4.8\pm 0.5}$$
(4)
(Ferrarese & Merritt, 2000), where $`\sigma _c`$ is the central velocity dispersion corrected to an aperture of $`r_e/8`$ with $`r_e`$ the effective radius. We used the aperture corrections of Jørgensen et al. (1995) to compute $`\sigma _c`$ from published values of $`\sigma `$ in Davies et al. (1987), Tonry & Davis (1981) and Di Nella et al. (1995).

## 3 Results

Results are given in Tables 1 and 2 and Figure 2; Table 1 gives mass deficits only for $`\gamma _0=2`$ while Figure 2 shows $`M_{\mathrm{def}}`$ computed using all three values of the fiducial slope, $`\gamma _0`$=(2, 1.75, 1.5). In the first panel of Figure 2 we have also plotted “dynamical” estimates of $`M_{\mathrm{def}}`$ for a sample of galaxies from Gebhardt et al. (1996), using equation (41) from Milosavljević & Merritt (2001), $`M_{\mathrm{dyn}}2(2\gamma )/(3\gamma )\sigma ^2R_b/G`$; $`M_{\mathrm{dyn}}`$ depends on the density profile only through $`R_b`$, the break radius of the surface bightness profile, and $`\gamma `$. When calculated for the deprojected galaxies in our sample, $`M_{\mathrm{dyn}}`$ was consistent within the scatter with $`M_{\mathrm{def}}`$.

For $`\gamma _0=2`$, the mass deficits are clustered about a linear relation defined by $`\mathrm{log}(M_{\mathrm{def}}/M_{})=0.92`$, $`1.0`$, and $`0.65`$, respectively, for Es and S0s; Es only; and S0s only, corresponding to $`M_{\mathrm{def}}(8.4,10,4.5)M_{}`$. Decreasing $`\gamma _0`$ decreases $`M_{\mathrm{def}}`$ (cf. Fig. 1b) and $`M_{\mathrm{def}}`$ becomes negative/undefined in galaxies when the minimum pointwise slope $`\gamma _{\mathrm{min}}`$ approaches $`\gamma _0`$. We do not cite values of $`\mathrm{log}(M_{\mathrm{def}}/M_{})`$ for these low values of $`\gamma _0`$ since the mean depends strongly on which galaxies are defined as having “cores.” However Figures 2b,c shows that $`M_{\mathrm{def}}`$ remains of order $`M_{}`$ or greater for many of the galaxies even when $`\gamma _0<2`$. We speculate that the lower mean value of $`M_{\mathrm{def}}`$ for the S0s may indicate a role for gaseous dissipation in the re-formation of cusps following mergers. We emphasize, however, that the mass deficit of lenticulars alone appears to be completely uncorrelated with the BH mass.

Table 2 lists linear regression fits to the data that were carried out using the routine of Akritas & Bershady (1996). For the Rest et al. (2001) ellipticals, as well as for the complete sample, the fitted power-law indices of the $`M_{\mathrm{def}}M_{}`$ relation calculated with $`\gamma _0=2.0`$ and $`1.75`$ are statistically consistent with unity. To test the sensitivity of the fitting parameters to the assumed power-law index $`\alpha `$ of the $`M_{}\sigma `$ relation (equation 4), we recomputed the fits with the Gebhardt et al. (2000) value $`\alpha =3.75\pm 0.3`$ in addition to the Ferrarese & Merritt (2000) value $`\alpha =4.85\pm 0.5`$ that is standard throughout the current paper. The effect of changing to a shallower version of the $`M_{}\sigma `$ relation is a steepening of the $`M_{\mathrm{def}}M_{}`$ relation, from $`d\mathrm{log}M_{\mathrm{def}}/d\mathrm{log}M_{}0.91`$ to $`d\mathrm{log}M_{\mathrm{def}}/d\mathrm{log}M_{}1.16`$.

We note that in most galaxies the break radius defined via $`\gamma (r_b)=\gamma _0=2`$ is close to the radius at which the density profile exhibits a visual break, i.e., where the curvature is greatest. The break radius, however, is not a good predictor of the mass deficit; in particular, some galaxies with large break radii, $`r_br_{}`$, with $`r_{}=GM_{}/\sigma ^2`$ the dynamical radius of the BH, have mass deficits below the mean. We speculate below that these large break radii may be produced by mechanisms other than BH mergers.

Although estimation of errors in Figure 2 is difficult, we believe that the scatter is at least partly intrinsic. Uncertainties in $`\mathrm{log}M_{}`$ are due primarily to uncertainties in $`\sigma `$ and are of order $`\mathrm{log}2`$. Uncertainties in $`\mathrm{log}M_{\mathrm{def}}`$ are also roughly $`\mathrm{log}2`$ based on variances between the redshift-based and SBF distances. A scatter greater than that due to measurement uncertainties would be reasonable given the different merger histories of galaxies with given $`M_{}`$. Similarly, if the progenitors of the galaxies exhibited a range of different $`\gamma _0`$, this could in itself explain the observed scatter in Figure 2.

## 4 Discussion

The mass ejected by a decaying BH binary is

$$M_{\mathrm{ej}}JM_{12}\mathrm{ln}\left(\frac{a_h}{a_{gr}}\right)$$
(5)
where $`M_{12}=M_1+M_2`$ is the binary mass, $`a_h`$ is the semi-major axis when the binary first becomes hard, and $`a_{gr}`$ is the separation at which the rate of energy loss to gravitational radiation equals the rate of energy loss to the stars (Quinlan, 1996). $`J`$ is a dimensionless mass-ejection rate; for equal-mass binaries, $`J0.5`$ (Milosavljević & Merritt, 2001). Quinlan (1996) claims that $`J`$ is nearly independent of $`M_2/M_1`$ even for extreme mass ratios, implying that a mass in stars of order $`M_1+M_2`$ is ejected during every accretion event, even when $`M_2`$ is tiny compared with $`M_1`$. This non-intuitive result is due to Quinlan’s ejection criterion, which for $`M_2M_1`$ includes stars that would not have gained sufficient energy to escape from the binary. If instead we equate the change in energy of the binary with the energy carried away by stars that are ejected with $`v>V_{bin}`$, we find $`M_{ej}M_2`$. This argument suggests a relation

$$M_{\mathrm{ej}}M_2\mathrm{ln}\left(\frac{a_h}{a_{gr}}\right)$$
(6)
which is consistent with equation (5) when $`M_1M_2`$.

Using equation (6), and adopting Merritt’s (2000) semi-analytic model for decay of a binary in a power-law cusp, we find

$$\frac{a_{gr}}{a_h}A\left|\mathrm{ln}A\right|^{0.4},A7.5\left(\frac{M_1}{M_2}\right)^{0.2}\frac{\sigma }{c}$$
(7)
and

$$M_{ej}4.6M_2\left[1+0.043\mathrm{ln}\left(\frac{M_2}{M_1}\right)\right],M_1M_2.$$
(8)
Thus $`M_{ej}/M_2`$ varies only negligibly with $`M_1/M_2`$. Henceforth we adopt $`M_{ej}5M_2`$.

If a BH grows by sequential accretion of smaller BHs, this result implies a mass deficit of order five times the final BH mass. However if the BH grows via a merger hierarchy of comparably-massive BHs, we expect $`M_{def}`$ to be larger. The idea here is that the damage done to cusps is cumulative: a merger of two galaxies whose cusps had previously been destroyed by binary BHs, will produce a shallower profile than a merger between two galaxies with initially steep cusps, even if the final BH mass is the same. Galaxies with masses $`M>\mathrm{\hspace{0.33em}10}^{11}M_{}`$, including most of the galaxies plotted on Figure 2, are believed to have undergone at least one major merger since a redshift of $`1`$ (e.g. Kauffmann, Charlot & Balogh 2001). Thus we predict $`M_{\mathrm{def}}>\mathrm{\hspace{0.33em}5}M_{}`$, consistent with Figure 2 if $`\gamma _0>\mathrm{\hspace{0.33em}1.5}`$.

Our interpretation of the mass deficit depends critically on the assumption that all of the change in $`\gamma `$ during a merger can be attributed to the BHs, i.e., that cusp slopes remain unchanged the absence of BHs. This is known to be the case in equal-mass mergers between galaxies with power-law cusps (Barnes, 1999; Milosavljević & Merritt, 2001), though in mergers with extreme mass ratios, features can appear in the density profile that are not due to BHs (Merritt & Cruz, 2001). We speculate that the break radii in some of the galaxies in our sample may be due to this process, particularly those galaxies (e.g., NGC 3640, 4168) where $`r_b`$ greatly exceeds the radius of gravitational influence of the BH. $`N`$-body simulations of cumulative mergers with unequal BH masses will be needed to assess this hypothesis.

Our model presents an interesting contrast to that of van der Marel (1999), who proposed that cores (in the sense of constant-density regions) were present in all galaxies ab initio, and that power-law cusps were generated by the growth of the BHs – roughly the opposite of our model in which BHs destroy pre-existing cusps. Van der Marel assumed that core mass correlated initially with bulge luminosity as $`M_{\mathrm{core}}L^{1.5}`$ and that $`M_{}L`$; hence $`M_{\mathrm{core}}M_{}^{1.5}`$, consistent with the correlation in Figure 2 if we identify $`M_{\mathrm{core}}`$ with $`M_{\mathrm{def}}`$. We believe that this agreement is coincidental. Van der Marel’s model relates core mass to BH mass via an ad hoc postulate, while the model discussed here contains a mechanism for core formation. Van der Marel also ignored the effects of mergers. Nevertheless, van der Marel’s model shows that our interpretation is not unique.

If stellar cusps are destroyed by binary BHs, the same should be true of dark-matter cusps, like those predicted in CDM theories of structure formation (e.g., Navarro, Frenk & White (1996); Moore et al. (1998); Bullock et al. (2001)). Destruction of dark matter cusps could be very efficient if supermassive BHs were present in dark matter halos at large redshifts (e.g. Fan et al. (2001); Haiman & Loeb (2001); Menou, Haiman & Narayanan (2001)) due to the cumulative effect mentioned above; furthermore, dark matter cusps would not be regenerated the way that stellar cusps might be via star formation. If our model for the formation of stellar cores is correct, we would predict tha cores of non-interacting CDM should be about as large as those observed in the stars, and perhaps larger.

We thank L. Ferrarese for her generous help with the data reduction and interpretation. The image of M87 used for calculating the mass deficit of that galaxy (Section 2) was kindly provided by P. Côté and A. Jordán in advance of publication. This work was supported by NSF grant AST 00-71099 and NASA grants NAG5-6037 and NAG5-9046 to DM.