text:
# 1 Introduction

## 1 Introduction

In this paper we investigate the constraints that can be posed on the existence of a possible fifth force by the analysis of the laser–ranged data to LAGEOS II and LAGEOS Earth satellites. We will consider a potential energy including a Yukawa term of the form \[Ohanian and Ruffini, 1994; Ciufolini and Wheeler, 1995; Nordvedt, 1998\]

$$U=U_0+U_\mathrm{Y}=\frac{GMm}{r}\left(1+\alpha e^{\frac{r}{\lambda }}\right),$$
($`1`$)
where $`G`$ is the Newtonian gravitational constant, $`M`$ is the mass of the central body, $`m`$ is the mass of the orbiting test particle, $`r`$ is the distance between the two bodies, $`\alpha =\frac{Kk}{GMm}`$ contains the couplings $`K`$ and $`k`$ of the new force to the two bodies111In general, $`K`$ and $`k`$ are not proportional to $`M`$ and $`m`$. and $`\lambda `$ is the finite range of the new force.

In checking the nature of the fifth force it is of the utmost importance to perform experiments spanning the widest range of length scales as possible \[Nordvedt, 1998\]: for experiments at laboratory scale see \[Krause and Fischbach, 2001\] and references therein. Using the orbits of LAGEOS and LAGEOS II satellites implies that we are testing the hypothesis of the fifth force at a length scale of almost two Earth radii, i.e. 104 km: at this scale the constraints on $`\alpha `$ are of the order of $`|\alpha |<10^510^8`$ (see Fig. 3.2 (a) of \[Ciufolini and Wheeler, 1995\]) and are derived from a data analysis of LAGEOS.

Our analysis includes also an evaluation of the error budget in order to account for various systematic errors induced by several classical aliasing forces. We will show that it is possible to improve sensibly the present limits by using suitably the data from the existing LAGEOS and LAGEOS II satellites and the present or near future knowledge of the terrestrial gravitational field whose uncertainties represent, as we will see later, the main sources of systematical errors.

The paper is organized as follows. In section 2 we derive the effects of the Yukawa gravitational potential on the orbit of a test body with the standard technique of the Gauss perturbative equations for the rates of change of the Keplerian orbital elements. In section 3 we apply the results obtained in section 2 to the Earth–LAGEOS system and discuss the constraints posed on $`\alpha `$ by using a suitable combination of the residuals of the perigee of LAGEOS II and the nodes of LAGEOS II and LAGEOS. The effects of various sources of systematical errors are investigated. The role of the proposed LARES satellite is considered as well. Section 4 is devoted to the conclusions.

## 2 The orbital effects of the Yukawa perturbation

The acceleration felt by a test body orbiting the mass $`M`$ in the potential energy given by eq.($`1`$) is

$$\text{a}=\text{a}_0+\text{a}_{\mathrm{pert}}=\frac{GM}{r^2}\widehat{r}\alpha GM\left(\frac{1}{r^2}\frac{1}{\lambda ^2}\right)\widehat{r}.$$
($`2`$)
In obtaining eq.($`2`$) it has been assumed that $`e^{\frac{r}{\lambda }}1\frac{r}{\lambda }`$, as it should be the case for an Earth orbiting satellite with $`r10^7`$ m. In eq.($`2`$) $`\widehat{r}`$ is the unit vector pointing from the central mass to the orbiting body.

The second term of the right–hand side of eq.($`2`$) can be considered a small, central perturbation of the Newtonian monopole acceleration. It may be interesting to note that, since, in general, $`\alpha `$ may vary from a body to another in the field of the same mass222This happens if $`k`$ is not proportional to $`m`$., eq.($`2`$) implies that different bodies may be accelerated differently in the field of $`M`$ \[Nordvedt, 1998\].

Let us work out explicitly the effects of the Yukawa perturbing acceleration on the orbit of an artificial satellite. We will adopt the standard approach based on the Gauss perturbative equations and the projections of the disturbing acceleration $`R,T,N`$ onto the radial, along–track and cross–track mutually orthogonal directions \[Milani et al., 1987\], respectively. The Gauss equations are

$`{\displaystyle \frac{da}{dt}}`$ $`=`$ $`{\displaystyle \frac{2}{n\sqrt{1e^2}}}\left[Re\mathrm{sin}f+T{\displaystyle \frac{p}{r}}\right],`$ ($`3`$)
$`{\displaystyle \frac{de}{dt}}`$ $`=`$ $`{\displaystyle \frac{\sqrt{1e^2}}{na}}\left[R\mathrm{sin}f+T\left(\mathrm{cos}f+{\displaystyle \frac{1}{e}}\left(1{\displaystyle \frac{r}{a}}\right)\right)\right],`$ ($`4`$)
$`{\displaystyle \frac{di}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{na\sqrt{1e^2}}}N{\displaystyle \frac{r}{a}}\mathrm{cos}(\omega +f),`$ ($`5`$)
$`{\displaystyle \frac{d\mathrm{\Omega }}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{na\mathrm{sin}i\sqrt{1e^2}}}N{\displaystyle \frac{r}{a}}\mathrm{sin}(\omega +f),`$ ($`6`$)
$`{\displaystyle \frac{d\omega }{dt}}`$ $`=`$ $`\mathrm{cos}i{\displaystyle \frac{d\mathrm{\Omega }}{dt}}+{\displaystyle \frac{\sqrt{1e^2}}{nae}}\left[R\mathrm{cos}f+T\left(1+{\displaystyle \frac{r}{p}}\right)\mathrm{sin}f\right],`$ ($`7`$)
$`{\displaystyle \frac{d}{dt}}`$ $`=`$ $`n{\displaystyle \frac{2}{na}}R{\displaystyle \frac{r}{a}}\sqrt{1e^2}\left({\displaystyle \frac{d\omega }{dt}}+\mathrm{cos}i{\displaystyle \frac{d\mathrm{\Omega }}{dt}}\right),`$ ($`8`$)
where $`a,e,i,\mathrm{\Omega },\omega `$ and $``$ are the satellite’s semimajor axis, eccentricity, inclination, longitude of the ascending node, argument of perigee and mean anomaly, respectively. Moreover, $`p=a(1e^2)`$, $`f`$ is the true anomaly and $`n=\sqrt{GMa^3}`$ is the Keplerian mean motion. As can be noticed from eq.($`2`$), the Yukawa disturbing acceleration has only the in–plane, radial component $`R`$, so that it can be straightforwardly inferred from eqs.($`5`$)-($`6`$) that the out–of–plane Keplerian orbital elements like the inclination $`i`$ and the longitude of the ascending node $`\mathrm{\Omega }`$ are not affected by it. By evaluating the Yukawa acceleration on the unperturbed Keplerian ellipse, for which

$$r=\frac{a(1e^2)}{1+e\mathrm{cos}f},$$
($`9`$)
inserting it in eqs.($`3`$)-($`8`$) and, then, averaging them over an orbital revolution, the long–period rates of change of the Keplerian orbital elements can be obtained. It turns out that, by neglecting terms of order $`𝒪(e^n)`$, with $`n2`$, in the satellite’s eccentricity, only the perigee $`\omega `$ and the mean anomaly $``$ are affected by long–term Yukawa perturbations. Indeed, from eq.($`2`$) and eqs.($`7`$)-($`8`$) it can be obtained

$`{\displaystyle \frac{d\omega }{dt}}`$ $`=`$ $`\alpha {\displaystyle \frac{n}{(1e^2)^{\frac{3}{2}}}},`$ ($`10`$)
$`{\displaystyle \frac{d}{dt}}`$ $`=`$ $`n+\alpha {\displaystyle \frac{n}{(1e^2)}}\alpha {\displaystyle \frac{2GM(1e^2)}{na\lambda ^2}}.`$ ($`11`$)
## 3 The constraint on $`\alpha `$

In order to constraint effectively the Yukawa coupling $`\alpha `$, let us focus on eq.($`10`$). It tells us that the Yukawa perturbation induces on the perigee of a near Earth satellite a secular rate which, for LAGEOS II, is proportional to $`3.06658517\times 10^{12}`$ milliarcseconds per year (mas/y) via the coupling constant $`\alpha `$. What is the constraint posed on $`\alpha `$ by the experimental accuracy with which it can be possible to measure the perigee rate? In the case of this Keplerian orbital element the observable quantity is $`r=ea\omega `$. So, by assuming an experimental error of, say, $`\delta r_{\mathrm{exp}}^{\mathrm{II}}`$ =1 cm over 1 year, for LAGEOS II, which has $`e=0.014`$ and $`a=1.2163\times 10^9`$ cm, we have $`\delta \omega _{\mathrm{exp}}^{\mathrm{II}}=12`$ mas333The other existing passive geodetic laser–ranged satellites are unsuitable because their eccentricities are smaller than LAGEOS II, except for Starlette whose orbit, however, is known less accurately than that of LAGEOS II for various reasons.. This yields a relative accuracy on $`\alpha `$

$$\left(\frac{\delta \alpha }{\alpha }\right)^{\mathrm{exp}}4\times 10^{12}.$$
($`12`$)
It is interesting to note that the same estimate for the proposed LARES satellite \[Ciufolini and Matzner, 1998\], which is planned to have a larger eccentricity, $`e_{\mathrm{LARES}}=0.04`$, would yield a relative accuracy of almost $`1\times 10^{12}`$.

However, this estimate does not include any systematic errors. As it is well known from various proposed or performed tests of General Relativity with LAGEOS satellites \[Ciufolini et al., 1997; 1998; Ciufolini, 2000; Iorio and Pavlis, 2001; Iorio, 2001; Iorio et al., 2002; Iorio, 2002\], in such kind of measurements there are lots of competing classical forces which may act as superimposed biases affecting sensibly the precision of the measurements. Then, the evaluation of the systematical errors induced by them is of the utmost importance.

The main source of systematic errors is represented by the mismodelled precessions of the even zonal harmonics of the static part of the geopotential. In particular, the first two even zonal harmonics $`J_2`$ and $`J_4`$ are the most insidious. In order to cancel their impact, as proposed in the PPN LAGEOS experiment \[Iorio et al., 2002\], the following combination of orbital residuals could be used444Notice that, since LAGEOS enters the combination of eq.($`13`$) only with its node $`\mathrm{\Omega }`$, and since the Yukawa perturbation does not affect such Keplerian orbital element, eq.($`13`$) is insensitive to the possible differential falling of LAGEOS and LAGEOS II in the terrestrial gravitational field.

$$\delta \dot{\omega }^{\mathrm{II}}+c_1\delta \dot{\mathrm{\Omega }}^{\mathrm{II}}+c_1\delta \dot{\mathrm{\Omega }}^\mathrm{I}=\alpha x_{\mathrm{Yuk}},$$
($`13`$)
where

$`c_1`$ $`=`$ $`0.86,`$ ($`14`$)
$`c_2`$ $`=`$ $`2.85,`$ ($`15`$)
$`x_{\mathrm{Yuk}}`$ $`=`$ $`3.06658517\times 10^{12}\mathrm{mas}/\mathrm{y}.`$ ($`16`$)
The coefficients of eq.($`13`$) depend on the orbital parameters of LAGEOS and LAGEOS II and are obtained in order to cancel the contributions of the first two even zonal harmonics of the geopotential to the measurement of $`\alpha `$. It is intended that the residuals would account for the Yukawa perturbation in the sense that it would be viewed as an unmodelled feature not included in the force models adopted in fitting the satellites’ orbits. According to the covariance matrix of the most recent available Earth gravity model EGM96 \[Lemoine et al., 1998\], the systematic error induced by the uncancelled even zonal harmonics of the geopotential amounts to555This estimate has been obtained by considering the geopotential harmonics up to degree $`l=20`$. This is well justified by the insensitivity of LAGEOS satellites to the higher degree terms. Moreover, this fact makes our estimate reliable because the higher degree terms of geopotential in EGM96 are not particularly well determined.

$$\left(\frac{\delta \alpha }{\alpha }\right)^{\mathrm{zonals}}7\times 10^{12}.$$
($`17`$)
Regarding the time–dependent part of the Earth gravitational field, the estimates of \[Iorio et al., 2002\], adapted to this context, yield for an observational time span of 5 years

$$\left(\frac{\delta \alpha }{\alpha }\right)^{\mathrm{harmonics}}1\times 10^{12}.$$
($`18`$)
The most relevant non–gravitational perturbations are the direct solar radiation pressure and the Earth’s albedo. Their impact can be evaluated from \[Lucchesi, 2001\] for 5 years as

$$\left(\frac{\delta \alpha }{\alpha }\right)^{\mathrm{non}\mathrm{grav}}8\times 10^{12}.$$
($`19`$)
The errors induced by the direct solar radiation pressure and the Earth’s albedo have been added quadratically, as suggested in \[Lucchesi, 2001\]. Regarding other subtle non–gravitational perturbations of thermal origin acting on the orbits of the LAGEOS satellites, their effects on the perigee of LAGEOS II are currently under accurate evaluation.

Then, a reliable estimate of the total systematic error induced by various classical perturbations on the measurement of $`\alpha `$ through eq.($`13`$) over a 5 years time span yields

$$\left(\frac{\delta \alpha }{\alpha }\right)^{\mathrm{total}}1\times 10^{11}.$$
($`20`$)
In obtaining eq.($`20`$) we have summed in a root–sum–square fashion the gravitational and non–gravitational errors assumed to be independent. Instead, the gravitational error due to the static and time–dependent parts of the Earth’s gravitational field have been simply summed up in view of a reciprocal correlation.

## 4 Conclusions

In this paper we have shown that by using a suitable combination of the orbital residuals of the perigee of LAGEOS II and the nodes of LAGEOS II and LAGEOS laser–ranged Earth satellites it would be possible to constrain effectively the Yukawa coupling $`\alpha `$ of a possible fifth force at a length scale of almost two Earth’s radii. The experimental sensitivity in measuring the perigee, which would be affected by the Yukawa force, of LAGEOS II would allow a precision on $`\alpha `$ of the order of almost $`4\times 10^{12}`$ over a time span of 1 year. According to the present knowledge of the terrestrial gravity field, the constraint posed by the systematic errors is of the order of $`1\times 10^{11}`$ over an observational time span of 5 years. These estimates should greatly improve when the new and more accurate data on the Earth gravity field from the CHAMP and GRACE missions will be available in the near future. The use of the proposed LARES satellite, with its larger eccentricity, would allow to improve the experimental constraint to $`1\times 10^{12}`$.


format: text

text:
# The runaway growth of intermediate-mass black holes in dense star clusters

## 1 Introduction

Using the Chandra X-ray observatory, Kaaret et al. (2000; 2001) and Matsumoto et al. (2000; 2001) recently discovered nine bright X-ray sources in the irregular galaxy M82. Their brightest source (No. 7 in Table 1 of Matsumoto et al. 2001) has a luminosity of $`9\times 10^{40}\mathrm{erg}\mathrm{s}^1`$ in the 0.2–10 KeV band, corresponding to the Eddington luminosity of a $`600`$$`\mathrm{M}_{}`$ compact object. The high luminosity and rather soft X-ray spectrum of the object indicates that it may be an intermediate-mass black hole (IMBH) with a mass of at least 600 $`\mathrm{M}_{}`$ (Kaaret et al. 2001; Matsumoto et al. 2001).

An optical follow-up in the infrared (J, H, and K-bands) with the CISCO instrument on the SUBARU telescope revealed a star cluster with an estimated mass of a few $`10^6`$$`\mathrm{M}_{}`$ at a position consistent with the X-ray location of the IMBH (Harashima et al. 2001). This star cluster appears to be very young ($`\stackrel{<}{}10`$ Myr), as it is extremely blue and expanding shells of molecular gas have been discovered in its vicinity (Matsushita et al. 2000), typical for a star-forming region of a few million years.

Matsushita et al. (2000) estimate that the environment has an age of only a few million years.

More unusually bright X-ray point sources have been discovered in the early spiral galaxies NGC 2403 (Kotoku et al. 2000) and NGC 4565 (Mizuno, et al. 1999). Most remarkable, however, is the discovery of many bright X-ray sources in the “Antennae” system (NGC 4038/4039) by Fabbiano et al. (2001), Zazas & Fabbiano (2002) and Zazas et al. (2002), also using Chandra. These authors conclude that many of these sources may be $`\stackrel{>}{}100`$$`\mathrm{M}_{}`$ accreting black holes (although alternative explanations exist—see e.g. Mizuno 1999; King et al. 2001). The Antennae contain many young star clusters with characteristics similar to those found in M82 (Mengel et al. 2001). However, it is not yet clear how many of the X-ray sources in the Antennae are associated with these clusters (Zazas & Fabbiano 2002). There may also be an example of an IMBH in our own Galaxy, as recent reverberation mapping of the globular cluster M15 by Gebhardt et al. (2000, 2001, and private communication) strongly suggests that the cluster may harbor a $`2500`$$`\mathrm{M}_{}`$ black hole at its center.

Several possible mechanisms for forming IMBHs in star clusters have recently been suggested. Miller & Hamilton (2001) have studied the possibility that an IMBH may form slowly (on a Hubble time scale) by occasionally encountering and devouring other cluster stars. Mouri & Taniguchi (2002) have proposed a much more rapid black-hole merger mechanism, operating in very high-density ($`10^6`$ black holes pc<sup>-3</sup>) environments on time scales as short as $`10^7`$ yr. In this paper we consider the possibility of forming a massive object in a young star cluster due to repeated collisions during an early phase of core collapse.

Sanders (1970), Lee (1987), and Quinlan & Shapiro (1990) have studied the possibility of collision runaways in spherical stellar systems of $`\stackrel{>}{}10^7`$ stars with high ($`>100`$$`\mathrm{km}\mathrm{s}^1`$) velocity dispersions. All studies began with stars of equal masses and found that, for sufficiently high densities and velocity dispersions, runaway mergers could indeed occur. Quinlan & Shapiro observed that the collision time scale for massive stars decreases faster with increasing mass than does the main-sequence lifetime, and concluded that clusters with initial relaxation times of 1–$`5\times 10^8`$ years could grow a massive $`\stackrel{>}{}100`$$`\mathrm{M}_{}`$ object by multiple mergers. Sanders’ (1970) Monte-Carlo calculations neglected the effects of mass segregation and found collision runaways only after mergers had driven the cluster into a state of high central density. However, in the self-consistent Fokker–Planck models of Lee and Quinlan & Shapiro, the runaway started well before core collapse occurred. All authors concluded that runaways would not occur in clusters containing less than $`10^610^7`$ stars because three-body binary heating in small $`N`$-systems provided sufficient energy to reverse core collapse before the runaway process could begin.

In contrast to the studies just described, the models discussed in this paper begin with a broad range of stellar masses. Vishniac (1978) demonstrated that a Salpeter (1955) initial mass function is Spitzer (1969) unstable. As a result, young star clusters may experience core collapse on the time scale on which the most massive stars segregate to the cluster center. This time scale may be much shorter than the main sequence lifetimes of the stars involved. Vishniac suggested that such a prompt collapse might lead to the formation of a massive compact object. We find that early core collapse in a relatively low-$`N`$ star cluster may result in a collision runaway, so long as the most massive stars remain on the main sequence while the collapse occurs.

The possibility of multiple collisions involving the same star in a dense star cluster was demonstrated convincingly by Portegies Zwart et al. (1999), using the special-purpose GRAPE-4 (Makino et al. 1997) to speed up their direct N-body calculations with up to 12288 stars. They concluded that, even in small clusters, runaway collisions may lead to the growth of a single massive star. The earlier arguments that three-body binary heating would drive the expansion of the cluster core appear to be unimportant in these simulations, as mergers between stars tend to destroy binaries before they can heat the cluster effectively. Indeed, in contrast to the underlying assumptions of previous collision studies, dynamically formed binaries in dense clusters act as a catalyst for stellar mergers, boosting the collision rate far beyond the simple two-body expressions used in earlier work. The N-body simulations covered a rather limited part of the available parameter space, but the initial conditions were selected to mimic known dense star clusters in the Galaxy and the Large Magellanic Cloud.

If the bright X-ray source in M82 does indeed correspond to a compact object of $`\stackrel{>}{}600`$$`\mathrm{M}_{}`$, this IMBH could have been formed by a collision runaway resulting from collapse of the cluster core early in the cluster. In fact, as we will see, it is quite natural to expect a $`10^3`$$`\mathrm{M}_{}`$ black hole in a million-solar-mass star cluster. We begin by deriving (§2) some simple analytic expressions describing the dynamical behavior observed in cluster simulations. In §3 we calibrate these relations using direct N-body simulations. We then (§4) extend these results to derive simple relations between black-hole formation and cluster parameters.

## 2 Runaway growth of a massive object in a dense star cluster

### 2.1 Core collapse and the first collision

A star cluster is a self-gravitating group of stars. So long as stellar evolution remains relatively unimportant, the cluster’s dynamical evolution is dominated by two-body relaxation, with characteristic time scale (Spitzer, 1987)

$$t_{\mathrm{rlx}}=\left(\frac{R_{\mathrm{c}}^{}{}_{}{}^{3}}{GM_\mathrm{c}}\right)^{1/2}\frac{N_\mathrm{c}}{8\mathrm{ln}\mathrm{\Lambda }_\mathrm{c}},$$
(1)
the half-mass relaxation time. Here $`G`$ is the gravitational constant, $`M_\mathrm{c}`$ is the total mass of the cluster, $`N_\mathrm{c}M_\mathrm{c}/m`$ is the number of stars and $`R_\mathrm{c}`$ is the characteristic (half-mass) radius of the cluster. The Coulomb logarithm $`\mathrm{ln}\mathrm{\Lambda }_\mathrm{c}\mathrm{ln}(0.1N_\mathrm{c})10`$ typically. In convenient units the two-body relaxation time becomes

$$t_{\mathrm{rlx}}1.9\mathrm{Myr}\left(\frac{R_\mathrm{c}}{1\mathrm{pc}}\right)^{3/2}\left(\frac{M_\mathrm{c}}{1\mathrm{M}_{}}\right)^{1/2}\left(\frac{1\mathrm{M}_{}}{m}\right)\left(\mathrm{ln}\mathrm{\Lambda }_\mathrm{c}\right)^1.$$
(2)
The dynamical evolution of the star cluster drives it toward core collapse (Antonov 1962; Spitzer & Hart, 1971) in which the central density runs away to a formally infinite value in a finite time. In an isolated cluster in which all stars have the same mass, core collapse occurs in a time $`t_{\mathrm{cc}}15t_{\mathrm{rlx}}`$ (Cohn 1980).

Realistic clusters have a broad range in initial stellar masses, generally from $`m_{\mathrm{min}}0.1`$$`\mathrm{M}_{}`$ to $`m_{\mathrm{max}}100`$$`\mathrm{M}_{}`$, with mean mass $`m`$ ranging from $`0.39`$$`\mathrm{M}_{}`$ (Salpeter 1955) to about 0.65 $`\mathrm{M}_{}`$ (Scalo 1986), depending on the specific mass function adopted. During the early evolution of the cluster, massive stars sink toward the cluster center via dynamical friction. Approximating the cluster structure as an isothermal sphere, we find (Binney & Tremaine 1987, Eq. 7-25) that a star of mass $`m`$ at distance $`r`$ from the cluster center drifts inward at a rate given by

$$r\frac{dr}{dt}=0.43\frac{Gm}{V_\mathrm{c}}\mathrm{ln}\mathrm{\Lambda }_\mathrm{c},$$
(3)
Here $`V_\mathrm{c}`$ is the cluster velocity dispersion. Using Eq. 2, we can integrate Eq. 3 with respect to time to obtain the dynamical friction inspiral time scale

$$t_\mathrm{f}=3.3\frac{m}{m}t_{\mathrm{rlx}}.$$
(4)
This is the time taken for a star of mass $`m`$ to sink to the cluster center from a circular orbit at initial distance $`rr_{\mathrm{core}}`$.

In a multi-mass system, core collapse is driven by the accumulation of the most massive stars in the cluster center. This process takes place on a dynamical friction time scale (Eq. 4). Empirically, we find, for initial mass functions of interest here, that core collapse (actually, the appearance of the first persistent dynamically formed binary systems) occurs at about

$$t_{\mathrm{cc}}0.20t_{\mathrm{rlx}}.$$
(5)
This core collapse time is taken in the limit where stellar evolution is unimportant, i.e. where stellar mass loss is negligible and the most massive stars survive until they reach the cluster center.

The collapse of the cluster core may initiate physical collisions between stars. The product of the first collision is likely to be among the most massive stars in the system, and to be in the core. This star is therefore likely to experience subsequent collisions, resulting in a collision runaway (see Portegies Zwart et al. 1999). The maximum mass that can be grown in a dense star cluster if all collisions involve the same star is $`m_\mathrm{r}`$, where

$$\frac{dm_\mathrm{r}}{dt}=𝒩_{\mathrm{coll}}\delta m_{\mathrm{coll}}.$$
(6)
Here $`𝒩_{\mathrm{coll}}`$ and $`\delta m_{\mathrm{coll}}`$ are the average collision rate and the average mass increase per collision (assumed independent). We now discuss these quantities in more detail. In §3 we present some N-body results that both motivated and calibrate the following discussion.

### 2.2 The collision rate $`𝒩_{\mathrm{coll}}`$

A key result from our simulations is the fact that collisions between stars generally occur in dynamically formed (“three-body”) binaries. The collision rate is therefore closely related to the binary formation rate, which we now estimate.

The flux of energy through the half-mass radius of a cluster during one half-mass relaxation time is on the order of 10% of the cluster potential energy, largely independent of the total number of stars or the details of the cluster’s internal structure (Goodman 1987). For a system without primordial binaries this flux is produced by heating due to dynamically formed binaries (Makino & Hut 1990). It is released partly in the form of scattering products which remain bound to the system, and partly in the form of potential energy removed from the system by escapers recoiling out of the cluster (Hut & Inagaki 1985). Makino & Hut argue, for an equal-mass system, that a binary generates an amount of energy on the order of $`10^2kT`$ via binary–single-star scattering (where the total kinetic energy of the stellar system is $`\frac{3}{2}N_\mathrm{c}kT`$). This quantity originates from the minimum binding energy of a binary that can eject itself following a strong encounter. Assuming that the large-scale energy flux in the cluster is ultimately powered by binary heating in the core. It follows that the required formation rate of binaries via three-body encounters is

$$n_{\mathrm{bf}}10^3\frac{N_\mathrm{c}}{t_{\mathrm{rlx}}}.$$
(7)
For systems containing significant numbers of primordial binaries, which segregate to the cluster core, equivalent energetic arguments (Goodman & Hut 1989) lead to a similar scaling for the net rate at which binary encounters occur in the core.

The above arguments apply to star clusters comprising identical point-mass stars. In a cluster with a range of stellar masses, three-body binaries generally form from stars which are more massive than average. After repeated exchange interactions, the binary will consist of two of the most massive stars in the cluster. Conservation of linear momentum during encounters with lower mass stars means that the binary receives a smaller recoil velocity, making it less likely to be ejected from the cluster. The binary must therefore be considerably harder—$`\stackrel{>}{}10^3kT`$—before it is ejected following a encounter with another star (see Portegies Zwart & McMillan 2000).

However, taking the finite sizes of real stars into account, it is quite likely that such a hard binary experiences a collision rather than being ejected. A strong encounter between a single star and a hard binary generally results in a resonant interaction. Three stars remain in resonance until at least one of them escapes, or a collision reduces the three-body system to a stable binary. For harder binaries it becomes increasingly likely that a collision occurs instead of ejection (McMillan 1986). In the calculations of Portegies Zwart et al. (1999) most binaries experience a collision at a binding energy of order $`10^2kT`$, considerably smaller than the binding energy required for ejection. Accordingly, we retain the above estimate of the binary formation rate (Eq. 7) and conclude that the collision rate per half-mass relaxation time is

$$𝒩_{\mathrm{coll}}10^3f_\mathrm{c}\frac{N_\mathrm{c}}{t_{\mathrm{rlx}}}.$$
(8)
Here we introduce $`f_\mathrm{c}1`$, the effective fraction of dynamically formed binaries that produce a collision. Note again that Eq. 8 is valid only in the limit where stellar evolution is unimportant.

The most massive star in the cluster is typically a member of the interacting binary and therefore dominates the collision rate. Subsequent collisions cause the runaway to grow in mass, making it progressively less likely to escape from the cluster. The star which experiences the first collision is therefore likely to participate in subsequent collisions. The majority of collisions thus involve one particular object—the runaway merger—generally selected by its high initial mass and proximity to the cluster center (see Portegies Zwart et al. 1999).

For systems containing many primordial binaries the above argument must be modified. Since dynamically formed binaries tend to be fairly soft—a few $`kT`$—we expect that the fraction of interactions with primordial binaries leading to collision is comparable to the value $`f_\mathrm{c}`$ above. However, a critical difference is that, in systems containing many binaries, the collisions involve many different pairs of stars, not just the binary containing the massive runaway. For our proposed runaway scenario to operate, we must assume that high-mass binaries are rapidly destroyed or merge following interactions in the core, in which case the above arguments apply. We note that, once a runaway begins, binaries have large interaction cross sections and hence are likely to participate in the runaway process. From the point of view of producing massive merger products, the worst-case scenario would be a substantial primordial population of wide high-mass binaries. We are currently carrying out N-body simulations to investigate the behavior of such systems.

### 2.3 Average mass increase per collision

Once begun, the collision runaway dominates the collision cross section. The average mass increase per collision depends on the characteristics of the mass function in the cluster core. A lower limit for stars which participate in collisions can be derived from the degree of segregation in the cluster. Inverting Eq. 4 results in an estimate (still assuming an isothermal sphere) of the minimum mass of a star that can reach the cluster core in time $`t`$ due to dynamical friction:

$$m_\mathrm{f}=1.9\mathrm{M}_{}\left(\frac{1\mathrm{Myr}}{t}\right)\left(\frac{R_\mathrm{c}}{1\mathrm{pc}}\right)^{3/2}\left(\frac{M_\mathrm{c}}{1\mathrm{M}_{}}\right)^{1/2}\left(\mathrm{ln}\mathrm{\Lambda }_\mathrm{c}\right)^1.$$
(9)
Thus, at time $`t`$ and for a given mass $`m`$, there is a maximum radius $`r(t)`$ inside of which stars of mass $`m`$ will have segregated to the core. The stars contributing to the growth of the runaway are likely to be among those more massive than $`m_\mathrm{f}`$, because their number density in the core is enhanced by mass segregation, their collision cross sections are larger, and they contribute more to $`\delta m_{\mathrm{coll}}`$ when they do collide.

The shape of the central mass function of a segregated cluster is not trivial to derive. In thermal equilibrium, the central number densities of stars of different masses would be expected to scale as

$$n_0(m)m^{3/2}\frac{dN}{dm},$$
(10)
where $`dN/dm`$ is the global (Scalo) initial mass function, which scales roughly as $`m^{2.7}`$ at the high-mass end ($`m\stackrel{>}{}10\mathrm{M}_{}`$). However, as discussed in §3, the distribution of secondary masses (i.e. the masses of the lighter stars participating in collisions) does not follow the above simple relation. Rather, we find that stars in the core do not reach thermal equilibrium (a result generally consistent with earlier findings by Chernoff and Weinberg 1990 and Joshi, Nave & Rasio 2001), and that the dynamical nature of the collisional processes involved mean that more massive stars tend to be consumed before lower-mass stars arrive in the core. In addition, most collisions involve three-body binary formation and binary interactions in a multi-mass environment, further complicating the connection between stellar densities and secondary masses.

Empirically, we find that the secondary mass distribution is quite well fit by a power-law, $`dN/dmm^{2.3}`$ (coincidentally very close to a Salpeter distribution). Integrating this expression from a minimum mass of $`m_\mathrm{f}`$ (and ignoring the upper limit) results in a mean mass increase per collision of

$$\delta m_{\mathrm{coll}}4m_\mathrm{f}.$$
(11)
We neglect stars with masses less than $`m_\mathrm{f}`$. Substitution of Eq. 1 into Eq. 9 and Eq. 11 then results in a mass increase per collision of

$$\delta m_{\mathrm{coll}}4\frac{t_{\mathrm{rlx}}}{t}m\mathrm{ln}\mathrm{\Lambda }_\mathrm{c}.$$
(12)
Taken over the entire “collisional” lifetime of the core, it is perhaps not surprising that the net distribution of secondary masses tends to follow the overall distribution of high-mass stars.

### 2.4 Lifetime of a cluster in a static tidal field

With simple expressions for $`𝒩_{\mathrm{coll}}`$ and $`\delta m_{\mathrm{coll}}`$ now in hand, we return to the determination of the runaway growth rate (Eq. 6). The evaporation of a star cluster which fills its Jacobi surface in an external potential is driven by tidal stripping. Portegies Zwart et al (2001a) have studied the evolution of young compact star clusters within $`200`$ pc of the Galactic center. Their calculations employed direct N-body integration, including the effects of both stellar and binary evolution and the (static) external influence of the Galaxy, and made extensive use of the GRAPE-4 (Makino et al. 1997) special-purpose computer. They found that the mass of a typical model cluster decreased almost linearly with time:

$$M_\mathrm{c}=M_{\mathrm{c0}}\left(1\frac{t}{t_{\mathrm{disr}}}\right).$$
(13)
Here $`M_{\mathrm{c0}}`$ is the mass of the cluster at birth and $`t_{\mathrm{disr}}`$ is the cluster’s disruption time. Portegies Zwart et al. (2001a) found that their model clusters dissolved within about 30% of the two-body relaxation time at the tidal radius (defined by substituting the tidal radius instead of the virial radius in Eq. 1). In terms of the half-mass relaxation time, we find $`t_{\mathrm{disr}}=1.6`$–5.4 $`t_{\mathrm{rlx}}`$, depending on the initial density profile (the range corresponds to King  dimensionless depths $`\mathrm{W}_0=3`$–7; more centrally condensed clusters live longer).

Substituting Eqs. 8 and 12 into Eq. 6, and defining $`M_{\mathrm{c0}}=N_\mathrm{c}m`$ to rewrite Eq. 13 in terms of the number of stars in the cluster, we find

$`{\displaystyle \frac{dm_\mathrm{r}}{dt}}`$ $`=`$ $`4\times 10^3f_\mathrm{c}{\displaystyle \frac{N_\mathrm{c}m\mathrm{ln}\mathrm{\Lambda }_\mathrm{c}}{t}}`$ (14)
$`=`$ $`4\times 10^3f_\mathrm{c}M_{\mathrm{c0}}\mathrm{ln}\mathrm{\Lambda }_\mathrm{c}\left({\displaystyle \frac{1}{t}}{\displaystyle \frac{1}{t_{\mathrm{disr}}}}\right).`$
Integrating from $`t=t_{\mathrm{cc}}`$ to $`t=t_{\mathrm{disr}}`$ results in

$$m_\mathrm{r}=m_{\mathrm{seed}}+4\times 10^3f_\mathrm{c}M_{\mathrm{c0}}\mathrm{ln}\mathrm{\Lambda }_\mathrm{c}\left[\mathrm{ln}\left(\frac{t_{\mathrm{disr}}}{t_{\mathrm{cc}}}\right)+\frac{t_{\mathrm{cc}}}{t_{\mathrm{disr}}}1\right].$$
(15)
Here $`m_{\mathrm{seed}}`$ is the seed mass of the star which initiates the runaway growth, most likely one of the most massive stars initially in the cluster. With $`t_{\mathrm{cc}}0.2t_{\mathrm{rlx}}`$, Eq 15 reduces to

$$m_\mathrm{r}=m_{\mathrm{seed}}+4\times 10^3f_\mathrm{c}M_{\mathrm{c0}}\gamma \mathrm{ln}\mathrm{\Lambda }_\mathrm{c},$$
(16)
where $`\gamma \mathrm{ln}t_{\mathrm{disr}}/t_{\mathrm{cc}}+t_{\mathrm{cc}}/t_{\mathrm{disr}}11`$.

## 3 Results of N-body simulations

The development of the GRAPE family of special-purpose computers makes it relatively straightforward to test and tune the above simple model using direct N-body calculations. Tab. 1 summarizes the results of an extensive series of detailed N-body simulations of core collapse and stellar collisions in dense star clusters containing up to 65536 stars. These simulations were performed using the “Starlab” software environment (see Portegies Zwart et al. 2001b) running on the GRAPE-6 (Makino 2000). The calculations were performed with initially single stars but, as just discussed, the presence of primordial binaries is not likely to change the picture qualitatively. To expand on these findings, we have performed an additional series of simulations with $`10^4`$ stars using the same software and hardware. Further simulations of systems containing substantial numbers of primordial binaries are in progress, but are much more time consuming, due to the complexity of following binary and multiple encounters in a large-$`N`$ context.

### 3.1 Core collapse

In our isolated star clusters (three calculations) with $`10^4`$ identical single point masses distributed as a Plummer model, core collapse occurs at $`t_{\mathrm{cc}}15.2\pm 0.1t_{\mathrm{rlx}}`$. This result is consistent with earlier calculations of e.g., Cohn (1980) and Makino (1996). Doubling the mass of 20% of the stars reduced the core collapse time to $`t_{\mathrm{cc}}7.2t_{\mathrm{rlx}}`$. Making 20% of the stars 10 or 100 times more massive reduced the time of core collapse further, to $`t_{\mathrm{cc}}1.4t_{\mathrm{rlx}}`$ and $`t_{\mathrm{cc}}0.16t_{\mathrm{rlx}}`$, respectively.

The more realistic models of Portegies Zwart et al. (1999) with 6144 and 12288 single stars drawn from a Scalo (1986) initial mass function also include mass loss from stellar evolution. The initial density distributions for these models were $`\mathrm{W}_0=6`$ King (1966) models. Core collapse in these models occurred at $`t_{\mathrm{cc}}0.19\pm 0.08t_{\mathrm{rlx}}`$. The slightly later collapse compared to the models just described, containing $`10^4`$ identical point masses and a heavy component, may be attributable to the rather different mass function, as well as to stellar mass loss, which tends to delay core collapse.

### 3.2 Collision rate

Relaxing the assumption of point masses to include finite stellar sizes introduces collisions into our models. In all calculations, the first collision occurred shortly after the formation of the first $`\stackrel{>}{}10kT`$ binary by a three-body encounter, i.e. close to the time of core collapse. When stars were given unrealistically large radii (100 times larger than normal), the first collisions occurred only slightly (about 5%) earlier.

As discussed earlier, the first star to experience a collision was generally one of the most massive stars in the cluster; this star then became the target for further collisions. In models with initial relaxation times greater than about 30 Myr the target star exploded in a supernova before experiencing runaway growth. The collision rates in these clusters were considerably smaller than for clusters with smaller relaxation times (see Fig. 1). As discussed in more detail in §4, the onset of stellar evolution terminates the collision process; premature disruption of the cluster also ends the period of runaway growth.

The 45 $`N`$-body calculations listed in Tab. 1 span a broad range of initial conditions. The number of stars varied from 1k (1024) to 64k (65536). Initial density profiles and velocity dispersion for the models were taken from Heggie-Ramamani models (Heggie & Ramamani 1995) with $`\mathrm{W}_0`$ ranging from 1 to 7. At birth, the clusters were assumed to fill their zero-velocity (Jacobi) surfaces in the Galactic tidal field. In most cases we adopted an initial mass function between 0.1 $`\mathrm{M}_{}`$ and 100 $`\mathrm{M}_{}`$ suggested for the Solar neighborhood by Scalo (1986). However, several calculations were performed using power-law initial mass functions with exponents of -2 or -2.35 (Salpeter) and lower mass limits of 1 $`\mathrm{M}_{}`$. The model with 64k stars (model N64R6r36) was initialized with a Scalo (1986) mass function, but with a lower mass limit of 0.3 $`\mathrm{M}_{}`$ instead of the 0.1 $`\mathrm{M}_{}`$ used in the other models. The characterization of the tidal field is discussed in Portegies Zwart et al. (2001b).

The number of collisions in these simulations ranged from 0 to 24. Fig. 1 shows the mean collision rate $`𝒩_{\mathrm{coll}}`$ per star per million years as a function of the initial half-mass relaxation time. The solid line in Fig. 1 is a fit to the simulation data, and has

$$𝒩_{\mathrm{coll}}=2.2\times 10^4\frac{N_\mathrm{c}}{t_{\mathrm{rlx}}},$$
(17)
for $`t_{\mathrm{rlx}}\stackrel{<}{}2030`$ Myr, consistent with our earlier estimate (Eq. 8) if $`f_\mathrm{c}=0.2`$. The quality of the fit in Fig. 1 is quite striking, especially when one bears in mind the rather large spread in initial conditions for the various models.

Figure 2 shows the cumulative mass distributions of the primary (more massive) and secondary (less massive) stars participating in collisions. We include only events in which the secondary experienced its first collision (that is, we omit secondaries which were themselves collision products). In addition, we distinguish between collisions early in the evolution of the cluster and those that happened later by subdividing our data based on the ratio $`\tau =t_{\mathrm{coll}}/t_\mathrm{f}`$, where $`t_{\mathrm{coll}}`$ is the time at which a collision occurred and $`t_\mathrm{f}`$ is the dynamical friction time scale of the secondary star (see Eq. 4). The solid lines in Figure 2 show cuts in the secondary masses at $`\tau \stackrel{<}{}1`$, $`\tau \stackrel{<}{}5`$ and $`\tau <\mathrm{}`$ (rightmost line). The mean secondary masses are $`m=4.0\pm 4.8`$$`\mathrm{M}_{}`$, $`8.2\pm 6.5`$ and $`m=13.5\pm 8.8`$$`\mathrm{M}_{}`$ for $`\tau \stackrel{<}{}1`$, 5 and $`\mathrm{}`$, respectively.

The distribution of primary masses in Figure 2 (dashed line) hardly changes as we vary the selection on $`\tau `$. We therefore show only the full ($`\tau \stackrel{<}{}\mathrm{}`$) data set for the primaries. In conrast, the distribution of secondary masses changes considerably with increasing $`\tau `$. For small $`\tau `$, secondaries are drawn primarily from low-mass stars. As $`\tau `$ increases, the secondary distribution shifts to higher masses while the low-mass part of the distribution remains largely unchanged. The shift from low-mass ($`\stackrel{<}{}8`$$`\mathrm{M}_{}`$) to high-mass collision secondaries ($`\stackrel{>}{}8`$$`\mathrm{M}_{}`$) occurs between $`\tau =1`$ and $`\tau =5`$. This is consistent with the theoretical arguments presented in Sec. 2.3. During the early evolution of the cluster ($`\tau \stackrel{<}{}1`$), collision partners are selected more or less randomly from the available (initial) population in the cluster core; at later times, most secondaries are drawn from the mass-segregated population.

Interestingly, although hard to see in Fig. 2, all the curves are well fit by power laws between $`8`$$`\mathrm{M}_{}`$ and $`80`$$`\mathrm{M}_{}`$(0.8 $`\mathrm{M}_{}`$ and 30 $`\mathrm{M}_{}`$for the leftmost curve). The power-law exponents are $`0.4`$, $`0.5`$ and $`2.3`$ for $`\tau \stackrel{<}{}1`$, 5, and $`\mathrm{}`$, and $`0.3`$ for the primary (dashed) curve. (Note that the Salpeter mass function has exponent $`2.35`$.)

Figure 3 shows the maximum mass of the runaway collision product as function of the initial mass of the star cluster. Only the left side ($`\mathrm{log}M/\mathrm{M}_{}\stackrel{<}{}7`$) of the figure is relevant here; we discuss the extrapolation to larger masses in Sec. 4.3. The N-body results are consistent with the theoretical model presented in Eq. 15.

## 4 Discussion

Early core collapse in dense star clusters may initiate a phase of runaway stellar growth, leading to an object containing up to $`0.1`$ % of the total cluster mass. We do not address here the state of this object, which could be a black hole or a star. If the object is a helium- or hydrogen-burning star, it may collapse into a compact object when it exhausts its central fuel. The amount of mass lost in the supernova explosion and whether the compact object receives a velocity kick are important considerations for the future evolution of the collision runaway. An extensive parameter study of the details of the supernova is beyond the scope of this paper. The basis of our analysis, however, is simple and robust to quite substantial perturbations.

We now consider the circumstances under which runaway growth may be prevented or terminated at an early stage. Premature termination of the runaway occurs when stellar mass loss starts to drive the expansion of the star cluster, or when the star cluster is disrupted by external influences. At the end of this section (Sec. 4.3) we briefly discuss the application of runaway growth to the possible formation of supermassive black holes in the bulges of galaxies.

### 4.1 Prevention of the collision runaway

Runaway growth in a star cluster can only occur when stellar evolution is relatively unimportant compared to the dynamical evolution of the cluster. Stellar mass loss tends to heat the cluster by loss of potential energy, and can easily reverse core collapse. This is particularly true for the most massive stars, which dominate the dynamics of the cluster core and are also the first to lose substantial amounts of mass in stellar winds and supernovae. The prevention of core collapse also prevents the first or any subsequent collisions, and a reversal of core collapse terminates the collision runaway. As a rule of thumb, we argue that runaway growth can be prevented when the time scale for the most massive stars to segregate to the cluster center exceeds the lifetimes of those stars.

The main-sequence lifetime ($`t_{\mathrm{ms}}`$) for stars more massive than $`25`$$`\mathrm{M}_{}`$ is a rather flat function of mass. For high metallicities, the most massive stars ($`m\stackrel{>}{}85`$$`\mathrm{M}_{}`$) actually live longer than stars with masses between 25 and $`65\mathrm{M}_{}`$ (Meynet et al. 1994). For $`Z=0.04`$ (recall that the solar metallicity is $`\mathrm{Z}_{}=0.02`$), the hydrogen plus helium burning lifetime varies from 6.27 Myr for a $`25\mathrm{M}_{}`$ star to 7.57 Myr for a $`120\mathrm{M}_{}`$ star. For low metallicities ($`Z=0.001`$), the range becomes 8.19–3.20 Myr for the same masses (see Meynet et al. for details).

For the star cluster to experience core collapse before the most massive stars evolve, we require $`t_{\mathrm{cc}}0.2t_{\mathrm{rlx}}\stackrel{<}{}t_{\mathrm{ms}}(120\mathrm{M}_{})3.20`$$`7.57`$ Myr for $`Z=0.001`$–0.04. Runaway growth therefore does not occur in star clusters with initial relaxation times $`t_{\mathrm{rlx}}\stackrel{>}{}16`$–38 Myr. For definiteness, we conclude that clusters with $`t_{\mathrm{rlx}}\stackrel{<}{}25`$ Myr experience core collapse before the most massive stars explode, and are therefore prone to runaway collision.

The half-mass radius of a tidally limited cluster expands during core collapse, causing the mean relaxation time to increase by about a factor of 4 (see Portegies Zwart et al. 2001). A cluster with an initial relaxation time of $`t_{\mathrm{rlx}}25`$ Myr will therefore have a relaxation time of about 100 Myr after core collapse. Such a cluster will not experience any further collision runaway, but may still contain the evidence of such a phase in the form of a central compact object with a mass $`\stackrel{<}{}0.1`$% of the initial cluster mass. The cluster may also be relatively depleted in low-mass compact objects (stellar mass black holes and neutron stars), as these are consumed during the runaway growth phase (Portegies Zwart et al. 1999).

### 4.2 Early termination of the runaway by tidal disruption

A star cluster in orbit around the Galactic center is subject to dynamical friction, in much the same way as dynamical friction drives massive stars toward the cluster center. This causes the cluster to spiral into the Galactic center, where it is destroyed (see Gerhard 2001). We derive here in some detail the dynamical friction time scale for a star cluster in the potential of the Galactic center. We assume constant cluster mass $`M_c`$, deferring the more realistic case of a time-dependent cluster mass (cf. Eq.13) to McMillan & Portegies Zwart (2002, in preparation).

The drag acceleration due to dynamical friction is (equation \[7-18\] in Binney & Tremaine, 1987)

$$a=\frac{4\pi \mathrm{ln}\mathrm{\Lambda }_\mathrm{G}G^2M_\mathrm{c}\rho _G(R_\mathrm{G})}{v_{\mathrm{o}}^{}{}_{}{}^{2}}\left[\mathrm{erf}(X)\frac{2X}{\sqrt{\pi }}e^{X^2}\right].$$
(18)
Here $`\mathrm{ln}\mathrm{\Lambda }_\mathrm{G}`$ is the Coulomb logarithm for the Galactic central region, for which we adopt $`\mathrm{ln}\mathrm{\Lambda }_\mathrm{G}R_\mathrm{G}/R_\mathrm{c}`$, $`\mathrm{erf}`$ is the error function and $`Xv_\mathrm{o}/\sqrt{2}V_\mathrm{G}`$, where $`V_\mathrm{G}`$ is the one-dimensional velocity dispersion of the stars at distance $`R_\mathrm{G}`$ from the Galactic center.

The mass of the Galaxy lying within the cluster’s orbit at distance $`R_\mathrm{G}`$ ($`\stackrel{<}{}500`$ pc) from the Galactic center is (Sanders & Lowinger 1972; Mezger et al. 1996)

$$M_\mathrm{G}(R_\mathrm{G})=4.25\times 10^6\left(\frac{R_\mathrm{G}}{1\mathrm{pc}}\right)^{1.2}\mathrm{M}_{}.$$
(19)
Its derivative, the local Galactic density (see Portegies Zwart et al. 2001) is

$$\rho _G(R)4.06\times 10^5\left(\frac{R_\mathrm{G}}{1\mathrm{pc}}\right)^{1.8}\mathrm{M}_{}\mathrm{pc}^3.$$
(20)
For inspiral through a sequence of nearly circular orbits, the function $`\mathrm{erf}(X)\frac{2X}{\sqrt{\pi }}\mathrm{exp}(X^2)`$ appearing in Eq. 18 may be determined as follows.

Following Binney & Tremaine (p. 226), we write the equation of dynamical equilibrium for stars near the Galactic center as

$$\frac{dP}{dR_\mathrm{G}}=\rho _G\frac{GM_\mathrm{G}(R_\mathrm{G})}{R_\mathrm{G}^2},$$
(21)
where $`P=kT\rho /m`$, $`\frac{3}{2}kT=\frac{1}{2}mv^2`$. Since $`\sigma ^2=\frac{1}{3}v^2`$, it follows that $`P=\sigma ^2\rho `$, and Eq. 21 becomes

$$\frac{d}{dr}(\sigma ^2\rho )=\frac{\rho }{r}v_\mathrm{o}^2,$$
(22)
where $`v_\mathrm{o}`$ is the circular orbital velocity at radius $`R`$: $`V_{\mathrm{G}}^{}{}_{}{}^{2}=GM_\mathrm{G}(R_\mathrm{G})/R_\mathrm{G}`$. For $`M_\mathrm{G}R_{\mathrm{G}}^{}{}_{}{}^{x}`$ (see Eq. 19), and assuming that $`V_{\mathrm{G}}^{}{}_{}{}^{2}v_\mathrm{o}^2R_{\mathrm{G}}^{}{}_{}{}^{x1}`$, we find $`V_{\mathrm{G}}^{}{}_{}{}^{2}\rho R_{\mathrm{G}}^{}{}_{}{}^{2x4}`$, so

$$r\frac{d}{dr}(V_{\mathrm{G}}^{}{}_{}{}^{2}\rho )=(2x4)V_{\mathrm{G}}^{}{}_{}{}^{2}\rho =\rho v_\mathrm{o}^2,$$
(23)
and hence $`X=\sqrt{2x}`$. Eq. 18 then becomes

$$a=1.2\mathrm{ln}\mathrm{\Lambda }_\mathrm{G}\frac{GM_\mathrm{c}}{R_\mathrm{G}^2}\left[\mathrm{erf}(X)\frac{2X}{\sqrt{\pi }}\mathrm{exp}(X^2)\right].$$
(24)
For $`x=1.2`$, $`X=0.89`$ and

$$a=0.41\mathrm{ln}\mathrm{\Lambda }_\mathrm{G}\frac{GM_\mathrm{c}}{R_\mathrm{G}^2}.$$
(25)
Again following Binney & Tremaine, defining $`L=R_\mathrm{G}v_\mathrm{o}`$ and setting $`dL/dt=aR_\mathrm{G}`$, we can integrate Eq. 25 with respect to time to find an inspiral time from initial radius $`R_i`$ of

$`T_\mathrm{f}`$ $``$ $`{\displaystyle \frac{1.28}{\mathrm{ln}\mathrm{\Lambda }_\mathrm{G}}}{\displaystyle \frac{M_\mathrm{G}(R_i)}{M_\mathrm{c}}}\left[{\displaystyle \frac{GM_\mathrm{G}(R_i)}{R_i^3}}\right]^{1/2}`$ (26)
$``$ $`1.4\left({\displaystyle \frac{R_i}{10\mathrm{p}\mathrm{c}}}\right)^{2.1}\left({\displaystyle \frac{10^6\mathrm{M}_{}}{M_\mathrm{c}}}\right)\mathrm{Myr}`$ (27)
For definiteness, we have assumed $`\mathrm{ln}\mathrm{\Lambda }_\mathrm{G}4`$ ($`\mathrm{\Lambda }_\mathrm{G}R_\mathrm{G}/R_\mathrm{c}100`$) in Eq. 27, corresponding to a distance of about 10–30 pc from the Galactic center.

The maximum mass of the runaway merger for clusters which are disrupted by inspiral (which of course always destroys the cluster before it reaches the center) may be calculated by replacing $`t_{\mathrm{disr}}`$ in Eq. 15 by $`T_\mathrm{f}`$. The right-hand side of that equation then becomes a function of

$$\frac{T_\mathrm{f}}{t_{\mathrm{cc}}}9.0\left(\frac{R_i}{10\mathrm{pc}}\right)^{2.1}\left(\frac{0.25\mathrm{pc}}{R_\mathrm{c}}\right)^{3/2}\left(\frac{10^5\mathrm{M}_{}}{M_\mathrm{c}}\right)^{3/2}$$
(28)
We can also estimate the maximum initial distance from the Galactic center for which core collapse occurs (and hence runaway merging may begin) before the cluster disrupts by setting $`T_\mathrm{f}=t_{\mathrm{cc}}`$. The result is $`R_i\stackrel{>}{}0.0025\mathrm{pc}\left(R_\mathrm{c}M_\mathrm{c}/[\mathrm{pc}\mathrm{M}_{}]\right)^{0.71}`$. For $`R_c=0.25`$ pc and $`M_c=10^5\mathrm{M}_{}`$, we find $`R_i\stackrel{>}{}3.3`$ pc.

### 4.3 Speculation on the formation of supermassive black holes

A million-solar-mass star cluster formed at a distance of $`\stackrel{<}{}30`$ pc from the Galactic center can spiral into the Galactic center by dynamical friction before being disrupted by the tidal field of the Galaxy (see Gerhard 2001). Only the densest star clusters survive to reach the center. These clusters are prone to runaway growth and produce massive compact objects at their centers. Upon arrival at the Galactic center, the star cluster dissolves, depositing its central black hole there. Black holes from in-spiraling star clusters may subsequently merge to form a supermassive black hole. Ebisuzaki et al. (2001) have proposed that such a scenario might explain the presence of the central black hole in the Milky Way galaxy.

If we simply assume that bulges and central supermassive black holes are formed from disrupted star clusters, this model predicts a relation between black hole and bulge masses in galaxies similar to the expression (Eq. 16) connecting the mass of an IMBH to that of its parent cluster. However, the ratio of stellar mass to black-hole mass might be expected to be smaller for galactic bulges than for star clusters, because not all star clusters produce a black hole and not all star clusters survive until the maximum black hole mass is reached. We would expect, however, that the general relation between the black hole mass and that of the bulge remains valid.

Figure 3 shows the relation between the black hole mass and the bulge mass for several Seyfert galaxies and quasars. The expression derived in Sec. 2 and the results of our N-body calculations (Sec. 3) are also indicated. The solid and dashed lines (Eq. 15) fit the N-body calculations and enclose the area of the measured black hole mass–bulge masses. On the way, the solid curve passes though two other black-hole mass estimates, for M82 and the globular cluster M15. We note that the observed relation between bulge and black hole masses has a spread of two orders of magnitude. If this bold extrapolation really does reflect the formation process of bulges and central black holes, this spread could be interpreted as a variation in the efficiency of the runaway merger process. In that case, only about one in a hundred star clusters reaches the galactic center, where its black hole is deposited.

### 4.4 Is the globular cluster M15 a special case?

The possible black hole in the globular cluster M15 may have been formed by a scenario different from the one described in this paper, as the cluster’s initial relaxation time probably exceeded our upper limit of 25 Myr. The current half-mass relaxation time of M15 is about 2.5 Gyr (Harris 1996), which is far more than our 100 Myr limit for forming a massive central object from a collision runaway.

An alternative is provided by Miller & Hamilton (2001), who describe the formation of massive ($`10^3`$$`\mathrm{M}_{}`$) black holes in star clusters with relatively long relaxation times. In their model the black hole grows very slowly over a Hubble time via occasional collisions with other stars, in contrast to the model described here, in which the runaway grows much more rapidly, reaching a characteristic mass of about 0.1% of the total birth mass of the cluster within a few megayears.

One possible way around M15’s long relaxation time may involve the cluster’s rotation. Gebhardt (2000; 2001; private communication) has measured radial velocities of individual stars in the crowded central field, down to two arcsec of the cluster center. He finds that, both in the central part of the cluster ($`r<0.1r_{\mathrm{hm}}`$) and outside the half mass radius, the average rotation velocity is substantial ($`v_{\mathrm{rot}}\stackrel{>}{}0.5v^2^{1/2}`$). Rotation is quickly lost in a cluster, so to explain a current rotation, M15’s initial rotation rate must probably have been even larger than observed today (see Einsel & Spurzem 1999). Hachisu (1978; 1982) found, using gaseous cluster models, that an initially rotating cluster tends to evolve into a ’gravo-gyro catastrophe’ which drives the cluster into core collapse far more rapidly than would occur in a non-rotating system. If the gravo-gyro-driven core collapse occurred within 25 Myr, a collision runaway might have initiated the growth of an intermediate mass black hole in the core of M15.

## 5 Conclusions

We study the runaway growth of a single star in a dense star cluster using a combination of complementary approaches. Our semi-analytic analysis is supported by detailed N-body calculations in which the effects of stellar evolution, stellar dynamics, binary evolution and the perturbing effect of a background Galactic potential are taken self-consistently into account.

Star clusters with initial half-mass relaxation times $`t_{\mathrm{rlx}}\stackrel{<}{}25`$ Myr experience a phase of runaway growth. In this phase a single seed star grows to a mass of about 0.1% of the total mass of the cluster. The first collision occurs at the moment the cluster core collapses. This happens at about 0.2 $`t_{\mathrm{rlx}}`$ but no later than about 5 Myr (the evolution time scale for a $`\stackrel{>}{}50\mathrm{M}_{}`$ star). The star which experiences the first collision becomes the target for further collisions, initiating runaway growth. The growth phase is terminated by (1) the disruption of the cluster in the tidal field of the Galaxy (at $`t\stackrel{<}{}5t_{\mathrm{rlx}}`$) or (2) the reversal of core collapse by mass loss from the evolving stellar population (after about 25 Myr).

A star cluster can survive for longer than $`5t_{\mathrm{rlx}}`$ if, for example, it did not initially fill its Jacobi surface (“Roche lobe”) in the Galactic tidal field. (Examples are NGC 3603 and R 136, the dense star cluster in the 30 Doradus region in the Large Magellanic cloud.) Such clusters go though a phase of runaway stellar growth, but recover after stellar mass loss drives the re-expansion of the cluster core.

¿From an observational point of view, a tidally limited cluster experiences three very distinct evolutionary phases: a pre-collapse phase until 0.2 $`t_{\mathrm{rlx}}`$, a phase of deep core collapse (from 0.2 $`t_{\mathrm{rlx}}`$ to about 25 Myr), followed by an expansion phase eventually leading to the disruption of the cluster. During the expansion phase the cluster half-mass radius expands, causing the mean relaxation time to increase by a factor of 4 (see Portegies Zwart et al. 2001). A cluster in this final phase will be observable with a current relaxation time less than $`4\times 25\mathrm{Myr}=100`$ Myr. The clearest indication of its previous phase of core collapse and runaway growth would be the presence of a central compact object with a mass $`\stackrel{<}{}0.1`$% of the initial cluster mass. The cluster may also be relatively depleted in low-mass compact objects (stellar mass black holes and neutron stars), as these are consumed during the runaway growth phase.

Star clusters with an initial relaxation time $`t_{\mathrm{rlx}}\stackrel{>}{}25`$ Myr do not experience a phase of runaway growth, as core collapse is prevented by mass loss from the most massive stars. These clusters may experience core collapse after $`100`$ Myr, when stellar evolution slows (Takahashi & Portegies Zwart 1999). This later core collapse, however, does not lead to a phase of runaway growth. In such old clusters multiple collisions are still likely to be common and may lead to blue stragglers with a mass more than twice the turn-off mass.

We are grateful to Karl Gebhardt, Ortwin Gerhard, Pranab Ghosh, Douglas Heggie, Piet Hut, Jun Makino, and Rainer Spurzem for discussions and to the anonymous referee for a careful reading of the manuscript. We also thank Tokyo University, the Institute for Advanced Study, the American Museum for Natural History and the University of Indiana for their hospitality and the use of their fabulous GRAPE hardware. This work was supported by NASA ATP grants NAG5-6964 and NAG5-9264, and by the Nederlands Organization for Scientific Research (NWO). SPZ was supported as Hubble Fellow (HF-01112.01-98A) and as fellow of the Royal Netherlands Academy of Arts and Sciences (KNAW).


text:
# Introduction to Diffraction and low 𝑥 Dynamics

## 1 Introduction

The geometrical picture of diffraction is one in which the target is (essentially) left intact. The forward scattering amplitude is exponentially suppressed with momentum transfer squared $`t`$. The slope $`b`$ of the distribution, $`e^{bt}`$ is directly related to the size of the scattering target.

In $`ep`$ (and $`pp`$) scattering diffraction inevitably is also related to low $`x`$ such that the inelasticity $`y`$ remains $`1`$. Low $`x`$ naturally implies high mass $`W`$ since in $`ep`$ scattering $`W^2=Q^2(1x)/x`$. Parton densities at low $`x`$ have been seen to be large. For large $`Q^2`$ there is hence a region of phase space which is calculable perturbatively in QCD even at low $`x`$. Figure 1 shows the regions of phase space in a pictorial manner. For very small $`x`$ initial states of high partonic density are formed. While in these regions methods of perturbative expansion have to be proven to be valid it is interesting in its own right to search for experimental evidence for such new partonic configurations. In particular the relation between low $`x`$ processes and diffraction comes into play.

In diffraction the perturbative exchange has to be colourless and involves at least two gluons. Theoretically the interaction is thus modelled phenomenologically. A particularly useful description is placed in the proton rest frame. At low $`x`$ the $`\gamma `$ emitted from the incoming electron may turn into a $`q\overline{q}`$-dipole long before the interaction with the proton. The associated formation time $`c\tau 1/xm_p`$ is typically much larger than the size of the proton target.

Various configurations of the $`q\overline{q}`$ states may be examined. In *vector meson* production the pair forms a bound state that leaves the interaction unchanged. In *deeply virtual compton scattering* (DVCS) the $`q\overline{q}`$ state even recombines to form a real photon. The study of *inclusive diffraction* enables a test of these dipole configurations under controlled experimental conditions.

Vector meson production has been studied in much detail and a wealth of data both from HERA and the Tevatron have been quantitatively examined in terms of QCD.

## 2 Geometrical Picture

From its origin in optics diffraction is intimately associated with a geometrical interpretation \- the size of the target which participates as a whole. The $`t`$-distribution is exponentially suppressed. The energy dependence arises when exciting higher angular momentum states that become accessible with increased phase space. This behaviour is beautifully supported by the observed energy dependence of the total $`\gamma p`$-cross section.

Vector mesons exhibit quite distinct variations in their $`t`$-behaviour reminiscent of the different sizes at work: from spatially extended $`\rho `$ mesons to much smaller $`J/\mathrm{\Psi }`$ mesons. It is however interesting to compare the production of vector mesons in elastic and dissociative processes. While the elastic process is only affected by the formation of a vector meson, the dissociative process is modified by the form factor of the proton. Their ratio should hence be independent of the vector meson species. Figure 2 depicts the ratio as a function of the momentum transfer $`t`$ at the proton vertex, which seems to be described by a universal behaviour in support of this picture.

## 3 Inclusive $`ep`$ Scattering

So far the probe has been regarded as a quasi-real photon (without restrictions on transverse dimensions). For $`Q^20`$ the probe becomes essentially point like and the spatial interplay between probe and target has to be considered. The transition between real photoproduction and deeply inelastic scattering is therefore particularly interesting.

### 3.1 Inclusive $`ep`$ Scattering and the low $`x`$ Behaviour

For $`Q^2`$ well below the $`Z`$-boson threshold ($`Q^2<M_Z^2`$) the cross section

$$\frac{\mathrm{d}\sigma }{\mathrm{d}x\mathrm{d}Q^2}=\left(\frac{2\pi \alpha ^2}{Q^4x}(1+(1y))^2F_2y^2F_L\right)$$
(1)
is dominated by the structure function $`F_2`$; $`F_L`$ becoming relevant only at large $`y`$ and at small $`Q^2`$. Figure 3a) shows the transition of $`F_2`$ towards photoproduction where $`F_2`$ has to vanish $`Q^2`$ to yield a finite $`\gamma ^{}p`$ cross section. While the region $`Q^2>1`$ GeV2 is well described by recent QCD fits, the region below $`1`$ GeV2 relies on Regge inspired fits typically involving vector meson formation. For the region explored the energy dependence is smooth. There are no explicit signs of alteration of the behaviour as a function of $`xQ^2/W^2`$.

It is therefore interesting to study the small $`x`$ behaviour of the structure function for not so small $`Q^2`$, i.e. the region related to the high energy behaviour of the scattering process and accessible to perturbative methods. Figure 3b) depicts the derivative

$$\lambda =\left(\frac{\mathrm{ln}F_2}{\mathrm{ln}x}\right)_{Q^2}.$$
(2)
The distribution is experimentally consistent with being constant for $`x<0.01`$ over a large range in $`Q^2`$. The implication $`F_2=c(Q^2)x^\lambda `$ then yields a value for $`\lambda `$ that rises linearly with $`Q^2`$ (cf. figure 4a)) quite different from the photoproduction value $`\lambda =0.08`$. In the Regge language of *Pomeron* exchange the exchange becomes harder the larger $`Q^2`$.

### 3.2 Inclusive diffraction

Characteristic of the diffractive interaction is the large region void of hadronic activity - a gap in rapidity between the forward going nucleon and the centrally detected final state $`X`$, cf. figure 4b). Such a selection procedure leads to an inclusive measurement. Figure 5 shows the fraction $`\rho ^{D(3)}(\beta ,Q^2,x_{IP})=M_X^2x/Q^2F_2^{D(3)}(\beta ,Q^2,x_{IP})/F_2(x,Q^2)`$ of diffractively produced events as a function of $`Q^2`$ and $`x`$. With $`\beta Q^2/(Q^2+M_X^2)`$ and for two intervals of $`x_{IP}(Q^2+M_X^2)/(Q^2+W^2)`$ the fraction of diffractive events is considerable and the production persists up to the largest values of $`Q^2`$ of $`60`$ GeV2 investigated.

When analyzing the diffractive exchange in terms of a partonic description in QCD one finds a large gluonic component. To understand the formation of the final state more explicit final state measurements have to be investigated.

### 3.3 Inclusive Diffraction of Charm Mesons and Jets

Charm mesons originate from gluonic interactions and are hence a sensitive probe of the gluonic content of the interaction. Figure 6 (left) demonstrates a large charm component in diffractive processes. Similarly, multi-jet formation involves several partons. From an analysis of such final states, cf. figure 6 (right) one may conclude on the dynamical characteristics of the diffractive process. - As becomes evident, the description thus achieved is satisfactory and indicative of the large gluon content.

## 4 Comparison to Tevatron Jets

The observations are quite at variance to the understanding of diffraction at the Tevatron. In $`pp`$ scattering the density of hadronic matter is so large that an initially formed gap is not likely to survive unaffectedly the remainder of the transition through hadronic matter. The application of the partonic description of diffractive $`ep`$ data thus leads to a gross overestimate of the production at the Tevatron, cf. figure 7. The environment of hadron-hadron scattering is less amenable to a theoretical understanding of diffraction.

## 5 Conclusion

As wealth of data has been amassed at HERA and the Tevatron on low $`x`$ physics and diffractive processes. The HERA data provide the clues on further understanding of high density (partonic) final states which to some extent are a prerequisite for diffractive interaction to take place. While at HERA diffractive production can be understood in terms of QCD the Tevatron data show that more complex mechanisms are at work in hadron-hadron scattering. The reach to low $`x`$ currently available at HERA does not yield experimental evidence for saturation of partonic matter density.

## Acknowledgments

I wish to thank organisers of the 31st ISMD for a very pleasant stay in Datong and the kind hospitality extended to the participants. Special thanks go to Liu Lianshou, Wu Yuanfang and Frans Verbeure for organizing the symposium and to Jim Crittenden for arranging a very fruitful session on diffraction.


text:
# About one forgotten mechanism of nonlinearity in the theory of hot electrons

## Abstract

It is shown that in general electron gas heating inevitably results in the change of the carrier concentration in the conduction band. It is proved, that this change, as a rule, leads to the kinetic coefficient nonlinearity of the same order, as the change of mobility does. The conditions are determined, when this change can be neglected.

It is well known, that the fundamental reason for nonlinearity of a current-voltage characteristic (CVC) of a homogeneous semiconductor in strong electrical fields is the change of the mean carrier energy (carrier heating). The classical theory of the hot carrier transport was developed a long time ago and rather explicitly . Thus, as a rule, it was considered, that the nonlinearity of the CVC is related to the carrier mobility alteration because of the change of carrier mean energy. Some number of works was devoted to nonlinearity caused by impact ionization , carrier lifetime change , inter-valley redistribution of carriers  or by non-parabolic form of the carrier dispersion law  in strong fields. In single-valley semiconductors, neglecting such processes as an impact ionization and carrier lifetime change in strong fields, it is usually considered  that only the carriers, which already exist in bands, are subjected to heating, i.e. during heating the concentration of carriers remains equal to its value at the state of thermodynamic equilibrium. However, as it was shown in Ref. , there exists one more mechanism of nonlinearity, connected with the fact that the violation of the energy equilibrium between electrons and holes (the difference between electron and hole temperatures) inevitably results in the violation of the concentration equilibrium between electrons of the conduction band and the valence band. This idea was developed in Ref. . However, due to the assumption made that the population of the impurity level does not depend on heating, the results turn out to be incorrect if the heating of electrons and holes is different (electron and holes temperatures are unequal). Besides, the question when this effect takes place was left open.

In general, the discrepancy between electron and hole temperatures should cause the change of the carrier concentrations in both bands and the impurity level population. Hence, the problem is reduced only to what is the magnitude of the contribution of the latter effect to kinetic coefficients in comparison with the change of mobility (relaxation time) at the same temperature discrepancy. Thus, a new origin of strong field nonlinear effects, related to the alteration of the energy level population in conduction and valence bands due to the difference between electron and hole temperatures in strong electrical fields, is discussed below. It is possible to neglect this phenomena in the theory of hot electrons only under special conditions indicated below.

Experimental verification of the magnitude of the considered effect is shown in Ref. . Unfortunately, in spite of all these facts, the considered mechanism of nonlinearity in the theory of hot electrons had been actually forgotten and did not obtain a further development.

Virtually, the carrier concentration change because of this mechanism is reduced to the carrier recombination rate change, which is owing to the alteration of carrier distribution functions.

Let’s consider the elementary model — a homogeneous single-valley semiconductor with one nondegenerate impurity level at an energy $`\epsilon _t`$ and concentration of impurity $`N_t`$. In the given model the kinetics of the carrier concentration change within bands due to heating is determined by the following processes: (1) the capture of electrons from conduction band to the impurity level, (2) the thermal emission of electrons from the impurity level into the conduction band, (3) the capture of holes to the impurity level, (4) the thermal emission of holes from the impurity level into the valence band. To simplify the calculations we neglect the influence of interband transitions on carrier concentration change at heating.

Let’s note, that the concentration change is caused by the alteration of the rate of recombination rather than thermal generation.

Suppose, that only the electron subsystem is subjected to heating. Then, if conditions for carrier concentrations indicated in  are fulfilled, the electron gas can be described by the Fermi distribution function with electron temperature $`T_e`$. Subsystems of holes and captured carriers have the lattice temperature $`T_0`$.

The capture rate of electrons onto the impurity level can be represented by following expression :

$$r_n=\alpha _n(T_e)N_t[1f_t(T_0)]n(T_e),$$
(1)
where $`f_t`$ is the distribution function of electrons on the impurity level, $`N_t[1f_t(T_0)]`$ represents the concentration of free impurity states, $`n(T_e)`$ is the concentration of electrons, $`\alpha _n(T_e)`$ is the capture factor of electrons by the trap. Let’s emphasize once again, that here we do not consider the explicit dependence of the capture factor on the magnitude of the electrical field applied, i.e. we do not take into account such processes as the change of the carrier life time in strong fields (see ). By definition,

$$\alpha _n(T_e)=\frac{_{\epsilon _c}^{\mathrm{}}c_n(\epsilon )\nu _n(\epsilon )f_n(\epsilon ,T_e)𝑑\epsilon }{_{\epsilon _c}^{\mathrm{}}\nu _n(\epsilon )f_n(\epsilon ,T_e)𝑑\epsilon }.$$
(2)
Here $`c_n(\epsilon )`$ is a probability of the electron transition from the impurity level into the state with an energy $`\epsilon `$, $`\nu _n(\epsilon )`$ is a density of states, $`f_n(\epsilon ,T_e)`$ is the Fermi distribution function with temperature $`T_e`$ for conductivity electrons, $`\epsilon _c`$ is the energy of an electron at the bottom of conduction band.

The rate of a thermal emission into the conduction band is assumed to be independent of the temperature of electrons in the conduction band (that is correct, at least, for wide-band semiconductors):

$$g_{nT}=\alpha _n(T_0)N_tf_t(T_0)n_1,$$
(3)
where $`n_1\nu _{n0}\mathrm{exp}\left(/T_0\right)`$ is a parameter describing the impurity level, $`\nu _{n0}`$ is the effective density of states in the conduction band, $`\epsilon _c\epsilon _t`$ is the ionization energy of the impurity level. The parameter $`n_1`$ represents concentration of electrons in the conduction band, which would take place, if the Fermi level would coincide with the impurity level.

Obviously, the recombination rate of electrons is equal to $`R_n=r_ng_{nT}`$.

Similar equations can be written for the hole subsystem as well:

$`r_p=\alpha _p(T_0)N_tf_t(T_0)p,`$ (4)
$`g_{pT}=\alpha _p(T_0)N_t[1f_t(T_0)]p_1,R_p=r_pg_{pT}.`$ (5)
Here $`p`$ is a hole concentration, the definition of quantities $`\alpha _p(T_0)`$ and $`p_1`$ is similar to the above mentioned one for electrons.

In a homogeneous semiconductor in steady state $`R_n=R_p=0`$, and the condition of the electroneutrality is fulfilled. It is convenient to express such a condition as:

$$\delta n+N_t\delta f_t=\delta p,$$
(6)
where $`\delta nnn_0`$, $`\delta ppp_0`$, $`\delta f_tf_tf_t^0`$, and $`n_0`$, $`p_0`$ and $`f_t^0`$ are respectively electron and hole concentrations and the impurity level population in the absence of heating ($`T_e=T_0`$).

After solving the set of equations $`R_n=R_p=0`$ together with the Eq. (6), one obtains the following expressions for the carrier concentration change caused by heating:

$`\delta n`$ $`=`$ $`\left[1+{\displaystyle \frac{p_1(n_0+n_1)^2}{n_1(p_0+p_1)^2+N_tn_1p_1}}\right]^1`$ (8)
$`\times {\displaystyle \frac{n_0}{\alpha _n(T_0)}}{\displaystyle \frac{\alpha _n(T_0)}{T}}\delta T,`$
$`\delta p`$ $`=`$ $`\delta n\left[1+{\displaystyle \frac{N_tp_1}{(p_0+p_1)^2}}\right]^1,\delta TT_eT_0T_0.`$ (9)
Assuming additionally, that the gas of carriers is nondegenerate, i.e. the relation $`n_0p_0=n_1p_1=n_i^2`$ holds, where $`n_i`$ is the intrinsic carrier concentration (i.e. at $`N_t=0`$) in the absence of heating, and, besides, bands are parabolic, we will analyze Eqs. (8)–(9) in two limiting cases.

1. An intrinsic semiconductor with low concentration of shallow traps for electrons .
In this case $`n_1n_0=p_0=n_ip_1`$ and Eqs. (8)–(9) acquire the form:

$`{\displaystyle \frac{\delta n}{n_i}}`$ $`=`$ $`{\displaystyle \frac{\alpha _n(T_0)}{T}}{\displaystyle \frac{\delta T}{2\alpha _n(T_0)}},`$ (10)
$`\delta p`$ $``$ $`\delta n.`$ (11)
Thus the magnitude of the carrier concentration change at heating in this case is determined only by the temperature dependence of the electron capture factor $`\alpha _n`$.

Let’s note also, that expressions (10)–(11) have the same form as in the case of only interband transition in an intrinsic semiconductor.

2. A $`n`$-type monopolar semiconductor with donor impurity ($`n_0n_ip_0`$).
The relative carrier concentration change in this case is described by the following expressions:

$`{\displaystyle \frac{\delta n}{n_0}}`$ $`=`$ $`\left[1+1/\left({\displaystyle \frac{n_i^2}{n_0^2}}+{\displaystyle \frac{N_t}{n_1}}\right)\right]^1{\displaystyle \frac{\alpha _n(T_0)}{T}}{\displaystyle \frac{\delta T}{\alpha _n(T_0)}},`$ (12)
$`\delta p`$ $`=`$ $`\delta n\left[1+{\displaystyle \frac{n_0^2}{n_i^2}}{\displaystyle \frac{N_t}{n_1}}\right]^1.`$ (13)
It is easy to verify, that in this case the carrier concentration change at heating is determined not only by the temperature dependence of the electron capture factor $`\alpha _n`$, but also depends on concentration of impurity and on the temperature $`T_0`$. In the range of high temperatures ($`T_0/\mathrm{ln}(\nu _{n0}/N_t)`$, i.e. $`N_tn_1`$) the functional dependence of $`|\delta n/n_0|`$ on $`N_t`$ has a deep minimum at $`N_t\sqrt[3]{n_i^2n_1}`$, being equal to

$`\left|{\displaystyle \frac{\delta n}{n_0}}\right|_{\mathrm{min}}`$ $`=`$ $`\left({\displaystyle \frac{2n_i}{n_1}}\right)^{2/3}{\displaystyle \frac{\alpha _n(T_0)}{T}}{\displaystyle \frac{\delta T}{\alpha _n(T_0)}}`$ (15)
$`\mathrm{exp}\left({\displaystyle \frac{\epsilon _g2}{3T}}\right)1,`$
where $`\epsilon _g`$ is the bandgap width. Thus, near the indicated concentration the additional contribution to the conductivity change is negligibly small and heating effects are described by existing theories .

At low temperatures or, that is equivalent, heavily doping ($`N_tn_1`$) we come back to expression (10). Let’s note, that in this case (in contrast to Eq. (11)) $`\delta p\delta n`$.

In monopolar case, as well as for an intrinsic semiconductor, the deviation from the conventional theories is determined by the temperature dependence of the electron capture factor $`\alpha _n`$. For simplicity we will analyze model dependence for two cases: attracting and repulsive potentials of impurity centers.

In the case of the electron capture by an attracting potential $`\alpha _n(T)T^m`$, and $`m`$ varies within the limits from $`m1`$ up to $`m5`$ depending on the nature of the semiconductors and impurities . Then

$$\frac{\delta T}{\alpha _n(T_0)}\frac{\alpha _n(T_0)}{T}=m\frac{\delta T}{T_0}.$$
(16)
In the case of electron capture by a repulsive potential of impurity (for example, ions of gold or copper in germanium) the temperature dependence of the capture factor is satisfactorily described by the expression :

$$\alpha _n(T)\mathrm{exp}\left[\left(\frac{T^{}}{T}\right)^{1/3}\right],$$
(17)
where $`T^{}`$ is a parameter depending on the specific kind of the semiconductor and the impurity in it.

Hence,

$$\frac{\delta T}{\alpha _n(T_0)}\frac{\alpha _n(T_0)}{T}=\frac{1}{3}\left(\frac{T^{}}{T_0}\right)^{1/3}\frac{\delta T}{T_0}.$$
(18)
As long as $`T^{}`$ usually lies in the range from $`10^4\mathrm{K}`$ to $`10^9\mathrm{K}`$ , in the case of a repulsive potential the contribution to the conductivity change in the heating field can be even more essential.

Thus, the electron concentration change caused by carrier heating alters kinetic coefficients as much as the mobility (relaxation time) change does. It should be mentioned that there exist only two specific, not very interesting, situation, namely: (1) the electron Fermi level lies far enough both from the middle of the gap and from the impurity level, and (2) the electron capture factor weakly depends on the temperature, when the traditional theory of hot electrons is correct.

We wish to thank Dr. F. Pérez-Rodríguez for helpful discussion. This work has been partially supported by CONACyT–México.


text:
# Primordial Helium Production in a Charged Universe

## Abstract

We use the constraints arising from primordial nucleosynthesis to bound a putative electric charge density $`|e|n_q`$ of the universe. We find $`|n_q/n_\gamma |10^{43}`$, four orders of magnitude more stringent than previous limits. We also work out the bounds on $`n_q`$ in models with a photon mass, that allows to have a charge density without large-scale electric fields.

preprint: UAB-FT-519

One of the fundamental parameters of the universe is its electric charge. It is usually assumed to be exactly vanishing, but it is of course desirable to have observational evidences of the hypothesis of a neutral universe. In this paper we will discuss constraints on $`n_q`$, defined in such a way that $`|e|n_q`$ is the electric charge density of the universe ($`e`$ is the electron charge). It is convenient to normalize this quantity to the photon number density, so we define the parameter

$$\eta _q\frac{n_q}{n_\gamma }$$
(1)
with

$$n_\gamma =\frac{2\xi (3)}{\pi ^2}T^3$$
(2)
where $`T`$ is the temperature, and $`\xi (3)1.20`$.

Speculations about a possible electric charge of the universe go back to the work of Lyttleton and Bondi lyttleton , who showed that an excess charge could account for the expansion of the universe. They assumed that the excess charge was due to a tiny charge imbalance among the proton and the electron, $`q_p+q_e0`$. Laboratory limits dylla ; Groom:2000in  on this imbalance rule out the assumption of Lyttleton and Bondi. Their argument can be somehow inverted Orito:1985cf : demanding that the gravitational attraction among cosmological objects is larger than the electromagnetic repulsion leads to the bound

$$|\eta _q|10^{28}$$
(3)
More recent speculations involving an excess charge are exotic suggestions of violation of electric charge conservation. Grand unified theories imply violations of baryon and lepton numbers, but conservation of electric charge is usually believed to be exact as a result of gauge invariance. Nonetheless one has the possibility of a (necessary small) breakdown of electric charge conservation at high energies voloshinetal . Also, we may have an excess charge in theories with extra dimensions, first introduced by Kaluza and Klein kaluzaklein  and subject to a recent strong revival extradim . The extra dimensions should be compactified to a small size, since our world appears four dimensional. To probe them one can consider the universe in its first instants; the energy is so high that one has a truly higher-than-four dimensional space. Subsequently, the universe evolves and cools down but some electric charge could remain in the ordinary space. Some authors have built classes of models where there is charge leakage into infinite extra dimensions dubovsky .

It is therefore interesting to examine the implications of a charge density $`n_q`$ and to find constraints on it. Orito and Yoshimura Orito:1985cf  examined some of these implications. They considered the anisotropies generated by the large-scale electric field that would be produced by a net charge density. Their best bound is obtained using observational limits on cosmic ray anisotropies:

$$|\eta _q|10^{39}$$
(4)
In this letter we show that a much more stringent bound can be obtained by considering the primordial nucleosynthesis period of the early universe and the change in the primordial helium production that would be originated by a charge density present at the nucleosynthesis era.

The physics involved in the primordial nucleosynthesis period is well understood and the theoretical predictions of the primordial yields of light elements are robust Sarkar:1996dd ; Olive:2000ij . The agreement with observation is considered one of the pillars of modern cosmology. For our purposes, it is useful to recall here the simple arguments that allow to understand the main physical features of helium production Sarkar:1996dd ; Kolb:1990vq . In the early universe, at time $`t1`$s and temperature $`T1\mathrm{MeV}`$, neutrons and protons are in kinetic and chemical equilibrium due to weak interactions. At a lower $`TT_f`$, the rate of these interactions becomes less than the expansion rate of the universe and they go out of equilibrium. Then, the relative neutron-to-proton density is given approximately by

$$\left(\frac{n}{p}\right)_fe^{Q/T_f}0.24$$
(5)
where

$$Q=m_nm_p=1.29\mathrm{MeV}$$
(6)
We have taken $`T_f0.9\mathrm{MeV}`$. After freeze-out, there is a decrease in the neutron number due to $`\beta `$-decay,

$$npe^{}\overline{\nu }$$
(7)
until the time $`t_{He}`$ where all the helium is produced. Then

$$\left(\frac{n}{p}\right)_{He}e^{t_{He}/\tau _n}\left(\frac{n}{p}\right)_f0.15$$
(8)
Here we have introduced the neutron lifetime Groom:2000in  $`\tau _n=887`$s and we have taken $`t_{He}400`$s.

Since nearly all neutrons are processed into helium, the expected mass fraction $`Y`$ of 4He is

$$Y\frac{2(n/p)_{He}}{1+(n/p)_{He}}0.26$$
(9)
The actual theoretical prediction for $`Y`$ has to be obtained by a numerical code Wagoner ; Kawano  that solves the relevant set of ordinary differential equations, and indeed gives a result not far from (9).

Primordial nucleosynthesis is known to be a useful tool to constrain non-standard physics Sarkar:1996dd . We now apply it to the issue of a charged universe. A charged particle in such universe has a potential $`V`$ (with the condition $`V=0`$ when $`n_q=0`$) and thus it has a momentum $`p`$ related to its total energy $`E`$ by

$$(EV)^2=p^2+m^2$$
(10)
With the introduction of an effective mass $`m^{}`$, the dispersion relation can be written as

$$E^2=p^2+(m^{})^2$$
(11)
We now assume $`VE`$ and that the charged particle is non-relativistic. We get

$$m^{}m+V$$
(12)
A similar phenomenon is at the basis of the MSW effect MSW ; eqs (10-12) in the relativistic case are discussed in bethe .

The horizon distance $`a_H`$ is defined as the propagation time of light since $`t=0`$, and then for $`r>a_H`$ the charge density could not interact with the charged particle. Thus, to the potential $`V`$ contributes all the charge inside a sphere of radius $`a_H`$

$$V=\pm _0^{a_H}(dr\mathrm{\hspace{0.17em}4}\pi r^2)\frac{\alpha n_q}{r}$$
(13)
To summarize, a net charge density $`|e|n_q`$ induces effective masses for the charged particles present in the universe at that time

$`m_pm_p^{}`$ $`=`$ $`m_p+\delta m`$
$`m_em_e^{}`$ $`=`$ $`m_e\delta m`$ (14)
(positrons are affected with the opposite sign than electrons), where from (12) and (13)

$$\delta m=2\pi \alpha a_H^2n_\gamma \eta _q$$
(15)
As we are in a radiation dominated universe we will put $`a_H=2t`$.

A non-vanishing $`\delta m`$ alters the predictions of primordial nucleosynthesis for the helium yield, and this leads to a bound on $`\delta m`$ and hence on $`\eta _q`$. Before presenting our numerical results let us show that we can get an approximate expression for the change in $`Y`$ due to the $`\delta m`$ shift, using the same type of simple analysis that we used in our discussion of the standard $`Y`$, from eq. (5) to (9). We work at first order in $`\delta m`$.

There are two main sources of change in $`Y`$. The first is through the dependence of $`Y`$ in eq. (5) on $`Q`$ at freeze-out,

$$\frac{\delta Y}{Y}\frac{1}{1+(n/p)_f}\frac{\delta m}{T_f}\frac{\delta m}{\mathrm{MeV}}$$
(16)
where we have introduced the numerical values for $`T_f`$ and $`(n/p)_f`$.

The second is due to the dependence of $`\tau _n`$, the lifetime of process (7), on the proton and electron masses. When decaying into a proton of mass $`m_p^{}`$ and an electron of mass $`m_e^{}`$, the neutron lifetime (for zero temperature) depends on masses as

$`{\displaystyle \frac{1}{\tau _n^{}}}`$ $``$ $`m_e^{\mathrm{\hspace{0.17em}5}}\lambda (q^{})`$ (17)
$`\lambda (q^{})`$ $`=`$ $`{\displaystyle _1^q^{}}𝑑xx(xq^{})^2\sqrt{x^21}`$ (18)
$`=`$ $`{\displaystyle \frac{1}{60}}\sqrt{q^{}1}(2q^49q^28)`$
$`+{\displaystyle \frac{1}{4}}q^{}\mathrm{ln}(\sqrt{q^21}+q^{})`$
with $`q^{}=(m_nm_p^{})/m_e^{}`$

After some algebra one finds the change in $`\tau _n`$

$$\frac{\delta \tau _n}{\tau _n}\left[5(q1)\left(\frac{1}{\lambda }\frac{d\lambda }{dq^{}}\right)_{q^{}=q}\right]\frac{\delta m}{m_e}$$
(19)
with $`q=(m_nm_p)/m_e`$. From (8) it is easy to see that a shift $`\delta \tau _n`$ induces a change $`\delta Y`$

$$\frac{\delta Y}{Y}\frac{1}{1+(n/p)_{He}}\frac{t_{He}}{\tau _n^2}\delta \tau _n$$
(20)
Introducing the numerical values for the parameters that appear in (19) and (20), we end up with

$$\frac{\delta Y}{Y}\frac{\delta m}{\mathrm{MeV}}$$
(21)
To find the expressions (16) and (21) we have made the assumption that $`\delta m`$ is independent of time, which is not true, since $`a_H=2t`$ and $`n_\gamma T^3t^{3/2}`$ in (15). Still, we are able to find the right order of magnitude for our bound on $`\eta _q`$ if we proceed as follows. Since (16) is approximately valid at freeze-out, we evaluate $`\delta m`$ in (15) using $`t1`$s, $`T1\mathrm{MeV}`$ and get

$$\delta m10^{41}\eta _q\mathrm{MeV}$$
(22)
Similarly, (21) is approximately valid during the period of neutron decay. Now we evaluate $`\delta m`$ in (15) with typical values for this period, $`t100`$s, $`T0.1\mathrm{MeV}`$, and get

$$\delta m10^{42}\eta _q\mathrm{MeV}$$
(23)
Comparing (22) and (23), we see that it is the change in $`Y`$ due to $`\delta \tau _n`$, eq. (23), that will dominate the whole effect. The origin of this dominance is that, as we see from (15), $`\delta mt^2T^3t^{1/2}`$, i.e., $`\delta m`$ increases with time. To estimate the order of magnitude of the bound on $`\eta _q`$, we put (23) in (21) and allow, for instance, a $`10\%`$ change in $`Y`$. We get $`|\eta _q|10^{43}`$

To find a precise bound, we have implemented the modification (14) and (15) in Kawano’s version Kawano  of the primordial nucleosynthesis code of Wagoner Wagoner  to obtain the effect of a density $`n_q`$ on the helium abundance. The abundance depends on the laboratory neutron lifetime $`\tau _n`$, the number of neutrinos $`N_\nu `$, and the baryon-to-photon ratio $`\eta _B`$. We use $`\tau _n=887\pm 2`$ s Groom:2000in , $`N_\nu =3`$ Groom:2000in , and for $`\eta _B`$ the value extracted from CMB anisotropy measurements CMB . Unfortunately, this value is still not completely settled but the hope is that in the future it will be with more refined experiments. For the time being, we take

$$10^{10}\eta _B10^9$$
(24)
a generous range that embraces the observations CMB .

The prediction for $`Y`$ as a function of $`\eta _q`$ is shown in Fig.1, for the two extremes of the range (24). (The experimental error in $`\tau _n`$ introduces a negligible change in our prediction.) The linear dependence is due to the expected dominance of the first-order in $`\eta _q`$. Our bound is stringent because $`Y`$ is measured with relatively small error. The observational status is discussed in Sarkar:1996dd ; Olive:2000ij . We adopt the range

$$0.228Y0.248$$
(25)
that the authors of Olive:2000ij  claim is “$`95\%`$ C.L.”. As we see from Fig.1, combining the experimental range (25) with the predictions for $`Y`$ in the range (24) constrains the charge density of the universe. The bound is

$$1.1\times 10^{43}\eta _q0.8\times 10^{43}$$
(26)
We notice that it is four orders of magnitude better than (4), and that it is not far from our estimate using (23).

A non-vanishing $`\eta _q`$ would also affect the yields of the other light elements D, 3He, and 7Li. However, since these yields are not as well measured as 4He, taking them into account could not significantly improve our bound (26).

Even if an hypothetical charge density of the universe $`n_q`$ has to be so tiny, one may worry about the induced large-scale electric fields that appear when $`n_q0`$, no matter how small. In fact, there are models having $`n_q0`$ but no large-scale fields. This can be achieved by endowing the photon with a small mass $`m_\gamma `$ Barnes\_Barry  and thus having an electromagnetic interaction of finite range $`\lambda `$. Bounds, like (4), that are based on the effects of a large-scale electric field are very much weakened in this kind of models (to get (4) the authors of Orito:1985cf  consider a scale $`500`$ pc for the electric field).

Let us show that our bound (26), when $`m_\gamma 0`$, is only modified by one order of magnitude. First, we notice that the mass shift $`\delta mt^{1/2}`$ shown in (15) is valid when $`a_H=2t<\lambda `$. However, for later times, $`\delta m`$ gets contribution only up to a radius $`\lambda `$,

$$\delta m=\frac{4\alpha }{\pi }\xi (3)\eta _q\lambda T^3t^{3/2}$$
(27)
It follows that $`\delta m`$ increases until $`t=\lambda /2`$ and afterwards it decreases.

Since helium production finishes when $`tt_{He}`$, it is clear that our bound (26) is still valid for $`\lambda >t_{He}`$. For smaller $`\lambda `$, we expect to find bounds on $`n_q`$ that are less severe. Obviously, for $`\lambda 0`$ the bound would disappear. However, $`\lambda `$ is subject to the experimental constraint

$$\lambda =\frac{1}{m_\gamma }10^9\mathrm{m}$$
(28)
coming from studies on torques on a toroid balance Lakes:1998mi . Similar limits are obtained from measurements of the Jovian magnetic field Davis:1975mn . We notice that the lower limit $`\lambda =10^9`$m is about the same as $`t_f1\mathrm{s}=3\times 10^9`$m. Thus, to estimate the bound on $`|\eta _q|`$ we can use (22) since it is valid for $`t1`$s. We expect to get $`|\eta _q|10^{42}`$, so that in models with a photon mass we expect to have a bound that is only about one order of magnitude worse than (26).

We have modified Kawano’s code Kawano  with the mass shift $`\delta m`$ shown in (27) to find the precise change in $`Y`$. By demanding that the predicted $`Y`$ is in the experimental interval (25), and allowing for $`\eta _B`$ in the range (24), we are able to limit $`\eta _q`$. In Fig.2, we show our bound as a function of $`\lambda `$ and also the experimental limit (28). We confirm that for $`\lambda 10^{11}`$m we get our previous limit (26) and we also see that the constraints relax for smaller $`\lambda `$. We can put a limit on $`\eta _q`$ in models with a photon mass, which clearly is the one for $`\lambda =10^9`$m. It reads

$$1.6\times 10^{42}\eta _q1.8\times 10^{42}$$
(29)
###### Acknowledgements.

Discussions with F. Ferrer, J. Garriga, J. A. Grifols, S. Sarkar, G. Senjanović and R. Toldrà are greatfully acknowledged. Work partially supported by the CICYT Research Project AEN99-0766, and by the EU network on Supersymmetry and the Early Universe (HPRN-CT-2000-00152)