text:
# Star Streams in the Milky Way Halo

## 1. Introduction

A revolution is underway in our understanding of how galaxies form: The idea of a monolithic collapse of a single system at earliest times (Eggen et al. 1962) has been challenged by observations at high redshift and by increasingly sophisticated models of the evolution of structures within the basic framework of the Hot Big Bang model (for example, Pearce et al. 1999, Steinmetz & Navarro 1999). These findings imply that we should see clear evidence of hierarchical formation processes within nearby galaxies. Indeed, some of the earliest work that suggested complex and possibly hierarchical evolutionary histories for galaxies such as ours came from studies of the stellar populations of the outer halo of the Milky Way (Searle and Zinn 1978). Studies of nearby and high-redshift systems represent complementary approaches.

We can learn much about accretion processes in galaxy formation using the Milky Way because we can view our halo in 3 dimensions. By adding the kinematics of halo stars to their 3-D positions, moving groups are given enhanced contrast (see, for example, Harding et al. 2001). This sort of study is impossible for any other spiral galaxy.

Recent studies of our Galaxy have begun to reveal tantalizing evidence of possible substructure in the halo, both locally (Majewski 1992, Côté et al. 1993, Arnold and Gilmore 1992, Helmi et al. 1999) and at larger distances (Ivezić et al. 2000, Yanny et al. 2000). And of course, the Sgr dwarf galaxy (Ibata et al. 1995) is injecting globulars and field stars that formed in a completely separate entity from the Milky Way into its halo. The question has now shifted from ‘Did accretion play any role in the formation of the Galaxy?’ to ‘How much of the Galaxy formed from accretion of hierarchical fragments?’ We will show below that the data that is needed to answer the second question definitively also has the potential to teach us about the dark matter halo of the Galaxy and may directly constrain cosmological formation models such as that of Klypin et al. (1999), Moore et al. (1999).

## 2. The “Spaghetti” Survey

Several years ago we began a search designed to find halo streams if they are present. Since then, the Sloan Digital Sky Survey (SDSS) team has found one stream in their first few hundred degrees of equatorial test data (Ivezić et al. 2000, Yanny et al. 2000). Our group has found evidence that has strengthened the case that this substructure is associated with the Sgr dwarf using velocities of halo giants (Dohm-Palmer et al.2001). Our survey is unique in using a set of halo tracers (giants) that can be observed spectroscopically out to distances of hundreds of kpc where it is likely that only dwarf galaxies are found.

Our photometric survey consists of several hundred pencil beams at high Galactic latitude, which will comprise 100 square degrees of data in the Washington system. Our goal is to use this survey, the spectroscopic followup, and our extensive modelling, to quantify the amount of halo substructure associated with completely destroyed satellites, with partially destroyed satellites like the Sgr dwarf, and with intact objects such as $`\omega `$ Cen. More details of the survey are available in Morrison et al. (2000), Dohm-Palmer et al.(2000), Morrison et al. (2001) and at http://www.astro.lsa.umich.edu/users/rdpalmer/spag/. In this paper we focus on the part of the survey which identifies distant halo giants.

Local studies of halo stars passing through the solar neighborhood are limited in numbers of stars, do not well sample all regions of the halo, and often have biased kinematics. Chiba and Yoshii(1998) have shown that the 250 red giants and RR Lyrae stars with space motions determined from Hipparcos occupy a region extending to at most 35 kpc from the Galactic center. Carney et al. (1996) derive the apogalacticon distances of approximately 1200 local stars. There are only 23 stars in this sample with calculated apogalacticon distances larger than 30 kpc. These particular stars, by definition, are on extremely elongated orbits. There are vast parts of the Milky Way, and vast parts of phase space, therefore left unsampled. Figure 1 illustrates this clearly: halo globular clusters are plotted with open symbols, and it can be seen that they extend out to almost 100 kpc from the Galactic center. Known field stars (taken from the compilation of Beers & Sommer-Larsen 1995) are shown with closed circles. Most halo field stars known are within a few kpc of the Sun.

Figure 2 shows the spectrum of one of the most distant giants we have identified so far, star l304.49b+60.51, which lies at a distance of $``$80 kpc from the Sun. We expect to identify several hundred giants from the extreme outer halo when our 100 square degree survey is complete: an order of magnitude more than are currently known.

## 3. Searching for Substructure in the Outer Halo

Halo red giants are useful probes because of their high luminosity, which makes spectroscopic followup much easier. However, they are rare: we average $``$1 giant in each 0.25 square degree field. This means that we cannot use many of the techniques which were developed to search for substructure in the larger local samples. These techniques are often one- and two-dimensional variations on the use of velocity histograms, whether to search for a spike (Côté et al. 1993), unusual outliers (Majewski 1992) or simply a non-gaussian velocity distribution (Harding et al. 2001).

We have used the extensive simulations of Harding et al. (2001), re-normalized so that one particle corresponds to one halo giant, to help us identify and test techniques to search for substructure. Our method uses the fact that that there are strong correlations between position and velocity along a disrupting stream, illustrated in Figure 3. Our sparse sampling reduces the numbers of stars detected in a given stream, but preserves the position/velocity correlations, as can be seen in Figure 3.

Since we can measure the distances to our giants with reasonable accuracy, we are able to derive all three spatial coordinates, but because they are too distant for proper motions until space missions such as GAIA and SIM are operating, we only have one of the three velocity coordinates, the line-of-sight radial velocity $`V_{los}`$. We define a scaled metric in these four dimensions to measure the correlations between distance and velocity: two stars with similar distances and velocities will have a very small 4-D “distance” apart.

For two stars with galactic latitude, galactic longitude, distance and velocity $`l_1,b_1,d_1,V_1`$ and $`l_2,b_2,d_2,V_2`$, we write the metric as:

$$D_{lbdV}=\left(w_1(l_1l_2)^2+w_2(b_1b_2)^2+w_3(d_1d_2)^2+w_4(V_1V_2)^2\right)^{1/2}$$

We choose the scalings $`w_i`$ for the l, b, distance and velocity terms to account for the different measurement accuracy of each quantity, and for the fact that streams can spread spatially but have a tight relation between position and velocity at a given location. We know l and b very accurately, velocity to $``$10% (20 km/s) and distance to $``$25%, so two values that differ by less than their errors will be given less weight than values that differ significantly.

We use two different simulated datasets to test this technique, with sample sizes equal to our expected number of giants when we have completed our survey. The first contains 500 points from an accreted halo simulation constructed with 23 disrupting tidal streams from the library of Harding et al. (2001). These points were chosen from a dataset containing 320,000 points in 23 streams in a way that simulates our real observations. We used our real field centers for the 50 square degrees where observations have been taken, and planned field centers plus reflections about the galactic plane to make 400 pencil-beam field centers at high galactic latitude. If the stream particle was in the 0.25 square degrees surrounding a field center, we applied observational errors of 20 km/s for velocity and 25% for distance. The second dataset contained points from a smooth halo constructed with an isotropic velocity ellipsoid, as described in Harding et al. (2001), chosen to simulate our real observations in the same way.

We calculate the value of the scaled 4-d metric described above for every pair of points in the smooth and accreted halo simulations. As expected, the accreted halo shows significantly more points with small values of the metric, as can be seen in Figure 4. We are working with Jiayang Sun (CWRU statistics) to develop sensitive statistical tests for this accretion signature.

## 4. First Results

Here we show some preliminary results for our sample of 58 halo giants with spectroscopic confirmation. The distances are calculated using the photometric metallicities only at this point, although we will soon have spectroscopic metallicities available, which will lead to more accurate distance measures. Figure 5 shows the histogram of metric values for this sample (calculated as described above) and the histogram of metric values for a smooth halo sample with the same size as our real halo giant sample, chosen from the simulated smooth halo described above. While the sample size is small, we are already seeing some of the signatures of an accreted halo.

Some of the smallest values of the 4-d metric in our giant sample come from a group of stars which are located in the same region of the sky as the SDSS overdensity (Ivezić et al. 2000, Yanny et al. 2000). This gives us the opportunity to test suggestions that this stream is associated with the Sgr dwarf by comparing the velocities of giant stars at similar distances to the SDSS BHB and RR Lyrae stars to the predictions of models of the disruption of Sgr. We have chosen to use the models of Helmi & White (2001), which, unlike earlier models of its disruption, can explain the fact that the galaxy was not completely disrupted on its first few passages without resort to unusual scenarios such as a very stiff and extended dark matter halo (Ibata & Lewis 1998) or a convenient recent collision with the Large Magellanic Cloud (Zhao 1998).

Figure 6 shows the predictions of the purely stellar model of Helmi & White (2001) for particles close to the minor axis of the Galaxy, compared with the actual positions and velocities of 21 confirmed metal-poor giants in this direction. It can be seen that all of the halo giants with distances greater than 30 kpc agree well with the predictions of the Sgr disruption model, confirming the suggestions of Ivezić et al. (2000) and Ibata et al. (2001a) that the SDSS stream is associated with Sgr. The right-hand panel of the plot shows, for comparison, one simulation taken from a smooth halo: it can be seen that the arrangement of distances and velocities that we observe is very unlikely to have occurred by chance. More details of these observations can be found in Dohm-Palmer et al. (2001).

The wrap at distances near 50 kpc is the one detected by the SDSS. It is interesting to note that we may also have detected stars from the wraps near 20 kpc (although the large velocity range predicted for these wraps makes this harder to test) and from the wrap at 80 kpc, which contains debris lost from the Sgr dwarf on one of its earliest passages. Such debris, if confirmed with a larger sample, will give us the opportunity to measure the M/L ratio of the progenitor to the Sgr dwarf, as the outer regions are the first to be lost to tidal disruption, and the presence of a massive dark matter halo in the outer regions will help retain more stars than in a galaxy with little or no dark matter (Helmi 2000).

Our sample of giants from 100 square degrees of imaging will be large enough to assign membership to different streams in the outer halo. This will allow their use in ways that were almost inconceivable ten years ago. For example, we can learn about the structure of the progenitor satellites: the red giants that are Sgr stream members have a mean metallicity of \[Fe/H\]=–1.5, significantly more metal-poor than the main body of Sgr (Layden and Sarajedini 2000). This suggests that the Sgr progenitor had a radial abundance gradient. As mentioned above, identification of stars lost on early passages of Sgr will allow us to constrain the mass-to-light ratio of its progenitor. The spread in longitude of the Sgr debris will constrain the flattening of the Galaxy’s dark matter halo (Helmi 2000, Ibata et al. 2001b). We can even use the stream velocity dispersion to constrain the lumpiness of the Galaxy’s dark halo, which will allow us to directly test the predictions of simulations of structure formation (eg Klypin et al. 1999, Moore et al. 1999).

### Acknowledgments.

The Galaxy’s halo and its formation is a subject which is close to Ken’s heart: he supervised two dissertations on the topic (Ratnatunga’s and Morrison’s). The technology available at the time (photographic plates) wasn’t equal to the challenge of studying the outer halo, so it’s a particular pleasure to be able to talk about a successful survey (and one that he is involved in) for his 60th year celebration. It looks like his early interest in accreted halos was spot on: the instrumentation just had to catch up!

HLM, EWO and MLM acknowledge support from the NSF on grants AST 96-19490, AST 00-98435, AST 96-19524, AST- 00-98518, AST 95-28367, AST 96-19632 and AST 00-98661.

## References

Arnold, R. and Gilmore, G. 1992, MNRAS, 257, 225

Beers, T. C. & Sommer-Larsen, J. 1995, VizieR On-line Data Catalog: J/ApJS/96/175. Originally published in: ApJS, 96, 175

Carney, B. W., Laird, J. B., Latham, D. W. and Aguilar, L. A. 1996, AJ, 112, 668

Chiba, M. and Yoshii, Y. 1998, AJ, 115, 168

Côté, P., Welch, D.L., Fischer, P. and Irwin, M.J. 1993, ApJ, 406, L59

Dohm-Palmer, R. C., Mateo, M., Olszewski, E., Morrison, H., Harding, P., Freeman, K. C., & Norris, J. 2000, AJ, 120, 2496

Dohm-Palmer, R. C. et al. 2001, ApJ, 555, L37

Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748

Harding, P., Morrison, H. L., Olszewski, E. W., Arabadjis, J., Mateo, M., Dohm-Palmer, R. C., Freeman, K. C., & Norris, J. E. 2001, AJ, 122, 1397

Helmi, A. 2000, doctoral dissertation, Leiden University

Helmi, A., White, S. D. M., de Zeeuw, P. T. and Zhao, H. 1999, Nature, 402, 53

Helmi, A. & White, S. D. M. 2001, MNRAS, 323, 529

Ibata, R. A., Gilmore, G. and Irwin, M. J. 1995, MNRAS, 277, 781

Ibata, R. A. & Lewis, G. F. 1998, ApJ, 500, 575

Ibata, R., Irwin, M., Lewis, G. F., & Stolte, A. 2001a, ApJ, 547, L133

Ibata, R., Lewis, G. F., Irwin, M., Totten, E., & Quinn, T. 2001b, ApJ, 551, 294

Ivezić, Željko et al. 2000, AJ, 120, 963

Klypin, A., Kravtsov, A. V., Valenzuela, O. and Prada, F. 1999, ApJ, 522, 82

Layden, A. C. & Sarajedini, A. 2000, AJ, 119, 1760

Majewski, S. R. 1992, ApJS, 78, 87

Moore, B., Ghigna, S., Governato, F., Lake, G., Quinn, T., Stadel, J. and Tozzi, P. 1999, ApJ, 524, L19

Morrison, H. L., Mateo, M., Olszewski, E. W., Harding, P., Dohm-Palmer, R. C., Freeman, K. C., Norris, J. E., & Morita, M. 2000, AJ, 119, 2254

Morrison, H. L., Olszewski, E. W., Mateo, M., Norris, J. E., Harding, P., Dohm-Palmer, R. C., & Freeman, K. C. 2001, AJ, 121, 283

Pearce F.R., Jenkins A., Frenk C.S., Colberg J.M., White S.D.M., Thomas P.A., Couchman H.M.P., Peacock J.A., Efstathiou G. 1999, ApJ, 521, 99

Searle, L. & Zinn, R. 1978, ApJ, 225, 357

Steinmetz M. & Navarro J.F. 1999, ApJ, 513, 555

Yanny, B., Newberg, H. et al. 2000, ApJ, 540, 825

Zhao, H. 1998, ApJ, 500, L149


format: text

text:
# Phenomenological analysis of experimental data on 𝜂-photoproduction on protons11footnote 1The contributed paper to MENU2001, Washington, USA, 26-31 July

## Abstract

The results of linear regresion of differential cross sections, $`\mathrm{\Sigma }`$\- and $`T`$-asymmetries of $`\eta `$-photoproduction on protons in energy region from threshold up to 1 GeV are presented. Serious contradictions between angular distributions measured in different laboratories are revealed. The energy dependance of regression coefficients may be due to transition from energy region of $`S_{11}(1535)`$ and $`D_{13}(1520)`$ to energy region of $`D_{15}(1675)`$ and $`F_{15}(1680)`$ resonances.

During the past years $`\eta `$-photoproduction on protons has attracted increasingly high interest. This is due not only because this item is a new physical phenomenon different from photoproduction of pions but mainly because $`\eta `$-photoproduction should proceed through the small number of nucleon resonances. Even in energy region up to 1 GeV there will be not too many overlapping resonances that permits to extract reliable information on resonance parameters from experimental data.

Complete phenomenological analysis of experimental data on photoproduction, as a rule, encounters a number of problems, e.g. solving of nonlinear equations, removal of continuous and discrete theoretical ambiguities, elimination of experimental ambiguities, etc. Analysis of experimental data may be naturally divided in two stages . The first is the linear regression which provides the information about the number of partial waves that contribute to the measured experimental characteristics of process and provides information on the resonances concerned. The linearity of the model used ensures that the estimates of regression coefficients are unbiased. The second is to determine the multipole amplitudes.

This paper is confined to the first stage of analysis. We have analysed all known experimental data on differential cross sections (angular distributions) of process $`\gamma p\eta p`$ , and also the data on polarization observables, i.e. angular distributions of asymmetry $`\mathrm{\Sigma },`$ measured with linear polarized beam  and angular distributions of asymmetry $`T,`$ measured on a polarized target . The energy independent analysis consist in expanding angular distribution of the observables at definite energy using Legendre polinomials. To find how many terms in this expansion provide the best description of data standart statistical procedures including the Fisher criterion were used. Unlike the energy dependent analysis that is based on parametric models and, generally, gives biased estimates, energy independent analysis relies on nonparametric model that provides unbiased estimates. Expansion of the observables and corresponding statistics are:

$$\frac{k}{q}\frac{d\sigma (\theta )}{d\mathrm{\Omega }}=a_nP_n(\mathrm{cos}\theta ),$$
(1a)
$$\frac{k}{q}\frac{d\sigma }{d\mathrm{\Omega }}\frac{1}{\mathrm{sin}^2\theta }\mathrm{\Sigma }=b_nP_n(\mathrm{cos}\theta ),$$
(1b)
$$\frac{k}{q}\frac{d\sigma }{d\mathrm{\Omega }}\frac{1}{\mathrm{sin}\theta }T=c_nP_n(\mathrm{cos}\theta ).$$
(1c)
Multipole decomposition of coefficients $`a_n,b_n,c_n`$ up to terms $`E_3`$ and $`M_3`$ may be found in . In all cases the best description of experimental data on $`d\sigma /d\mathrm{\Omega },`$ were obtained with three terms of the expansion. The dominance of $`s`$-wave, the coefficient $`a_0,`$ was already pointed out . However, the coefficients $`a_1`$ and $`a_2`$ connected, correspondingly, with the $`sp`$\- and $`sd`$-interferences demonstrate the existence of serious contradictions between results in . This is also diplayed by Fig. 1.

Since the observables $`d\sigma /d\mathrm{\Omega },`$ $`\mathrm{\Sigma },`$ $`T`$ were measured at different energies and angles to form the statistics with $`\mathrm{\Sigma }`$ and $`T`$ we used interpolated values of $`d\sigma /d\mathrm{\Omega }.`$ The polarization statistics $`\mathrm{\Sigma }`$ and $`T`$ were analysed with both $`d\sigma /d\mathrm{\Omega }`$ obtained in the same laboratory and $`d\sigma /d\mathrm{\Omega }`$ from another laboratories. To get the description of $`\mathrm{\Sigma }(\theta )`$ it was necessary to keep three terms in expantion. For the $`T(\theta )`$ it was sufficient to keep two terms. The energy dependence of $`b_n`$ and $`c_n`$ is shown on Fig. 2 and 3.

Contradictions between angular distributions obtaned in different laboratories are not reflected in the general behavior of the coefficients $`b_n`$ and $`c_n.`$ In other words, these contradictions between $`d\sigma /d\mathrm{\Omega }`$ do not appear in polarizational observables.

It seems to be instructive to consider the energy behaviour of coefficients $`a_1,b_1,b_2`$ and $`c_0,c_1.`$ The change of energy dependence of this coefficients at 0.9 GeV might indicate the change of regime of the process. For instance, the decrease of $`a_1`$  from the threshold to 0.9 GeV may be due to the damping of $`s`$-wave and to weakening of $`sp`$-interference. The further rise of $`a_1`$ may be related with the contribution of higher partial waves. The decrease of $`b_2`$ at energies below 0.9 GeV may be related with resonance $`D_{13}(1520)`$; the growth of $`b_2`$ at energies 0.9–1.1 GeV may be due to influence of resonances $`D_{15}(1675)`$ and $`F_{15}(1680).`$ The interference of $`d`$\- and $`f`$-waves should lead to the shift of angular distribution $`\mathrm{\Sigma }(\theta )`$ to the smaller angles in the CM system as really seen in experiment . The behaviour of $`b_1`$ at energies higher than 0.9 GeV can be attributed to $`sf`$-interference and so on.

Thus, energy dependance of regression coefficients found in our analysis may be due to the transition from energy region of $`S_{11}(1535)`$ and $`D_{13}(1520)`$ to energy region of $`D_{15}(1675)`$ and $`F_{15}(1680).`$


text:
# New Insights from a study of the Globular Cluster Systems of 60 Early-Type Galaxies

## 1. Introduction

Globular clusters (GC) are among the most pristine objects in galaxies and hence provide invaluable insights into the chemical and kinematic conditions prevalent during the early epochs of galaxy formation. Just how early is presently a matter of some debate. Initial studies of the GCSs of the Milky Way and other nearby galaxies led to the conclusion that all globular clusters are old, metal-poor systems that were formed during the initial collapse phase of a protogalactic gas cloud. The subsequent discovery of metal-rich GCSs and bimodal color (metallicity) distributions in some galaxies has led to much discussion about alternate mechanisms of globular cluster formation, including at this very symposium. Schweizer (1987) and Ashman $`\&`$ Zepf (1992) proposed that some GCs may be formed during the interaction or merger of galaxies. While Zepf $`\&`$ Ashman (1993) suggest that the bimodality in the color distribution is a natural consequence of the merger scenario, Forbes, Brodie & Grillmair (1997) favor multiple phases of cluster formation during the collapse of a gas cloud into a galaxy.

Only recently has it become possible to test these and other theories of globular cluster formation and evolution, due to the explosion in the number of observations of GCSs of external galaxies. The advent of the HST with its superior angular resolution has given a major boost to the field because of the ease with which cluster candidates can be identified. While most HST-based studies concentrate on a single galaxy, or just a handful, there are few integrated analyses that study a large number of galaxies simultaneously (Forbes et al. 1996; Gebhardt & Kissler-Patig 1999; Brodie & Larsen elsewhere in this volume) which helps eliminate the possible systematic differences that make comparison of the results of different authors a somewhat tricky task. We report here the results from our recent analysis of the GCS of over 60 S0 and elliptical galaxies using WFPC2 V and I images (Kundu & Whitmore 2001a; 2001b). This is the largest and most comprehensive analysis of GCSs to date.

## 2. Color Distributions and Implications on Formation Models

Our analysis of the WFPC2 images of the inner regions of the 29 elliptical galaxies in our sample revealed a measurable GCS in each of them (Kundu & Whitmore 2001a). However, only 30 of the 35 S0s (Kundu & Whitmore 1998; 2001b) showed evidence of a GCS. This may be partly due to the fact that S0s have intrinsically fewer clusters than ellipticals (see §4), and partially because for 34 of the 35 S0 candidates we analyzed were short exposure ’snapshot’ images. The majority of GC candidates lie in a narrow range of color between 0.5$`<`$V-I$`<`$1.5 with a mean color near V-I $``$ 1.0 mag, which is fairly typical for old globular cluster systems. Hence, the color of the cluster is predominantly a function of the metallicity, and throughout this analysis we use the linear color-metallicity relationship for Milky Way GCs (Kundu & Whitmore 1998) to quote metallicities (see Kissler-Patig in this volume for an alternate view).

The mean color of the GCSs of our entire elliptical sample is 1.04$`\pm `$0.04 (0.01) which is indistinguishable from the value of 1.00$`\pm `$0.07 (0.01) for the S0s. 111We use the following convention in quoting uncertainties: The number following the $`\pm `$ sign is the standard deviation, while a number in parentheses refers to the uncertainty in the mean. Therefore, inasmuch as the broad-band colors trace metallicity, we can conclude that the metallicity of the GCS of early type galaxies is not a function of Hubble type. What then drives this metallicity/color spread? It has been suggested that that the mean metallicities of elliptical galaxy GCSs are correlated to the mass of the host. To check the veracity of this claim we plot the mean color of the elliptical and S0 GCSs vs the absolute magnitude (mass) of the host galaxies in Fig 1. We find that \[Fe/H\] generally increases with galaxy luminosity, and that there is no offset between the the S0 and elliptical samples. While there are no low luminosity S0s or ellipticals with high mean metallicities, there appear to be a number of luminous elliptical galaxies that have metallicities significantly lower than that expected for a linear \[Fe/H\]-M$`{}_{}{}^{T}{}_{V}{}^{}`$ relationship. This difference in the mean metallicities of the cluster systems of hosts with similar luminosities may either be due to an offset in the entire metallicity distribution of the two cluster systems, or differing fractions of two or more sub-populations having different metallicities. To this end it is critical to study the detailed color distributions of individual GCSs.

### 2.1. Bimodality in Color Distributions

Ashman & Zepf (1992) predicted bimodality in GCS color distributions based on their merger model, and subsequently observed it (Zepf & Ashman 1993). While bimodality has been confirmed in a large number of elliptical galaxies since then, the frequency of bimodality and it’s correlation with various host properties, which is central to the various competing theories of GCS formation, has not been very well studied.

We searched for bimodality in our galaxy sample, using the KMM mixture modelling algorithm (Ashman, Bird & Zepf 1994), fitting Gaussians to binned data, and visual inspection. The galaxies were sorted into three broad categories, confirmed bimodal, very likely bimodal, and those with no statistically measurable evidence of bimodality. Our analysis revealed that 18 of the 29 ellipticals can be better described by two Gaussian sub-populations than a single Gaussian, of which 9 are confirmed as bimodal without a doubt. The other 11 showed no reasonable partitions for two or more groups. As showed in Fig 2, roughly 60% of ellipticals at all host masses are likely to be bimodal. This number is probably an underestimate as incidences of multiple populations may be difficult to discern due to photometric uncertainties or age-metallicity degeneracies. None of the GCSs in our sample show meaningful partitions for three or more groups.

Due to the short exposure times of the ’snapshot’ S0 images and the smaller number of clusters per galaxy (§4) we were unable to conclusively test for multiple populations of GCs in most of our S0s. Although only 2 of the S0s revealed definitive evidence for bimodality, statistical tests suggest that at least 15-20% of these galaxies are bimodal at the present level of photometric accuracy.

### 2.2. Implications on Formation Models

There are three competing models/scenarios that attempt to explain bimodality in GCS colors. While the merger model (Schweizer 1987; Ashman $`\&`$ Zepf 1992) and the multiple collapse model (Forbes, Brodie, & Grillmair 1997) both suggest that this is the consequence of globular clusters forming during two distinct epochs in the metal enrichment histories of these galaxies, the Côté et al. (1998) model proposes that all cluster systems formed at roughly the same time and that the bimodal distribution is a result of the merger (without creation of new clusters) or cannibalism of a metal-poor GCS by a larger galaxy with a pre-existing metal-rich GCS.

In Fig 3 we plot the mean GCS colors of the galaxies with confirmed bimodal distributions and the ones with confirmed lack of bimodality as a function of host luminosity. Both distributions appear indistinguishable with no evidence that the unimodal galaxies follow any kind of well defined color(metallicity)-luminosity relation. To recall, the Côté et al. (1998) model assumes that the input unimodal galaxies that are cannibalized to form giant ellipticals have a second order luminosity-cluster metallicity relationship. A set of merging/coalescing cluster systems with widely varying metallicity distributions, as seen in Fig 3, is unlikely to produce a bimodal distribution in the absence of such a relation, which suggests that the Côté et al. (1998) model is unviable. Moreover, within the range of host luminosities of our sample we find no unimodal cluster system with a mean GCS color of V-I$``$1.2 (corresponding to the red peak of giant galaxies) which could serve as the progenitor host for the Côté et al. cannibalism model. We suspect that the reason for the shallower input color-luminosity relationship used by Côté et al. (1998) is that they calculated it on the basis of ground-based data. It has been shown that in galaxies with bimodal color distributions the blue clusters are more spatially extended than the red ones (e.g. Geisler, Lee & Kim 1996), so it is very likely that a relationship derived on the basis of ground-based data, that mostly samples the outer regions of a cluster system, returns a shallower color-luminosity relationship. The low luminosity elliptical galaxies present a further problem for the Côté et al model. In order to create bimodal low mass systems (See Fig 2) through the ’cannibalism’ model one would require a population of low mass, unimodal metal-rich systems, which could then accrete other low mass metal-poor systems. There are two problems with this hypothesis: We observe no low mass metal-rich systems in our sample. And such a requirement would violate the smooth second order metallicity(color)-host luminosity relationship required for the formation of bimodal systems in giant galaxies through the Côté et al model. Thus the observed bimodality in low luminosity galaxy could not have been caused by the mergers of two or more systems of different metallicities via the ’cannibalism’ model, and some other mechanism is at work. If such an independent method for creating bimodality exists in low luminosity galaxies, one could turn this argument around and claim that the same process could occur in giant galaxies, and that the accreting low mass galaxies bring bimodal populations with them. While it is very likely that giant galaxies cannibalize dwarfs and their predominantly metal-poor cluster systems to some extent, it is unlikely to be the primary mechanism of creating bimodality in a majority of galaxies. Thus, our observations do not seem to support the Côté et al. (1998) model and suggests that the bimodality in globular cluster colors is evidence for the formation of clusters in two different epochs of the metal-enrichment process. By extension this implies that for old cluster systems the red clusters are younger than the blue ones.

There are no obvious ways to distinguish between the multiple collapse and merger models based on the colors of clusters. They both predict bimodal color distributions with the red clusters situated closer to the center of the galaxy than the blue ones (which is observed in our sample). While the merger model provides a mechanism for the second episode of cluster formation, one of the drawbacks of the multiple collapse model is that at present there is no known viable trigger. It is possible that the second burst of cluster formation in the collapse model is triggered by a merger. The massive gas budget required to form the many thousands of metal-rich GCs in giant ellipticals, coupled with the roughly similar inferred ages of the red and blue clusters observed in galaxies like M87, suggests that for the merger model to work the event must have taken place very early in the history of the progenitors when they were still gaseous. The difference between the merger of two largely gaseous bodies and the collapse of one large gaseous entity (with possible fragments within this body) may largely be one of semantics.

In NGC 3115, the only S0 with deep observations and sufficient number of clusters for independent analysis of the red and blue GCS, we find that the red clusters are associated with the thick disk while the blue ones are members of a spheroidal population. This is consistent with the minor merger scenario of Schweizer (1990). If the metal-poor GCs were formed during the initial collapse phase of the galaxy, the merger of a smaller galaxy (mass ratio $``$0.3) could cause a starburst in the disk that heats the gas and induces star formation in giant molecular clouds, most of which are likely to be associated with the disk of the larger of the two merging galaxies.

Based on our analysis of early type galaxies we conclude that a majority of the GCSs of ellipticals are formed in two episodes. The second episode is likely triggered by a merger; possibly a major merger in the case of cD galaxies like M87 and a minor merger in the case of S0 galaxies like NGC 3115 that preserves the disk component.

## 3. The Globular Cluster Luminosity Function

The turnover luminosity of the globular cluster luminosity function (GCLF) has been found to be remarkably constant over a wide range of galaxies and environments. So constant in fact that it has been used as a secondary distance indicator (Harris 1991, Jacoby 1992). The theoretical basis for this rather remarkable result - which implies that for any reasonable range of M/L ratio the underlying mass distribution of globular clusters is the same in all galaxies - is not well understood.

After detailed corrections for completeness effects, contamination etc. we fitted Gaussians curves to the GCLFs of our sample. While we found weak evidence that the width of the GCLF i.e. the dispersion increases with host mass, the variation in the dispersions is smaller than that reported in the literature. Overall a Gaussian with $`\sigma `$=1.3 mag describes both V and I band GCLFs quite well. The large scatter in previous observations can probably be attributed to the larger uncertainties induced by incompletenes and contamination in ground-based analyses.

In Fig 4 we compare the GCLF turnover luminosities of the GCLFs of the 11 systems with the deepest data in the V and I bands with the weighted distance moduli of three other distance indicators from Ferrarese et al. (2000) and the surface brightness fluctuation distances from Neilsen (1999). We also compare the Neilsen distances to the weighted means. The GCLF turnover, in both the V and I-band tracks the distance to these galaxies very well throughout the entire distance range and appears to be an excellent distance indicator. Comparing our turnover magnitudes with the Ferrarese et al. distances we calculate a turnover luminosity of M$`{}_{}{}^{0}{}_{V}{}^{}`$ = -7.41(0.03) in V, and M$`{}_{}{}^{0}{}_{I}{}^{}`$=-8.46 (0.03) in the I-band. We also find evidence that the difference in the turnover luminosities in V and I increases with metallicity, as predicted by Ashman, Conti & Zepf (1995). In sum we conclude that the GCLF is an excellent distance indicator with an accuracy that is as good as the surface brightness fluctuation method.

## 4. Specific Frequency

The specific frequency (SN) is a measure of the luminosity normalized globular cluster density in a host galaxy, and is typically measured for the entire cluster system. Due to the limited field of view of the WFPC2 we could only measure the local value of SN in the inner regions of our galaxies. Fig 5 plots SN(local) for the S0s and ellipticals with reliable estimates as a function of the host galaxy luminosity. The reduction methods used to analyze the ’snapshot’ S0 images were different from those used for the deeper ellipticals; hence the uncertainties are underestimated for the former by a factor of $``$2. However, it is clear that at all host masses ellipticals are more cluster rich than S0s. The mean SN(Local) for S0s, 1.0$`\pm `$0.6 (0.1) is significantly smaller than the value of SN(Local) = 2.4$`\pm `$1.8 (0.4) derived for ellipticals. Furthermore, there is weak evidence that SN(local) of early type S0s is slightly larger than that of the late types. Thus we conclude that specific frequency is a function of the Hubble type of the host galaxy.

For a more detailed description of our results we refer the reader to the recently published and forthcoming accepted paper Kundu & Whitmore (2001a) and (2001b) respectively.

## References

Ashman, K. M., Bird, C. M., & Zepf, S. E. 1994, AJ, 108, 2348

Ashman, K. M., Conti, A., & Zepf, S. E. 1995, AJ, 110, 1164

Ashman, K. M., & Zepf, S. E. 1992, ApJ, 384, 50

Côté, P., Marzke, R. O., & West, M. J. 1998, ApJ, 501, 544

Ferrarese, L. et al. 2000, ApJ, 529, 745

Forbes, D. A., Brodie, J. P., Grillmair, C. J. 1997, AJ, 113, 1652

Forbes, D. A., Franx, M., Illingworth, G. D., & Carollo, C. M. 1996, ApJ, 467, 126

Gebhardt, K. & Kissler-Patig, M. 1999, AJ, 118, 1526

Geisler, D., Lee, M. G., & Kim, E. 1996, AJ, 111, 1529

Harris, W. E. 1991, ARA&A, 29, 543

Jacoby, G. H. et al. 1992, PASP, 104, 599

Kundu, A., & Whitmore, B. C. 1998, AJ, 116, 2841

Kundu, A., & Whitmore, B. C. 2001a, AJ, 121, 2950

Kundu, A., & Whitmore, B. C. 2001b, AJ, in press

Neilsen, E. H. 1999, PhD Thesis, Johns Hopkins University

Schweizer, F. 1987, in Nearly Normal Galaxies, ed. S. Faber (New York: Springer)

Schweizer, F. 1990, in Dynamics and Interaction of Galaxies, ed. R. Wielen (Berlin: Springer-Verlag)

Zepf, S. E., & Ashman, K. M. 1993, MNRAS, 264, 611


text:
# 1 Introduction

## 1 Introduction

Recently Verlinde have made two interesting things . The related issues appeared in . First he proposed that using the AdS/CFT correspondence , the entropy of a conformal field theory (CFT) in any dimension can be expressed in terms of a generalized form of the Cardy formula . We consider a CFT residing in an $`(n+1)`$-dimensional spacetime with the static metric for the Einstein space

$$ds_{CFT}^2=d\tau ^2+R^2d\mathrm{\Omega }_n^2,$$
(1.1)
where $`d\mathrm{\Omega }_n^2`$ denotes an unit $`n`$-dimensional sphere. Especially for a strongly coupled CFT with the anti de Sitter (AdS) dual, one obtains the Cardy-Verlinde formula which states a relation between entropy ($`S`$) and energy ($`E`$)

$$S=\frac{2\pi R}{n}\sqrt{E_c(2EE_c)}$$
(1.2)
with the Casimir energy ($`E_c`$) for a finite system. Indeed, this formula was checked to hold for various kinds of AdS-bulk spacetimes: AdS-Schwarzschild black holes ; AdS-Kerr black holes ; AdS-charged black holes ; AdS-Taub-Bolt spacetimes .

The other is to connect the above Cardy-Verlinde formula with the $`(n+1)`$-dimensional Friedmann equation based on the Friedman-Robertson-Walker (FRW) metric for the closed universe

$$ds_{FRW}^2=d\tau ^2+^2(t)d\mathrm{\Omega }_n^2.$$
(1.3)
Although two metrices (1.1) and (1.3) have different natures, these are conformally equivalent. Hence it is possible to make connection between these. For a radiation-dominated closed universe, two Friedman equations are given by

$`H^2={\displaystyle \frac{16\pi G_{n+1}}{n(n1)}}{\displaystyle \frac{E}{V}}{\displaystyle \frac{1}{^2}},`$ (1.4)
$`\dot{H}={\displaystyle \frac{8\pi G_{n+1}}{n1}}\left({\displaystyle \frac{E}{V}}+p\right)+{\displaystyle \frac{1}{^2}},`$ (1.5)
where $`H=\dot{}/`$ is the Hubble parameter, the dot stands for the differentiation with respect to the proper time $`\tau `$, $`E`$ is the total energy of matter filling in the universe, $`p`$ denotes the pressure and $`V=^n\text{Vol}(S^n)`$ is the volume of the universe. In addition, $`G_{n+1}`$ is the $`(n+1)`$-dimensional newtonian constant. Verlinde pointed out that the Friedmann equation (1.4) can be expressed in terms of three cosmological entropy bounds:111In Ref. , the first bound is called the Bekenstein bound. In fact, this bound is slightly different from the original Bekenstein entropy bound proposed in . So we call this the Bekenstein-Verlinde bound.

$$\begin{array}{cc}\hfill \text{Bekenstein-Verlinde bound}:& S_{BV}=\frac{2\pi }{n}E,\hfill \\ \hfill \text{Bekenstein-Hawking bound}:& S_{BH}=(n1)\frac{V}{4G_{n+1}},\hfill \\ \hfill \text{Hubble bound}:& S_H=(n1)\frac{HV}{4G_{n+1}}.\hfill \end{array}$$
(1.6)
Then the Friedmann equation (1.4) can be rewritten as the Verlinde’s entropy relation

$$S_H^2+(S_{BV}S_{BH})^2=S_{BV}^2.$$
(1.7)
The above equation can be solved by introducing the conformal time coordinate $`\eta `$ as

$$S_H=S_{BV}\mathrm{sin}\eta ,S_{BH}=S_{BV}(1\mathrm{cos}\eta ).$$
(1.8)
This means that $`S_{BV}`$ is constant with respect to the cosmic time $`\tau `$, while $`S_H`$ and $`S_{BH}`$ depend on the cosmic time. Actually $`S_{BV}`$ is constant throughout the entire evolution, because $`E^1`$ for a radiation-dominated universe. We note that the Bekenstein-Verlinde bound is valid for the weakly self-gravitating universe ($`H1`$), while the Hubble bound holds for the strongly self-gravitating universe ($`H1`$). To decide whether a system is strongly or weakly gravitating, we have to introduce another quantity like $`S_{BH}`$. When $`S_{BV}S_{BH}`$, the system is weakly gravitating, while for $`S_{BV}S_{BH}`$ the self-gravity is strong. This is identified with the holographic Bekenstein-Hawking entropy of a black hole with the size of the universe. It grows like an area instead of the volume. Also the maximal entropy inside the universe is bounded by the black holes of the size of the Hubble horizon. This is the Hubble entropy bound $`S_H`$. It is clear from the Friedmann equation (1.4) that at the critical point of $`H=1`$, three entropy bounds coincide exactly with each other.

Further let us propose $`E_{BH}`$ corresponding to the Bekenstein-Hawking energy by using the Bekenstein-Verlinde bound such a way that $`S_{BH}=(n1)V/4G_{n+1}2\pi E_{BH}/n`$. Equation (1.7) then takes the form

$$S_H=\frac{2\pi }{n}\sqrt{E_{BH}(2EE_{BH})}.$$
(1.9)
It is very important to note here that this relation is the same form of the Cardy-Verlinde formula (1.2) except that the roles of the entropy $`(S)`$ and Casimir energy $`(E_c)`$ are taken over by the Hubble entropy bound $`(S_H)`$ and Bekenstein-Hawking energy $`(E_{BH})`$. This connection between the Cardy-Verlinde formula and the Friedmann equation can be interpreted as a consequence of the holographic principle . This implies that two (Friedman equation and Cardy-Verlinde formula) can be derived from the same first principle.

In this direction, Savonije and Verlinde  have studied a concrete model by using the one-side brane cosmology in the background of $`(n+2)`$-dimensional AdS-Schwarzschild spacetime

$$ds_{AdSS_{n+2}}^2=g_{MN}dx^Mdx^N=h(r)dt^2+\frac{1}{h(r)}dr^2+r^2\left[d\chi ^2+\mathrm{sin}^2\chi (d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)\right],$$
(1.10)
where $`h(r)`$ is given by

$$h(r)=1\frac{m}{r^{n1}}+\frac{r^2}{\mathrm{}^2}$$
(1.11)
with $`\omega _{n+1}M=\frac{16\pi G_{n+2}}{n\mathrm{Vol}(s^n)}M`$. Here $`\mathrm{}`$ is the curvature radius of AdSn+2 space and $`M`$ is the ADM mass of the black hole as measured by an observer who uses $`t`$ as his time coordinate. $`G_{n+2}`$ is the bulk newtonian constant. In the one-side brane world scenario, we have the relation of $`G_{n+2}=\mathrm{}G_{n+1}/(n1)`$ between the bulk and boundary constants. In the case of $`m=0`$, we have an exact AdSn+2-space. However, $`m0`$ generates the electric part of the Weyl tensor $`E_{00}=C_{0N0Q}n^Nn^Qm/r^{n1}`$ . This corresponds to the nonlocal effects arisen from the free gravitational field in the bulk , transmitted through the projection $`E_{MP}`$ of the bulk Weyl tensor. This nonlocal Weyl term will contributes corrections to the Friedmann equations on the brane. Actually we focus on the role of this term in cosmology. It was argued that the energy $`(E)`$, entropy $`(S)`$, and temperature $`(T)`$ of a CFT at high temperature can be related to the mass $`(M)`$, entropy $`(S_{BH})`$, and Hawking temperature $`(T_H)`$ of the AdS black hole. According to the GKPW prescription of the AdS/CFT correspondence , the conformal class of the boundary CFT metric is not fixed. Let us introduce its boundary metric from the bulk one in Eq.(1.10)

$$ds_{BCFT}^2=\underset{r\mathrm{}}{lim}\left[\frac{\mathrm{}^2}{r^2}ds_{AdSS}^2\right]=d\tau ^2+\mathrm{}^2d\mathrm{\Omega }_n^2.$$
(1.12)
From this we deduce the static relation between the quantities on the CFT-boundary and those in the AdS bulk : $`t\tau =t\mathrm{}/r;T_HT=T_H\mathrm{}/r;ME=M\mathrm{}/r.`$ But we have the same entropy $`S=S_{BH}=\frac{r_+^n\mathrm{Vol}(S^n)}{4G_{n+2}},`$ where $`r_+`$ is the event horizon of the AdS black hole. It is well-known that the equation governing the motion of the brane (Moving Domain Wall: MDW) are exactly given by the $`(n+1)`$-dimensional Friedmann equation with radiation matter222 There also exists the other brane cosmology: BDL approach. . In this case the radiation-matter which comes from the nonlocal Weyl term can be identified as a strongly coupled CFT, by making use of the AdS/CFT correspondence. Importantly, it turned out that the Friedmann equation is exactly matched with the Cardy-Verlinde formula for the CFT when the brane crosses the black hole horizon.

On the other hand authors in  pointed out that in general, the Cardy-Verlinde formula is not valid in weakly coupled CFTs. Further we wish to comment that the entropy relation of Eq.(1.7) is suitable only for the restricted case such as a radiation-dominated CFT.

In this article we will clarify again that a deep connection between the Cardy-Verlinde formula and Friedmann equation is just a peculiar property of the dynamic brane moving under the “5D” anti de Sitter Schwarzschild black hole spacetime. Here we provide a counter example where this connection fails for a matter-dominated universe. Actually all moving domain walls in the AdSS5 black hole can always take a kind of radiation-dominated matter $`(\rho E/V,V=^n\mathrm{Vol}(S^n))`$ which arises originally from the nonlocal term $`M/^4`$ of the Schwarzschild black hole through the relation $`E=M\mathrm{}/`$. However, if we consider the brane which moves in the ordinary bulk matter, this connection is no longer satisfied. For example, if one considers the moving domain wall (brane) in the bulk matter with $`\rho _B0,P_B=0`$ which does not include any black hole, then one finds equation of state for the matter-dominated universe of $`\rho _B1/a^3`$ on the brane. There is also a way to obtain a radiation-dominated CFT from the bulk spacetime in the framework of the BDL brane cosmology .

## 2 Brane Cosmology with MDW approach

For definiteness we choose $`n=3`$ (five-dimensional AdS-Schwarzscild balck hole spacetime). Now we introduce the radial location of a MDW in the form of $`r=(\tau ),t=t(\tau )`$ parametrized by the proper time $`\tau `$ to define a cosmic embedding: $`(t,r,\chi ,\theta ,\varphi )(t(\tau ),(\tau ),\chi ,\theta ,\varphi )`$. Then we expect that the induced metric of dynamical domain wall will be given by the FRW-type. Hence $`\tau `$ and $`(\tau )`$ will imply the cosmic time and scale factor of the FRW-universe, respectively. A tangent vector (proper velocity) of this MDW

$$u=\dot{t}\frac{}{t}+\dot{}\frac{}{},$$
(2.1)
is introduced to define this embedding properly. Here overdots mean differentiation with respect to $`\tau `$. This is normalized to satisfy

$$u^Mu^Ng_{MN}=1.$$
(2.2)
Given a tangent vector $`u_M`$, we need a normal 1-form directed toward to the bulk. Here we choose this as

$$n=\dot{}dt\dot{t}d,n_Mn_Ng^{MN}=1.$$
(2.3)
This convention is consistent with the Randall-Sundrum case in the limit of $`m=0`$ . Using Eq.(2.1) either with (2.2) or with Eq.(2.3), we can express the proper time rate of the AdSS5 time $`\dot{t}`$ in terms of $`\dot{}`$ as

$$\dot{t}=\frac{\sqrt{\dot{}^2+h()}}{h()}.$$
(2.4)
From the above, we worry about that $`\dot{t}`$ is not defined at $`=r_+`$ because $`h(r_+)=0`$. This also happens in the study of static black hole. Usually one introduces a tortoise coordinate $`r^{}=h^1𝑑r`$ to resolve it. Then Eq.(1.10) takes a form of $`ds_{AdSS_5}^2=h(dt^2dr^2)\mathrm{}`$ and one finds the Kruscal extension. This means that $`r=r_+`$ is just a coordinate singularity. We confirm this from the computation of $`R_{MNPQ}R^{MNPQ}=40/\mathrm{}^4+72m^2/r^8`$, which shows that $`r=0`$ $`(r=r_+)`$ are true (coordinate) singularities. It is found that there does not exits such a problem even for the dynamic case . Using Eq.(2.4), Eq.(1.10) leads to the 4D induced metric for the brane

$`ds_{FRW}^2`$ $`=d\tau ^2+(\tau )^2\left[d\chi ^2+\mathrm{sin}^2\chi (d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)\right]`$ (2.5)
$`h_{\mu \nu }dx^\mu dx^\nu ,`$
where we use the Greek indices only for the brane. Actually the embedding of the FRW-universe into AdSS5 space is a $`2(t,r)1(\tau )`$-mapping. The projection tensor is given by $`h_{MN}=g_{MN}n_Mn_N`$ and its determinant is zero. Hence its inverse $`h^{MN}`$ cannot be defined. This means that the above embedding belongs to a peculiar mapping to obtain the induced metric $`h_{\mu \nu }`$ in the AdSS5 black hole spacetime $`g_{MN}`$ with $`n_M`$. Now we have to calculate the scale factor $`(\tau )`$ from the Israel junction condition by introducing the extrinsic curvature . For this case the extrinsic curvature is defined by

$`K_{\tau \tau }=K_{MN}u^Mu^N=(h()\dot{t})^1(\ddot{}+h^{}()/2)={\displaystyle \frac{\ddot{}+h^{}()/2}{\sqrt{\dot{}^2+h()}}},`$ (2.6)
$`K_{\chi \chi }=K_{\theta \theta }=K_{\varphi \varphi }=h()\dot{t}=\sqrt{\dot{}^2+h()},`$ (2.7)
where prime stands for derivative with respect to $``$. A localized matter on the brane implies that the extrinsic curvature jumps across the brane. This jump is described by the Israel junction condition for the one-side brane world scenario

$$K_{\mu \nu }=\kappa _5^2\left(T_{\mu \nu }\frac{1}{3}T_\lambda ^\lambda h_{\mu \nu }\right)$$
(2.8)
with $`\kappa ^2=8\pi G_5`$. For cosmological purpose, we may introduce a localized stress-energy tensor on the brane as the 4D perfect fluid

$$T_{\mu \nu }=(\varrho +p)u_\mu u_\nu +ph_{\mu \nu }.$$
(2.9)
Here $`\varrho =\rho _m+\sigma `$ $`(p=p_m\sigma )`$, where $`\rho _m`$ $`(p_m)`$ are the energy density (pressure) of the localized matter and $`\sigma `$ is the brane tension. Here we consider the cosmological evolution without any localized matter on the brane. This means that there exists only the matter from the bulk configuration on the brane. In this case of $`\rho _m=p_m=0`$, the r.h.s. of Eq.(2.8) leads to a form of the RS case as $`\frac{\sigma \kappa ^2}{3}h_{\mu \nu }`$. From Eq.(2.8), one finds the space component of the junction condition

$$\sqrt{h()+\dot{}^2}=\frac{\kappa ^2}{3}\sigma .$$
(2.10)
For the one-side AdSS5 space, we have the brane tension $`\sigma =3/(\kappa ^2\mathrm{})`$ for the fine-tuning. The above equation leads to

$$H^2=\frac{1}{^2}+\frac{m}{^4},$$
(2.11)
where $`m/^4`$ originates from the electric (Coulomb) part of the 5D Weyl tensor, $`E_{00}m/r^2`$ . For $`n=3`$, we have $`m=\frac{16\pi G_5M}{3V(S^3)},M=\frac{}{\mathrm{}}E,V=^3\mathrm{Vol}(S^3),G_5=\frac{\mathrm{}}{2}G_4`$. Then one finds a CFT-radiation dominated universe

$$H^2=\frac{1}{^2}+\frac{8\pi G_4}{3}\rho _{CFT},\rho _{CFT}=\frac{E}{V}.$$
(2.12)
It seems that the equation (2.10) is well-defined even at $`=r_+`$. Thus this leads to $`H=\pm 1/\mathrm{}`$ at the horizon, which was the case mentioned first in ref. . At this moment when the brane crosses the horizon of the AdSS5 black hole, we find that the entropy density $`s=S/V`$ and the temperature of the CFT can be expressed in terms of the Hubble parameter $`H`$ and its time derivative $`\dot{H}`$ only

$$s=\frac{H}{2G_4},T=\frac{\dot{H}}{2\pi H},\mathrm{at}=r_+.$$
(2.13)
Now let us discuss thermodynamics of the CFT itself. Furthermore from the first law of thermodynamics ($`TdS=dE+PdV`$) and the CFT-radiation matter ($`\rho _{CFT}=M\mathrm{}/(V),P_{CFT}=\rho _{CFT}/3`$), we derive

$$\frac{3}{2}(\rho _{CFT}+P_{CFT}sT)=\frac{\gamma }{^2}$$
(2.14)
with

$$T=\frac{1}{4\pi }\left(\frac{4r_+}{\mathrm{}}+\frac{2\mathrm{}}{r_+}\right),\gamma =\frac{3}{8\pi G_4}\frac{r_+^2}{^2}.$$
(2.15)
Here the r.h.s of Eq.(2.14) represents the geometric Casimir part of the energy density. Also from the $`(3+1)`$-dimensional Cardy-Verlinde formula Eq.(1.2), we find

$$s^2=\left(\frac{4\pi }{3}\right)^2\gamma \left(\rho _{CFT}\frac{\gamma }{}\right).$$
(2.16)
Two of Eqs.(2.14) and (2.16) are valid at all times. Let us check what happens for these at the moment when the MDW crosses the horizon. In this case we have $`\gamma =\frac{3}{8\pi G_4}`$. Using this and Eq.(2.13), we can recover the first Friedmann equation (1.4) from Eq.(2.16). Also the second Friedmann equation (1.5) can be derived from (2.14). These imply that the Friedmann equations know about thermodynamics of the CFT.

For a charged background, we find the same result except the appearance of the negative energy density . The similar result was found for the dilatonic black hole background .

## 3 Entropy bound relation revisited

In this section we focus on the $`n=3`$ case. The relation for the entropy bounds Eq.(1.7) corresponds to just an algebraic version of the radiation-dominated FRW equation (1.4). In other words, this relation reflects partly the nature of a kind of newtonian equation :

$$\dot{}^2+V_r()=1$$
(3.1)
with

$$V_r()=\frac{_r^2}{^2},_r^2=\frac{4\pi G_4E}{3\mathrm{V}\mathrm{o}\mathrm{l}(S^3)}.$$
(3.2)
The solution of this differential equation can be solved parametrically in terms of an arc parameter $`\eta `$ ,

$$=_r\mathrm{sin}\eta ,\tau =_r(1\mathrm{cos}\eta ).$$
(3.3)
The range of $`\eta `$ from start of expansion to end of recontraction is just $`\pi `$ and the curve relating radius $``$ to time $`\tau `$ is a semicircle. From Eq.(3.3), we find a relation representing the diagram for cosmic evolution

$$^2+(_r\tau )^2=_r^2$$
(3.4)
which is the same form as in the entropy relation Eq.(1.7). Of course this is valid for the strongly self-gravitating universe with $`H>1`$. In this case we expect a naive correspondence between $`(,\tau ,_r)`$ and $`(S_H,S_{BH},S_{BV})`$:

$$S_H,\tau S_{BH},_rS_{BV}.$$
(3.5)
Also, we observe that for a radiation-dominated universe, the solution (3.3) to the differential equation takes exactly the same form as in Eq.(1.8) to the algebraic equation for the entropy bounds. This property can be regarded as an important factor for establishing the connection between the Cardy-Verlinde formula and the FRW-equation in the radiation-dominated CFT universe. In the next we introduce a counter example where this connection is no longer satisfied.

## 4 Matter-dominated universe

In this section we consider the moving domain wall in the the AdS5 spacetime with the negative cosmological constant $`\mathrm{\Lambda }`$ including the other bulk matter $`(\stackrel{~}{M})`$. We assume that the AdS space does not include any object like black hole. We introduce the matter-dominated Friedmann equation

$$\dot{}^2+V_m()=1$$
(4.1)
with

$$V_m()=\frac{_m}{},_m=\frac{8\pi G_4\stackrel{~}{M}}{3\mathrm{V}\mathrm{o}\mathrm{l}(S^3)}.$$
(4.2)
This situation may be figured out from the brane world scenario . This can be derived from the newtonian cosmology . For example, we consider a spherical ball $`(S^3)`$ on the brane with matter $`m`$ in it. Of course this arises from a kind of the bulk matter $`\rho _B\stackrel{~}{M},p_B=0`$. In this case the potential of the ball on its surface is given by $`\varphi _{ball}(r=(\tau ))=G_4m/`$. Let us introduce a point-like probe with unit mass on the ball. The energy conservation condition for the bound-motion with negative total energy $`k/2`$ is given by

$$\frac{\dot{}^2}{2}+\varphi _{ball}=\frac{k}{2}$$
(4.3)
which can be led to the matter-dominated FRW-equation Eq.(4.1) if one chooses $`m=4\pi \stackrel{~}{M}/3\mathrm{V}\mathrm{o}\mathrm{l}(S^3)`$ and $`k=1`$. However, one handicap of the newtonian cosmology is that this choice of $`m,k`$ is unclear.

The solution to the equation (4.1) can be expressed parametrically in terms of an arc parameter $`\eta `$ ,

$$=\frac{_m}{2}(1\mathrm{cos}\eta ),\tau =\frac{_m}{2}(\eta \mathrm{sin}\eta ).$$
(4.4)
Here the range of $`\eta `$ from start of expansion to end of recontraction is $`2\pi `$ and the curve relating radius $``$ to time $`\tau `$ is not a semicircle. These are different points when comparing with the radiation-dominated case. From Eq.(4.4), we find a relation between $``$ and $`\tau `$ expressed as cycloid, compared with a semicircle for the radiation-dominated universe Eq.(3.4)

$$(_m2)^2+(_m\eta 2\tau )^2=_m^2.$$
(4.5)
This is valid for the weakly self-gravitating universe with $`H<1`$. The newtonian cosmology also belongs to this category. Assuming a naive correspondence between $`(,\tau ,_m)`$ and $`(S_H,S_{BH},S_{BV})`$ as in the radiation-dominated case, then one finds one entropy relation from Eq.(4.5) :

$$(S_{BV}2S_H)^2+(S_{BV}\eta 2S_{BH})^2=S_{BV}^2,$$
(4.6)
where $`S_{BV}`$ is constant.

On the other hand, if one follows closely the definition of each entropy bound, we may propose the new entropy bounds : Bekenstein-Verlinde, Bekenstein-Hawking, and Hubble bounds which may be useful for describing the matter-dominated universe

$$\begin{array}{cc}\hfill \text{Bekenstein-Verlinde bound}:& \stackrel{~}{S}_{BV}=\frac{2\pi }{n}\stackrel{~}{M},\hfill \\ \hfill \text{Bekenstein-Hawking bound}:& \stackrel{~}{S}_{BH}=(n1)\frac{V}{4G_n},\hfill \\ \hfill \text{Hubble bound}:& \stackrel{~}{S}_H=(n1)\frac{HV}{4G_n}.\hfill \end{array}$$
(4.7)
These are the same forms as in Eq.(1.6) except for the definition of the Bekenstein-Verlinde bound. In the radiation-dominated CFT case, $`S_{BV}`$ is constant because of $`E=M\mathrm{}`$(=constant), whereas for the matter-dominated universe, $`\stackrel{~}{S}_{BV}`$ scales as $``$ as the universe evolves. This means that $`\stackrel{~}{S}_{BV}`$ is not constant here. But we suggest a familiar relation between new entropy bounds

$$\stackrel{~}{S}_H^2+(\stackrel{~}{S}_{BV}\stackrel{~}{S}_{BH})^2=\stackrel{~}{S}_{BV}^2.$$
(4.8)
As it stands, this is different from Eq.(1.7). Eq.(4.8) is not obviously an equation for the semicircle because $`\stackrel{~}{S}_{BV}`$ is not constant.

## 5 Conclusions

At this time we do not know exactly which one among Eqs.(4.6) and (4.8) is appropriate for describing the entropy relation for the matter-dominated universe. Also we note that the Verlinde’s entropy relation Eq.(1.7) which derives from the first Friedmann equation Eq.(1.4) for the radiation-dominated case is not suitable for the matter-dominated case. Finally we wish to comment that a close relationship between the Cardy-Verlinde formula and Friedmann equation is realized only for a special circumstance such as a radiation-dominated CFT within the AdS/CFT correspondence through the brane cosmology.

## Acknowledgments

We thank to R.G. Cai and H. W. Lee for helpful discussions. This was supported in part by the Brain Korea 21 Program, Ministry of Education, Project No. D-0025 and KOSEF, Project No. 2000-1-11200-001-3.


text:
# Level Set Approach to Reversible Epitaxial Growth

## I Introduction

At the present time, there is no practical approach to epitaxial growth modelling that bridges the gap between microscopic and macroscopic length scales. Rate equations offer some hope,  but their reliance on uncontrolled mean-field approximations remains a serious obstacle. Kinetic Monte Carlo (KMC) simulations are very popular,  but scale-up to the micron range is very doubtful, even with future supercomputers.

One approach to the multiple-scale problem uses an atomic description in the vertical (growth) direction and a continuum description in the lateral directions. Specifically, the random walk of individual atoms on a flat terrace is replaced by the solution of a diffusion equation for the monomer density on each terrace. This is not a new idea , but its recent rebirth in the context of the level set (LVST) method  is particularly promising in light of the relatively low computational cost needed to treat arbitrarily complicated surface morphologies. So far, good success has been achieved for irreversible epitaxial growth where LVST calculations quantitatively reproduce the results of KMC calculations for the distribution of two-dimensional islands in the submonolayer regime.

The purpose of this paper is to extend the LVST method to the case of reversible epitaxial growth where thermal detachment of atoms from island edges is allowed. This step is necessary if one hopes to produce a model that is relevant to growth at elevated temperatures. Moreover, a reversible LVST growth model has significant computational advantages over a reversible KMC model. This is so because KMC keeps track of every detaching atom, including those that eventually return to the island from whence they came. Such events leave the system unchanged overall  and slow down the simulation significantly. By contrast, the reversible LVST scheme we develop below replaces these events by their time average and so includes only those detachments that do not lead to subsequent reattachments. Moreover, because of the mean field approach, a large number of detachment events can be treated within a single simulation time step.

## II Method

### A Level Sets

The level set method  models the time evolution of arbitrarily shaped objects in $`n`$ dimensions that can undergo topological changes. In this paper, the relevant objects are two-dimensional islands and topological changes occur due to nucleation, dissociation, and coalescence. The key idea is to represent a curve or interface $`\mathrm{\Gamma }`$ in $`\mathrm{}^n`$ by the level $`k`$ of a function $`\varphi (𝐱,t)`$

$`\mathrm{\Gamma }_k`$ $`=`$ $`\{𝐱:\varphi (𝐱,t)=k\}𝐱\mathrm{}^n.`$ (1)
Here, $`\mathrm{\Gamma }_k`$ is the set of closed curves that constitute the perimeters of the islands with height $`(k+1)a`$ ($`a`$ is the lattice parameter).

Fig. 1 illustrates the level set description of a typical epitaxial growth scenario. The left panel is a side view of two islands on a terrace (a) that grow to a precoalesence state (b) and subsequently merge (c). Later, a new island nucleates on top (d). The right panel shows the corresponding level set functions $`\varphi `$. Note that it is not $`\varphi `$ that represents the surface morphology, but only the level sets ($`\varphi `$ = 0 and $`\varphi `$ = 1).

The motion of $`\mathrm{\Gamma }`$ is partly deterministic and partly stochastic. There is a deterministic part because a mean field theory is sufficient to model the time-average of many of the physical processes that contribute to growth. These effects are built into a velocity function $`v_n(𝐱,t)`$ that evolves the function $`\varphi (𝐱,t)`$ in time according to the partial differential equation

$`\dot{\varphi }(𝐱,t)+v_n(𝐱,t)|\varphi (𝐱,t)|`$ $`=`$ $`0.`$ (2)
As the notation suggests, $`v_n`$ is the component of the growth velocity in the direction of the local surface normal $`𝐧=\varphi /|\varphi |`$. The stochastic motion of $`\mathrm{\Gamma }`$ is associated with nucleation events and small-island dynamics. There is no unique algorithm to incorporate these effects into $`\dot{\varphi }(𝐱,t)`$. The particular choice we make is explained in detail below.

### B Deterministic Evolution

It is convenient to write the velocity $`v_n`$ in Eq. (2) in the form

$`v_n`$ $`=`$ $`v^{\mathrm{att}}v^{\mathrm{det}}`$ (3)
where $`v^{\mathrm{att}}`$ accounts for attachment processes that grow islands and $`v^{\mathrm{det}}`$ accounts for detachment processes that shrink islands. The first of these is proportional to the diffusive flux of atoms that approach an island edge from its bounding terraces. Therefore, if $`D`$ is the surface diffusion constant and $`\rho (𝐱,t)`$ is the adatom density, mass conservation gives

$`v^{\mathrm{att}}=a^2D\left({\displaystyle \frac{\rho }{n}}|_{\text{terrace}}{\displaystyle \frac{\rho }{n}}|_{\text{top of island}}\right).`$ (4)
We compute the required density from the mean-field, driven, diffusion equation

$`\dot{\rho }(𝐱,t)`$ $`=`$ $`D^2\rho (𝐱,t)+F2{\displaystyle \frac{dN_{\mathrm{nuc}}}{dt}}.`$ (5)
The loss term in Eq. (5),

$`{\displaystyle \frac{dN_{\mathrm{nuc}}(t)}{dt}}`$ $`=`$ $`D\sigma _1{\displaystyle _\mathrm{\Omega }}\rho (𝐱,t)^2d^2x,`$ (6)
accounts for the dimers that nucleate as a result of binary collision between monomers. The total simulation area is $`\mathrm{\Omega }`$ and $`\sigma _1`$ is the so-called ”capture number” for an adatom. We solve Eq. (5) subject to the boundary condition that $`\rho =0`$ at every point on $`\mathrm{\Gamma }`$. This differs from the boundary condition usually used for reversible aggregation  because we have elected to incorporate all detachment effects into the velocity $`v^{\mathrm{det}}`$ (see Appendix A).

To find an explicit expression for $`v^{\mathrm{det}}`$, we note first that most particles that detach from an island are driven back to that island by the diffusion field.  Our interest here is those particles that detach without subsequent reattachment, i.e., those that escape from the ”capture zone” of the island. This is so because, by definition, the adatoms in the capture zone of a given island are guaranteed to attach to that island eventually. A relevant quantity is thus $`p_{\mathrm{esc}}`$, the probability that a detached particle reaches the border of the capture zone. The positions of these borders are easily calculated (the locus of points where $`\rho =0`$) and so is $`p_{\mathrm{esc}}`$. We find (Appendix B) that

$`p_{\mathrm{esc}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{log}((R_{\mathrm{is}}+a)/R_{\mathrm{is}})}{\mathrm{log}(R_{\mathrm{cz}}/R_{\mathrm{is}})}}`$ (7)
where $`R_{\mathrm{is}}`$ and $`R_{\mathrm{cz}}`$ are the radii of the circularly averaged island and capture zone.

We now define an effective escape rate per unit length of island perimeter as

$`R_{\mathrm{det}}`$ $`=`$ $`D_{\mathrm{det}}p_{\mathrm{esc}}\lambda `$ (8)
where $`D_{\mathrm{det}}`$ is an effective detachment rate and $`\lambda `$ is the linear density of detaching particles (singly coordinated edge atoms). We use $`\lambda =2`$ for “small” islands and the expression

$`\lambda `$ $`=`$ $`{\displaystyle \frac{16}{3}}\left({\displaystyle \frac{16}{15}}v^{\mathrm{att}}{\displaystyle \frac{1}{D_{\mathrm{edg}}}}\right)^{\frac{2}{3}}`$ (9)
for “large” islands. The distinction between “large” and “small” islands will be made clear below. $`D_{\mathrm{edg}}`$ is the edge diffusion constant. From $`R_{\mathrm{det}}`$ we get the desired expression for the detachment velocity that enters Eq. (3):

$`𝐯^{\mathrm{det}}`$ $`=`$ $`a^2R_{\mathrm{det}}=a^2D_{\mathrm{det}}p_{\mathrm{esc}}\lambda `$ (10)
Significantly (see below), there is negligible extra computational overhead needed to incorporate detachment in this way.

It remains only to find a home for the atoms that escape from the islands. In the spirit of the mean-field approximation, we return them to the adatom pool by simply augmenting the external flux. That is, the variable $`F`$ in Eq. (5) is

$`F`$ $`=`$ $`F_0+F_{\mathrm{rev}}`$ (11)
where $`F_0`$ is the deposition flux and

$`F_{\mathrm{rev}}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }}}{\displaystyle _\mathrm{\Gamma }}R_{\mathrm{det}}𝑑\mathrm{\Gamma }`$ (12)
is the escape rate of atoms from all island edges. In this integral, $`\mathrm{\Gamma }`$ runs over all level sets of $`\varphi (𝐱,t)`$ (see Eq. (1)).

### C Stochastic Evolution

A nucleation event occurs when the variable $`N_{\mathrm{nuc}}`$ in Eq. (6) becomes larger than the next integer. This implies that $`\varphi (𝐱,t)`$ increases by a discrete amount at a discrete point. In the interest of numerical stabilty, we smooth out this increase over several points on the numerical gid used to solve Eq. (5). The exact position where the dimer nucleates is chosen randomly with the integrand of Eq. (6) as a weight factor.

Randomness is also important for detachment from “small” islands that consist of only a few atoms. Our approach is to choose an island area $`A_{\mathrm{cut}}`$ and treat all islands smaller than this size statistically. Thus, in a given time interval, we use Eq. (8) to calculate the total number of adatoms that detach from all islands smaller than $`A_{\mathrm{cut}}`$. This corresponds to a total area-loss $`A_{\mathrm{loss}}`$. If $`A_{\mathrm{loss}}>A_1`$, the area occupied by one atom, we detach an atom from one of the islands smaller than $`A_{\mathrm{cut}}`$. The specific island that loses an atom is chosen randomly with $`p_{\mathrm{esc}}`$ as a weight factor. We then decrement $`A_{\mathrm{loss}}`$ by $`A_1`$ and repeat the process until $`A_{\mathrm{loss}}<A_1`$. This value is stored and added to the loss that occurs in the next time interval. If a detachment process leads to an island smaller than a dimer, we dissociate the dimer and decrement $`A_{\mathrm{loss}}`$ accordingly.

We emphasize that $`A_{\mathrm{cut}}`$ is not the area of a critical nucleus, i.e., an island that is absolutely stable against breakup. Instead, $`A_{\mathrm{cut}}`$ is merely the size of the smallest “large” island which we use as a parameter to switch between statistical and continuous detachment. In the context of our approach, the critical nucleus is defined by the condition $`v^{\mathrm{att}}=v^{\mathrm{det}}`$ (see Eq. (3)).

## III Results

### A Parameter Choices & Systematics

All the LVST simulations reported here use the value $`D/F_0=10^6`$. $`D_{\mathrm{edg}}`$ in Eq. (9) was chosen as 104 and calculations where carried out on a 200$`\times `$200 lattice represented on a numerical grid of 568$`\times `$568 points. To determine $`A_{\mathrm{cut}}`$, we compared runs with different choices for this parameter and looked for stabilization of the physical results. Thus, Fig 2 compares island and adatom densities as a function of coverage $`\mathrm{\Theta }`$ obtained for $`A_{\mathrm{cut}}`$ = 4,6, and 8 at a detachment rate of $`D_{\mathrm{det}}`$/$`D`$ = 0.001. Based on this data and related statistical tests, we find that our results are independent of $`A_{\mathrm{cut}}`$ if $`A_{\mathrm{cut}}6`$. Therefore we use $`A_{\mathrm{cut}}=6`$ in all subsequent simulations.

Fig. 3 illustrates the role of detachment during growth. It shows $`F_{\mathrm{rev}}`$ from Eq. (11) as a function of $`F_0t`$ for different values of $`D_{\mathrm{det}}`$/$`D`$. The fact that $`F_{\mathrm{rev}}`$ has its highest values at the earliest stages of growth (for the higher values of $`D_{\mathrm{det}}`$/$`D`$) shows that most of the effective detachment takes place when the islands are quite small. Indeed, as islands grow bigger, fewer particles escape from the island boundaries since the number of escaping particles relative to the perimeter of the growing islands becomes smaller and smaller.

### B Comparison to KMC

As a critical test, we compared our LVST results directly with KMC simulations. In our KMC simulations adatoms are allowd to hop to a nearest-neighbor site at a rate $`r_n=D\mathrm{exp}(nE_N/k_\mathrm{B}T)`$, where $`n`$ is the number of nearest-neighbours, $`k_\mathrm{B}`$ the Boltzmann-constant, and $`T`$ the temperature. Adatoms are deposited at a rate $`F_0`$. In all simulations presented in this paper the ratio $`D`$/$`F_0`$ is set to $`10^6`$. The energy barrier $`E_N`$ is chosen such that $`D_{\mathrm{det}}=D\mathrm{exp}(E_N/k_\mathrm{B}T)`$. In addition, singly coordinated edge-atoms are allowed to diffuse along the step edge at a rate $`D_{\mathrm{edg}}`$. We chose $`D_{\mathrm{edg}}`$=$`D_{\mathrm{det}}`$, but the results were not sensitive to this parameter as long as the islands were compact (as our LVST model assumes).

Fig. 4, Fig. 5, and Fig. 6 respectively show the adatom density, island density, and island size distributon obtained by both methods as a function of detachment rate. These curves agree with previous reversible KMC simulations  that discuss, e.g., the physical origin of the observed saturation of the island densities and the sharpening of the island size distributions. Evidently, there is semi-quantitative agreement between LVST and KMC. The size distributions results are particularly notable because they reflect information about spatial correlations that are averaged over to get $`\rho `$ and $`N`$.

The only disagreement we find between LVST and KMC is for the saturation value of the island density for the smallest value of the detachment rate. We understand this based on recent research  with irreversible growth where a corresponding disagreement vanishes when the $`\rho =0`$ boundary condition for Eq. (5) is applied not at the true perimeter of each island (as we do) but instead at a closed boundary that exceeds the perimeter everywhere by one lattice constant. The disagreement disappears altogether when the island density is small. This is consistent with our observations because we get good agreement with KMC when the detachment rate is large (and the island density is small).

Fig. 7 shows some systematic features of reversible growth as a function of the effective detachment rate. The LVST data (collected at 0.25 ML coverage) show that the average adatom density increases linearly with $`D_{\mathrm{det}}`$ while the average island density decreases exponentially with $`D_{\mathrm{det}}`$. Unfortunately, we have been unable to derive theses interesting results analytically using rate equations.

The data in Fig. 6 was easy to obtain because, beginning with an irreversible growth simulation, the extra computational cost to include detachment is very small for LVST compared to KMC. This is so because LVST precisely suppresses the time-consuming detachment/attachment fluctuations that occupy a KMC simulation. Moreover, the LVST has essentially no restriction on the number of detachment events simulated during a specific time step. Quantitatively, Fig. 8 shows that the LVST method requires only negligibly more run time to include detachment whereas the KMC simulation cost increases sharply as the detachment rate increases. The LVST results depend very weakly on the rate of detachment because the increased cost is associated wholly with reductions in the step advance rate. Specifically, the adatom density rises with increasing detachment rate so the gradients evaluated in Eq. (4) become larger. But the step advance rate is limited by the condition that the boundary of the level set function can advance only one grid point in each simulation step. This implies that the scaling should be even better for simulations of, say, annealing processes, where there is no deposition flux and the adatom density is very low.

## IV Conclusion

In summary, we have developed a method to model epitaxial growth including atomic detachment from island edges within the context of the level set method. By all reasonable measures, the results are in excellent agreement with KMC simulations. Moreover, the LVST simulations scale significantly in CPU-time demand than KMC simulations when the effective detachment rate is large. This is so because our mean field method eliminates the many atomic detachment events (each processed separately in KMC) that do not lead to successful escape from an island.

## V Acknowledgments

We thank M. Kang, C. Andersen, S. Osher, and D. Vvedensky for helpful discussions. M.P. and A.Z. gratefully acknowledge the support of NSF DMR-953-1115. This work was also supported by DARPA through cooperative agreement DMS-9615854 as part of the Virtual Integrated Prototyping Initiative.

## A Boundary Conditions

The usual boundary condition for reversible growth  sets $`\rho =\rho _{\mathrm{eq}}`$ at every island perimeter. In this Appendix, we relate $`\rho _{\mathrm{eq}}`$ to the detachment velocity $`v^{\mathrm{det}}`$ used in this paper. The idea is to consider two adatom densities: $`\overline{\rho }`$ and $`\rho `$ as schematically shown in Fig. 9. $`\overline{\rho }`$ ($`\rho `$) is the system where the adatom density drops to $`\rho _{\mathrm{eq}}`$ ($`\rho `$=0) at the island boundary. We then write down diffusion for both adatom densities and derive from those the respective reversible growth velocities.

Let $`\overline{\rho }`$ (the solid line in Fig. 9) be the exact solution of

$`0`$ $`=`$ $`D^2\overline{\rho }(𝐱)+F_0𝐱\mathrm{\Omega }`$ (A1)
$`\overline{\rho }(𝐱)`$ $`=`$ $`\rho _{\mathrm{eq}}𝐱B(t)`$ (A2)
where $`\mathrm{\Omega }`$ is the domain and $`B(t)`$ the island boundaries. $`\rho _{\mathrm{eq}}`$ is the (finite) equilibrium value of $`\overline{\rho }(𝐱)`$ at the island boundaries due to detachment of atoms.

If $`C(t)`$ is the capture zone and $`C(t)`$ its boundary (dashed line in Fig. 9), then with the boundary condition

$`{\displaystyle \frac{\overline{\rho }}{n}}=0𝐱C(t).`$ (A3)
we can uniquely solve the diffusion equation (A1).

In this case the reversible growth velocity is given by Eq. (4) (labeled $`𝐯_\mathrm{I}`$ in Fig. 9).

Now we want to replace $`\overline{\rho }`$ by a adatom density $`\rho `$ (dashed line in Fig. 9) such that

$`0`$ $`=`$ $`D^2\rho (𝐱)+F_0+F_{\mathrm{rev}}𝐱\mathrm{\Omega }`$ (A4)
$`\rho (𝐱)`$ $`=`$ $`0𝐱B(t)`$ (A5)
were $`F_{\mathrm{rev}}`$ is as explained in Eq. (11) the diffusive flux away from the island boundary due to detachment.

It follows from (A1) and (A4) that

$`^2\rho (𝐱)`$ $`=`$ $`{\displaystyle \frac{F_0+F_{\mathrm{rev}}}{F_0}}^2\overline{\rho }(𝐱)`$ (A6)
and therefore

$`\rho (𝐱)`$ $`=`$ $`{\displaystyle \frac{F_0+F_{\mathrm{rev}}}{F_0}}(\overline{\rho }(𝐱)\rho _{\mathrm{eq}}).`$ (A7)
If we want the two systems $`\overline{\rho }`$ and $`\rho `$ to be equivalent, then we must require

$`{\displaystyle _{C(t)}}\rho (𝐱)d^2x`$ $`=`$ $`{\displaystyle _{C(t)}}\overline{\rho }(𝐱)d^2x`$ (A8)
Combining (A7) and (A8) we obtain

$`{\displaystyle \frac{F_0+F_{\mathrm{rev}}}{F_0}}{\displaystyle _{C(t)}}\left[\overline{\rho }(𝐱)\rho _{\mathrm{eq}}\right]d^2x`$ $`=`$ $`{\displaystyle _{C(t)}}\overline{\rho }(𝐱)d^2x`$ (A9)
which is equivalent to

$`F_{\mathrm{rev}}{\displaystyle _{C(t)}}\overline{\rho }(𝐱)d^2x`$ $`=`$ $`(F_0+F_{\mathrm{rev}}){\displaystyle _{C(t)}}\rho _{\mathrm{eq}}d^2x.`$ (A10)
If $`A(C)`$ is the area of the capture zone, the relation between $`F_{\mathrm{rev}}`$ and $`\rho _{\mathrm{eq}}`$ is

$`F_{\mathrm{rev}}`$ $`=`$ $`{\displaystyle \frac{F_0}{\frac{1}{A(C)}_{C(t)}\frac{\overline{\rho }(𝐱)}{\rho _{\mathrm{eq}}}d^2x1}}.`$ (A11)
Assuming that $`R_{\mathrm{det}}`$ is independent of $`d\mathrm{\Gamma }`$ we can rewrite Eq. (12) as

$`F_{\mathrm{rev}}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }}}R_{\mathrm{det}}L.`$ (A12)
Then, Eq. (10) reduces to

$`v^{\mathrm{det}}`$ $`=`$ $`a^2{\displaystyle \frac{\mathrm{\Omega }}{L}}F_{\mathrm{rev}}`$ (A13)
$`=`$ $`a^2{\displaystyle \frac{\mathrm{\Omega }}{L}}{\displaystyle \frac{F_0}{\frac{1}{A(C)}_{C(t)}\frac{\overline{\rho }(𝐱)}{\rho _{\mathrm{eq}}}d^2x1}}.`$ (A14)
## B Derivation of $`p_{\mathrm{esc}}`$

Any atom that detaches from an island but does not reach the capture zone boundary is driven by the diffusion field back to the island. In such a case, $`p_{\mathrm{esc}}`$ is zero. But if the detached adatom does reach the border of the capture zone, $`p_{\mathrm{esc}}`$ is one. Therefore, to determine the probability distribution $`P(𝐱)`$ for all possible paths of an detaching adatom, we solve the diffusion equation

$`^2P(𝐱)`$ $`=`$ $`0`$ (B1)
with the boundary condition

$`P(𝐱)`$ $`=`$ $`\{\begin{array}{cc}0:\hfill & 𝐱S_{\mathrm{is}}\hfill \\ 1:\hfill & 𝐱S_{\mathrm{cz}}\hfill \end{array}`$ (B4)
where $`S_{\mathrm{is}}`$ and $`S_{\mathrm{cz}}`$ are the border of the island and of the capture zone respectively. These boundary conditions represent the either successfully escape (probability 1 to reach $`S_{\mathrm{cz}}`$) or the reattachment to the island (0 at $`S_{\mathrm{is}}`$).

Eq. (B1) and (B4) can in principle be solved for any island and capture zone geometry. But in our calculations we assume for simplicity a spherical average for both, the islands and the capture zones. For this case the general solution to Eq. (B1,B4) is

$`P(𝐱)={\displaystyle \frac{\mathrm{log}(|𝐱|/R_{\mathrm{is}})}{\mathrm{log}(R_{\mathrm{cz}}/R_{\mathrm{is}})}}`$ (B5)
where $`R_{\mathrm{is}}`$ and $`R_{\mathrm{cz}}`$ are the radii of the island and of the capture zone. From this we obtain the escape probability of a detached atom by taking $`|𝐱|=R_{\mathrm{is}}+a`$ ($`a`$ is the lattice parameter) because when an adatom detaches, it is roughly at distance $`a`$ from the island boundary. Therefore

$`p_{\mathrm{esc}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{log}((R_{\mathrm{is}}+a)/R_{\mathrm{is}})}{\mathrm{log}(R_{\mathrm{cz}}/R_{\mathrm{is}})}}.`$ (B6)