text:
# Observation of a Narrow Resonance of Mass 2.46 GeV/c2 Decaying to 𝑫_𝒔^{∗+}𝝅^𝟎 and Confirmation of the 𝑫_{𝒔𝑱}^∗(𝟐𝟑𝟏𝟕) State

## I Introduction

The BaBar Collaboration has recently reported babar  evidence for a new narrow resonance with a mass near 2.32 GeV/c2, which decays to $`D_s^+\pi ^0`$. The BaBar data are consistent with the identification of this state as one of the four lowest-lying members of the $`c\overline{s}`$ system with orbital angular momentum $`L=1`$, and provisionally it has been named the $`D_{sJ}^{}(2317)`$ meson. A natural candidate would be the $`{}_{}{}^{3}P_{0}^{}`$ $`c\overline{s}`$ meson with spin-parity $`J^P=0^+`$, but other possibilities, including exotic states, are not ruled out. In this paper, we report on a search for the $`D_{sJ}^{}(2317)`$ meson, as well as other, possibly related states, in data collected with the CLEO II detector in symmetric $`e^+e^{}`$ collisions at the Cornell Electron Storage Ring, at center-of-mass energies $`\sqrt{s}10.6`$ GeV.

The spectroscopy of $`P`$-wave $`c\overline{s}`$ mesons is summarized in Ref. bartelt . Prior to the observation of the $`D_{sJ}^{}(2317)`$, theoretical expectations rgg ; gi ; iw ; gk ; dipe  were that: (1) all four states with $`L=1`$ are massive enough that their dominant strong decays would be to the isospin-conserving $`DK`$ and/or $`D^{}K`$ final states, (2) the singlet and triplet $`J^P=1^+`$ states could mix, and (3) in the heavy quark limit, the two states with $`j=3/2`$ would be narrow while the two with $`j=1/2`$ would be broad, where $`j`$ is the sum of the strange quark spin and the orbital angular momentum. Existing experimental evidence pdg ; cleods2  for the narrow $`D_{s1}(2536)`$ and $`D_{sJ}^{}(2573)`$ mesons which decay dominantly to $`D^{}K`$ and $`DK`$ respectively, and the compatibility of the $`D_{sJ}^{}(2573)`$ with the $`J^P`$ assignment as $`2^+`$ support this picture.

The observation by BaBar babar  of the new state at a mass of 2.32 GeV is surprising because: (1) it is narrow (with intrinsic width $`\mathrm{\Gamma }<10`$MeV), (2) it has been observed in the isospin-violating $`D_s\pi ^0`$ channel, and (3) its mass ($`2316.8\pm 0.4\text{[stat.]}`$MeV/c2) is smaller than most theoretical predictions for a $`0^+`$ $`c\overline{s}`$ state that could decay via this channel. However, points (1) and (2) would be obvious consequences of the low mass, since the $`D^{()}K`$ decay modes are not allowed kinematically. We also note that at least two theoretical calculations nrz ; bh  prior to the $`D_{sJ}^{}(2317)^+`$ observation had suggested that, in the heavy quark limit, the $`j=1/2`$ states with $`J^P=0^+`$ and $`1^+`$ could be thought of as chiral partners of the $`D_s`$ and $`D_s^{}`$ mesons, and thus would be relatively light. In one model bh  it was proposed that the mass splittings between the $`0^+`$ and $`0^{}`$ states of heavy flavored mesons could be as small as 338 MeV/c2, which is near the $`D_{sJ}^{}(2317)^+D_s^+`$ mass splitting of $`348.3`$MeV/c2 measured by BaBar.

Since the initial observation, a number of explanations have appeared cahnjackson ; barnescloselipkin ; beverenrupp ; bardeeneichtenhill ; szczepaniak ; godfrey ; colangelodefazio ; bali . Cahn and Jackson cahnjackson  apply non-relativistic vector and scalar exchange forces to the constituent quarks. Barnes, Close and Lipkin barnescloselipkin  consider a quark model explanation unlikely and propose a $`DK`$ molecular state. Similarly, Szczepaniak szczepaniak  suggests a $`D\pi `$ atom. Also going beyond a simple quark model description, Van Beveren and Rupp beverenrupp  present arguments for a low mass $`0^+`$ $`c\overline{s}`$ state based on a unitarized meson model, by analogy with members of the light scalar meson nonet. Bali bali  reports on lattice QCD calculations that predict signficantly larger $`0^+0^{}`$ meson mass splittings than what has been observed for the $`D_{sJ}^{}(2317)D_s`$ splitting.

On the contrary, Bardeen, Eichten and Hill bardeeneichtenhill  couple chiral perturbation theory with a quark model representation in heavy quark effective theory, building on the model described in Ref. bh . They infer that the $`D_{sJ}^{}(2317)`$ is indeed the $`0^+`$ $`c\overline{s}`$ state expected in the quark model, predict the existence of the $`1^+`$ partner of this state with a $`1^+1^{}`$ mass splitting equal to the $`0^+0^{}`$ mass splitting, and compute the partial widths for decays to allowed final states. Godfrey godfrey  and Colangelo and De Fazio colangelodefazio  find that the radiative transistion of the $`D_{sJ}^{}(2317)`$ should be significant if it is indeed a $`c\overline{s}`$ state.

The goals of the analysis presented here are to use CLEO data to provide independent evidence regarding the existence of the $`D_{sJ}^{}(2317)`$, to shed additional light on its properties, and to search for decays of other new, possibly related states. In particular, we address the following questions. Are the electromagnetic decays $`D_s\gamma `$ or $`D_s^{}\gamma `$ observable in light of the isospin suppression of the strong decay to $`D_s\pi ^0`$? Are other strong decays observable such as $`D_s^{}\pi ^0`$, or the isospin-conserving but Okubo-Zweig-Iizuka (OZI) suppressed ozi  decay $`D_s\pi ^+\pi ^{}`$? If the $`D_{sJ}^{}(2317)`$ is the expected $`0^+`$ $`c\overline{s}`$ state, might the remaining $`1^+`$ state also be below threshold for decay to $`D^{}K`$, as suggested in Ref. bardeeneichtenhill , and thus be narrow enough to be observable in its decays to $`D_s^{}\pi ^0`$, $`D_s\gamma `$ or $`D_s^{}\gamma `$?

This article cleohepex  is organized as follows. After describing the detector and data set in Section II, we summarize the reconstruction of the $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ decay channel in Section III, including efforts to understand and exclude contributions from known background processes. We then report in Section IV on searches for other possible decay channels as described in the preceding paragraph. In Section V, we report on the appearance of a statistically significant signal in the $`D_s^+\pi ^0`$ channel at a mass of 2.463 GeV/c2, not compatible with a kinematic reflection of the $`D_{sJ}^{}(2317)^+`$. We describe a quantitative analysis of the signals in the $`D_s^+\pi ^0`$ and $`D_s^+\pi ^0`$ channels, leading us to infer the existence of two distinct states. Based on this conclusion, we discuss the properties of these two states in Section VI, after which we summarize the principal results of the analysis.

## II Detector and Data Set

The analysis described here is based on 13.5 fb-1 of $`e^+e^{}`$ collision data collected between 1990 and 1998. CLEO II is a general purpose, large solid angle, cylindrical detector featuring precision charged particle tracking and electromagnetic calorimetry, and is described in detail in Refs. Kubota:1992ww ; Hill:1998ea . In its initial configuration, the tracking system was comprised of a six-layer straw tube chamber just outside of a 3.2 cm radius beryllium beam pipe, followed by a 10 layer hexagonal cell drift chamber and a 51 layer square cell drift chamber, immersed in a $`1.5`$T magnetic field generated by a superconducting solenoid. In 1995, the beam pipe and straw tubes were replaced by a 2.0 cm radius beam pipe plus three layers of silicon strip detectors each with double-sided readout, and a helium-propane gas mixture replaced the argon-ethane mixture previously used in the main drift chamber.

Beyond the tracking system, but within the solenoid, were also located a 5 cm thick plastic scintillation counter system for time-of-flight measurement and triggering, and a barrel calorimeter consisting of 6144 tapered CsI(Tl) crystals 30 cm in length, arrayed in a projective geometry, with their long axis oriented radially with respect to the $`e^+e^{}`$ interaction point. An additional 1656 crystals were deployed in two end caps to complete the solid angle coverage. The excellent energy and angular resolution of the calorimeter is critical for the reconstruction of $`\pi ^0\gamma \gamma `$ decays as well as single low-energy photons such as those emitted in the $`D_s^+D_s\gamma `$ transition.

## III Confirmation of $`𝑫_{𝒔𝑱}^{\mathbf{}}\mathbf{(}\mathrm{𝟐𝟑𝟏𝟕}\mathbf{)}^\mathbf{+}\mathbf{}𝑫_𝒔^\mathbf{+}𝝅^\mathrm{𝟎}`$

The search for the $`D_{sJ}^{}(2317)`$ was carried out by reconstructing the $`D_s^+\pi ^0`$ state, using the $`D_s^+\varphi \pi ^+`$ channel with $`\varphi K^+K^{}`$. Charge conjugation is implied throughout this article. Pairs of oppositely charged tracks were considered as candidates for the decay products of the $`\varphi `$ if the specific ionization ($`dE/dx`$) was measured in the main drift chamber to be within 2.5 standard deviations of the expectation for a kaon, and if the invariant mass of the $`K^+K^{}`$ system was within $`\pm 10`$MeV/c2 of the $`\varphi `$ mass. A third track with $`dE/dx`$ consistent with the expectation for a pion was combined with the $`K^+K^{}`$ system to form a $`D_s^+`$ candidate with mass $`M(KK\pi )`$. To improve resolution we adjust the momenta of the three particles subject to the constraint that their trajectories intersect at a common point corresponding to the decay point of a $`D_s`$ meson. When fitted to a Gaussian, the observed $`D_s^+`$ mass peak has a standard deviation ($`\sigma `$) of $`6.5\pm 0.4`$MeV/c2 in our data, consistent with CLEO Monte Carlo simulations of $`D_s`$ production and decay plus a GEANT-3 geant  based simulation of particle propagation and detector response.

Clusters of energy deposition in the calorimeter unassociated with charged particle interactions were identified as potential photon candidates. To be considered as candidates for the photons from $`\pi ^0\gamma \gamma `$ decay, clusters with energy greater than 100 MeV located in the central region of the calorimeter ($`|\mathrm{cos}\theta |<0.71`$, where $`\theta `$ is measured with respect to the beam axis) were selected. Pairs of photons were required to satisfy $`3.0<[M(\gamma \gamma )M_{\pi ^0}]/\sigma (\gamma \gamma )<2.5`$, where $`M(\gamma \gamma )`$ is the invariant mass of the two photons and $`\sigma (\gamma \gamma )`$ is the expected resolution on this mass. For each cluster being considered as a photon candidate, we additionally required that the lateral profile of energy deposition in the calorimeter be consistent, at the $`99\%`$ confidence level, with expectations for photons. This requirement removes spurious photon candidates that are mainly due to inelastic interactions of charged hadrons or long-lived neutral hadrons. The peak in the $`M(\gamma \gamma )`$ distribution for photon-pairs accompanying a $`D_s^+`$ candidate with $`M(D_s)=M(KK\pi )`$ between 1.9565 and 1.9805 GeV/c2 has $`\sigma =5.8\pm 0.4`$MeV/c2 in our data, consistent with expectations from the Monte Carlo simulations. Once identified as a $`\pi ^0`$ candidate, the directions and energies of the two photons are adjusted with a kinematic fit to reconstruct to the known value pdg  for the $`\pi ^0`$ mass $`M_{\pi ^0}`$.

To suppress combinatoric backgrounds, we further required that the momentum of the $`D_s^+\pi ^0`$ candidate be greater than $`3.5`$GeV/c. We also required that the helicity angle of the $`\varphi K^+K^{}`$ decay satisfy the requirement $`|\mathrm{cos}\theta _h|>0.3`$, where $`\theta _h`$ is the angle between the $`K^+`$ momentum vector measured in the $`\varphi `$ rest frame, and the $`\varphi `$ momentum vector measured in the $`D_s`$ rest frame. The expected distribution from real $`\varphi `$ decays varies as $`\mathrm{cos}^2\theta _h`$, whereas combinatoric backgrounds tend to be flat. For $`D_s\pi ^0`$ combinations satisfying these requirements, we plot the mass $`M(D_s\pi ^0)=M(KK\pi \pi ^0)`$ and the mass difference $`\mathrm{\Delta }M(D_s\pi ^0)=M(D_s\pi ^0)M(D_s)`$ in Fig. 1(a) and (b), respectively.

To improve the experimental resolution on $`M(D_s\pi ^0)`$, the known value of the $`D_s`$ mass, $`M_{D_s}=1968.5\pm 0.6`$MeV/c2 pdg , has been used to determine the energy of the $`KK\pi `$ system from its measured momentum in Fig. 1(a); this substitution is not done for $`\mathrm{\Delta }M(D_s\pi ^0)`$ in Fig. 1(b), or for the calculation of other mass differences entering this analysis.

The narrow peaks in Fig. 1 at a mass near 2.32 GeV/c2 and a $`\mathrm{\Delta }M(D_s\pi ^0)`$ near 350 MeV/c2 are in qualitative agreement with the BaBar observation. We note that there are no peaks in this region when $`KK\pi `$ combinations with $`M(KK\pi )`$ lying in $`D_s`$ side band regions are combined with a $`\pi ^0`$. The other feature of note in the spectra is the sharp signal from $`D_s^+D_s^+\pi ^0`$ cleodspi0  near the kinematic threshold. In addition, Monte Carlo simulations of inclusive multi-hadron production via $`e^+e^{}q\overline{q}`$ ($`q=u,d,s,c`$) give $`M(D_s\pi ^0)`$ and $`\mathrm{\Delta }M(D_s\pi ^0)`$ spectra that reproduce the features observed in the data, except for the peaks near 2.32 GeV/c2 and 350 MeV/c2 in the respective plots. This is also illustrated in Fig. 1, where the normalization for the $`q\overline{q}`$ Monte Carlo spectra is fixed by the ratio of the luminosity of the data sample to the equivalent luminosity of the Monte Carlo sample. This normalization is known to a precision of approximately $`\pm 5\%`$.

The agreement between the Monte Carlo and data distributions in Fig. 1 in normalization as well as shape demonstrates that the simulation of ‘random’ photons accompanying $`D_s`$ decays is accurate. The accuracy of this simulation is important for our detailed analysis of this signal, described in Section V.1.

We have investigated mechanisms by which a peak at 2.32 GeV/c2 could be generated from decays involving known particles, either through the addition, omission or substitution of a pion or photon, or through the mis-assignment of particle masses to the observed charged particles. In no cases were narrow enhancements in the $`M(D_s\pi ^0)`$ spectrum near 2.32 GeV/c2 observed. We will discuss the issue of backgrounds from a new resonance at 2.46 GeV/c2 when we describe our studies of the $`D_s^+\pi ^0`$ final state.

From a binned maximum likelihood fit of the $`\mathrm{\Delta }M(D_s\pi ^0)`$ distribution to a Gaussian signal shape and second-order polynomial background function, we obtain a yield of $`165\pm 20`$ events in the peak near 350 MeV/c2. In this fit, the mean and Gaussian width of the peak are allowed to float. These parameters are determined to be $`\mathrm{\Delta }M(D_s\pi ^0)=349.4\pm 1.0`$ MeV/c2 and $`\sigma =8.0_{1.1}^{+1.3}`$MeV/c2, where the errors are due to statistics only. The peak is somewhat broader than the expected mass resolution of $`6.0\pm 0.3`$MeV/c2, determined from Monte Carlo simulations. The detection efficiency associated with the reconstruction of the full $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$, $`D_s^+\varphi \pi ^+`$, $`\varphi K^+K^{}`$ decay chain is $`(9.73\pm 0.57)\%`$ for the portion of the $`D_{sJ}^{}(2317)^+`$ momentum spectrum above $`3.5`$GeV/c, where this efficiency does not include the $`D_s`$ and $`\varphi `$ decay branching fractions.

Thus, we confirm the existence of a peak in the $`D_s\pi ^0`$ mass spectrum that cannot be explained as reflections from decays of known particles. Our measurements of the mean mass difference and width of the peak are consistent with the values obtained by BaBar babar  for the $`D_{sJ}^{}(2317)^+`$. Further discussions of the width, as well as of systematic errors in the measurements of the mass and width of the $`D_{sJ}^{}(2317)`$ appear later in this article.

## IV Searches for $`𝑫_{𝒔𝑱}^{\mathbf{}}\mathbf{(}\mathrm{𝟐𝟑𝟏𝟕}\mathbf{)}`$ in other channels

The conclusion that the $`D_{sJ}^{}(2317)`$ is a new narrow resonance decaying to $`D_s\pi ^0`$ leads to two questions: (1) are there other observable decay modes, and (2) might additional new $`c\overline{s}`$ resonances also exist in which normally suppressed decay modes such as $`D_s^{()}\pi ^0`$ are dominant? To answer these questions we have searched in the channels $`D_s\gamma `$, $`D_s^{}\gamma `$, $`D_s^{}\pi ^0`$ and $`D_s\pi ^+\pi ^{}`$.

If the $`D_{sJ}^{}(2317)`$ is a $`0^+`$ $`L=1`$ $`c\overline{s}`$ meson, as has been suggested bardeeneichtenhill , it could decay via $`S`$\- or $`D`$-wave to $`D_s^{}\gamma `$, but would not be able to decay to $`D_s\gamma `$ due to parity and angular momentum conservation. Consequently, observation of one or both of these channels would be interesting. On the other hand, if neither channel is seen, this would not be too surprising since these are electromagnetic decays, and the $`D_s\pi ^0`$ decay, while isospin violating, is not as severely phase-space suppressed as in the case of the corresponding decay of the $`D_s^{}`$ where the electromagnetic decay dominates. The BaBar data show no evidence for either channel, however no upper limits were reported on the branching ratios for these channels.

With regard to strong decays, the $`D_s\pi ^+\pi ^{}`$ final state is kinematically allowed and isospin conserving, but would be suppressed by the OZI rule. This is in contrast to the $`D_s\pi ^0`$ channel for which one mechanism would be decay to a $`D_s`$ plus a virtual $`\eta `$, with production of the $`\pi ^0`$ via $`\eta `$-$`\pi ^0`$ mixing cho . However, angular momentum and parity conservation forbid the decay of a $`0^+`$ state to three pseudoscalars. Thus, observation of the $`D_s\pi ^+\pi ^{}`$ channel would be strong evidence against the interpretation of the $`D_{sJ}^{}(2317)`$ as a $`0^+`$ meson.

Finally, it is possible that the remaining $`L=1`$ $`c\overline{s}`$ state with $`J^P=1^+`$ could also be light enough that decays to $`D^{}K`$ would be kinematically forbidden. In this case, the strong isospin-violating decay of this $`1^+`$ state to $`D_s^{}\pi ^0`$ could occur via $`S`$-wave (the electromagnetic decays to $`D_s\gamma `$ or $`D_s^{}\gamma `$ would also be possible), and thus a narrow peak in the $`\mathrm{\Delta }M(D_s^{}\pi ^0)=M(D_s^{}\pi ^0)M(D_s^{})`$ spectrum would be a signature of such a state.

### IV.1 Searches for $`𝑫_{𝒔𝑱}^{\mathbf{}}\mathbf{(}\mathrm{𝟐𝟑𝟏𝟕}\mathbf{)}^\mathbf{+}`$ decays to $`𝑫_𝒔^\mathbf{+}𝝅^\mathbf{+}𝝅^{\mathbf{}}`$, $`𝑫_𝒔^\mathbf{+}𝜸`$ and $`𝑫_𝒔^\mathbf{}\mathbf{+}𝜸`$

To look for these channels we select events containing $`D_s^+\varphi \pi ^+`$ candidates as in the $`D_s\pi ^0`$ analysis. For the $`D_s\pi ^+\pi ^{}`$ channel, we combine the $`D_s`$ candidates with two oppositely charged tracks, and plot the mass difference $`\mathrm{\Delta }M(D_s\pi \pi )=M(D_s\pi \pi )M(D_s)`$. As shown in Fig. 2, no signal is evident in the vicinity of 350 MeV/c2.

To search for states decaying to $`D_s^+\gamma `$, we have formed $`D_s^+\gamma `$ combinations by selecting photons of energy greater than 150 MeV. To select $`D_s^+`$ candidates for use in other searches, we relax this to include photon candidates with energy above 50 MeV. We ignore photons that can be paired with another photon such that $`M(\gamma \gamma )`$ is consistent with $`\pi ^0`$ decay. The inclusive $`\mathrm{\Delta }M(D_s\gamma )=M(D_s\gamma )M(D_s)`$ spectrum for this sample is plotted in Fig. 3(a), illustrating that a large $`D_s^{}`$ sample can be obtained. For decay modes with a $`D_s^{}`$ in the final state, we select $`D_s\gamma `$ combinations where the mass difference $`\mathrm{\Delta }M(D_s\gamma )`$ is reconstructed to be between 0.1308 and 0.1568 GeV/c2.

Also visible in Fig. 3(a), are regions of the $`\mathrm{\Delta }M(D_s\gamma )`$ spectrum where decays of the $`D_{sJ}^{}(2317)`$ (or of a possible higher mass state) into $`D_s\gamma `$ would appear. There is no evidence for a signal near 350 MeV/c2 corresponding to a $`M(D_s\gamma )`$ in the vicinity of 2.32 GeV/c2.

The same conclusion holds for the $`D_s^{}\gamma `$ final state, shown in Fig. 3(b), where we combine selected $`D_s^{}`$ candidates with photons of energy above 150 MeV. The peak in the $`\mathrm{\Delta }M(D_s^{}\gamma )`$ spectrum in Fig. 3(b) near 150 MeV/c2 is due to real $`D_s^+D_s^+\gamma `$ decays in which a random photon has been combined with the $`D_s^+`$ candidate to form the $`D_s^{}`$ candidate, and the actual photon from this transition is combined with this system to form the $`D_{sJ}^{}`$ candidate. There is no sign of any structure in this spectrum near 205 MeV/c2, where a signal from $`D_{sJ}^{}(2317)`$ decay would be expected.

### IV.2 Search for $`𝑫_{𝒔𝑱}^{\mathbf{}}\mathbf{(}\mathrm{𝟐𝟑𝟏𝟕}\mathbf{)}^\mathbf{+}`$ decays to $`𝑫_𝒔^\mathbf{}\mathbf{+}𝝅^\mathrm{𝟎}`$

We have also searched in the $`D_s^+\pi ^0`$ channel for $`D_{sJ}^{}`$ states. To maintain efficiency for this final state, we do not veto $`D_s^+`$ candidates where the photon used in the $`D_s^+`$ reconstruction can be combined with an extra photon to form a $`\pi ^0`$ decay candidate. We also applied slightly less restrictive track quality and shower shape criteria than in the $`D_s\pi ^0`$ analysis. As with the modes involving $`D_s^{}`$ candidates described the preceding section, the energy of photons selected for reconstruction of the $`D_s^{}D_s\gamma `$ decay is required to satisfy $`E_\gamma >50`$MeV. The $`D_s^{}\pi ^0`$ candidates are required to have momenta above 3.5 GeV/c. Fig. 4(a) shows the mass difference plot for events with candidate $`D_s^+\varphi \pi ^+`$, $`D_s^+\gamma D_s^+`$ decays plus di-photon combinations consistent with $`\pi ^0`$ decay.

If the $`D_{sJ}^{}(2317)^+`$ were to decay to the $`D_s^+\pi ^0`$ final state, a peak would be expected at a $`\mathrm{\Delta }M(D_s^{}\pi ^0)205`$MeV/c2. Although we see no evidence for such a peak, there is a significant excess in a narrow region near 350 MeV/c2. We discuss the properties of this new peak in the following section.

## V Observation of a new state at $`\mathbf{2.463}\text{GeV/c}^\mathrm{𝟐}`$

From a fit to a signal Gaussian signal function plus a polynomial background function, we observe a peak in Fig. 4(a), comprised of $`55\pm 10`$ $`D_s^{}\pi ^0`$ combinations, at $`\mathrm{\Delta }M(D_s^{}\pi ^0)=349.8\pm 1.3`$ MeV/c2. The fit yields a Gaussian width of $`6.1\pm 1.0`$MeV/c2 for the peak, consistent with our mass resolution of $`6.6\pm 0.5`$MeV/c2. The existence of this peak leads us to investigate the possibility of a second narrow resonance with a mass near 2.46 GeV/c2 that decays to $`D_s^+\pi ^0`$. We note that a similar peak is also present in the $`M(D_s^{}\pi ^0)`$ spectrum observed by BaBar babar , although BaBar does not claim this as evidence for a new state. For ease of notation, we refer to the postulated particle as the $`D_{sJ}(2463)^+`$.

### V.1 Analysis of cross feed between $`𝑫_𝒔^\mathbf{+}𝝅^\mathrm{𝟎}`$ and $`𝑫_𝒔^\mathbf{}\mathbf{+}𝝅^\mathrm{𝟎}`$ samples

The kinematics of the $`D_s^+\pi ^0`$ and $`D_s^+\pi ^0`$ decays are quite similar, and it is possible that they can reflect into one another. For example, by ignoring the photon from the $`D_s^{}`$ decay in $`D_{sJ}(2463)^+D_s^+\pi ^0`$ decays, nearly all the putative signal combinations form a peak in the $`\mathrm{\Delta }M(D_s\pi ^0)`$ spectrum in the same region as the $`D_{sJ}^{}(2317)`$ signal described in previous sections of this article. We refer to the background entering via this scenario in the $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ sample as ‘feed down’ from the $`D_{sJ}(2463)`$ state. The impact of neglecting the photon is that this peak in $`\mathrm{\Delta }M(D_s\pi ^0)`$ is broader than that for real $`D_{sJ}^{}(2317)`$ decays. From Monte Carlo simulations, we determine the width of this smeared peak to be $`\sigma =14.9\pm 0.4`$ MeV/c2.

It is also possible that a $`D_s^+\pi ^0`$ candidate can be combined with a random photon such that the $`D_s^+\gamma `$ combination accidentally falls in the $`D_s^+`$ signal region defined earlier. In this case, $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ decays would reflect or ‘feed up’ into the $`D_{sJ}(2463)^+D_s^+\pi ^0`$ signal region. A Monte Carlo simulation of $`D_{sJ}^{}(2317)^+`$ production and decay to $`D_s^+\pi ^0`$ shows that this does happen, but only for approximately $`9\%`$ of the reconstructed decays. The peak in the $`\mathrm{\Delta }M(D_s^{}\pi ^0)`$ distribution generated by this feed up sample is also broadened relative to the expectation for real $`D_{sJ}(2463)`$ decays, analogous to the smearing of the feed down kinematics discussed in the preceding paragraph.

We can extract the number of real $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ decays reconstructed in our data, denoted as $`R_0`$, as well as the number of real $`D_{sJ}(2463)^+D_s^+\pi ^0`$ decays, denoted as $`R_1`$, taking into account that the corresponding real signal decays in one channel can enter the candidate sample for the other channel as described above. The following linear equations relate the real to observed numbers:

$$N_0=R_0+f_1R_1$$
(1)
$$N_1=R_1+f_0R_0,$$
(2)
where $`N_0`$ and $`N_1`$ are the numbers of observed decays in the $`D_s\pi ^0`$ and $`D_s^{}\pi ^0`$ channels respectively, and $`R_0`$ and $`R_1`$ are the number of real decays produced times the efficiency to observe them in the corresponding signal decay channels. The coefficients $`f_0`$ and $`f_1`$ are the feed up and feed down probabilities relative to the reconstruction efficiency for the respective signal modes. We note that these relations represent first-order approximations; higher-order corrections, such as that due to the scenario where the $`D_s\pi ^0`$ system from a real $`D_{sJ}(2463)`$ decay is combined with an unrelated photon to form a feed up $`D_{sJ}(2463)`$ candidate, are negligible in the present case.

The observed number of decays in the $`D_s^{}\pi ^0`$ channel is $`N_1=55\pm 10`$, obtained from the fit to the peak in Fig. 4(a) described above. For $`N_0`$, it is desirable to obtain a $`D_s\pi ^0`$ sample selected with criteria that most closely match those used to select $`D_s^{}\pi ^0`$ combinations, and that is enriched in $`D_{sJ}^{}(2317)`$ decays relative to feed down from $`D_{sJ}(2463)`$ decays. Thus we apply the same selection criteria that were used for the $`D_s^{}\pi ^0`$ sample, but without selecting the photon from the $`D_s^{}D_s\gamma `$ transition. To measure the event yield in this sample, we fit the peak in the $`\mathrm{\Delta }M(D_s\pi ^0)`$ distribution to a Gaussian with its width fixed to the Monte Carlo expectation for $`D_{sJ}^{}(2317)`$ decays. In this fit, a significant fraction of feed down combinations is counted as part of the combinatoric background rather than as signal. We obtain $`N_0=190\pm 19`$ candidates. This sample effectively constitutes the source of potential feed up candidates. The difference between this yield and the $`165\pm 20`$ events reported in Section III is consistent with the different acceptances for the two sets of selection criteria.

From the Monte Carlo simulations we measure $`f_0=0.091\pm 0.007\pm 0.015`$ for the probability that a reconstructed $`D_{sJ}^{}(2317)D_s\pi ^0`$ can be combined with a random photon to mimic a $`D_{sJ}(2463)D_s^{}\pi ^0`$ decay. The first error is due to limited Monte Carlo statistics while the second is due to systematic uncertainties associated with: (1) the modeling of extra photons in the simulations, and (2) the fraction of such combinations that are counted by the fit to the $`\mathrm{\Delta }M(D_s^{}\pi ^0)`$ distribution as contributing to the Gaussian signal. As indicated above, this fraction counted by the fit is less than one due to the smearing of $`\mathrm{\Delta }M(D_s^{}\pi ^0)`$ that results when an unrelated photon is added to the $`D_s\pi ^0`$ system. The agreement between the data and Monte Carlo distributions in Fig. 1 lends confidence in the modeling of extra photons. We assign a relative systematic uncertainty of $`5\%`$ based on this and on studies of combinations entering $`D_s^{}`$ side bands described in the following section. To study the second source of systematic uncertainty, we have carried out fits to the $`\mathrm{\Delta }M(D_s^{}\pi ^0)`$ distribution in which the width of the Gaussian signal function was fixed to $`\pm 1\sigma `$ relative to the central value obtained from the nominal fit. Based on the resulting variation in event yields, we have estimated a relative uncertainty on $`f_0`$ of $`16\%`$ from this source.

We also obtain $`f_1=0.84\pm 0.04\pm 0.10`$ from Monte Carlo simulations, where the first error is statistical and the second is due to systematic errors. This includes the probability of feed down as well as the photon finding efficiency. If all $`D_{sJ}(2463)D_s^{}\pi ^0`$ decays with a reconstructed $`D_s`$ plus $`\pi ^0`$ combination were to be counted as $`D_{sJ}^{}(2317)`$ decays, $`f_1`$ would simply be one divided by the acceptance for finding the photon from the $`D_s^{}D_s\gamma `$ transition. However, because the $`\mathrm{\Delta }M(D_s\pi ^0)`$ distribution for the feed down background is broadened, a significant fraction of these combinations are not counted as part of the Gaussian signal, instead being absorbed into the polynomial background. The contributions to the relative systematic error on $`f_1`$ are estimated to be $`5\%`$ from the uncertainty on the photon-finding efficiency and $`11\%`$ from the uncertainty on the probability of feed down, obtained by performing alternate fits to the $`\mathrm{\Delta }M(D_s\pi ^0)`$ distribution.

Inverting Equations 1 and 2, we find that $`R_0=155\pm 23`$ decays and $`R_1=41\pm 12`$ decays, where the uncertainties include both statistical and systematic sources. The result for $`R_1`$ demonstrates the existence of a state at 2463 MeV/c2. The significance of the signal for this state, accounting for statistical and systematic errors, is determined to be in excess of $`5\sigma `$ by computing the probability for the combinatoric background plus the feed up background to fluctuate up to give the observed yield in the signal region in Fig. 4(a).

### V.2 Further evidence for the $`𝑫_{𝒔𝑱}\mathbf{(}\mathrm{𝟐𝟒𝟔𝟑}\mathbf{)}^\mathbf{+}\mathbf{}𝑫_𝒔^{\mathbf{}}𝝅^\mathrm{𝟎}`$ decay

We conclude from the analysis described in the preceding section that a new state, the $`D_{sJ}(2463)`$, exists in addition to the $`D_{sJ}^{}(2317)`$ state reported by BaBar, because feed up from the $`D_{sJ}^{}(2317)`$ is only a minor background component ($`25\%`$) of the narrow peak observed in Fig. 4(a). To provide further support for this conclusion, we have directly measured the feed up background in Fig. 4(a) due to $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ plus random photon combinations, by selecting combinations in $`D_s^{}`$ side band regions in the $`D_s\gamma \pi ^0`$ sample. The $`M(D_s\gamma \pi ^0)M(D_s\gamma )`$ distribution for this sample, plotted in Fig. 4(b), shows only a small enhancement in the region of the $`D_{sJ}(2463)`$, demonstrating that the background from $`D_{sJ}^{}(2317)`$ decays indeed constitutes only a small fraction of the entries in the $`D_{sJ}(2463)`$ peak.

We performed a binned likelihood fit of the spectrum in Fig. 4(a) to a Gaussian signal shape plus a second-order polynomial plus the spectrum from the $`D_s^{}`$ side band region in Fig. 4(b) with its normalization fixed. From this fit, we obtain $`R_1=45.7\pm 11.6`$ decays, consistent with the value of $`R_1`$ obtained from Equations 1 and 2. From the change in the likelihood of fits performed with and without the $`D_{sJ}(2463)`$ signal contribution, we infer that the statistical significance of the signal is $`5.7\sigma `$.

Finally we note that the width of the peak in Fig. 4(a), $`\sigma =6.1\pm 1.0`$MeV/c2, is consistent with the detector resolution. If the origin of this peak was feed up from $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ decays, then the effect of including unrelated photons to form $`D_s^{}\pi ^0`$ candidates would be to smear out the $`\mathrm{\Delta }M(D_s^{}\pi ^0)`$ distribution, in the same way that the feed down background to the $`D_{sJ}^{}(2317)`$ state is broadened as described in the preceding section. From fits to Monte Carlo simulations of this feed up process, the expectation for the width is determined to be $`\sigma =14.9\pm 0.6`$ MeV/c2. Thus, the narrowness of the peak in Fig. 4(a) also rules out the possibility that the peak is dominantly due to feed up from $`D_{sJ}^{}(2317)^+`$ decays.

## VI Properties of the $`𝑫_{𝒔𝑱}^{\mathbf{}}\mathbf{(}\mathrm{𝟐𝟑𝟏𝟕}\mathbf{)}^\mathbf{+}`$ and $`𝑫_{𝒔𝑱}\mathbf{(}\mathrm{𝟐𝟒𝟔𝟑}\mathbf{)}^\mathbf{+}`$ States

### VI.1 Mass and width of the $`𝑫_{𝒔𝑱}^{\mathbf{}}\mathbf{(}\mathrm{𝟐𝟑𝟏𝟕}\mathbf{)}^\mathbf{+}`$

Having obtained evidence for the $`D_{sJ}(2463)`$ state, and having characterized the background that it contributes in the $`\mathrm{\Delta }M(D_s\pi ^0)`$ mass difference spectrum, we are now able to further address properties of the $`D_{sJ}^{}(2317)`$ state. We recall that our measurement of the width of the peak in Fig. 1 is $`\sigma =8.0_{1.1}^{+1.3}`$MeV/c2, somewhat larger than our mass difference resolution, $`\sigma =6.0\pm 0.3`$MeV/c2. This difference is consistent with predictions from Monte Carlo simulations where we include both $`D_{sJ}(2463)`$ and $`D_{sJ}^{}(2317)`$ production, since roughly $`18\%`$ of the observed $`D_s^+\pi ^0`$ decays in the $`D_{sJ}^{}(2317)`$ signal region enter as feed down from the $`D_{sJ}(2463)`$ state, this ‘background’ peak having an expected width of $`\sigma =14.9\pm 0.4`$MeV/c2.

To better determine the mass and natural width of the $`D_{sJ}^{}(2317)`$, we carry out a binned likelihood fit of the peak in the $`\mathrm{\Delta }M(D_s\pi ^0)`$ spectrum in Fig. 1(b) to a sum of two Gaussians, one for the $`D_{sJ}^{}(2317)`$ signal, and one to account for the feed down from the $`D_{sJ}(2463)`$. Allowing the means and widths of both Gaussians to float, we obtain $`\mathrm{\Delta }M(D_s\pi ^0)=350.0\pm 1.2`$MeV/c2 with $`\sigma =6.0\pm 1.2`$MeV/c2 for the $`D_{sJ}^{}(2317)`$ component. The mean mass difference and width for the feed down component are $`344.9\pm 6.1`$MeV/c2 and $`16.5\pm 6.3`$MeV/c2, respectively. The errors in the above values are due to statistics only; systematic errors are discussed below. Both widths are consistent with predictions from Monte Carlo simulations in which the two states are modeled with a natural width of zero.

We have also carried out fits in which one or both of the widths of the Gaussians were fixed to values determined by the Monte Carlo. In all cases the results were consistent with the results from the fit described above. We have also tried to obtain a purer $`D_{sJ}^{}(2317)`$ sample by vetoing combinations with photons that can be combined with the $`D_s`$ candidate to form a $`D_s^{}`$, thereby removing some of the feed down background from the $`D_{sJ}(2463)`$. This veto marginally improves the $`D_s\pi ^0`$ signal when we fit with two Gaussians, and the mass and width change by only a small fraction of the statistical uncertainty. The systematic uncertainty for $`\mathrm{\Delta }M(D_s\pi ^0)`$ receives contributions from uncertainties in the characterization of the $`D_{sJ}(2463)`$ feed down and from uncertainties in the modeling of the energy resolution of the calorimeter. We estimate the total systematic error on the mass difference to be 1.0 MeV/c2. Based on these studies, we limit the natural width of the $`D_{sJ}^{}(2317)`$ to be $`\mathrm{\Gamma }<7`$ MeV at the $`90\%`$ confidence level (C.L.).

### VI.2 Mass and width of the $`𝑫_{𝒔𝑱}\mathbf{(}\mathrm{𝟐𝟒𝟔𝟑}\mathbf{)}^\mathbf{+}`$

From the fit to the distribution resulting from the subtraction of Fig. 4(b) from Fig. 4(a) reported in Section V.2, we obtain $`\mathrm{\Delta }M(D_s^{}\pi ^0)=351.2\pm 1.7\pm 1.0`$MeV/c2 for the mass difference between the $`D_{sJ}(2463)`$ and the $`D_s^{}`$. The first error is statistical and the second is the systematic uncertainty which is the same as that presented in the previous section for the $`D_{sJ}^{}(2317)D_s`$ mass difference. From our fits to data and Monte Carlo $`\mathrm{\Delta }M(D_s^{}\pi ^0)`$ distributions, we also infer a $`90\%`$ C.L. upper limit on the natural width ($`\mathrm{\Gamma }`$) of the $`D_{sJ}(2463)^+`$ state to be $`7`$MeV.

### VI.3 Production Properties

We now give a measure of the production rates of $`D_{sJ}^{}(2317)`$ and $`D_{sJ}(2463)`$ mesons. A full understanding would require the determination of the fragmentation functions of both particles and their branching ratios into the final states we observe. To minimize systematic errors, we report the relative yields with respect to $`D_s^+`$ production, where all putative charmed-antistrange systems have momenta greater than 3.5 GeV/c. We use all observed events for each channel, which includes direct production and any contributions from decays of higher mass objects. Then

$`{\displaystyle \frac{\sigma \left(D_{sJ}^{}(2317)D_s^+\pi ^o\right)}{\sigma \left(D_s^+\right)}}`$ $`=`$ $`(7.9\pm 1.2\pm 0.4)\times 10^2,`$ (3)
$`{\displaystyle \frac{\sigma \left(D_{sJ}(2463)D_s^+\pi ^o\right)}{\sigma \left(D_s^+\right)}}`$ $`=`$ $`(3.5\pm 0.9\pm 0.2)\times 10^2,`$ (4)
We also note that

$$\frac{\sigma \left(D_s^+(2112)D_s^+\gamma \right)}{\sigma \left(D_s^+\right)}=0.59\pm 0.03\pm 0.01,$$
(5)
here and above, the first error includes the statistical and systematic errors on the event yields while the second includes the systematic errors for photon detection ($`2\%`$), and for $`\pi ^0`$ detection ($`5\%`$).

### VI.4 Decays of $`𝑫_{𝒔𝑱}^{\mathbf{}}\mathbf{(}\mathrm{𝟐𝟑𝟏𝟕}\mathbf{)}`$ to other final states

With regard to the alternate $`D_{sJ}^{}(2317)`$ decay channels described earlier, in which no signals were observed, we summarize the limits on the branching fractions relative to the $`D_s^+\pi ^0`$ mode in Table 1. The normalization for these limits is based on the determination that $`(81.7\pm 5.7)\%`$ of the observed yield of $`165\pm 20`$ entries in the peak of the $`\mathrm{\Delta }M(D_s\pi ^0)`$ spectrum in Fig. 1(b) are attributable to $`D_{sJ}^{}(2317)D_s\pi ^0`$ decay after accounting for the feed down from decays of the $`D_{sJ}(2463)`$ state to $`D_s^{}\pi ^0`$. We have estimated the systematic error on this yield to be $`\pm 16`$ entries by varying selection criteria and the parameterization of signal and background shapes used in the fit to Fig. 1.

The event yields for the various final states are obtained by fitting the mass difference distributions to Gaussians with each mean fixed to the result from the $`D_s^+\pi ^0`$ channel and each width given by the resolution determined from the simulation of the corresponding decay mode. Uncertainties are dominated by the statistical error on the fitted yields and limits on the relative rates are calculated assuming a Gaussian distribution with negative values not allowed.

### VI.5 Decays of $`𝑫_{𝒔𝑱}\mathbf{(}\mathrm{𝟐𝟒𝟔𝟑}\mathbf{)}`$ to other final states

Unlike the case of a $`0^+`$ state, the $`D_s\pi ^+\pi ^{}`$ decay mode, as well as both radiative decay modes $`D_s\gamma `$ and $`D_s^{}\gamma `$ are allowed for a state with $`J^P=1^+`$. From fits to the mass difference distributions displayed in Figs. 2 and 3 for peaks in the regions where a contribution from the $`D_{sJ}(2463)`$ would appear, we find no evidence of decays to any of these final states. We summarize the limits obtained on these decays, relative to $`D_s^{}\pi ^0`$, in Table 2.

Despite a high relative efficiency, the limit on the decay $`D_{sJ}(2463)^+D_s^+\gamma `$ is less stringent than those on the decays to $`D_s^^+\gamma `$ and $`D_s^+\pi ^+\pi ^{}`$. This is due to an excess of combinations in the signal region. From fits performed with and without the signal Gaussian, we determine that the statistical significance of this excess is 2.4 standard deviations.

If the $`D_{sJ}(2463)^+`$ is a $`1^+`$ state, then it is also possible for it to undergo a P-wave radiative decay to $`D_{sJ}^{}(2317)^+\gamma `$ noteonxfeed . We have looked for this transition in our $`D_s\gamma \pi ^0`$ sample. To reduce backgrounds from $`D_{sJ}(2463)^+D_s^+\gamma `$, we required that the $`D_s\pi ^0`$ system be consistent with the decay of the $`D_{sJ}^{}(2317)`$, namely that $`|\mathrm{\Delta }M(D_s\pi ^0)350.0\text{MeV/c}^2|<13.4`$MeV/c2 ($`2\sigma `$ based on Monte Carlo simulations). We also required that the $`D_s\gamma `$ system be inconsistent with $`D_s^{}`$ decay at the $`1\sigma `$ level (the corresponding $`\mathrm{\Delta }M(D_s\gamma )`$ must deviate from the expected value for this decay by more than $`4.4`$MeV/c2), and that the momentum of the $`\pi ^0`$ be inconsistent with the $`D_{sJ}(2463)D_s^{}\pi ^0`$ transistion, also at the $`1\sigma `$ level. The $`M(D_s\pi ^0\gamma )M(D_s\pi ^0)`$ distribution, plotted in Fig. 5, provides no evidence for a signal in the vicinity of 150 MeV/c2.

Due to the tightness of these cuts, the efficiency for detecting this decay is roughly a factor of three smaller than for the $`D_s^{}\pi ^0`$ decay mode. The $`90\%`$ C.L. upper limit for this channel is reported in the bottom row of Table 2.

## VII Conclusions and Discussion

In summary, data from the CLEO II detector have provided confirming evidence for the existence of a new narrow resonance decaying to $`D_s^+\pi ^0`$, with a mass near 2.32 GeV/c2. This state is consistent with being the $`0^+`$ member of the lowest-lying $`P`$-wave $`c\overline{s}`$ multiplet. As summarized in Table 1, we have set upper limits on other decay modes of this state. We have measured the mass splitting of this state with respect to the $`D_s`$ meson to be $`350.0\pm 1.2\text{[stat.]}\pm 1.0\text{[syst.]}`$MeV/c2, and we find its natural width to be $`\mathrm{\Gamma }<7`$MeV at $`90\%`$ C.L.

We have observed and established the existence of a new narrow state with a mass near 2.46 GeV/c2 in its decay to $`D_s^+\pi ^0`$, which we have denoted $`D_{sJ}(2463)`$. We have demonstrated that the signal for this decay cannot be interpreted as a reflection from the $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ decay. The measured properties of this state are consistent with its interpretation as the $`1^+`$ partner of the $`0^+`$ state in the spin multiplet with light quark angular momentum of $`j=1/2`$. We have measured the mass splitting of this state with respect to the $`D_s^{}`$ meson to be $`351.2\pm 1.7\text{[stat.]}\pm 1.0\text{[syst.]}`$MeV/c2. The natural width of this state is found to be $`\mathrm{\Gamma }<7`$MeV at $`90\%`$ C.L. Since the $`D_{sJ}(2463)`$ mass lies above the kinematic threshold for decay to $`DK`$ (but not for $`D^{}K`$), the narrow width suggests this decay does not occur. Since angular momentum and parity conservation laws forbid a $`1^+`$ state from decaying to two pseudoscalars, this provides additional evidence for the compatibility of the $`D_{sJ}(2463)`$ with the $`J^P=1^+`$ hypothesis.

In the model of Bardeen, Eichten and Hill bardeeneichtenhill , a $`J^P=1^+`$ state is predicted with the same mass splitting $`\mathrm{\Delta }M`$ with respect to the $`1^{}`$ state as that between the $`0^+`$ and $`0^{}`$ states. Taking the difference between the two mean mass differences reported above, we obtain $`\delta (\mathrm{\Delta }M)=(351.2\pm 1.7)(350.0\pm 1.2)=1.2\pm 2.1`$MeV/c2 for the difference between the $`1^+1^{}`$ and $`0^+0^{}`$ mass splittings, where the dominant uncertainty is due to statistics. Thus our observations are consistent with these predictions.

We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. We thank W. Bardeen, E. Eichten, S. Godfrey, C. Hill and J. Rosner for useful discussions. M. Selen thanks the Research Corporation, and A.H. Mahmood thanks the Texas Advanced Research Program. This work was supported by the National Science Foundation, and the U.S. Department of Energy.


format: text

text:
# Spontaneous Magnetization and Electron Momentum Density in 3D Quantum Dots

## Abstract

We discuss an exactly solvable model Hamiltonian for describing the interacting electron gas in a quantum dot. Results for a spherical square well confining potential are presented. The ground state is found to exhibit striking oscillations in spin polarization with dot radius at a fixed electron density. These oscillations are shown to induce characteristic signatures in the momentum density of the electron gas, providing a novel route for direct experimental observation of the dot magnetization via spectroscopies sensitive to the electron momentum density.

Quantum dots, which may be thought of as artificial “atoms”, not only provide a nanoscale laboratory for fundamental quantum mechanical studies, but also hold the promise of novel spintronic and other applications ashoori96 ; wolf01 ; ciorga02 ; smet02 . Spin polarization in nanosystems is therefore drawing great attention in view of its intrinsic importance tarucha96 ; zabala98 ; fujisawa02 , and its possible relationship to other phenomena, such as kinks in the positions of Coulomb blockade peaks in 2D quantum dots baranger00 , or the so-called 0.7 conduction anomaly in quantum point contacts wang98 . In 2D dots, experiments point to shell-like filling of levels in accord with Hund’s first rule tarucha96 , and a recent spin density-functional calculation indicates the presence of a ground state with high spin in relatively large dots jiang02 . In 3D dots, shell structure of energy levels has been reported bakkers01 , and results on the fine structure in the energy spectra of gold nanoparticles have been interpreted in terms of a ground state spin magnetization larger than 1/2 davidovic00 . However, Ref. jacquod00,  argues that in finite disordered systems fluctuations in off-diagonal interaction matrix elements may suppress ground state magnetization. It is clear that spin polarization of quantum dots is a subject of considerable current interest and controversy.

In this work, we develop an exactly solvable many-body model Hamiltonian of wide applicability for describing correlated electrons confined to a quantum dot. The ground state is clearly found to exhibit spontaneous magnetization and strong oscillations in spin polarization as a function of the dot radius. In order to identify observable signatures of the remarkable aforementioned oscillations in polarization, we consider the momentum density of the electron gas in the dot at some length. The key quantity in this connection turns out to be the width $`\mathrm{\Delta }p`$ over which the electron momentum density falls to zero around the Fermi cut-off momentum. In particular, oscillations in magnetization are shown to induce similar oscillations in $`\mathrm{\Delta }p`$, indicating that this effect would be amenable to direct experimental observation via electron momentum density spectroscopies such as angular correlation of positron annihilation radiation (ACAR) and Compton scattering. Our results also bear more broadly on the applications of quantum dots in that they point to the possibility of obtaining a system in a ferromagnetic ground state without requiring low temperatures and/or external magnetic fields.

We consider the Hamiltonian

$$\widehat{H}=\underset{\nu \sigma }{}ϵ_\nu ^0a_{\nu \sigma }^{}a_{\nu \sigma }+\frac{1}{2}U\underset{\nu \nu ^{}\sigma }{}a_{\nu \sigma }^{}a_{\nu ^{}\sigma }^{}a_{\nu ^{}\sigma }^{}a_{\nu \sigma }^{}.$$
(1)
The first term ($`\widehat{H}_0`$) describes the noninteracting system. It corresponds to the spherical square well: $`V(r)=V_0`$, for $`rR`$, and $`V(r)=0`$ otherwise, and the associated one-particle eigenfunctions $`\phi _\nu `$ with eigenvalues $`ϵ_\nu ^0`$; $`\nu nlm`$ is a composite index and $`\sigma `$ labels spin. The second term ($`\widehat{H}_I`$) describes the interaction and allows Coulomb repulsion via the parameter $`U`$ only between electrons of opposite spin—electrons of like spin being kept apart by the Pauli exclusion principle. In this connection, we have carried out extensive numerical estimates, and find that direct Coulomb matrix elements are larger, on the average, by an order of magnitude or more, than exchange or other terms for dot radii up to 30 Å and electron densities ranging from $`r_s=2`$ to 5 (Bohr radii) comment0 . Spin-orbit effects here are also expected to be small, since the confining potential is relatively weak due to the screening of the nuclei by the core electrons. For these reasons, form (1) should provide a reasonable model for the electronic structure of 3D quantum dots.

The standard one-particle Green function for the many electron system can be expanded as $`𝒢_{\sigma \sigma ^{}}(\mathrm{𝐫𝐫}^{},\tau )=_{\nu \nu ^{}}\phi _\nu ^{}(𝐫)\phi _\nu ^{}(𝐫^{})𝒢_{\sigma \sigma ^{}}(\nu \nu ^{},\tau )`$, with $`𝒢_{\sigma \sigma ^{}}(\nu \nu ^{},\tau )=T_\tau a_{\nu \sigma }(\tau )a_{\nu ^{}\sigma ^{}}^{}(0)`$. It is readily shown that $`\mathrm{}a_\nu /\tau =\xi _\nu a_{\nu \sigma }U\widehat{N}_\sigma a_{\nu \sigma }`$, where $`\xi _\nu =ϵ_\nu ^0\mu `$ ($`\mu `$ is the chemical potential), and $`\widehat{N}_\sigma =_\nu ^{}a_{\nu ^{}\sigma }^{}a_{\nu ^{}\sigma }^{}`$. Then, the equation of motion for the Green function $`𝒢_{\sigma \sigma ^{}}(\nu \nu ^{},\tau )`$ is

$`\mathrm{}{\displaystyle \frac{}{\tau }}𝒢_{\sigma \sigma ^{}}(\nu \nu ^{},\tau )`$ $`=`$ $`\mathrm{}\delta _{\sigma \sigma ^{}}\delta _{\nu \nu ^{}}\delta (\tau )+\xi _\nu 𝒢_{\sigma \sigma ^{}}(\nu \nu ^{},\tau )`$ (2)
$``$ $`UT_\tau \widehat{N}_\sigma (\tau )a_{\nu \sigma }^{}(\tau )a_{\nu ^{}\sigma ^{}}^{}(0).`$
Since our Hamiltonian possesses the property that $`[\widehat{H}_I,\widehat{H}_0\mu \widehat{N}]=0`$, Wick’s theorem allows us to factorize the preceding ground state average, i.e., $`T_\tau \widehat{N}_\sigma (\tau )a_{\nu \sigma }^{}(\tau )a_{\nu ^{}\sigma ^{}}^{}(0)=N_\sigma 𝒢_{\sigma \sigma ^{}}(\nu \nu ^{},\tau )`$, where $`N_\sigma =_\nu ^{}f_{\nu ^{}\sigma }`$, with $`f_{\nu ^{}\sigma }`$ denoting the Fermi occupation function. A Fourier transform in $`\tau `$ immediately yields the exact Green function

$$𝒢_{\sigma \sigma ^{}}(\nu \nu ^{},i\omega _n)=\delta _{\sigma \sigma ^{}}\delta _{\nu \nu ^{}}\frac{1}{i\omega _n\xi _\nu /\mathrm{}UN_\sigma /\mathrm{}}.$$
(3)
The self-energy thus is $`\mathrm{\Sigma }_\sigma ^{}=UN_\sigma /\mathrm{}`$. The interacting energy levels, $`ϵ_{\nu \sigma }`$, are given by the solutions of $`\omega \xi _\nu /\mathrm{}\mathrm{\Sigma }_\sigma ^{}=0`$ (frequencies measured with respect to $`\mu /\mathrm{}`$) minimizing the ground state energy, $`E=_{\nu \sigma }(ϵ_\nu ^0+\frac{1}{2}\mathrm{}\mathrm{\Sigma }_\sigma ^{})f_{\nu \sigma }`$. Given the number of particles $`N`$, this defines a set of nonlinear equations for the populations, $`N_{}`$ and $`N_{}`$, of up and down spin states, respectively. The resulting splitting in energy for states of opposite spin, $`\mathrm{\Delta }=U(N_{}N_{})`$, is uniform, i.e., it does not depend on the quantum numbers $`\nu `$. Note that temperature enters our formulae only through the Fermi function, $`f_{\nu \sigma }`$; all results in the remainder of this article refer to the zero temperature limit.

We comment briefly on the nature of the four essential parameters which describe our model: $`r_s`$, the electron density; $`R`$, the dot radius or, equivalently, the number $`N`$ of electrons \[$`N=(R/r_s)^3`$\]; $`V_0`$, the well-depth; and, $`U`$, the Coulomb energy. The values of $`r_s`$ and $`R`$ are related to the type and size of the systems considered. In this study we use $`r_s=5`$, a relatively low density, enabling us to consider a wide range of dot radii. $`V_0`$ may then be obtained reasonably in terms of the work function of the system studied and the Fermi energy of the free electron gas with the same $`r_s`$. For $`r_s=5`$, following Ref. saniz02, , we take $`V_0=8.62`$ eV. The choice of the remaining parameter $`U`$ is somewhat tricky, mainly because correlations in quantum dots are not well understood. A handle on $`U`$ can be obtained by considering the average value $`U_{\mathrm{TF}}`$ of the Thomas-Fermi (TF) screened Coulomb interaction, $`𝑑𝐫𝑑𝐫^{}|\phi _\nu (𝐫)|^2v_{\mathrm{TF}}(|𝐫𝐫^{}|)|\phi _\mu (𝐫^{})|^2`$, averaged over all occupied states in the noninteracting system, with $`v_{\mathrm{TF}}(r)=e^2\mathrm{exp}(r/l_{\mathrm{TF}})/r`$, where $`l_{\mathrm{TF}}=\sqrt{r_s}/1.56`$ is the TF screening length. However, actual screening in real materials is likely to be weaker. For example, in noble metals, Ref. delfatti00, , reports an effective screening length $`l_{\mathrm{TF}}/0.73`$. Using this value, we find an enhanced Coulomb interaction given roughly by $`1.75U_{\mathrm{TF}}`$, which is the value we have used throughout this work.

Fig. 1 illustrates the nature of the energy levels $`ϵ_{\nu \sigma }`$ for a 12.11 Å radius dot, for which $`N=96`$, and $`U=27.13`$ meV ($`U_{\mathrm{TF}}=15.50`$ meV). The self-consistent splitting between up and down spin electrons is found to be, $`\mathrm{\Delta }=1.03`$ eV. The highest occupied molecular orbital (HOMO) energy level is, $`E_{\mathrm{HOMO}}=5.86`$ eV, which may be thought of as the dot “Fermi energy”. We have $`N_{}=67`$ and $`N_{}=29`$, so that the ground state spin polarization per electron is $`\zeta =(N_{}N_{})/N=0.4`$.

We discuss now $`\zeta `$ as a function of dot radius with reference to Fig. 2. In the weakly interacting limit ($`U0`$), $`\mathrm{\Delta }0`$, still the degeneracy between up and down spin in the $`nl`$ shells is lifted, and as Fig. 2(a) clearly shows, shells fill in accord with Hund’s first rule. Indeed, as the number of particles increases with dot radius, the spin polarization reaches a peak each time a shell is half filled with up-spin electrons, and falls to zero when the shell is completed with down-spin electrons. The sequence of these so-called “magic numbers”, i.e., $`N`$ values for which all occupied shells are completely filled, will of course differ here from that in atoms or 2D quantum dots due to differences in the details of the underlying spectrum. The magic numbers up to 106 are indicated on the upper axis in Fig. 2(a). Turning next to the strongly interacting case, Fig. 2(b) shows large deviations in $`\zeta `$ from a simple Hund’s rule filling comment1 . This is because $`\mathrm{\Delta }`$ changes with each added electron in order to minimize the total energy. The first 4 electrons enter the up spin $`1s`$ and $`1p`$ levels, so that $`\zeta =1`$, while the next 4 enter the corresponding down spin levels until $`\zeta `$ vanishes for $`N=8`$. The 9th electron induces a strong change in $`\mathrm{\Delta }`$, so that the $`1s`$ and $`1p`$ down spin levels are pushed above the $`1d`$ up spin level, causing the former to empty in favor of the latter and the system to be become completely polarized again. Similar level reorderings are involved in the other “kinks” seen in Fig. 2(b).

In order to identify observable signatures of the remarkable changes in the polarization depicted in Fig. 2, we consider the electron momentum density (EMD), defined by $`n(𝐩)=(2\pi )^1_{\mathrm{}}^{\mathrm{}}𝑑\omega f(\omega )A(𝐩,\omega )`$, where $`f`$ is the Fermi function and $`A`$ is the spectral function $`A(𝐩,\omega )=2\mathrm{Im}G^R(𝐩,\omega )`$, with $`G^R(𝐩,\omega )=𝒢(𝐩,i\omega _n\omega +i\delta )`$. The inset in Fig. 3 illustrates $`n(p)`$ and the magnitude of its first derivative, $`|n^{}(p)|`$, for the same parameter values as in Fig 1. The region of rapid variation in $`n(p)`$ can be characterized via the position, $`p_\mathrm{F}`$, of the peak in $`|n^{}(p)|`$ and the associated full-width-at-half-maximum (FWHM), $`\mathrm{\Delta }p`$. In the bulk limit in a metallic system $`p_\mathrm{F}`$ will tend to the Fermi momentum, where $`n(p)`$ will develop a break and $`|n^{}(p)|`$ a $`\delta `$-function peak, correspondingly. For simplicity, therefore, we may refer to $`p_\mathrm{F}`$ loosely as the dot “Fermi momentum”, even though there are no breaks in the EMD in a finite system comment2 ; comment3 . Fig. 3 considers $`\mathrm{\Delta }p`$ systematically and shows the presence of a dramatic series of peaks (solid line) spaced more or less regularly as a function of $`R`$. A comparison with corresponding changes in $`\zeta `$ (dotted curve), makes it evident that peaks in $`\mathrm{\Delta }p`$ are well correlated with those in $`\zeta `$, some differences notwithstanding. This correlation arises because, at a peak in $`\zeta `$, the up and down spin EMDs become substantially shifted with respect to each other, reflecting the polarization of the system, so that the EMD presents a much larger value of $`\mathrm{\Delta }p`$ overall. In short, peaks in $`\mathrm{\Delta }p`$ clearly are a signature of peaks in $`\zeta `$.

Motivated by the preceding considerations, we have obtained an approximate expression for $`\mathrm{\Delta }p`$ as

$$\mathrm{\Delta }p/p_\mathrm{F}=\mathrm{\Delta }p_0/p_\mathrm{F}+c\mathrm{\Delta }(\zeta \zeta _0)R,$$
(4)
where $`\zeta _0`$ is the spin polarization in the weakly interacting case ($`U0`$), and $`c`$ is a constant, which depends on various dot parameters, with a fitted value in the present case of, $`c=4.29\times 10^2`$ eV-1Å-1. The first term in Eq. (4) describes the $`U0`$ limit and can be shown to follow the scaling law $`0.93\mathrm{\AA }(r_s/a_0)/R`$ saniz02 . The second term incorporates the effect of the spin polarization induced through the interaction $`U`$. Fig. 4 shows that Eq. (4) provides an excellent description of the exact $`\mathrm{\Delta }p`$ data as a function of $`R`$. The approximation continues to be equally good up to $`R=30`$ Å, although for clarity results over a shorter $`R`$-range are shown in the figure.

We emphasize that the form of our solution in Eq. (3) does not depend explicitly on that of the noninteracting Hamiltonian ($`\widehat{H}_0`$). Although details of $`\zeta `$ will vary with those of the spectrum of $`\widehat{H}_0`$, we expect our prediction of oscillations in $`\zeta `$ with dot radius to be generally robust to changes in the shape and dimensionality of the confining potential comment5 ; comment6 . Thus, we expect our results to be relevant to almost one dimensional systems. Indeed, in metallic nanowires comment6b , Zabala et al. zabala98  report a series of peaks in the magnetic moment per electron as a function of wire radius, correlated strongly with another observable of the system (the elongation force of the nanowire). The physics driving these results is similar to that in our work comment7 , despite the fact that Ref. zabala98,  considers the stabilized jellium model within the framework of the spin-dependent density-functional formalism, while we treat correlation effects differently, in terms of the direct interaction $`U`$ and neglect various off-diagonal exchange and Coulomb contributions comment8 .

The relationship between $`\zeta `$ and $`\mathrm{\Delta }p`$ provides a powerful handle for a direct experimental verification of the polarization oscillations in quantum dots predicted in this study via the use of solid state spectroscopies sensitive to the EMD, namely, angular-correlation of positron-annihilation radiation (ACAR) \[positrons have the advantage of being local probes, see Weber et al. in Ref. comment2, \], inelastic X-ray scattering (IXS) in the deeply inelastic (Compton) regime and scanning tunneling microscopy (STM). We have carried out further computations at other $`r_s`$ values and find that, generally speaking, polarization effects become less prominent with decreasing $`r_s`$, as the kinetic energy dominates, as well as with increasing $`R`$, as electrons are less confined and $`U`$ becomes weaker. Nevertheless, oscillations in $`\mathrm{\Delta }p`$ are observable even for $`r_s3`$, particularly for dot sizes $`R10`$ Å.

In summary, we have presented an exactly solvable model Hamiltonian for discussing the properties of the interacting electron gas in the confined geometry of a quantum dot. Although the results presented in this article are based on a confining potential in the shape of a spherical square well, computations for other geometries and confining potentials would be quite straightforward. The ground state is found to exhibit spontaneous magnetization and striking oscillations in spin polarization $`\zeta `$ as a function of the dot radius $`R`$, at the fixed electron density considered. The oscillations in $`\zeta `$ are shown to induce similar oscillations in the width $`\mathrm{\Delta }p`$ around the Fermi cut-off in the momentum density, which we refer to for simplicity as the dot “Fermi momentum”, providing a novel route for direct experimental observation of the dot magnetization via spectroscopies sensitive to momentum density (specially, positron annihilation and Compton scattering). A simple expression for $`\mathrm{\Delta }p`$ is discussed, which gives an excellent approximation to the exact numerical data as a function of radius. We expect the results of the present study to be robust to the details of the confining potential and to the approximations inherent in the treatment of correlations in our model Hamiltonian.

This work was supported by the U. S. Department of Energy contracts DE-AC03-76SF00098 and W-7405-ENG-48 and the U. S. Department of Energy Office of Basic Energy Science, Division of Materials Science, and benefitted from the allocation of supercomputer time at the NERSC and the Northeastern University Advanced Scientific Computation Center (ASCC). R.S. acknowledges the hospitality and support of Northeastern University’s Physics Department during a research visit.


text:
# Modeling the dynamical evolution of the M87 globular cluster system

## 1 Introduction

A number of recent observational studies have provided detailed information on the properties of the M87 GCS. Most recently, Kundu, Zepf & Ashman (2002) have used multiple HST pointings to show that the turnover of the luminosity function of the GCS (hereafter we will indicate the luminosity function and the mass function of the globular cluster system by GCLF and GCMF respectively) is constant over a large range of galactocentric distances. The data analyzed include clusters up to a projected galactocentric distance equal to about 70 kpc and confirm and extend the lack of radial gradient in the GCLF turnover already found by Kundu et al. (1999) for clusters in the inner regions of M87. Because of the very large number of clusters in the M87 GCS and the high quality of the HST data, the uncertainty in the determination of the GCLF turnover is small and a stringent constraint on its possible radial variation is thus imposed by the new data.

The lack of a radial gradient in the GCS mean luminosity is sometimes interpreted as an indication that evolutionary processes, whose efficiency is known to depend on the galactocentric distance, did not play an important role in determining the current properties of GCSs. However, the studies of the dynamical processes and their effects on the properties of globular clusters are based on fundamental stellar dynamics and there is no way to keep dynamical evolution from occurring. The results of all the theoretical investigations of the dynamical evolution of globular clusters and GCSs (see e.g. Fall & Rees 1977, Gnedin & Ostriker 1997, Murali & Weinberg 1997a,b, Baumgardt 1998, Vesperini 1997, 1998, 2000, 2001 and references therein) show that evolutionary processes lead to the disruption of a significant number of clusters (the exact fraction of the initial number of clusters depending on the properties of the host galaxy; see e.g. Murali & Weinberg 1997a, Vesperini 2000, 2001) and to the modification of the properties of the surviving clusters.

Several recent papers on models of the dynamical evolution of GCSs have addressed the issue of the radial behavior of the GCMF. Vesperini (1998, 2000), in studies of the evolution of the Galactic GCS and of GCSs of elliptical galaxies, has shown that for a log-normal initial GCMF (or more in general for GCMFs which are bell-shaped in $`\mathrm{log}M`$; see Vesperini 2002) with mean mass and dispersion similar to those currently observed, evolutionary processes are very efficient in disrupting a large fraction of clusters but they do not always produce a significant radial variation of the GCMF parameters. On the other hand, several theoretical investigations (see e.g. Baumgardt 1998, Murali & Weinberg 1997a,b, Vesperini 1998,2001) have shown that for a power-law initial GCMF as steep as that observed in young GCS in merging galaxies, evolutionary processes lead to the development of mean mass radial gradients not consistent with observations.

In a recent study of the evolution of the Galactic GCS, Fall & Zhang (2001) have shown that, adopting an initial orbital distribution characterized by a strong radial anisotropy increasing with the galactocentric distance, the final GCMF has a modest radial variation that is in fair agreement with the observed one, even for a steep power-law for the initial GCMF. Nevertheless, the radial anisotropy necessary to obtain the desired results is very large and, although the initial anisotropy can be reduced by disruption of clusters on more eccentric orbits, observational estimates of the current anisotropy of the Galactic GCS (see e.g Dinescu, Girard & van Altena 1999) indicate a much smaller anisotropy.

In this paper we will model the evolution of the M87 GCS and explore the range of initial conditions leading to final properties consistent with those observed. The GCLF properties determined by Kundu et al. (2002) and the kinematical properties determined by Cote et al. (2001), Romanowsky & Kochanek (2001), Cohen (2000) and Cohen & Ryzov (1997) will be used to constrain viable initial conditions. The outline of this paper is the following. In section 2 we describe the method used for our study and the initial conditions considered; in section 3 we show the results obtained, in section 4 we discuss our results and in section 5 we summarize our conclusions.

## 2 Method and initial conditions

The method we use to follow the evolution of a model GCS is the same as that adopted in Vesperini (2000,2001) to study the evolution of GCS in elliptical galaxies. The evolution of the masses of individual clusters is calculated using the results of N-body simulations by Vesperini & Heggie (1997). The effects of mass loss due to stellar evolution (assuming a stellar IMF with a slope $`\alpha =2.3`$; similar results are obtained if the 2-slope power-law stellar IMF suggested by Kroupa (2001) is adopted), two-body relaxation, dynamical friction and of the tidal field of the host galaxy have been included. The mass loss rate due to two-body relaxation has been calculated assuming the tidal radius of each cluster is determined by the pericentric distance. We refer the reader to Vesperini (2000) for further details.

We adopt a mass model for M87 proposed by McLaughlin (1999a) on the basis of various observational constraints available on the stellar and dark matter distribution in M87 and the Virgo Cluster; in this model, adopting a distance of 15 Mpc to M87, the total mass enclosed within a radius $`r`$ is equal to

$`M_{tot}(r)`$ $`=`$ $`8.1\times 10^{11}M_{}\left[(r/5.1\text{kpc})/(1+r/5.1\text{kpc})\right]^{1.67}+`$ (1)
$`7.06\times 10^{14}M_{}[\mathrm{ln}(1+r/560\text{kpc})(r/560\text{kpc})/(1+r/560\text{kpc})].`$
For the initial spatial distribution of the GCS we have considered the following profiles:

* a Navarro, Frenk & White (1996) (hereafter NFW) profile with scale radius, $`r_s`$, equal to 9.1 kpc. This model has been shown (McLaughlin 1999b) to fit the observed surface number density profile of the M87 GCS.
* An NFW profile with $`r_s=9.1`$ kpc for $`r<95`$ kpc and an NFW profile with $`r_s=560`$ kpc for $`r>95`$ kpc. The choice of the NFW profile with $`r_s=560`$ kpc for $`r>95`$ kpc, where there are no observational constraints on the GCS spatial distribution, has been proposed by McLaughlin (1999b) and it relies on the assumption that the density profile of the GCS becomes the same as that of gas and stars at $`r95`$ kpc.
* It has been often claimed that the inner flattening in the spatial distribution of the M87 GCS can not be the result of the effects of dynamical disruption of clusters (Lauer & Kormendy 1986). In order to further explore if this is indeed the case we have also considered two profiles equal to those considered in A) and B) above but with $`r_s=0.1`$ kpc instead of the current best fit value of 9.1 kpc, giving much more concentrated initial GCS profile.

As for the initial GCMF, we have restricted our choice to the functional forms and the parameters indicated by observational studies of young and old cluster systems: our fiducial choice was a power-law function, $`f(M)M^\alpha `$ for $`10^4M_{}<M<10^7M_{}`$ , with index $`\alpha =1.8`$ as observed in young clusters systems in merging galaxies (see e.g. Zepf et al. 1999 and Whitmore 1999 for a review). Since, with an initial power-law GCMF, none of our models satisfied all the observational constraints, we have explored the evolution of GCSs with a two-slope power-law with $`\alpha =1.8`$ for $`10^{5.25}M_{}<M<10^7M_{}`$ and $`\alpha =0.2`$ for $`10^4M_{}<M<10^{5.25}M_{}`$. The latter GCMF is bell-shaped with a turnover at $`\mathrm{log}M=5.25`$ if binned in $`\mathrm{log}M`$ (see e.g. McLaughlin 1994) and is similar to that currently observed in old cluster systems (see e.g. Ashman & Zepf 1998, Harris 2001). As discussed in Vesperini (2002), the evolution of this GCMF is not significantly different from that of other bell-shaped GCMFs, such as a log-normal or a t-Student distribution, with the same turnover.

The initial density profiles and GCMFs considered are summarized in Table 1.

As for the GCS velocity distribution, we adopted an initial distribution characterized by a radial anisotropy increasing with the galactocentric distance as in the Osipkov-Merrit models (Osipkov 1979, Merritt 1985); for these models the ratio of the radial to the tangential velocity dispersion increases with galactocentric distance, $`r`$, as

$$\frac{\sigma _R^2}{\sigma _T^2}=1+r^2/r_a^2.$$
(2)
Beyond the anisotropy radius, $`r_a`$, the distribution is dominated by radial orbits whereas for $`r<r_a`$, it is essentially isotropic. According to the calculations of Fall & Zhang (2001), for the Galactic GCS, this distribution along with a particular choice of the anisotropy radius ($`r_a=5`$ kpc) appears to be the only one leading to a GCS mean mass radial gradient in fair agreement with observations if a power-law initial GCMF as steep as that of young cluster systems in merging galaxies is adopted.

In order to fully explore the dependence of the GCS evolution on the initial anisotropy in the velocity distribution, for each of the initial density profiles mentioned above, simulations with initial values of $`r_a`$ ranging from 1 kpc to 250 kpc have been carried out.

Each GCS investigated was initially made of 400000 clusters with positions randomly drawn from the spatial profiles described above; given the current position of a cluster, $`r`$, the local radial velocity dispersion was calculated using the Jeans equation and the local tangential velocity dispersion was then calculated using the value of $`\sigma _R^2/\sigma _T^2`$ at $`r`$; the radial and the tangential velocities of each cluster were then drawn from Maxwellian velocity distributions with radial and tangential dispersions calculated above. This method, introduced by Hernquist (1993), has been used in many studies to set initial conditions for N-body simulations of both spherical and disk galaxies.

The evolution of each GCS has been followed until $`t=15`$ Gyr and the GCS properties at several intermediate ages have also been studied.

## 3 Results

### 3.1 Power-law initial GCMF

In Fig. 1, we plot the mean mass of the globular clusters against the projected galactocentric distance for theoretical models with different initial anisotropy radii. We also plot the observational data from Kundu et al. (2002) analysis. The absence of any significant radial gradient in the mean mass observed for the M87 GCS requires small values of the initial anisotropy radius ($`r_a\text{ }<\text{ }3`$ kpc). Initial conditions with anisotropy radii larger than 3 kpc produce final mean masses smaller than those observed and steep radial gradients inconsistent with data.

We have also considered possible radial variations in $`M/L_V`$ for the M87 GCS. The only observed property of the globular clusters that varies with galactocentric radius which is relevant for $`M/L_V`$ is the ratio of blue to red globular clusters, which increases with galactocentric distance for M87 (Kundu et al. 2002) as it does for most other galaxies (see Ashman & Zepf 1998). This color difference tends to produce a difference in $`M/L_V`$ for typical models of stellar populations, but this difference is much smaller than the mass variation predicted by models of dynamical evolution with a power-law initial GCMF and anisotropy radii larger than a few kpc. Specifically, using current stellar populations models (Bruzual & Charlot 2000, Maraston et al. 2002) we find that $`(M/L_V)_{red}1.6(M/L_V)_{blue}`$ for a range of possible age and metallicity combinations that are consistent with the observed colors of the red and blue populations of globular clusters. Applying this estimate to the M87 system, for which we approximate the color gradient by adopting a ratio of the number of red to blue clusters equal to 1 for $`R_g<20`$ kpc and equal to 0 beyond 20 kpc, we find the inner clusters would have a mass-to-light ratio 1.3 times larger than the mass-to-light ratio of outer clusters leading to a $`\mathrm{\Delta }\overline{\mathrm{log}M}_f0.1`$. This is significantly smaller than the radial gradient obtained in models with large initial anisotropy radii. Moreover, the mean mass of outer clusters from simulations with large anisotropy radii is much smaller than the mass obtained from the observed mean luminosity with any reasonable value of the mass-to-light ratio.

This comparison with the observational constraints on the radial variation of the mean mass of clusters has significantly narrowed the range of possible initial values of $`r_a`$ for the M87 GCS with a power-law initial GCMF.

We now turn our attention to the observational constraints imposed by the kinematics of the M87 GCS. In fig.2 we show the observed projected velocity dispersion profile of the M87 GCS as determined by Cote et al. (2001) along with the dispersion profiles of various models with different orbital anisotropy. Density profile A (fig.2a), density profile B (fig.2b) and an Osipkov-Merritt radial profile for the anisotropy (see eq.2 above) are used to calculate the theoretical profiles. Values of $`r_a`$ larger than 50 kpc are necessary to fit the projected velocity dispersion profile if density profile A is adopted and larger than 200 kpc for density profile B. These values are much larger than the initial anisotropy radius required to produce a radial gradient of the mean mass of clusters consistent with observations. Fig. 3 shows the final projected velocity dispersion profiles from our simulations for initial values of $`r_a`$ leading to final GCMF properties consistent with observations: it is clear that these profiles do not fit the observational data. For the small initial values of $`r_a`$ necessary to produce GCMF properties consistent with observations, the number density profile does not evolve significantly (see discussion below) and, therefore, simulations starting with profiles C and D lead to both kinematical properties and final density profiles inconsistent with observations.

Fig.4 shows the radial profile of $`\sigma _R^2/\sigma _T^2`$ at the end of the simulations with initial $`r_a=2`$ and $`r_a=3`$ kpc and with initial density profile A, and with initial $`r_a=1`$ and $`r_a=2`$ kpc with initial density profile B: the preferential disruption of clusters on high eccentricity orbits is not sufficient to reduce the strong initial anisotropy to a level consistent with observations; the $`\sigma _R^2/\sigma _T^2`$ profile from the most and the least anisotropic models used to fit the observed velocity dispersion profile in fig. 2 along with the value of $`\sigma _R^2/\sigma _T^2`$ of the different models proposed by Romanowsky & Kochanek (2001) calculated at the outermost distance constrained by observations are also shown in fig.4.

It is thus clear from these results that, adopting a strong initial radial anisotropy increasing with galactocentric distance, although necessary to produce final GCMF properties consistent with observations from a steep initial power-law GCMF, is not a viable solution as the final anisotropy is stronger than shown by observations of the radial velocities of the globular clusters. It is also important to point out that, although numerical simulations of galaxy formation predict a radial anisotropy increasing with galactocentric distance, the initial anisotropy required to obtain GCMF properties consistent with data is much stronger than that found in simulations. For example, numerical simulations of the formation of giant galaxies at the center of groups and clusters carried out by Dubinski (1998) predict a $`\sigma _R^2/\sigma _T^2`$ growing from 0 to 2 from the center of the galaxy to 100 kpc; the initial anisotropy required to fit the GCMF properties in our models is such that $`\sigma _R^2/\sigma _T^210^3`$ at 100 kpc.

Fig. 5 shows the time evolution of $`\overline{\mathrm{log}M}_f(R_g)`$ for two models with $`r_a`$ equal to 2 and 150 kpc. In the model with initial anisotropy radius equal to 2 kpc the mean mass profile is always approximately flat and, as a result of the preferential disruption of low-mass clusters, the GCS mean mass increases with time. In the model with $`r_a=150`$ kpc, a strong radial gradient quickly develops and it becomes steeper and more extended as the system evolves; the quick formation of the mean mass radial gradient implies that a significant radial gradient would have to be expected in models with large anisotropy radii even if the time needed for the potential of M87 to reach its equilibrium state were considered and the GCS evolved in the current M87 potential for a time shorter than that (15 Gyr) considered in our simulations.

As shown above, none of the initial conditions considered leads to final properties consistent with both the observed GCLF properties and the observed kinematical properties. Although it is therefore not necessary to explore the evolution of the GCS number density profile to further constrain the range of viable initial conditions, we conclude this subsection with a few remarks on the evolution of the GCS spatial distribution as this can shed further light on the GCS evolution.

Fig.6a and 6b show the final GCS surface density for a few different initial values of $`r_a`$ along with the initial profile (profile A for fig.6a and profile C for fig.6b). The initial velocity distribution adopted is such that the anisotropy increases with galactocentric distance and, in particular, radial orbits dominate beyond $`r_a`$. For small values of $`r_a`$, most clusters, regardless of their current galactocentric distances, have small pericentric distances (see fig.8a) and, since we have assumed that the pericentric distance determines the timescale of evolution of clusters, the disruption rate is approximately constant and independent of the current galactocentric distances of clusters (see fig.7). As a consequence of that, the shape of the final spatial profile is very similar to the initial one; this is also the reason why small values of $`r_a`$ do not produce a gradient in the mean mass: for small values of $`r_a`$, clusters at any distance from the galactic center are disrupted and lose mass at approximately the same rate. More in general, no radial gradient in any other cluster property will be produced as a result of evolution for GCS with initial small values of $`r_a`$. Note that if the disruption rate does not depend on galactocentric distance, the decrease of specific frequency in the inner regions of M87 reported by McLaughlin (1999b; see Rhode & Zepf 2001 for a similar radial dependence of the specific frequency in NGC 4472) can not be ascribed to dynamical evolution and it should be explained in terms of a radial dependency of cluster formation efficiency.

As $`r_a`$ increases and the velocity distribution becomes more isotropic, the range of pericentric distances broadens (see fig.8b) and clusters closer to the center of the host galaxy are preferentially disrupted (see fig.7); the disruption rate is therefore higher near the center of the galaxy and a central flattening of the density profile ensues. Therefore we find that evolutionary processes, in particular disruption of inner low-mass clusters, can be responsible for the observed central flattening in the spatial distribution of the M87 GCS (see also Murali & Weinberg 1997a, Capuzzo Dolcetta & Tesseri 1997) and, contrary to previous claims, this does not have to due to initial conditions.

### 3.2 Two-slope power-law initial GCMF

As discussed in the previous subsection, starting with a power-law initial GCMF similar to that of young GCSs observed in merging galaxies, we could not find any model satisfying all the observational constraints. We have therefore continued our investigation by exploring the evolution of GCSs with a two-slope power-law initial GCMF. The parameters of the initial GCMF adopted are summarized in Table 1. We note that this GCMF is bell-shaped if binned in $`\mathrm{log}M`$ and, as discussed in Vesperini (2002), the results obtained using this functional form for the initial GCMF do not differ significantly from those obtained adopting other bell-shaped functions often used in the literature to fit the GCLF of old GCS, such as a log-normal, or a t-Student distribution (see e.g. Secker 1992). As shown in the previous subsection, an initial number density profile similar to that observed preserves its shape only for an initial velocity distribution characterized by a strong radial anisotropy inconsistent with the kinematical data available. For this reason, here, we focus our attention only on number density profiles C and D.

Fig.9 shows the $`\overline{\mathrm{log}M}_f`$ radial profile for our simulations of the M87 GCS with different initial values of $`r_a`$: the GCS final $`\overline{\mathrm{log}M}_f`$ and its radial profile are consistent with the observed mean mass and lack of a radial gradient in the mean mass for the entire range of values of $`r_a`$ considered and no constraint on the anisotropy of the initial velocity distribution is imposed by the final GCMF properties.

In fig.10 we compare the final projected velocity dispersion radial profile from our simulations of M87 GCS with that determined by Cote et al’s observational data. Fig. 10a shows the results for initial number density profile C and $`r_a=30,60,200`$ kpc; initial values of the anisotropy radius around 60 kpc ($`50<r_a<100`$ kpc) appear to be those which better fit the observed velocity dispersion profile, but all the profiles resulting from simulations with initial anisotropy radii ranging from 30 kpc to 200 kpc fall within the $`90\%`$ confidence limits and increase with radius as in the observed profile. As expected (see fig.2 above), if initial number density profile D is adopted (fig.10b), larger initial values of $`r_a`$ ($`r_a>200`$ kpc) are required to obtain final velocity dispersion profiles fitting the data.

In fig. 11 we compare the theoretical $`\sigma _R^2/\sigma _T^2`$ profile obtained in the simulations leading to final GCMF properties and velocity dispersion profile consistent with observations and the $`\sigma _R^2/\sigma _T^2`$ profile derived from our fits in fig.2; fig.11 shows that the results from our simulations are consistent with the constraints on the anisotropy imposed by observational data.

Finally, in fig.12 we compare the observed surface density profile with the final surface density profile from our simulation with initial density profile C and initial $`r_a`$ equal to 60 kpc. Final surface density profiles from simulations with $`r_a>30`$ kpc are essentially indistinguishable from that shown in fig.12. This figure shows that a density profile initially more peaked than the observed one can evolve, as a result of the preferential disruption of inner clusters, toward a final shape in good agreement with the observed one.

Fig.13 shows the radial dependence of the fraction of surviving clusters; assuming a similar initial density profile for stars and for the GCS, this plot clearly shows that disruption can produce the observed central flattening of the GCS density profile and explain a rising local specific frequency up to about 20-25 kpc.

A large fraction of the total initial GCS mass is lost because of cluster disruption and mass loss from surviving clusters. Although the majority of the mass lost is from clusters in the inner regions of the galaxy, a full understanding of the relevance of this mass to the build-up of the central black hole would require a careful study of the orbital evolution of the mass lost from clusters which is beyond the scope of this paper.

The dependence of the mass loss on the initial mass of clusters is shown in Fig. 14. This figure shows the ratio of the final to the initial mass of surviving clusters, $`M_f/M_i`$, versus the initial mass of clusters from one of our simulations: as expected, most of the surviving low-mass clusters have lost more than half of their initial mass but it is interesting to point out that this is also the case for a large number of clusters with high initial mass. Although, as expected, low-mass clusters are preferentially disrupted by dynamical processes, a significant fraction of clusters with large initial masses are disrupted as well; for example, for a model with initial density profile C and $`r_a=60`$ kpc, about 40 per cent of the initial population of clusters with $`\mathrm{log}M_i5`$ is disrupted by evaporation in 15 Gyr.

The results of the simulations discussed in this subsection show that, starting from a two-slope power-law initial GCMF with parameters similar to those currently observed in old GCSs, it is possible to produce models satisfying all the observational constraints (the two-slope power-law GCMF adopted is bell-shaped if binned in log M and the results would be similar if other bell-shaped GCMFs with the same initial turnover were adopted). Additional observational data, particularly on the kinematics of the M87 GCS, would be valuable to further constrain the range of viable initial conditions.

## 4 Discussion

The simulations carried out in this investigation to explore the range of viable initial conditions for the M87 GCS and the role of evolutionary processes have shown that models with a power-law initial GCMF require a very strong initial radial anisotropy to avoid the development of a steep mean mass radial gradient inconsistent with the roughly constant $`\overline{\mathrm{log}M}_f`$ radial profile observed. On the other hand, adopting the strong radial anisotropy required to obtain a flat $`\overline{\mathrm{log}M}_f`$ radial profile leads to final kinematical properties inconsistent with the kinematical data available, which are best-fit by models with a much weaker radial anisotropy. We have explored the role of several factors which could have affected this conclusion but we found it to be very robust. In particular we have shown that a steep $`\overline{\mathrm{log}M}_f`$ radial gradient in models with kinematical properties matching the observed kinematical data develops in a few gigayears and even if the current M87 potential has not been in place for one Hubble time, a $`\overline{\mathrm{log}M}_f`$ radial gradient inconsistent with observational data is still expected. We have also considered the effect of a radial variation of the mass-to-light ratio, but the $`\overline{\mathrm{log}M}_f`$ radial gradient obtained from the observed radial profile of the mean luminosity is still much smaller than that produced in theoretical models. In our simulations we started with an initial GCMF independent on the galactocentric distance; although it would be possible to assume a GCMF with an ad hoc initial radial variation to counteract the gradient produced by evolutionary processes, there is neither theoretical nor observational support for such a choice.

As already shown in a number of previous theoretical investigations, a power-law initial GCMF can be turned by the evolutionary processes considered in our study into a bell-shaped (in $`\mathrm{log}M`$) GCMF resembling the GCMF of old GCSs. However, if the evolution is driven only by disruption and evolution due to two-body relaxation, interaction with the tidal field of the host galaxy and dynamical friction, none of the models with a power-law initial GCMF we have studied could match both the observed kinematical properties and the observed $`\overline{\mathrm{log}M}_f`$ radial profile.

Our study shows that in models starting with a bell-shaped initial GCMF with parameters similar to those observed, evolutionary processes lead to final properties satisfying all the observational constraints ($`\overline{\mathrm{log}M}_f`$ radial profile, spatial distribution and kinematical data). If the initial GCMF is indeed a power-law, as both observational studies of young cluster systems and theoretical studies of cluster formation suggest, the details of the evolution of an initial power-law GCMF into a bell-shaped GCMF are still to be clarified.

We are currently exploring this issue by focussing our attention on the effects of clusters disruption (see e.g. Chernoff & Weinberg 1990, Fukushige & Heggie 1995) caused by the early mass loss due to the evolution of massive stars (Vesperini & Zepf 2003). In particular we are studying to what extent early disruption of clusters can alter the initial GCMF and how the evolution of the initial GCMF depends on the distribution of cluster concentrations and on the correlation of cluster concentrations with other cluster properties.

## 5 Conclusions

In this paper we have carried out a number of simulations of the M87 GCS evolution. We have considered a number of different initial conditions and determined which can lead, driven by the effects of internal relaxation, interaction with the tidal field of the host galaxy and dynamical friction, to final properties consistent with the observed GCS properties (radial profile of the GCS mean mass, GCS spatial distribution and GCS kinematical properties).

We have adopted two different initial GCMF and focussed our attention only on the functional forms and the parameters indicated by observational studies of young and old cluster systems: a power-law similar to that observed in young cluster systems in merging galaxies and a two-slope power-law with parameters similar to that currently observed in old cluster systems. The initial velocity distribution considered is such that the ratio of the radial to the tangential velocity dispersion increases with the galactocentric distance as in the Osipkov-Merritt models ($`\sigma _R^2/\sigma _T^2=1+r^2/r_a^2`$); a wide range of values for the anisotropy radius have been considered ($`1<r_a<250`$ kpc). NFW number density profiles with different values of the scale radius have been studied (see Table 1 for a summary of the initial conditions considered).

For a power-law initial GCMF we could not find any initial condition leading to final properties satisfying all the observational constraints. For the two-slope power-law initial GCMF considered here (but the results are similar for other initial GCMF bell-shaped in $`\mathrm{log}M`$ such as a log-normal or a t-Student distribution) we found a range of initial conditions whose evolution lead to final properties consistent with all the observational constraints.

In particular we find that:

1. for a power-law initial GCMF, only models with a strong initial radial anisotropy fit the flat $`\overline{\mathrm{log}M}_f`$ profile observed. Specifically, only models with an initial value of the anisotropy radius $`r_a\text{ }<\text{ }3`$ kpc have a final mean mass radial profile fairly matching the observed $`\overline{\mathrm{log}M}_f(r)`$;
2. the small initial values of $`r_a`$ necessary to obtain a flat mean mass radial profile from a power-law initial GCMF lead to final kinematical properties characterized by a strong radial anisotropy. This is inconsistent with observational kinematical data which, instead, favor models characterized by more isotropic velocity distributions.
3. For a two-slope power-law initial GCMF with parameters similar to those observed in old GCSs, the final $`\overline{\mathrm{log}M}_f`$ radial profile is approximately flat and consistent with observations for all the values of $`r_a`$ considered ($`1<r_a<250`$ kpc) and for all the initial number density profiles considered.
4. For the two-slope power-law initial GCMF adopted in this study, a number of the models studied also have final kinematical properties and final spatial distributions consistent with observations and they, therefore, satisfy all the observational constraints available. In particular for the initial number density profile C (single NFW profile with initial scale radius equal to 0.1 kpc; see Table 1) values of the anisotropy radius larger 30 kpc lead to final kinematical properties consistent with observations with the best-fitting models being those with $`50\text{ }<\text{ }r_a\text{ }<\text{ }100`$ kpc. For this range of values of $`r_a`$, the initially peaked spatial distribution significantly flattens as a results of the preferential disruption of inner clusters and the final surface number density is perfectly consistent with observations. For the initial number density profile D (double NFW profile with scale radii equal to 0.1 and 560 kpc; see Table 1) more isotropic initial velocity distributions are necessary to produce final kinematical properties consistent with observations ($`r_a\text{ }>\text{ }200`$ kpc).
5. Although low-mass clusters are those preferentially affected by evolutionary processes, we have shown that for a significant fraction of clusters with initial masses larger than $`10^5M_{}`$ the current mass is less than 50 per cent of the initial mass.
6. We have explored the dependence of the fraction of disrupted clusters on the galactocentric distance and shown that, although clusters are disrupted more efficiently at small galactocentric distances, disruption can be significant over a large radial range and explain a rising local specific frequency out to large radii.

In a future work (Vesperini & Zepf 2003) we will study the role of early disruption of clusters with different initial concentrations caused by the mass loss due to the evolution massive stars and we will explore the effect of this process on the evolution of the GCMF.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge support from NASA via the ATP grant NAG5-11320, the LTSA grant NAG5-11319 and grant AR-09208.01 from Space Telescope Science Institute, which is operated by AURA under NASA contract NAS 5-26555. We thank D.E. McLaughlin for kindly providing us with the data of the M87 GCS velocity dispersion profile and an anonymous referee for useful comments.

## REFERENCES

Ashman K.M., Zepf S.E., 1998, Globular Cluster Systems, Cambridge University Press
Baumgardt, H., 1998, A&A, 330, 480
Bruzual, A.G., Charlot, S. 2000, private communications
Capuzzo Dolcetta, R., Tesseri, A., 1997, MNRAS, 292, 808
Chernoff D.F., Weinberg M.D., 1990, ApJ, 351, 121
Cohen, J., 2000, AJ, 119, 162
Cohen, J., Ryzhov, A., 1997, ApJ, 486, 230
Cote, P., et al., 2001, ApJ, 559, 828
Dinescu, D.I., Girard, T.M., van Altena, W.F., 1999, AJ, 117, 1792
Dubinski, J., 1998, ApJ, 502, 141
Fall, S.M., Rees, M.J., 1977, MNRAS, 181, 37P
Fall, S.M., Zhang, Q., 2001, ApJ, 561, 751
Fukushige T., Heggie, D.C., 1995, MNRAS, 276, 206
Gnedin O.Y., Ostriker J.P., 1997, ApJ, 474, 223
Harris, W.E., 2001, in Star Clusters, Lectures for 1998 Saas-Fee Advanced School, L. Labhardt, B. Binggeli eds., (Springer-Verlag, Berlin, 2001)
Hernquist, L., 1993, ApJS, 86, 389
Kroupa, P., 2001, MNRAS, 322, 231
Kundu A., Whitmore B.C., Sparks W.B., Macchetto F.D., Zepf S.E., Ashman K.M., 1999, ApJ 513, 733
Kundu, A., Zepf, S.E., Ashman, K.M., 2002, in preparation
Lauer, T. R., Kormendy, J., 1986, ApJ, 303, L1
Maraston, C. et al., 2002, A&A submitted
McLaughlin D.E. 1994, PASP, 106, 47
McLaughlin D.E. 1999a, ApJ, 512, L9
McLaughlin D.E. 1999b, AJ, 117, 2398
Merritt, D., 1985, AJ, 90, 1027
Murali C., Weinberg M.D., 1997a, MNRAS, 288, 767
Murali C., Weinberg M.D., 1997b, MNRAS, 291, 717
Navarro, J., Frenk, C., White, S., 1996, ApJ, 462, 563
Osipkov, L.P., 1979, PAZh, 5, 77
Rhode,K.L., Zepf, S.E., 2001, AJ, 121, 210
Romanowsky, A.J., Kochanek, C.S., 2001, ApJ, 553, 722
Secker,J., 1992, AJ, 104, 1472
Vesperini, E., 1997, MNRAS, 287, 915
Vesperini, E., 1998, MNRAS, 299, 1019
Vesperini, E., 2000, MNRAS, 318, 841
Vesperini, E., 2001, MNRAS, 322, 247
Vesperini, E., 2002, in Grebel E.K., Geisler D., Minniti D. eds., Extragalactic Star Clusters, p.664
Vesperini, E., Heggie, D.C., 1997, MNRAS, 289, 898
Vesperini, E., Zepf, S.E., 2003, ApJ, 587, L97
Whitmore, B.C., 1999, in Barnes J.E., Sanders D.B., eds., Proc. IAU Symp. 186, Galaxy interactions at low and high redshift. Kluwer, Dordrecht, p.251
Zepf, S.E., Ashman, K.M., English, J., Freeman, K.C., Sharples, R.M., 1999, AJ, 118, 752