text:
# Cygnus OB2 – a young globular cluster in the Milky Way

## 1 Introduction

The study of galactic OB associations provides the key to a number of astrophysical questions, such as the star formation process and efficiency, the interaction of massive stars with the interstellar medium, the characterisation of the initial mass function at the high-mass end, the study of stellar nucleosynthesis, chemical evolution, and galactic recycling processes, and the evolution of binary systems. The Cygnus OB2 association is a particularly good region to address such questions, since it is extremely rich (e.g. Reddish et al. 1966, hereafter RLP), and contains some of the most luminous stars known in our Galaxy (e.g. Torres-Dodgen et al. 1991).

The most comprehensive study of the size and shape of Cyg OB2 has been performed by RLP who inferred an elliptical shape with major and minor axes of $`48^{}`$ and $`28^{}`$, respectively (see also Fig. 5). They estimate more than $`3000`$ members of which at least $`300`$ are of OB spectral type, resulting in a total stellar mass between $`(0.62.7)\times 10^4`$ $`M_{\mathrm{}}`$. For their analysis, RLP performed star counts on the blue and red plates of the Digitized Sky Survey (DSS), reaching limiting magnitudes around 20m. Although this limit assures a reasonable complete census for unobscured associations, the extreme reddening in and around Cyg OB2 hampers the detection of even OB stars. With an estimated distance of 1.7 kpc (e.g. Massey & Thompson 1991) Cyg OB2 is located behind the Great Cygnus Rift, causing visual extinction $`A_V`$ from 4m to at least 10m. A number of observations suggest that Cyg OB2 could indeed be larger than the RLP estimate, and that the observed morphology is rather an artifact of the particular extinction pattern in the field. The association boundary determined by RLP fits suspiciously well in a region of low CO column density (cf. Fig. 1) and low visual extinction (Dickel & Wendker 1978), indicating that the visual star densities are probably biased by the extinction pattern. There is a considerable number of early-type stars in the obscured area south and south-east of Cyg OB2 that are estimated to lie at the same distance as Cyg OB2, and that could indeed be bright member stars of the association. Examples are the Wolf-Rayet stars WR 145 and WR 146 (Niemela et al. 1998), the potential Luminous Blue Variable star G79.29+0.46 (Higgs et al. 1994), the massive binary system MWC 349 (Cohen et al. 1985), or the recently discovered group of massive stars around the H II region DR 18 (Comerón & Torra 1999).

The availability of the Two Micron All Sky Survey (2MASS) provides now an excellent opportunity to re-address the question on the morphology and stellar content of Cyg OB2. This survey covers the infrared bands $`J`$, $`H`$, and $`K`$ which have proven to be an excellent tool for unveiling embedded star clusters due to the reduced impact of dust extinction at longer wavelengths. In the following I will use these data to determine the morphology and stellar content of Cyg OB2. It will turn out that the association is much larger and much richer than previously thought, making it the most massive young stellar association known in our Galaxy.

## 2 Sample selection and colour-magnitude diagrams

The analysis presented in this work is based on the second incremental release of the 2MASS point source catalogue (PSC) that contains positional and photometric data for 162 million objects in the $`J`$, $`H`$, and $`K`$ bands. For the investigation of Cyg OB2, all stars within a square field of $`2^{}45^{}`$ in size, centred on $`\alpha =20^\mathrm{h}33^\mathrm{m}6^\mathrm{s}`$ and $`\delta =+41^{}22.5^{}`$ (J2000)111All coordinates given in this paper are for the epoch J2000. have been extracted from the catalogue, resulting in a total of 203318 selected objects. Differential star count distributions for these stars in the 3 photometric bands are shown in Fig. 2. The positions of the turn-over in these distributions suggest that the selected sample is complete to magnitudes of 16.6, 15.5, and 14.8 in the $`J`$, $`H`$, and $`K`$ bands, respectively. This is 0.5-0.8 magnitudes fainter than the nominal completeness limits of the PSC, as expected for less crowded areas in the galactic plane (Cutri et al. 2000). However, the limits have to be understood as an average for the investigated field, and enhanced crowding towards the centre of Cyg OB2 may lower the actual completeness magnitude for the association.

Colour-magnitude diagrams (CMDs) of the selected stars are presented in Fig. 3 as grey-scale plots, which have been obtained by counting the number of stars within pixels of 0.1 magnitudes in $`JK`$ and $`K`$. For illustration, the sample has been divided into two distinct subsets:

* the Cyg OB2 sample, containing all stars within a circular region of $`1.05^{}`$ in radius, centred on $`\alpha =20^\mathrm{h}33^\mathrm{m}10^\mathrm{s}`$ and $`\delta =+41^{}12^{}`$
* the field star sample, containing all remaining stars outside the circular region.

The choice of the separation was motivated by the location and size of the Cyg OB2 association (cf. Sect. 4). Geometrically, the Cyg OB2 sample covers a field of 3.46 degrees squared while the field star sample covers 4.10 degrees squared. The number of stars in the CMDs is 100754 for the Cyg OB2 and 102167 for the field stars sample; 397 stars lie outside the plotted range. Normalising the number of stars in the field star sample to the geometrical area of the Cyg OB2 sample and subtraction from the Cyg OB2 sample results in an excess of 14420 stars within the association field. The CMD of these excess stars is shown in Fig. 4. This figure has been obtained by subtracting the field star CMD, multiplied by 0.85 to account for the different geometrical areas, from the Cyg OB2 sample CMD. To reduce statistical fluctuations the resulting CMD was smoothed with a boxcar average of 0.3 magnitudes in $`JK`$ and $`K`$.

Figure 4 clearly illustrates the wide spread of visual extinctions in Cyg OB2, covering the range from $`A_V5`$ to $`20`$ magnitudes. Between $`A_V510`$ magnitudes, a clear main sequence distribution is visible which spatially coincides with regions of low visual obscuration, as depicted by the light areas in the CO map (cf. Fig. 1). Stars more heavily reddened are related to regions of higher CO column density, located in particular south of the classical Cyg OB2 boundaries.

In the CMDs the PSC sample completeness limits manifest as a colour dependent $`K`$ magnitude limit. For $`JK<1.8`$, the $`K`$ magnitude limit of 14.8 is driving the sample completeness, while for $`JK>1.8`$ the $`J`$ magnitude limit of 16.6 becomes dominant, resulting in a colour dependent completeness limit of $`K=16.6(JK)`$. The resulting limit is shown in Figs. 3 and 4 as solid line. For illustration, the unreddened main sequence for O to F type stars with a distance modulus of 11.2, as suggested by Massey & Thompson (1991) for Cyg OB2, is also shown (see Appendix A for a reference of the photometric calibrations). Reddened main sequences are also shown for $`A_V=5`$ to $`30`$ in steps of 5 magnitudes, assuming $`A_V/A_K=10`$ and a reddening slope of $`A_K/E(JK)`$ of 0.66 (Rieke & Lebofsky 1985).

Comparison of the reddened main sequences with the completeness limit suggests that the PSC is complete for OB stars in Cyg OB2 up to visual extinctions of $`A_V20`$ magnitudes. This is comparable to the maximum reddening found for the association stars (cf. Fig. 4), hence the OB stars sample in the PSC should be complete.

## 3 Morphology from star counts

To study the morphology of the Cyg OB2 association, the stellar density distribution of the field has been determined by binning the PSC stars in pixels of size $`30^{\prime \prime }\times 30^{\prime \prime }`$. Only stars with $`K`$ magnitude brighter than $`13.2+0.66\times (JK)`$ have been selected in order to minimise any bias through variable extinction, and to maximise at the same time the number of stars in the stellar density distribution. This selection is equivalent to choosing stars intrinsically brighter than $`M_K2.2^\mathrm{m}`$ at the distance of Cyg OB2, corresponding to main sequence stars of spectral type F3V and earlier (cf. Appendix A). A comparison to the completeness limit in Fig. 4 shows that stars later than spectral type B9 and suffering from visual extinction in excess of $`A_V12^\mathrm{m}`$ could be missed in the sample, but the CMD of Cyg OB2 suggests that only few association stars show these properties. Hence, the stellar density distribution should be relatively unaffected by absorption.

The presence of several bright stars in the field results in a couple of blind areas in the PSC where no information is available. These blind areas were filled by estimating their stellar densities from adjacent fields. A small ($`10\%`$) gradient in the field star density distribution has been removed by subtracting a 2-dimensional linear function that was estimated by fitting a lower envelope to the data. The resulting stellar density distribution is shown in the left panel of Fig. 5. For display, the distribution has been smoothed with a boxcar average of $`15^{}`$ to reduce statistical density fluctuations and to emphasise the global morphology of the association. Still, most of the fine structure in the image is probably of statistical nature and should be interpreted with care.

The image reveals a rather regular and almost circular stellar density profile for Cyg OB2 with a pronounced maximum at $`\alpha =20^\mathrm{h}33^\mathrm{m}10^\mathrm{s}`$ and $`\delta =+41^{}15.7^{}`$. The maximum is slightly offset from the centre, which is determined to $`\alpha =20^\mathrm{h}33^\mathrm{m}10^\mathrm{s}`$ and $`\delta =+41^{}12^{}`$ (see Sect. 4). Star counts in concentric radial annuli around the centre show that $`50\%`$ of the members are located within a radius of $`21^{}`$, and $`90\%`$ within a radius of $`45^{}`$ around the centre. The association merges with the field stars at a radius of $`1^{}`$. Assuming a total diameter of $`2^{}`$ results in a physical diameter of 60 pc at a distance of 1.7 kpc.

A slight density enhancement at $`\alpha =20^\mathrm{h}32^\mathrm{m}12^\mathrm{s}`$ and $`\delta =+40^{}21^{}`$ at the southern edge of Cyg OB2 possibly indicates some substructure in the association. The stars in this area show $`JK`$ colours $`>2^\mathrm{m}`$, and are compatible with heavily reddened ($`A_V>12^\mathrm{m}`$) main sequence objects. The feature coincides spatially with the H II region DR 15 which has been identified as a nursery of embedded B stars (Odenwald et al. 1990). However, the distance of DR 15 is quite uncertain ($`14.2`$ kpc; Wendker et al. 1991), making the physical association to Cyg OB2 highly speculative.

Two other features of localised stellar density enhancements are visible in Fig. 5. At $`\alpha =20^\mathrm{h}35^\mathrm{m}30^\mathrm{s}`$ and $`\delta =+42^{}25^{}`$ a group of red ($`JK>2^\mathrm{m}`$) stars coincide spatially with the extended H II region DR 17 (Wendker et al. 1991). Comparison with the CO map (cf. Fig. 1) suggests that also this feature could be an embedded star nursery, but its angular separation and estimated distance of 800 pc - 1.5 kpc (Odenwald & Schwartz 1993) make a physical relation to Cyg OB2 highly improbable. Another group of stars is found at $`\alpha =20^\mathrm{h}29^\mathrm{m}5^\mathrm{s}`$ and $`\delta =+40^{}25^{}`$ which spans a wide range of colours, and which does not coincide with any known H II region or stellar cluster. The angular separation from Cyg OB2 suggests that this group is also not related to the association.

The morphology of Cyg OB2 as derived from the 2MASS catalogue is quite different to that inferred by RLP from DSS plates. While RLP found an elliptically shaped association of $`29\times 17`$ pc in size, the PSC analysis presented in this work suggests a spherical association of 60 pc in diameter. To investigate the origin of this discrepancy, the stellar density distribution around Cyg OB2 has also been derived from the DSS red plates, reproducing Fig. 23 of RLP. The result of this analysis is shown in the right panel of Fig. 5, again smoothed by a boxcar average of $`15^{}`$ to enhance the morphology of the association. Indeed, the 2MASS and the DSS stellar density distributions are markedly different. Comparison with the matter distribution in the region (cf. Fig. 1) clearly shows that the morphology of the DSS density distribution is largely influenced by the absorption pattern. The elliptical shape of Cyg OB2 in the DSS data is produced by an elongated region of rather low extinction that is enclosed by two dust lanes running from south-east to north-west. Several local maxima that appear like extensions from the Cyg OB2 association are simply areas of reduced absorption, where the number of background stars is locally enhanced. A massive system of molecular clouds below $`\delta =+40^{}30^{}`$ provides an efficient barrier for visible light, pushing the apparent centroid of Cyg OB2 towards the north. Hence DSS star counts are apparently not very suitable for the determination of the morphology of Cyg OB2. In contrast, a comparison with Fig. 1, and the above discussion of the completeness of the PSC sample, suggest that the 2MASS stellar density distribution is basically unaffected by absorption, hence it reveals indeed the true morphology of Cygnus OB2.

## 4 Size, mass, and stellar content

In order to determine the size, centre, and stellar content of Cyg OB2, radial star density profiles have been extracted from the 2MASS data by counting the number of stars intrinsically brighter than $`M_K2.2^\mathrm{m}`$ in concentric radial annuli around an assumed centre, divided by the annuli surface. Geometrical corrections have been applied for the outer annuli that partially fall outside the rectangular survey region. The profiles were fitted by a King law (King 1962)

$$f(r)=k\{[1+(r/r_\mathrm{c})^2]^{1/2}[1+(r_\mathrm{t}/r_\mathrm{c})^2]^{1/2}\}^2$$
(1)
on top of a constant to determine the core radius $`r_\mathrm{c}`$, the tidal radius $`r_\mathrm{t}`$, and the field star density. By searching the central position that minimises the radial extent of the profile the centre of Cyg OB2 has been determined to $`\alpha =20^\mathrm{h}33^\mathrm{m}10^\mathrm{s}`$ and $`\delta =+41^{}12^{}`$. The corresponding density profile is shown in Fig. 6. The central stellar density reaches $`4.5`$ stars arcmin-2 above the field star density, and drops to $`50\%`$ at a radius of $`13^{}`$, resulting in a half light radius of $`R_h=6.4`$ pc at a distance of 1.7 kpc.

Best fitting King parameters for the profile are $`r_c=29^{}\pm 5^{}`$ and $`r_t=93^{}\pm 20^{}`$, leading to a concentration parameter $`\mathrm{log}r_t/r_c`$ of $`0.5`$. The reduced $`\chi ^2`$ of the fit is only $`10.9`$, indicating that the King law is not a very accurate description of the radial density profile. Indeed, there is no physical reason to believe that Cyg OB2 should follow a King profile. The basic aim of using King profiles was the estimation of the field star density, and comparison with Fig. 6 convinces that at least this goal was reached.

By integrating the radial profile after subtraction of the fitted field star density, the total number of association stars with $`M_K2.2^\mathrm{m}`$ (corresponding to spectral types of F3V and earlier) amounts to $`8600\pm 1300`$. The error (as all following error quotations) includes a possible systematic uncertainty from the field star subtraction. This systematic uncertainty was estimated by applying different methods for the extraction of the association stars, such as analysing star density profiles along constant Right Ascension or declination, estimating the field star density from a circular region around the association, or by performing the analysis without removing the gradient in the field star density distribution.

Selecting only $`K`$ magnitudes brighter than $`11.55+0.66\times (JK)`$ limits the sample to stars intrinsically brighter than $`M_K0.37^\mathrm{m}`$, corresponding roughly to spectral type B9V and earlier. Repeating the radial density profile analysis for this limited sample allows the determination of the total number of OB stars in Cyg OB2 to $`2600\pm 400`$. This number is at the upper end of the range quoted by RLP (300 - 3000), confirming their suggestion that many of the highly reddened stars in the DSS survey are indeed OB stars. Further, restricting $`M_K3.55^\mathrm{m}`$, by requiring $`K7.76+0.66\times (JK)`$, selects only O stars from the PSC, resulting in a total number of $`120\pm 20`$ objects. In their survey, Massey & Thompson (1991) find 40 O stars within a central field of 0.35 degrees squared of the Cyg OB2 association. Integrating over the same region in the 2MASS data gives $`39\pm 6`$ O stars, demonstrating that the present analysis is in excellent agreement with the Massey & Thompson observations. Hence, the large number of O stars found in this analysis is mainly due to the large extent of the association, which previously was missed due to the high extinction in the area.

Taking an initial mass of 1.5 $`M_{\mathrm{}}`$ for a F3V star (Schmidt-Kaler 1982), and assuming the slope $`\mathrm{\Gamma }`$ of the initial mass function (IMF)222The Salpeter IMF has a slope of $`\mathrm{\Gamma }=1.35`$ in this prescription. to be comprised between $`1.1`$ and $`1.7`$, the number of association stars converts into a total stellar mass of $`(35)\times 10^4`$ $`M_{\mathrm{}}`$ above 1.5 $`M_{\mathrm{}}`$. Extrapolation of the IMF to lower masses using the prescription of Kroupa et al. (1993) results in a total association mass of $`(410)\times 10^4`$ $`M_{\mathrm{}}`$, where the boundaries correspond to a lower mass cut-off of $`1.0`$ and $`0.08`$ $`M_{\mathrm{}}`$, respectively. Since the actual value of the mass cut-off in Cyg OB2 is unknown, the total association mass is uncertain to within the quoted limits. However, even at the lower mass limit, Cyg OB2 would be the most massive OB association known in the Galaxy, being comparable in mass to a small globular cluster.

From the radial density profile and the total mass estimates, the central mass density of Cyg OB2 can be estimated. Assuming a total mass of $`4\times 10^4`$ $`M_{\mathrm{}}`$ results in a mass density of $`\rho _0=4060`$ $`M_{\mathrm{}}`$ pc-3, where the lower value comes from the extrapolation of the best fitting King profile into the centre, and the upper value comes from the observed central stellar density. Correspondingly, the central mass density estimate for a total mass of $`1\times 10^5`$ $`M_{\mathrm{}}`$ amounts to $`\rho _0=100150`$ $`M_{\mathrm{}}`$ pc-3. As before, the distance to Cyg OB2 has been assumed to 1.7 kpc.

## 5 Initial mass function

To derive the mass spectrum of Cyg OB2, sub-samples of the PSC data containing all member stars within a given interval $`[M_K^{\mathrm{low}},M_K^{\mathrm{up}}]`$ of intrinsic $`K`$ magnitudes have been defined. Additionally, the $`K`$ band extinctions $`A_K`$ were restricted to the range $`0A_K2`$ that has been found to enclose the association members (cf. Sect. 2). To define the intrinsic $`K`$ magnitude intervals, photometric cuts along the reddening lines have been applied using

$$\stackrel{~}{K}_{\mathrm{low}}KDM0.66\times (JK)\stackrel{~}{K}_{\mathrm{up}}$$
(2)
where $`M_K=0.0037+1.036\times \stackrel{~}{K}`$. The extinction cut was realised using the relation

$$A_K\frac{K+0.05717.835\times (JK)DM}{117.835/0.66}$$
(3)
which is an approximation of the extinction at the distance of Cyg OB2 as function of the position in the CMD for $`M_K2^\mathrm{m}`$ (see Appendix B for details about the derivation of these relations). In Eqs. (2) and (3), $`DM=11.2`$ is the assumed distance modulus of Cyg OB2.

Using these equations, the $`K`$,$`JK`$ plane has been divided into 7 areas of identical size which define 7 stellar sub-samples in the PSC data (cf. Fig. 7). The corresponding $`[\stackrel{~}{K}_{\mathrm{low}},\stackrel{~}{K}_{\mathrm{up}}]`$ and $`[M_K^{\mathrm{low}},M_K^{\mathrm{up}}]`$ intervals are summarised in Table 1. The $`[M_K^{\mathrm{low}},M_K^{\mathrm{up}}]`$ intervals have then been converted into initial mass intervals using the mass-luminosity relation derived in Appendix A. The corresponding mass boundaries are quoted in the last two columns of Table 1.

Two methods have been employed to extract the IMF from the PSC data. First, radial density profiles have been derived for the 7 sub-samples which were then fitted by King profiles in order to estimate the field star density in each of the samples. The field star density estimate was then subtracted from the radial density profile, and integration over the residual profile provided an estimate of the number of member stars in the sub-sample. Note that this method properly accounts for possible variations in the size of the association for different intrinsic magnitude intervals, which could result from mass segregation. However, the analysis revealed similar radial extents for all 7 sub-samples, indicating that no mass segregation has taken place in Cyg OB2. The second method consisted simply in integrating the number of stars in the field star subtracted CMD in the 7 areas indicated in Fig. 7. Here, the field star density is estimated from the area beyond $`1.05^{}`$ from the association centre, as discussed in Sect. 2. Note that in this case the gradient in the field stars density distribution has not been removed.

The results of both analyses are summarised in Fig. 8. The points with error bars were derived using radial density profiles (i.e. the first method), while the results from integrating the CMD are show as diamonds. Apparently, both methods lead to consistent IMFs. The data points are reasonably well fitted by a power law IMF of slope $`\mathrm{\Gamma }=1.6\pm 0.1`$ for both methods. The slope is compatible with the Kroupa et al. (1993) analysis of the solar neighbourhood, but considerably steeper than $`\mathrm{\Gamma }=1.0\pm 0.1`$ determined by Massey & Thompson (1991) for the central part of Cyg OB2. Note, however, that photometric mass determinations, as performed in this work, are less reliable than a spectroscopic study, as done by Massey & Thompson (1991), and a change in the employed mass-luminosity relation could easily flatten the IMF. On the other hand, the Massey & Thompson (1991) analysis did not extend to such low masses as the present analysis does, and the flat slope they determined could also reflect an incompleteness in the lower mass domain.

Finally, Fig. 7 suggests that the lowest mass bin could be subject to incompleteness since the integration region extends below the completeness limit. In particular, the completeness limit for the central part of Cyg OB2 may lie at even smaller magnitudes due to enhanced crowding, making $`\mathrm{\Gamma }=1.6`$ formally to a lower limit for the slope of the IMF. Nevertheless, the lowest mass bins in Fig. 8 fall not particularly below the IMF extrapolation from the higher masses, indicating that there are probably not many member stars which are missed. The number of $`8600\pm 1300`$ members earlier than spectral type F3V, derived in Sect. 4, is therefore probably a solid value.

## 6 Conclusions

The analysis of 2MASS data from the Cygnus region has revealed that Cygnus OB2 is considerably larger and more massive than previously thought. Previous surveys in the visible wavelength range, such as the extensive work of RLP, were substantially biased by the local extinction pattern in the field, which is caused by a foreground molecular cloud structure, known as the Great Cygnus Rift. The infrared data reveal now a spherically shaped association, with stellar extinctions reaching $`A_V20^\mathrm{m}`$ in the southern part of the field. Several highly reddened early-type objects that formerly were situated outside the association boundary are now lying within the association, indicating that they possibly are members of Cyg OB2. Examples are the massive binary system MWC349 (Cohen et al. 1985), the potential Luminous Blue Variable star G79.29+0.46 (Higgs et al. 1994), the Wolf-Rayet stars WR 145 and WR 146 (Niemela et al. 1998), and the recently discovered group of massive stars around the H II region DR 18 (Comerón & Torra 1999). Deep spectroscopic surveys of the stars in the highly obscured regions should help to consolidate the stellar population of Cyg OB2, and will improve our knowledge about the evolutionary state of the association.

RLP already recognised from their analysis that Cyg OB2 is unusually massive and compact, and hence suggested that it should be regarded as a new globular cluster of stars, similar to the blue globular clusters discovered by Hodge (1961) in the Large Magellanic Cloud (LMC). Indeed, the comparison of the association parameters, summarised in Table 2, to typical parameters for galactic OB associations, young open clusters, and globular clusters leads inevitably to this conclusion. For an OB association, Cyg OB2 is simply too compact. The typical mass density of OB associations is well below $`0.1`$ $`M_{\mathrm{}}`$ pc-3 (Blaauw 1964) and the mean size amounts to 137 pc (Garmany 1994). Cyg OB2, however, shows an average mass density of $`512`$ $`M_{\mathrm{}}`$ pc-3 within the inner 10 pc and a diameter of only 60 pc. Hence, in contrast to classical OB association, Cyg OB2 should be a gravitationally bound system (Blaauw 1964). Garmany & Stencel (1992) suggest that Cyg OB2 might be the nucleus of an OB association, hence it could be considered as an open cluster. However, for an open cluster (as well as for an OB association) Cyg OB2 is by far too massive. Typical masses for open clusters are at most $`10^3`$ $`M_{\mathrm{}}`$ and do not exceed $`10^4`$ $`M_{\mathrm{}}`$ (Bruch & Sanders 1983) while the total mass of Cyg OB2 amounts to $`410\times 10^4`$ $`M_{\mathrm{}}`$. Such a high mass is more representative for a small globular cluster.

In particular, the parameters of Cyg OB2 (as listed in Table 2) correspond fairly well to the typical parameters derived for young globular clusters in the LMC (e.g. Elson et al. 1987; Fischer et al. 1992). They have masses $`10^410^5`$ $`M_{\mathrm{}}`$, central densities $`10^2`$ $`M_{\mathrm{}}`$ pc-3, and extend to radii $`80`$ pc. Little or no mass segregation is observed in these objects and they show only a small age spread (Elson et al. 1989). The same is true for Cyg OB2. The determination of radial profiles for different intrinsic luminosities did not reveal any mass segregation, and the analysis of a sub-sample of massive stars in Cyg OB2 by Massey & Thompson (1991) did not suggest a considerable age spread. Thus I believe it is reasonable to conclude that Cygnus OB2 shows the same properties as the system of young globular clusters in the LMC – hence Cygnus OB2 should be considered as a member of this population.

Cyg OB2 is the first object of this class that has been identified within our own Galaxy. Young globular clusters (or young populous clusters as termed by Hodge (1961)) seem not to be rare objects. They were first identified in the LMC but are now also found using the Hubble Space Telescope in a variety of extragalactic star forming regions (Ho 1997). Although young globulars are often observed in association with starburst phenomena, they also appear to form in more quiescent environments, such as circumnuclear rings in relatively undisturbed galaxies (e.g. Maoz et al. 1996). It seems therefore plausible that such a system can also form in our own Galaxy, although the formation conditions are yet to be explored.

In a review about young globulars in the LMC, Freeman (1980) asked why these systems “form in the LMC and not in the Galaxy”? Apparently, they do also form in our Galaxy, but heavy obscuration through the interstellar gas and dust make them difficult to detect – or may leave them unrecognised. Infrared surveys may help to improve the detectibility of such systems, although the sensitivity and confusion limits restrict such studies to within a few kpc. Alternatively, one may search for the fingerprints that these systems leave on the surrounding interstellar medium (ISM). An example is the thermal radio emission related to the ionisation of the ISM by the numerous massive stars. Indeed, the 120 O stars in Cyg OB2 produce a Lyman continuum flux of roughly $`10^{51}`$ photons s-1, making it probably the most important source of ionisation in the Cygnus X complex. It has been suggested that the Cygnus X complex could be a large Strömgren sphere powered by Cyg OB2 (Landecker 1984 and references therein), hence the search for similar complexes in surveys of thermal radio emission could provide hints for further young globulars in the Galaxy.

Another tracer may come from the emerging field of gamma-ray line spectroscopy. The massive stars in Cyg OB2 are an important source of radioisotopes that are expelled through stellar winds and supernova explosions into the interstellar medium. These isotopes eventually emit monoenergetic gamma-ray photons that arise from nuclear de-excitations following the radioactive decay. An example is the 1.809 MeV gamma-ray line arising from the decay of radioactive 26Al, an isotope with a lifetime of $`10^6`$ years (see Prantzos & Diehl 1996 for a review). Indeed, besides the central galactic radian, the Cygnus X region is the most intense source of 1.809 MeV photons known in the Galaxy, which can be well explained by 26Al production in massive stars within Cyg OB2 (Knödlseder et al., in preparation). A more detailed study of the galactic plane in the 1.809 MeV gamma-ray line by the upcoming INTEGRAL mission (Vedrenne et al. 1999) could lead to the identification of further regions of high star forming activity, and hence provide a unique tool to unveil young globular clusters in our own Galaxy.

###### Acknowledgements.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. The author wants to thank T. Dame for kindly providing the CO data of the Cygnus region. The author is supported by an external ESA fellowship.

## Appendix A Stellar calibrations

The analysis presented in this work makes use of calibrations between intrinsic $`K`$ magnitudes and $`JK`$ colours on the one side, and spectral types and initial masses on the other side. The employed calibrations, valid for luminosity class V, are summarised in Table 3. The use of luminosity class V calibrations for all stars is justified by the young age of Cyg OB2 (Massey & Thompson 1991).

The relation between spectral type and absolute visual magnitudes $`M_V`$ has been compiled from data published by Vacca et al. (1996) for O3V - B0.5V stars, Humphreys & McElroy (1984) for B1V - B3V stars, and Schmidt-Kaler (1982) for B5V - G0V stars. The intrinsic colour calibrations $`(VK)_0`$ and $`(JK)_0`$ were taken from Wegner (1994) for spectral types O6V - B9V and from Koornneef (1983) for later spectral types. The $`(VK)_0`$ calibration of Wegner (1994) was extrapolated in the $`M_V`$,$`(VK)_0`$ plane to spectral type O3V using the relation

$$(VK)_0=0.296+0.114\times M_V$$
(4)
which has been obtained form a fit to the data for spectral types B2V and earlier. The same relation was also used to interpolate the $`(VK)_0`$ calibration for some O subtypes that were not covered by the analysis of Wegner (1994).

Based on the $`M_V`$ and $`(VK)_0`$ calibrations, the intrinsic $`K`$ magnitudes were derived using

$$M_K=M_V(VK)_0.$$
(5)
To derive the mass-luminosity relation, the initial stellar mass to spectral type relation has been taken from Schaerer & de Koter (1997) for O and Schmidt-Kaler (1982) for later type stars. Using a fit in the $`\mathrm{log}M`$,$`M_K`$ plane, the following relation has been established

$$\mathrm{log}M=\{\begin{array}{cc}\hfill 0.6050.201M_K+0.0055M_K^2& M_K>3.70\hfill \\ \hfill 0.6300.557M_K& \mathrm{else}\hfill \end{array}$$
(6)
## Appendix B Photometric sample selection

For the photometric sample selection two constraints have been applied in the $`K`$,$`JK`$ plane. First, all association stars with identical intrinsic magnitude $`M_K`$ are collected along the reddening line. To understand the selection, recall that the apparent $`K`$ magnitude of a Cyg OB2 member star is related to the intrinsic magnitude $`M_K`$, the extinction $`A_K`$, and the distance modulus $`DM`$ via

$$K=M_K+DM+A_K.$$
(7)
The extinction of a star can be estimated from the colour excess

$$E(JK)=(JK)(JK)_0$$
(8)
using

$$A_K=0.66\times E(JK)$$
(9)
(Rieke & Lebofsky 1985). Thus all member stars of identical intrinsic magnitude but different extinction lie in the $`K`$,$`JK`$ CMD on a line of constant slope, given by

$$K=M_K0.66\times (JK)_0+DM+0.66\times (JK).$$
(10)
Hence, selecting stars within a band

$$\stackrel{~}{K}_{\mathrm{low}}KDM0.66\times (JK)\stackrel{~}{K}_{\mathrm{up}}$$
(11)
results in a reddening independent selection of stars, where

$$\stackrel{~}{K}M_K0.66\times (JK)_0.$$
(12)
The values of $`\stackrel{~}{K}`$ are listed in column 6 of Table 3. Fitting a linear relation to these data results in

$$M_K=0.0037+1.036\times \stackrel{~}{K}.$$
(13)
Second, the $`K`$ band extinction $`A_K`$ of an association star can be estimated by calculating its displacement along the reddening line. To understand this, notice that for main sequence stars, $`M_K`$ is a function of intrinsic colour $`(JK)_0`$, which for massive stars can be roughly approximated by a linear function. From the calibration table (cf. Table 3) the relation

$$M_K0.057+17.835\times (JK)_0$$
(14)
has been determined for $`M_K<2^\mathrm{m}`$. Using Eqs. (8) and (9) the relation can be rewritten as

$$M_K0.057+17.835\times \left((JK)\frac{A_K}{0.66}\right).$$
(15)
Replacing this approximation in Eq. (7) and solving for $`A_K`$ result in an approximative relation between apparent magnitude and colour and extinction in the $`K`$ band:

$$A_K\frac{K+0.05717.835\times (JK)DM}{117.835/0.66}.$$
(16)


format: text

text:
# Brane world from IIB matrices

## Introduction.—

The idea that our four-dimensional space time is realized as a brane in a non-compact higher-dimensional space time has recently attracted much attention. Through many works during the last few years, it is expected to provide natural resolutions to many long-standing problems in the Standard Model. The hierarchy problem is transmuted into a geometrical one , and it was further argued that the exponential dependence of the “warp” factor in the extra directions reduces the problem to a fine-tuning of order 50 . The cosmological constant problem may also be resolved in such a setup . It has been argued that any nonzero higher-dimensional cosmological constant is absorbed into the warp factor, and that the four-dimensional cosmological constant is automatically tuned to zero (or to a very small number). A possible obstruction to the idea (as opposed to a more conventional idea using Kaluza-Klein compactifications) is that gravity may propagate in higher dimensions and thereby contradicts the 4D Newton’s law observed in the low energy scale. However, the particular AdS-type background metric that arises naturally in such a setup allows a normalizable zero mode of the graviton bound to the brane . A small correction to the 4D Newton’s law due to the continuum spectrum of massive modes is argued to be small enough to be compatible with the experimental bound. All these attractive features of the idea lead us to hope that there is a natural string-theory realization of the brane world scenario.

In Ref. , we have proposed a dynamical mechanism which may realize a flat four-dimensional space time as a brane in type IIB superstring theory. Obviously, such a mechanism should inevitably be of nonperturbative nature. Indeed, our mechanism was based on the IKKT version  of the Matrix Theory , namely the IIB matrix model, which is conjectured to be a nonperturbative definition of type IIB superstring theory. The model is a supersymmetric matrix model composed by ten $`N\times N`$ hermitian bosonic matrices and sixteen $`N\times N`$ hermitian fermionic matrices. The space time is represented by the eigenvalues of the bosonic matrices. The model has manifest ten-dimensional Lorentz invariance, where the bosonic and fermionic matrix elements transform as a vector and a Majorana-Weyl spinor, respectively. The integral over the fermionic matrices yields a pfaffian which is complex in general. This poses a technical difficulty known as the ‘complex action’ problem in studying the IIB matrix model by Monte Carlo simulation. Monte Carlo studies incorporating only the modulus of the pfaffian (and omitting the phase by hand) showed that the space-time becomes isotropic in ten dimensions in the large-$`N`$ limit . This result suggests that the phase of the pfaffian must play a crucial rôle if a brane world naturally arises in the type IIB superstring theory. The effect of the phase is to favour configurations for which the phase becomes stationary. Such an effect has been studied within a saddle-point approximation and found to enhance lower-dimensional brane-like configurations considerably .

In this Letter, we demonstrate our mechanism more explicitly by studying a simplified model using Monte Carlo simulation. In this case, we find that the dominant saddle-points are given by configurations with only three-dimensional extent.

## The mechanism.—

The IIB matrix model  is formally a zero-volume limit of $`D=10`$, $`𝒩=1`$, pure super Yang-Mills theory. The partition function of the IIB matrix model (and its generalizations to $`D=4`$ and $`D=6`$) can be written as

$$Z_{\mathrm{IIB}}=\text{d}A\text{e}^{S_\mathrm{b}}Z_f[A],$$
(1)
where $`S_\mathrm{b}=\text{Tr }([A_\mu ,A_\nu ]^2)/4`$, and $`\mathrm{\Gamma }[A]=\mathrm{ln}Z_f[A]`$ represents the effective action induced by integration over the fermionic matrices. The dynamical variables $`A_\mu `$ ($`\mu =1,\mathrm{},D`$) are $`D`$ bosonic $`N\times N`$ traceless hermitian matrices. Expanding $`A_\mu `$ as $`A_\mu =_{a=1}^{N^21}A_\mu ^at^a`$ in terms of the generators $`t^a`$ $`(a=1,\mathrm{},N^21)`$ of SU($`N`$), the integration measure $`\text{d}A`$ is given as $`\text{d}A_{\mu =1}^D_{a=1}^{N^21}dA_\mu ^a`$. The generators $`t^a`$ are normalized as $`\text{Tr }(t^at^b)=2\delta _{ab}`$.

The fermion integral $`Z_f[A]`$ is complex in general for $`D=10`$, $`N4`$ and for $`D=6`$, $`N3`$ . We restrict ourselves to these cases in what follows. In the $`D=10`$ case, the fermion integral $`Z_f[A]`$ is given by the pfaffian $`\text{Pf}`$, where $``$ is a $`16(N^21)\times 16(N^21)`$ complex antisymmetric matrix defined by

$$_{a\alpha ,b\beta }\text{Tr }\left(t^a(𝒞\mathrm{\Gamma }_\mu )_{\alpha \beta }[A_\mu ,t^b]\right)$$
(2)
regarding each of $`(a\alpha )`$ and $`(b\beta )`$ as a single index. Here, $`\mathrm{\Gamma }_\mu `$ $`(\mu =1,\mathrm{},10)`$ are ten-dimensional Weyl-projected $`16\times 16`$ gamma matrices, satisfying $`𝒞\mathrm{\Gamma }_\mu 𝒞^{}=(\mathrm{\Gamma }_\mu )^{}`$ with $`𝒞=𝒞^{}`$ being the unitary charge conjugation matrix. Similarly in the $`D=6`$ case, the fermion integral $`Z_f[A]`$ is given by the determinant $`det^{(6)}`$, where $`^{(6)}`$ is a $`4(N^21)\times 4(N^21)`$ complex matrix defined by

$$_{a\alpha ,b\beta }^{(6)}\text{Tr }\left(t^a(\mathrm{\Gamma }_\mu ^{(6)})_{\alpha \beta }[A_\mu ,t^b]\right)$$
(3)
regarding each of $`(a\alpha )`$ and $`(b\beta )`$ as a single index. Here, $`\mathrm{\Gamma }_\mu ^{(6)}`$ $`(\mu =1,\mathrm{},6)`$ are six-dimensional Weyl-projected $`4\times 4`$ gamma matrices.

Since the fermion integral $`Z_f[A]`$ is a complex quantity for the cases under consideration, let us write it as $`Z_f[A]=\mathrm{exp}(\mathrm{\Gamma }^{(\mathrm{r})}i\mathrm{\Gamma }^{(\mathrm{i})})`$. In Ref. , the effect of the phase $`\mathrm{\Gamma }^{(\mathrm{i})}`$ in the path integral (1) has been studied using a saddle-point approximation, whose validity has been also discussed. The saddle-point equation for $`\mathrm{\Gamma }^{(\mathrm{i})}`$ is given by

$$\frac{\mathrm{\Gamma }^{(\mathrm{i})}}{A_\mu ^a}=0a,\mu .$$
(4)
It is useful to introduce the following classification of “brane” configurations

$$\mathrm{\Omega }_d=\{\{A_\mu \}n_\mu ^{(i)}(i=1,\mathrm{},Dd),n_\mu ^{(i)}A_\mu =0\},$$
(5)
where $`n_\mu ^{(i)}`$ ($`i=1,\mathrm{},Dd`$) are $`(Dd)`$ linearly independent $`D`$-dimensional real vectors. Namely, $`\mathrm{\Omega }_d`$ represents a set of configurations with less than $`d`$-dimensional extent. Note that $`\mathrm{\Omega }_1\mathrm{\Omega }_2\mathrm{}\mathrm{\Omega }_D`$, where $`\mathrm{\Omega }_D`$ is nothing but the whole configuration space of the model. In Ref.  we proved that all configurations in $`\mathrm{\Omega }_{D2}`$ are solutions to the saddle-point equation (4). Assuming that the configurations in $`\mathrm{\Omega }_{D2}`$ are the dominant saddle-point configurations, we still have to integrate over those configurations to determine the actual dimensionality of the space-time. In fact, the gaussian fluctuation of the phase $`\mathrm{\Gamma }^{(\mathrm{i})}`$ around the saddle-points gives a huge enhancement to the brane configurations with lower dimensionality, and this enhancement cancels exactly the entropical barrier against having such configurations. In the $`D=10`$ case, this provides a dynamical mechanism for the possible appearance of four-dimensional space time as a brane in ten-dimensional space time.

## A simplified model.—

In order to investigate how our mechanism works, we consider a simplified model describing an integration over the saddle-point configurations. Specifically, we consider the integral

$$Z=\text{d}A\text{e}^{\beta F[A]\gamma f[A]},$$
(6)
where the functions $`F[A]`$ and $`f[A]`$ are defined as

$`F[A]`$ $``$ $`{\displaystyle \underset{\mu =1}{\overset{D}{}}}{\displaystyle \underset{a=1}{\overset{N^21}{}}}\left({\displaystyle \frac{\mathrm{\Gamma }^{(\mathrm{i})}}{A_\mu ^a}}\right)^2,`$ (7)
$`f[A]`$ $``$ $`{\displaystyle \underset{\mu =1}{\overset{D}{}}}{\displaystyle \underset{a=1}{\overset{N^21}{}}}(A_\mu ^a)^2.`$ (8)
Since the function $`F[A]`$ vanishes if and only if the configuration $`\{A_\mu \}`$ satisfies the saddle-point equation (4), the integral (6) is dominated by the saddle-point configurations in the large-$`\beta `$ limit. The function $`f[A]`$ makes the integral (6) convergent as long as the coefficient $`\gamma `$ is fixed to be a real positive number. In fact, the parameter $`\gamma `$ can be absorbed by an appropriate rescaling of $`A_\mu ^a`$ and $`\beta `$. Therefore we take $`\gamma =1/2`$ in what follows without loss of generality. Note also that the simplified model (6) is invariant under a SO($`D`$) transformation $`A_\mu \mathrm{\Lambda }_{\mu \nu }A_\nu `$, where $`\mathrm{\Lambda }_{\mu \nu }\text{SO}(D)`$, and a SU($`N`$) transformation $`A_\mu gA_\mu g^{}`$, where $`g\text{SU}(N)`$, which are the symmetries of the original model (1).

Using the invariance of the partition function (6) under the change of variables $`A_\mu ^a\lambda A_\mu ^a`$, one can obtain an exact relation

$$\frac{1}{2}f_\beta \beta F_\beta =\frac{D}{2}(N^21),$$
(9)
where the ensemble average $`_\beta `$ is defined with the partition function (6). Assuming that $`f_\beta `$ goes to a constant $`c`$ for $`\beta \mathrm{}`$, we obtain the asymptotic behavior of $`F_\beta `$ for large $`\beta `$ as

$$F_\beta \frac{b}{\beta },$$
(10)
where the coefficient $`b`$ is given as $`b=c/2D(N^21)/2`$. This confirms the above claim that the integral (6) is dominated by the saddle-point configurations in the large-$`\beta `$ limit.

A quantity which fully characterizes the dimensionality $`d`$ of a given configuration can be given by the moment of inertia tensor $`T`$ defined by the $`D\times D`$ real symmetric matrix

$$T_{\mu \nu }=\underset{a=1}{\overset{N^21}{}}A_\mu ^aA_\nu ^a.$$
(11)
A configuration $`\{A_\mu \}`$ belongs to $`\mathrm{\Omega }_d`$, if and only if the number of zero eigenvalues of the matrix $`T`$ is more than or equal to $`(Dd)`$. Let us denote the eigenvalues of the matrix $`T`$ as $`\lambda _i`$ ($`i=1,\mathrm{},D`$), where $`\lambda _1\lambda _2\mathrm{}\lambda _D0`$. We can determine the dimensionality of the dominant saddle-point configurations from the ensemble average of the eigenvalues $`\lambda _i_\beta `$ in the $`\beta \mathrm{}`$ limit.

We address this issue by performing Monte Carlo simulation using a Metropolis algorithm. We create a trial configuration $`\{A_\mu ^{}\}`$ by replacing an element $`A_\mu ^a`$ of the previous configuration $`\{A_\mu \}`$ with a new one generated with the probability distribution $`\frac{1}{\sqrt{2\pi }}\mathrm{exp}[\frac{1}{2}(A_\mu ^a)^2]`$. The trial configuration is accepted with the probability $`\mathrm{min}(1,\mathrm{exp}(\beta \mathrm{\Delta }F))`$, where $`\mathrm{\Delta }F=F[A^{}]F[A]`$. This procedure is repeated for all the elements of the configuration. The computational effort required for the above algorithm is of order O$`(N^8)`$ per one sweep, which is much larger than that for the simulation encountered in Refs. . Due to this, results with high statistics are obtained only for the case of $`D=6`$ and $`N=3`$ (we have made 192,000 accepted updates for each $`\beta 384`$ and 768,000 for $`\beta =512`$).

## Results.—

In the upper part of Fig. 1 we plot the left hand side of (9), which demonstrates the validity of our simulation. In the lower part of Fig. 1 we plot the average $`F[A]_\beta `$ against $`\beta `$ in a log-log scale. The straight line represents the fit of the data for $`\beta 16`$ to the predicted large-$`\beta `$ behavior (10) with $`b=1.7(1)`$.

Fig. 2 shows the six eigenvalues of the moment of inertia tensor $`T`$ as a function of $`\beta `$. We find that the three smallest eigenvalues $`\lambda _i_\beta `$ ($`i=4,5,6`$) are monotonously decreasing with a pronounced power law behavior. Fitting the data for $`\beta 16`$ to the power law behavior, we extract the powers $`0.040(4)`$,$`0.199(2)`$ and $`0.450(6)`$, respectively. Similarly, the powers for $`N=4`$ are extracted to be $`0.11(1)`$,$`0.26(2)`$ and $`0.36(2)`$. Thus we conclude that the dominant saddle-point configurations of the simplified model (6) has only three-dimensional extent. Preliminary results for $`D=10`$, $`N=4`$ suggest that this is the case also for $`D=10`$.

## Discussion.—

The results presented in the previous section shows clearly that the stationarity of the phase $`\mathrm{\Gamma }^{(i)}`$ indeed enhances lower-dimensional brane configurations considerably, thus demonstrating our mechanism. In particular, while the existence of saddle-point configurations other than the configurations in $`\mathrm{\Omega }_{D2}`$ is not excluded, our results suggest that such configurations, even if they exist, can safely be neglected on statistical grounds. Given this observation, the reason why we obtain the dimensionality ‘three’ from the dominant saddle-point configurations of the simplified model (6) can be understood analytically. As has been done in Ref.  for the IIB matrix model, we can rewrite the $`\beta \mathrm{}`$ limit of the simplified model (6) as an integral over the configurations in $`\mathrm{\Omega }_{D2}`$. The gaussian fluctuation of the phase $`\mathrm{\Gamma }^{(\mathrm{i})}`$ should be taken into account by the corresponding Hesse matrix, which is, in the present case, just the square of the one for the IIB matrix model (or its $`D=6`$ version). Due to this, the gaussian fluctuation enhances lower-dimensional brane configurations even more strongly than in the IIB matrix model and overwhelms the entropical barrier against having those configurations.

The fact that the lowering of the dimensionality stops at three instead of continuing further down can be understood as follows. We first note that the fermion integral $`Z_f[A]`$ vanishes for configurations in $`\mathrm{\Omega }_2`$ . Therefore, the phase $`\mathrm{\Gamma }^{(i)}`$ is actually ill-defined for configurations in $`\mathrm{\Omega }_2`$. Still, we can consider configurations with $`A_1`$, $`A_2`$ being generic and $`A_i(i=3,\mathrm{},D)`$ being of order $`ϵ`$. One can easily see that the function $`F[A]`$ is diverging as $`ϵ^2`$ for $`ϵ0`$. Therefore, the $`\beta `$ term in (6) suppresses configurations in $`\mathrm{\Omega }_2`$ strongly. (In other words, configurations in $`\mathrm{\Omega }_2`$ are not saddle-point configurations, although $`\mathrm{\Omega }_2\mathrm{\Omega }_{D2}`$.)

As is clearly shown in the present work, the enhancement occurs precisely when the space-time becomes a flat lower-dimensional hyperplane, which we consider as a very attractive feature of our mechanism. Namely, our mechanism has a built-in structure in which the brane that appears as a result of the nonperturbative string dynamics is very likely to be flat. According to our mechanism, scenarios with two separated branes (i.e. our world and the so-called ‘Planck’ brane as in Ref. ), their extensions to many branes , and scenarios with mutually intersecting branes  seem to be unnatural.

Let us also comment on a connection of our mechanism to the brane world scenario. In Ref. , the IIB matrix model is expanded around a D3-brane configuration perturbatively and four-dimensional noncommutative Yang-Mills theory has been obtained . The (perturbatively stable) theory, which is obtained in this way from the IIB matrix model, has been recently identified  with a type IIB superstring theory in $`\text{AdS}_5\times \text{S}^5`$ with an infinite $`B`$-field background. Remarkably the metric that appears in the corresponding supergravity solution takes the form of Randall-Sundrum’s type , and thereby allows for a four-dimensional Newton’s law. We expect that brane configurations similar to the D3-brane configuration considered above as a background in the IIB matrix model should appear dynamically as a result of our mechanism. Thus our mechanism is rather directly related to the brane world scenario.

In the IIB matrix model, the enhancement and the entropical barrier are exactly balanced and the actual dimensionality of the brane is expected to be determined as a result of large $`N`$ dynamics. In this regard, we recall that the low-energy effective theory of the IIB matrix model has been shown to be described by a branched-polymer like system in Ref. . There it was further argued that a typical double-tree structure that appears in the effective theory might cause a collapse of the configuration to a lower-dimensional manifold. Whether the actual dimensionality of the brane turns out to be four or not can be investigated directly by performing the integration over the saddle-point configurations as formulated in Ref. . We hope that Monte Carlo techniques developed in Ref.  will enable us to address such an issue in near future.

## Acknowledgments.—

We thank S. Iso, H. Kawai, F.R. Klinkhamer, C. Schmidhuber, R.J. Szabo and K. Takenaga for helpful discussions. J.N. is supported by the Japan Society for the Promotion of Science as a Postdoctoral Fellow for Research Abroad. The work of G.V. is supported by MURST (Italy) within the project “Fisica Teorica delle Interazioni Fondamentali”.


text:
# REFERENCES

Comment on “A Convergent Series for the
Effective Action of QED”

In a recent letter, Cho and Pak claim to have found an additional contribution to the quantum electrodynamic one-loop effective action: a “logarithmic correction term” \[see the final result in Eq. (2) and the remark in Ref.  of \]. However, the “logarithmic correction term” found by Cho and Pak vanishes when the final result is written in terms of the finite, renormalized, physical electron charge.

Because current determinations of fundamental constants  rely on renormalized QED perturbation theory – without “logarithmic correction terms” of the kind advocated by Cho and Pak –, it is of prime general interest to point out that these terms do not appear if on-mass shell renormalization is used. In the on-mass shell scheme, the renormalized QED effective Lagrangian (see e.g. \[3, Eq. (3.43)\]) reads

$`\mathrm{\Delta }={\displaystyle \frac{e^2}{8\pi ^2}}\underset{ϵ,\eta 0^+}{lim}{\displaystyle _\eta ^{\mathrm{i}\mathrm{}+\eta }}{\displaystyle \frac{\mathrm{d}s}{s}}\mathrm{e}^{(m^2\mathrm{i}ϵ)s}`$ (2)
$`\times \left[ab\mathrm{coth}(eas)\mathrm{cot}(ebs){\displaystyle \frac{a^2b^2}{3}}{\displaystyle \frac{1}{(es)^2}}\right].`$
The two latter terms in the integrand are counterterms. The last term simply removes a divergent constant from the Lagrangian, while the term $`(a^2b^2)/3`$ – if it were not removed –, would lead to a logarithmic divergence at small eigentime $`s`$. This logarithmically divergent term, however, is proportional to the leading-order Maxwell Lagrangian $`_{\mathrm{cl}}=(b^2a^2)/2`$ and leads to a $`Z_3`$–renormalization \[see Eq. (8-97) of \]. Specifically, the introduction of the cut-off parameter $`\mu `$ in  leads to a logarithmic term $`\left[1(e^2/12\pi ^2)\mathrm{ln}(m^2/\mu ^2)\right]`$ which multiplies $`_{\mathrm{cl}}`$. In order to insure compliance with the renormalization conditions of on-mass shell renormalization \[see Eqs. (8-96d) and (8-96e) of \], a further counterterm $`+(e^2/12\pi ^2)\mathrm{ln}(m^2/\mu ^2)_{\mathrm{cl}}`$ has to be added to the Lagrangian. As a consequence, the logarithmic correction term is absent \[see Eqs. (2) – (6) of \]. For the particular problem at hand, the on-mass shell scheme is well motivated even from a purely mathematical point of view, as it is evident from the partial fraction theorem discussed in Sec. 3 of .

If Cho and Pak use different renormalization conditions, then the logarithmic correction term has to be reabsorbed into the physical charge of the electron, by considering the effect that the term has on matrix elements of transition currents \[see the elucidating discussion on page 325 of \]. In this case, we are forced to interpret $`e_{\mathrm{ph}}^2(\mu )=e^2[1(e^2/12\pi ^2)\mathrm{ln}(m^2/\mu ^2)]+𝒪(e^4)`$ as the physical charge, in which case the “logarithmic correction term”  is reabsorbed in a renormalization of charge. When expressing $`\mathrm{\Delta }`$ in terms of $`e_{\mathrm{ph}}^2(\mu )`$ instead of $`e^2`$, the resulting further modification of $`\mathrm{\Delta }`$ is of the same order as the two-loop effective Lagragian and therefore beyond the validity of the one-loop approximation inherent to Eq. (2). Finally, we would like to remark here that a renormalization-group (RG) improved running of the electron charge, based on the RG invariance of the effective action, has been discussed by Dittrich and Reuter (Ch. 8 of ) and Ritus , and that, in the latter case, two-loop effects are consistently taken into account in the analysis of the evolution of the electromagnetic charge.

Finally, we stress that the potentially important remark in Ref.  of  falsely suggests that the “usual” result for $`\mathrm{\Delta }`$ given in Eq. (2) is incomplete without the “logarithmic correction term”. Helpful conversations with H. Gies, B. R. Holstein, D. G. C. McKeon, C. Schubert, V. M. Shabaev and G. Soff are gratefully acknowledged.
D. R. Lamma), S. R. Vallurib), U. D. Jentschurac), and E. J. Wenigerd)
a)Georgia Institute of Technology (GTRI/EOEML),
Atlanta, GA 30332-0834
b)University of Western Ontario (Appl. Maths., Physics & Astronomy), London N6A 3K7, Canada
c)TU Dresden, D-01062 Dresden, Germany
d)Universität Regensburg (Theoretische Chemie),
D-93040 Regensburg, Germany
PACS: 12.20.Ds, 12.20.-m, 11.10.Jj, 11.15.Tk


text:
# Direct calculation of the hard-sphere crystal/melt interfacial free energy

## Abstract

We present a direct calculation by molecular-dynamics computer simulation of the crystal/melt interfacial free energy, $`\gamma `$, for a system of hard spheres of diameter $`\sigma `$. The calculation is performed by thermodynamic integration along a reversible path defined by cleaving, using specially constructed movable hard-sphere walls, separate bulk crystal and fluid systems, which are then merged to form an interface. We find the interfacial free energy to be slightly anisotropic with $`\gamma `$ = 0.62$`\pm 0.01`$, 0.64$`\pm 0.01`$ and 0.58$`\pm 0.01k_BT/\sigma ^2`$ for the (100), (110) and (111) fcc crystal/fluid interfaces, respectively. These values are consistent with earlier density functional calculations and recent experiments measuring the crystal nucleation rates from colloidal fluids of polystyrene spheres that have been interpreted \[Marr and Gast, Langmuir 10, 1348 (1994)\] to give an estimate of $`\gamma `$ for the hard-sphere system of $`0.55\pm 0.02k_BT/\sigma ^2`$, slightly lower than the directly determined value reported here.

A detailed microscopic description of the interface between a crystal and its melt is necessary for a full understanding of such important phenomena as homogeneous nucleation and crystal growth. Computer simulation studies of model materials have had some success in elucidating the phenomenology of such systems, the importance of such work being enhanced by the near absence of reliable experimental studies. These efforts, however, have been primarily focused on structural and dynamical properties, since the central thermodynamic property, the crystal/melt interfacial free energy, is difficult to measure by simulation or experiment. In this work we report the results of a direct calculation via molecular-dynamics (MD) simulation of the crystal/melt surface free energy of the hard-sphere system, one of the most important reference models for simple materials.

The crystal/melt surface free energy, $`\gamma `$, is defined as the (reversible) work required to form a unit area of interface between a crystal and its coexisting melt. Experimentally, $`\gamma `$ can be measured either indirectly from measurements of crystal nucleation rates interpreted through classical nucleation theory, or directly by contact angle measurements . Using the former method, Turnbull  estimated $`\gamma `$ for a number of materials and found a strong empirical correlation between the values obtained and the latent heat of fusion for each material given by the relation $`\gamma C_T\mathrm{\Delta }_fH\rho ^{2/3}`$, where $`\rho `$ is the number density of the crystal and with $`C_T`$ (the Turnbull coefficient) taking on the value 0.45 for most metals and 0.32 for other mostly nonmetallic materials. For the hard-sphere system considered in this work, recent experiments of the crystallization kinetics of a colloidal suspension of uncharged polystyrene spheres, which closely mimic hard spheres, have been interpreted within a classical nucleation model to yield an estimate for $`\gamma `$ of the hard-sphere system of 0.55$`\pm 0.02`$$`k_BT/\sigma ^2`$ . This value is in agreement  with that predicted using the empirical relationship above assuming a $`C_T`$ of 0.45 and values of $`\mathrm{\Delta }_fH`$ and coexistence densities for hard spheres as determined by MD simulation . Unfortunately, the accuracy of the values of $`\gamma `$ obtained from nucleation rates is severely limited by the approximations inherent in classical nucleation theory. More accurate values can be obtained directly using contact angles, but such experiments are difficult and only a few materials have been studied to date. One notable example is a series of grain boundary contact angle experiments on bismuth  that determined $`\gamma `$ to be relatively independent of crystal orientation at 61.3$`\times 10^3`$ J/m2, which is about 10% higher than Turnbull’s estimate from nucleation rates of 54.4$`\times 10^3`$ J/m2.

In recent years, the primary theoretical approach to studying the structure and thermodynamics of the crystal/melt interface has been density-functional theory (DFT) . For these studies, the hard-sphere system has been the benchmark calculation, due to the simplicity of the interaction and the availability of accurate, analytical formulas for the properties of the fluid. However, as discussed by Marr  the value of $`\gamma `$ obtained is highly dependent on the approximations used to construct the DFT and the reported values range from 0.25 to 4.00$`k_BT/\sigma ^2`$. The DFT studies also disagree dramatically in the degree of orientation dependence of the interfacial free energy. Unfortunately, in the absence of simulation results it has been difficult to assess the validity of the individual approaches, although only the DFT approach of Curtin ($`\gamma =0.62k_BT/\sigma ^2`$) and the related one of Marr and Gast ($`\gamma =0.60k_BT/\sigma ^2`$) come close to the nucleation result of $`0.55\pm 0.02k_BT/\sigma ^2`$.

To date, the only reliable calculation of the crystal/melt interfacial free energy via simulation is that of a system of particles interacting via a truncated Lennard-Jones potential by Broughton and Gilmer . In that work, a series of continuous, external “cleaving” potentials are used to separate (cleave) separate samples of bulk liquid and fcc crystal, prepared at the calculated coexistence temperature and densities. The solid and liquid slabs thus produced are then placed next to one another and the cleaving potentials removed to merge them into a coexisting interface. The reversible work to perform these steps can be calculated by thermodynamic integration, giving a direct calculation of $`\gamma `$ for this system. The values of $`\gamma `$ at the triple point were found to be statistically isotropic with $`\gamma \sigma ^2/ϵ=0.35\pm 0.02,0.34\pm 0.02`$ and $`0.36\pm 0.02`$ for the (111), (100) and (110) interfaces, respectively.

For the hard-sphere system, the Broughton-Gilmer procedure must be modified since the algorithm for MD simulation of discontinuous hard-core potentials is conceptually very different from the algorithm for continuous potentials. The latter are performed by integrating the system of ordinary differential equations, while in the former the dynamical algorithm proceeds on a collision by collision basis. Therefore, incorporating continuous cleaving potentials into the collisional algorithm would result in lost efficiency and substantial modification of the structure of the algorithm.

In this Letter we introduce an approach, which uses only hard-sphere interactions in order to cleave the bulk hard-sphere systems. This allows us to apply the Broughton-Gilmer cleaving procedure to the hard-sphere system with only minor changes to the algorithm structure. The idea of our approach is illustrated in Fig. 1. To cleave the bulk system at a cleaving plane (shown by the dashed line in Fig. 1), the spheres are assigned types 1 or 2 based on their position relative to the plane. Next, two walls, which are also assigned types 1 and 2, are placed on the opposite sides of the cleaving plane. The two types are introduced in order to specify interaction between the spheres and the walls; namely, the walls interact only with the spheres of similar type. Therefore, when the walls are placed as shown in Fig. 1, and the distance from the walls to the cleaving plane is larger than the sphere radius, the walls do not interact with the spheres. It is important that, during a simulation run, a sphere changes its type whenever it crosses the cleaving plane. Because of the periodic boundary conditions in the $`z`$ direction, another plane must be defined sufficiently far away from the cleaving plane, where the spheres also change type.

The cleaving of the system is achieved by slowly moving the walls towards each other (as shown by the arrows in Fig. 1), starting from the initial position $`z_i`$, where the walls do not interact with the system, till $`z_f`$, where the spheres of different types no longer collide with each other at the cleaving plane.

During the process, the average pressure on the walls, $`P(z)`$, is measured as a function of wall position. The work per unit area required to perform the cleaving is given by the integral

$$w=_{z_i}^{z_f}P(z)𝑑z.$$
(1)
Thus the crystal-fluid interfacial free energy, $`\gamma `$, can be measured in the reversible process involving the following four steps: (1) Cleave the bulk crystal by inserting two walls at the cleaving plane and moving them from $`z_i`$ to $`z_f`$; (2) Cleave the bulk fluid in a similar way; (3) Juxtapose the cleaved crystal and fluid systems by changing the periodic boundary conditions while retaining the crystal and fluid systems restricted by the respective cleaving walls; (4) Slowly move the walls back to their initial positions with respect to the cleaving planes. This series of steps is the same as that used by Broughton and Gilmer, except that, in our case, no work is done on the system in step 3. The interfacial free energy is given by the sum $`\gamma =w_1+w_2+w_4`$, where $`w_4`$ is negative and consists of the work done by the walls on the crystal and fluid parts of the system during the process of removing the cleaving walls.

The structure of the cleaving walls is crucial to the success of the procedure. The main requirement is that the insertion of the walls perturbs the systems as little as possible. Our approach is to use walls made of layers of ideal crystal oriented in correspondence with the orientation of the crystal system. For the (100) and (111) orientations it is sufficient to use one layer, while for the (110) orientation we use two layers. Such a choice of the wall structure ensures minimal perturbation of the cleaved crystal, while the cleaved fluid is expected to form properly oriented interfacial layers. On a technical side, the implementation of such a wall structure is quite simple, since we can treat collisions with the walls in the same manner we treat collisions between all the spheres in the system, except that the spheres forming the walls are assigned infinite mass.

The position of the walls, $`z`$, is measured by the distance of the centers of the spheres forming the walls to the cleaving plane. Obviously, the walls do not interact with the system when $`z>\sigma `$. Therefore, we set the initial position of the walls at $`z_i=\sigma `$. The pressure $`P(z)`$ is obtained by moving the walls from $`z_i`$ to $`z_f`$ in steps of 0.01$`\sigma `$. In order to move the walls to a new position, we assign a small velocity (about 0.1% of the average particle velocity) to the spheres forming the walls, and run the simulation until the walls reach the new position. At that moment the wall velocity is set to zero, and the velocities of the particles are rescaled in order to restore the initial value of the total kinetic energy of the system. Before measuring the pressure, we allow the system to relax in an equilibration run. In order to verify the reversibility of the cleaving process, we have also simulated the reverse process and measured the pressure while moving the walls from $`z_f`$ back to $`z_i`$. The details of the simulation process and obtained results follow.

Step 1: Cleaving the crystal For each of the three orientations, we start with a crystal at a density $`\rho _c=1.037\sigma ^3`$, which corresponds to the crystal-fluid coexistence pressure of 11.55$`k_BT\sigma ^3`$. In order to minimize the size effects and the amount of stress in the crystal introduced by the cleaving walls, we use large systems of about 8000 spheres and approximate dimensions of $`14\sigma \times 14\sigma \times 40\sigma `$. (To perform the simulations efficiently for such large systems, we use the algorithm of Rappaport .) The cleaving plane is placed in the middle between two crystal layers. The dependence of the pressure on the wall position is shown in Fig. 2. The walls do not interact with the crystal until they move sufficiently close to the crystal layers (about 0.7$`\sigma `$ for all orientations). Then the pressure on the walls quickly rises to slightly above the bulk crystal pressure. The steepness of the rise is directly related to the compactness of the layers for each orientation. The final positions, $`z_f`$, where the spheres of different types no longer collide across the cleaving plane, were determined to be 0.31$`\sigma `$, 0.16$`\sigma `$, and 0.35$`\sigma `$ for the (100), (110), and (111) system orientations, respectively. No hysteresis was observed in the reverse process.

Step 2: Cleaving the fluid The fluid systems consisting of about 7400 particles are prepared at the coexistence density $`\rho _f=0.939\sigma ^3`$ using box dimensions nearly identical to the crystal systems. Unlike in step 1, the cleaving walls begin to interact with the fluid system as long as $`z<\sigma `$. As can be seen in Fig. 3, the pressure on the walls increases approximately linearly until the fluid near the cleaving walls begins to develop significant crystal-like ordering commensurate with the wall structure. At that point the pressure in the bulk fluid decreases to about 11.2$`k_BT\sigma ^3`$, after which the dependence of pressure on the wall position follows essentially the same curve as during the cleaving of the crystal system, which leads to the same values of $`z_f`$ as in step 1.

Simulation of the reverse process shows that the ordering of the fluid against the walls is the source of some hysteresis. However, we have found that the magnitude of the hysteresis can be always reduced to within the statistical error by increasing the duration of the equilibration run. In other words, at every position of the cleaving walls, the pressure on the walls eventually converges to the same value (within the statistical error bounds) in both forward and reverse processes.

Step 3: Changing boundary conditions The combined system has two cleaving planes and four walls. The crystal part of the system and the two walls restraining it are assigned type 1, while the fluid part and the other two walls are assigned type 2.

Step 4: Removing the cleaving walls As can be seen in Fig. 4, the pressure on the walls restraining crystal and fluid parts of the system is essentially the same as in steps 1 and 2, respectively, except that the fluid part retains its structure in the interfacial region. Thus the main contribution to the interfacial free energy comes from the pressure of the fluid on the cleaving walls before significant crystal-like ordering at the wall develops.

The work done during each step and the resulting interfacial free energy for each of the three orientations is given in Table I. The average values of $`\gamma `$ for the three orientations is about 0.61$`k_BT/\sigma ^2`$, which corresponds to a Turnbull coefficient of 0.51. This average value is about 10% higher than the value determined from nucleation rates on colloidal crystals, consistent with the differences found in other materials, such as bismuth (discussed above). Note that the hard-sphere $`\gamma `$ values are slightly anisotropic and increase in the order of (111), (100) and (110). That (111) has the lowest interfacial free energy is perhaps not surprising, since the (111) crystal face resembles most the structure adopted by the fluid against a structureless wall .

It is generally accepted that the structure and thermodynamics of dense simple fluids is dominated by the repulsive part of the potential, which can often be well approximated as a hard sphere, and the effect of the attractive part of the potential can be viewed as a small perturbation. If one considers the truncated Lennard-Jones system studied by Broughton and Gilmer and calculates an effective hard-sphere diameter at the triple point ($`T=0.617ϵ/k_B`$), using the Barker-Henderson approach, one obtains a value of ($`0.39ϵ/\sigma ^2`$) simply by rescaling the hard-sphere result calculated here. Thus, the attractive part of this potential accounts for only about 10% of the total $`\gamma `$, which is consistent with a similar observation by Curtin on the basis of DFT calculations .

To summarize, we have determined the crystal/melt interfacial free energy, $`\gamma `$ for the hard-sphere system directly from simulation using a method that is similar to that Broughton and Gilmer  used for the truncated Lennard-Jones system except that we have replaced their external cleaving potentials with specially constructed cleaving walls. Although the method of cleaving walls is especially advantageous for the hard-sphere system, it could also be easily applied in modified form to continuous potentials. The hard-sphere $`\gamma `$ obtained is only slightly dependent upon orientation and averages about 0.61$`k_BT/\sigma `$, consistent with some previous theoretical predictions from density-functional theory. This value is also about 10% higher than that determined from nucleation experiments on colloidal suspensions of polystyrene spheres, giving a rare comparison between direct evaluations of $`\gamma `$ and less accurate indirect determinations via nucleation theory.

The authors gratefully acknowledge support from NSF under Grant CHE9970903, as well as the Kansas Center for Advanced Scientific Computing for computational support.


text:
# Charged particle multiplicity near mid-rapidity in central Au+Au collisions at √𝑠= 56 and 130 AGeV

## Abstract

We present the first measurement of pseudorapidity densities of primary charged particles near mid-rapidity in Au+Au collisions at $`\sqrt{s}=`$ 56 and 130 AGeV. For the most central collisions, we find the charged particle pseudorapidity density to be $`dN/d\eta |_{|\eta |<1}=408\pm 12\text{(stat)}\pm 30\text{(syst)}`$ at 56 AGeV and $`555\pm 12\text{(stat)}\pm 35\text{(syst)}`$ at 130 AGeV, values that are higher than any previously observed in nuclear collisions. Compared to proton-antiproton collisions, our data show an increase in the pseudorapidity density per participant by more than 40% at the higher energy.

PACS numbers: 25.75.-q

In June 2000, the Relativistic Heavy-Ion Collider (RHIC) at Brookhaven National Laboratory delivered the first collisions between Au nuclei at the highest center of mass energies achieved in the laboratory to date. In this paper we present data taken with the PHOBOS detector during the first collider run at energies of $`\sqrt{s}`$ = 56 and 130 AGeV. The ultimate goal of our work is to understand the behavior of strongly interacting matter at conditions of extreme density and temperature. Quantum chromodynamics (QCD), the fundamental theory of strong interactions, predicts that for sufficiently high energy density a new state of matter will be formed, the so-called quark-gluon plasma (QGP) . The measurements shown here represent the first step toward the development of a full picture of the dynamical evolution of nucleus-nucleus collisions at RHIC energies.

Studying the dependence of charged particle densities on energy and system size provides information on the interplay between hard parton-parton scattering processes, which can be calculated using perturbative QCD, and soft processes, which are treated by phenomenological models that describe the non-perturbative sector of QCD. Predictions for multi-particle production in high-energy heavy-ion collisions, obtained from a variety of models, typically vary by up to a factor of two .

In this letter we report data for the most central Au+Au collisions detected in our apparatus. We have determined the energy dependence of the density of primary charged particles emitted near $`90^{}`$ to the beam axis, characterized by the pseudorapidity density $`dN/d\eta |_{|\eta |<1}`$, where $`\eta =\mathrm{ln}\mathrm{tan}(\theta /2)`$ and $`\theta `$ is the polar angle from the beam axis. These data provide the first means to constrain models of heavy-ion collisions at RHIC energies. They will allow the extraction of basic information about the initial conditions in these collisions, in particular the energy density, and thus form an essential element for the proper prediction or description of other observables.

The PHOBOS detector employs silicon pad detectors to perform tracking, vertex detection and multiplicity measurements. Details of the setup and the layout of the silicon sensors can be found elsewhere. For the initial running period of the accelerator only a small fraction of the full setup was installed. It included the first 6 layers of the silicon spectrometer (SPEC), part of the two-layer silicon vertex detector (VTX) and one ladder of the large acceptance Octagon multiplicity detector (see Fig. 1). In total, the installed sensors had 20000 readout channels, of which less than 2% were non-functional. The detector setup also included two sets of 16 scintillator counters (“paddle counters”) located at distances of $`3.21`$ m (PN) and $`3.21`$ m (PP) from the nominal interaction point along the beam ($`z`$) axis. These counters subtended pseudorapidities between $`3<|\eta |<4.5`$. They served as the primary event trigger and were used for event selection. Two zero-degree calorimeters (ZDCP,ZDCN) at $`z=\pm 18.5`$ m provided additional information for event selection by measuring the energy deposited by spectator neutrons.

Monte Carlo (MC) simulations of the detector performance were based on the HIJING event generator  and the GEANT 3.21 simulation package, folding in the signal response for scintillator counters and silicon sensors.

In this analysis, we used the data from the SPEC subdetector to reconstruct tracks and to obtain the position, $`z_{vtx}`$, of the collision vertex. The SPEC and VTX subdetectors were then used to count the multiplicity of charged particles near $`\eta =0`$. Their acceptance, covering $`1<\eta <2.5`$, is shown in the insert in Fig. 1. The strong dependence on $`z_{vtx}`$ allows us to check our understanding of acceptance and backgrounds by testing the stability of our results as a function of $`z_{vtx}`$.

The readout trigger selected events based on the coincidence of at least one hit in both PP and PN within a time window of 38 ns. This condition selected most collision events as well as various backgrounds. Offline event selection was achieved by requiring the paddle time difference to be less than 8 ns, corresponding to a maximum displacement of roughly $`\pm 120`$ cm relative to $`z=0`$. Additional rejection of residual background can be achieved by requiring the ZDC time difference to be less than 20 ns. Since the ZDCs only detect spectator neutrons, they are slightly inefficient for the most central events. Thus, we accepted events with good ZDC timing or high multiplicity in the paddles or both. Double beam backgrounds were studied by using experimental runs where the beams did not collide. We found that the overall rate was less than 1% of all events and that there was no background mis-identified as central events, which are characterized by large signals in both paddle counters.

To select the most central events we have used an estimator based on the mean of gain-normalized ADC values in the 16 scintillator counters in PP and PN. Primary charged particles each leave approximately 1.8 MeV in the scintillator. Slow secondary particles that traverse the counters at large angles may deposit larger amounts of energy, and thus mimic a larger multiplicity. To reduce this effect, for each event we discarded the four scintillators in each set of paddle counters with the largest signals. The average of the pulse heights of the remaining 12 scintillators, PP12 and PN12, was then calculated for each set. Finally, we used the paddle mean $`\frac{1}{2}`$ (PP12 +PN12) as an observable proportional to the number of particles hitting the paddle counters.

Assuming that the number of particles hitting the paddle counters increases monotonically with increasing number of participants, a cut on large paddle mean selects the events with the largest number of participants, in both the data and simulation. This has been confirmed with the ZDCs for our more central events. We observe a good anti-correlation between the paddle mean and the total energy deposited in the ZDCs. For this analysis central events were selected by choosing the 6% of events with the largest paddle mean (see Fig. 2). We obtained a sample of 382 events for $`\sqrt{s}=56`$ AGeV and 724 events for $`\sqrt{s}=130`$ AGeV.

We estimated the average number of participants $`N_{part}`$ chosen by the 6% cut using simulations . MC studies showed that up to 10% of the inelastic cross section fails to produce enough particles to satisfy our trigger conditions. We also studied the distribution of $`PP_{12}`$ after a 6% cut on $`PN_{12}`$ and vice-versa. This showed the effect of additional fluctuations introduced by slow secondary particles. Considering both of these effects, we estimate that the systematic uncertainty in $`N_{part}`$ is 5%. Applying the 6% centrality cut to MC events we deduce that the number of participants in the data is $`N_{part}=330\pm 4`$ (stat) for $`\sqrt{s}=`$56 AGeV and $`343\pm 4`$ (stat) for 130 AGeV. The RMS of the distribution of $`N_{part}`$ for the MC events is 27 at 56 AGeV and 25 at 130 AGeV.

The analysis proceeded as follows: Using the information from the SPEC sensor planes a road-following tracking algorithm reconstructed straight line tracks passing through at least 4 layers of the detector and having a $`\chi ^2`$-probability of better than 1% for a non-vertex constrained straight line fit. MC studies show that in this step the tracking efficiency is better than 90%, with less than 10% ghost tracks. For central events, on average 13 tracks were reconstructed at $`\sqrt{s}=56`$ AGeV and 18 tracks at 130 AGeV.

The position of the primary collision vertex was determined as the point of closest approach for the found tracks. The distribution of collision vertices, which define the beam orbit, was found to be very stable within each data set, while an offset of 0.7 mm between the two beam energies was observed. To optimize the stability of the vertex finding algorithm, which is essential for the multiplicity determination, we performed a second track finding analysis of the data by confining the vertex position to within 3 mm distance to the previously determined beam orbit in the transverse direction and within $`25<z_{vtx}<15`$ cm longitudinally. For central events in this region, MC simulations show that the vertex finding algorithm is more than 99% efficient. As a cross-check of the vertex finding procedure we compared the pointing accuracy of the found tracks to the fitted vertex for data and MC and found a most probable value of 1 mm in both cases.

Within the vertex fiducial volume we found a total of 103 central events for the low energy and 151 central events in the high energy data set. For the selected events, the multiplicity of charged particles was determined using the SPEC and VTX subdetectors independently. The measurement was done by counting tracklets, which are two hit combinations in consecutive layers of the SPEC or VTX subdetector consistent with a track originating at the primary vertex. For every hit the coordinates in $`\eta `$ and $`\varphi `$ were calculated relative to the fitted vertex position, where $`\varphi `$ is the azimuthal angle in radians. Then we determined for each hit in the first layer the closest hit in the second layer. Finally, hit combinations with a distance of $`\sqrt{\delta \eta ^2+\delta \varphi ^2}<0.015`$ in the SPEC subdetector and $`|\delta \eta |<0.1`$ in the VTX subdetector were retained as tracklets. To obtain a constraint on possible random channel-by-channel inefficiencies in the sensors and readout-chain, we compared the multiplicities of single hits in VTX and SPEC, two-hit tracklets in VTX and SPEC and tracks in SPEC. A detection efficiency of less than 100% would enter in different powers into the different measurements. We find no evidence for inefficiencies beyond the $`<2`$% of dead channels identified in our standard calibration procedure.

MC studies show that the number of tracklets per event near $`\eta =0`$ is approximately proportional to the number of primary charged particles per event. Here we define as primary particles all charged hadrons produced in the collision, including the products of strong and electromagnetic decays, but excluding feed-down products from weak decays and hadrons produced in secondary interactions.

The pseudorapidity density of primary charged particles, $`dN/d\eta |_{|\eta |<1}`$, was obtained from the observed distribution of tracklets, $`dN/d\eta (\text{tracklets})`$, by accounting for several factors: The efficiency of tracklet finding, background due to noise hits, secondary and feed-down particles and the acceptance. The correction was done using the following procedure: Events from a standard event generator were propagated through the GEANT simulation of our setup. The simulated data were then processed with the identical reconstruction chain as the real data, including the event selection based on the paddle counter signals. For each simulated event we determined the number of tracklets, $`N_{tracklets}`$, found in the region $`|\eta |<1`$ and the number of primary charged particles $`N_{primaries}`$ in the same pseudorapidity region. We then determined the vertex dependent proportionality factor $`\alpha (z_{vtx})`$ by calculating $`N_{tracklets}/N_{primaries}`$ as a function of $`z_{vtx}`$.

We repeated this procedure using both the HIJING and VENUS  event generators, which predict multiplicities that differ by a factor 2-3. We also compared $`\alpha (z_{vtx})`$ obtained for running the event generators at various different energies. For all MC runs, $`\alpha `$ agreed to better than 5%, demonstrating that the tracklet counting procedure provides a robust measure of the primary charged particle pseudorapidity density.

MC studies showed that less than 5% of the particles emitted into the angular acceptance of the SPEC and VTX detectors stop in the material of the beam pipe or the first detector layer. The multiplicities we report are corrected for this missing fraction based on the HIJING momentum distributions.

We further studied the contamination of the tracklet distribution by feed-down products from weak decays of strange particles. Due to the proximity of our detectors to the beamline and the good pointing accuracy in the tracklet reconstruction, the contribution was found to be small ($`<4\%`$). Again, the multiplicities reported here are corrected based on the HIJING distributions.

Fig. 3 shows a direct comparison of the SPEC (left) and VTX (right) tracklet $`dN/d\eta `$ distributions for data (symbols) and MC events (solid lines), normalized per event. Scaling factors of 1.15 for $`\sqrt{s}=56`$ AGeV and 0.98 for $`\sqrt{s}=130`$ AGeV were applied to the MC distribution to match the integrals to the data. The shape of the distribution agrees well between simulation and data. SPEC and VTX distributions are both consistent with the same scaling factor. Based on the MC studies, the comparison of data and MC tracklet distributions and the comparison of results from the SPEC and VTX tracklet analysis, we conclude that the proportionality factors can be applied to the measured tracklet distributions with an overall systematic uncertainty of less than 8%.

As the result of this procedure we obtain a primary charged particle density of $`dN/d\eta |_{|\eta |<1}=408\pm 12\text{(stat)}\pm 30\text{(syst)}`$ for $`\sqrt{s}=56`$ AGeV and $`555\pm 12\text{(stat)}\pm 35\text{(syst)}`$ at $`\sqrt{s}=130`$ AGeV. From the simulation of the paddle counter selection we obtain for the mean number of participating nucleons $`N_{part}=330\pm 4\text{(stat)}_{15}^{+10}\text{(syst)}`$ for $`\sqrt{s}=56`$ AGeV and $`343\pm 4\text{(stat)}_{14}^{+7}\text{(syst)}`$ for $`\sqrt{s}=130`$ AGeV.

Normalizing per participant pair, we deduce $`dN/d\eta |_{|\eta |<1}/0.5N_{part}=2.47\pm 0.1\text{(stat)}\pm 0.25\text{(syst)}`$ and $`3.24\pm 0.1\text{(stat)}\pm 0.25\text{(syst)}`$, respectively.

Finally, taking the strong correlation between the systematic errors at the two energies into account, we obtain an increase in the charged particle density per participant by a factor of $`1.31\pm 0.04\text{(stat)}\pm 0.05\text{(syst)}`$.

In Fig. 4 we show the normalized yield per participant obtained for Au+Au collisions, proton-antiproton ($`p\overline{p}`$) collisions  and central Pb+Pb collisions at the CERN SPS . The $`dN/d\eta `$ value for the Pb+Pb data was obtained by numerically integrating the momentum distributions shown in .

Several important features of the data emerge: First, the central Au+Au collisions show a significantly larger charged particle density per participant than for example non-single diffractive (NSD) $`p\overline{p}`$ collisions at comparable energies. This rules out simple superposition models such as the wounded nucleon model  and is compatible with predictions of models like HIJING that include particle production via hard-scattering processes.

Secondly, the observed increase by 31% from 56 to 130 AGeV in central Au+Au collisions is significantly steeper than the increase shown by a $`p\overline{p}`$ parametrization (see Fig. 4) for the same energy interval . Finally, comparing our data to those obtained at the CERN SPS for Pb+Pb collisions at $`\sqrt{s}=17.8`$ AGeV, we find a 70% higher particle density per participant near $`\eta =0`$ at $`\sqrt{s}=130`$ AGeV. General arguments (c.f. Bjorken’s estimate ) suggest that this increase should correspond to a similar increase in the maximal energy density achieved in the collision.

Acknowledgements: We acknowledge the generous support of the entire RHIC project personnel, C-A and Chemistry Departments at BNL. We thank Fermilab and CERN for help in silicon detector assembly. We thank the MIT School of Science and LNS for financial support. This work was partially supported by US DoE grants DE-AC02-98CH10886, DE-FG02-93ER40802, DE-FC02-94ER40818, DE-FG02-94ER40865, DE-FG02-99ER41099, W-31-109-ENG-38. NSF grants 9603486, 9722606 and 0072204. The Polish groups were partially supported by KBN grant 2 P03B 04916. The NCU group was partially supported by NSC of Taiwan under contract NSC 89-2112-M-008-024.