text:
# The first determination of Generalized Polarizabilities of the proton by a Virtual Compton Scattering experiment

## Abstract

Absolute differential cross sections for the reaction $`epep\gamma `$ have been measured at a four-momentum transfer with virtuality $`Q^2=0.33`$ GeV2 and polarization $`ϵ=0.62`$ in the range 33.6 to 111.5 MeV/c for the momentum of the outgoing photon in the photon-proton center of mass frame. The experiment has been performed with the high resolution spectrometers at the Mainz Microtron MAMI. From the photon angular distributions, two structure functions which are a linear combination of the generalized polarizabilities have been determined for the first time.

PACS numbers: 13.60.Fz, 14.20.Dh, 13.10.+q, 25.30.Rw

Polarizabilities are fundamental quantities that characterize the response of a composite system to static or slowly varying external electric or magnetic fields. For the proton, Compton scattering of real photons has been used to determine the electric and magnetic polarizabilities . In this letter, we report measurements of Compton scattering of virtual photons, leading to the first determination of generalized polarizabilities (GPs). Virtual Compton scattering (VCS) off the proton refers to the reaction $`\gamma ^{}pp\gamma `$ where $`\gamma ^{}`$ stands for an incoming virtual photon of four-momentum squared $`Q^2`$. This reaction is experimentally accessed through photon electroproduction $`epep\gamma `$. We study this reaction in the regime where the produced photon has a small enough energy that its electric (E) and magnetic (M) fields look constant over the size of the nucleon. In this regime the reaction can be interpreted as electron scattering on a nucleon placed in a quasi-constant applied EM field . The induced motion of the nucleon as a whole can be eliminated thanks to a low energy theorem , so one is left with the deformation, due to the applied field, of the nucleon internal currents $`\delta J^\mu (r)`$ and the electron scattering measures its Fourier transform $`\delta J^\mu (Q)`$. To lowest order in $`\alpha _{QED}`$, $`\delta J^\mu (Q)`$ is linear in the applied field and the 6 coefficents of proportionality are the GPs . When $`Q^2`$=0 two of them reduce to the usual polarizabilities $`\alpha _E`$ and $`\beta _M`$ measured in real Compton scattering . A measurement of how the nucleon’s internal current distribution is deformed by an external EM field will yield valuable information about its non perturbative structure. This is illustrated by the quite different results obtained for the GPs according to the used model . As an example of new information provided by VCS and which could not be obtained either from ordinary electron scattering or real Compton scattering, consider the case of Coulomb scattering on a nucleon placed in a magnetic field. The action of the latter is mainly to flip the spin of the quarks, which hardly modifies the charge density. Therefore the corresponding GPs (namely $`P^{(11,00)1}`$ and $`P^{(11,02)1}`$ in Eq. 8) receive no contribution from the quarks. This is confirmed by actual calculations. As a consequence the value of these GPs should be mainly determined by the pion cloud.

First indications of the feasibility of a VCS experiment were reported in . However only recently, with the advent of high-intensity continuous beam electron accelerators, a thorough experimental investigation of the VCS process became possible. This letter reports on the first VCS experiment dedicated to the GPs which has been performed at the Mainz Microtron MAMI.

The general theoretical framework for VCS has been extensively described by Guichon et al.  and by Drechsel et al.  and only the relevant issues will be discussed here. In the reaction $`epep\gamma `$, the final photon can be emitted either by the electron, which is called the Bethe-Heitler (BH) process, or by the proton. The latter process is the VCS process which consists of Born and Non-Born amplitudes; while the Born amplitude depends only on static properties of the proton (charge, mass, form factors), the Non-Born amplitude contains dynamical internal structure information in terms of GPs. Despite of the strong dominance of the BH process, the VCS process can be measured under appropriately chosen kinematics through the interference with the dominant BH amplitude. In the zero-energy limit of the final photon, the cross section is independent of the dynamical nucleon structure , and can be evaluated using only the known BH and Born amplitudes. This is summarized by the following expression:

$`d^5\sigma ^{exp}(\mathrm{q},\mathrm{q}^{},ϵ,\theta ,\phi )`$ $`=d^5\sigma ^{BH+Born}(\mathrm{q},\mathrm{q}^{},ϵ,\theta ,\phi )`$ (1)
$`+\varphi \mathrm{q}^{}\mathrm{\Psi }_0(\mathrm{q},ϵ,\theta ,\phi )+𝒪(\mathrm{q}^2)`$ (2)
where $`d^5\sigma `$ is a notation for the differential cross section $`d^5\sigma /d\mathrm{k}_{lab}^{}[d\mathrm{\Omega }_e]_{lab}[d\mathrm{\Omega }_p]_{CM}`$ wherein $`\mathrm{k}_{lab}^{}`$ is the absolute value of the outgoing electron momentum in the laboratory, $`\mathrm{q}`$ and $`\mathrm{q}^{}`$ are the absolute values of the three-momenta of the virtual and real photons (in the photon-proton center of mass (CM) system) respectively and $`ϵ`$ is the virtual photon polarization. $`\theta `$ is the angle between the real and virtual photon in the CM system and $`\phi `$ is the angle between the electron and the photon-proton plane. $`\varphi `$ stands for a phase space factor. $`\mathrm{\Psi }_0(\mathrm{q},ϵ,\theta ,\phi )`$ is the leading term in the expansion in powers of the real photon momentum $`\mathrm{q}^{}`$. It contains the dynamical internal structure information of the proton, expressed by 6 GPs. The latter are denoted by $`P^{(\rho ^{}L^{},\rho L)S}(\mathrm{q})`$ where $`L(L^{})`$ are the initial (final) photon orbital angular momentum, $`\rho (\rho ^{})`$ the type of multipole transition (0 for Coulomb, 1 for magnetic) , and $`S`$ distinguishes between non-spin-flip ($`S=0`$) and spin-flip ($`S=1`$) transitions at the nucleon side. One has 2 non-spin-flip GPs, $`P^{(01,01)0}`$, $`P^{(11,11)0}`$, proportional to $`\alpha _E`$ and $`\beta _M`$ at $`Q^2=0`$ respectively, and 4 spin-flip GPs, $`P^{(11,11)1}`$, $`P^{(11,00)1}`$, $`P^{(11,02)1}`$, $`P^{(01,12)1}`$. The GPs are functions of $`\mathrm{q}`$ (or $`\stackrel{~}{Q}^2=Q^2|_{\mathrm{q}^{}=0}`$). In an unpolarized measurement, $`\mathrm{\Psi }_0(\mathrm{q},ϵ,\theta ,\phi )`$ can be written as:

$`\mathrm{\Psi }_0(\mathrm{q},ϵ,\theta ,\phi )`$ $`=v_1(\theta ,\phi ,ϵ)(P_{LL}(\mathrm{q})P_{TT}(\mathrm{q})/ϵ)`$ (3)
$`+v_2(\theta ,\phi ,ϵ)P_{LT}(\mathrm{q})`$ (4)
where $`v_1(\theta ,\phi ,ϵ)`$, $`v_2(\theta ,\phi ,ϵ)`$ are known kinematical factors and $`P_{LL}(\mathrm{q})`$, $`P_{TT}(\mathrm{q})`$ and $`P_{LT}(\mathrm{q})`$ are structure functions defined by:

$`P_{LL}=`$ $`2\sqrt{6}mG_EP^{(01,01)0}`$ (5)
$`P_{TT}=`$ $`3G_M\text{q}^2\left(\sqrt{2}P^{(01,12)1}P^{(11,11)1}/\stackrel{~}{q}_0\right)`$ (6)
$`P_{LT}=`$ $`\sqrt{{\displaystyle \frac{3}{2}}}{\displaystyle \frac{m\text{q}}{\stackrel{~}{Q}}}G_EP^{(11,11)0}`$ (7)
$`+{\displaystyle \frac{\sqrt{3}}{2}}{\displaystyle \frac{\stackrel{~}{Q}}{\text{q}}}G_M(P^{(11,00)1}+{\displaystyle \frac{\mathrm{q}^2}{\sqrt{2}}}P^{(11,02)1})`$ (8)
where $`m`$ stands for the proton mass, $`G_E`$ and $`G_M`$ denote the form factors evaluated at $`\stackrel{~}{Q}^2`$ and $`\stackrel{~}{q}_0`$ is the CM virtual photon energy at $`\mathrm{q}^{}=0`$.

Absolute cross sections $`d^5\sigma ^{exp}`$  have been measured using the three-spectrometer facility  of the A1 collaboration at the 855 MeV Mainz Microtron MAMI. The scattered electron and the recoiling proton were detected in coincidence with two of the high-resolution magnetic spectrometers. The photon production process was selected by a cut on the missing mass around zero, which was possible thanks to the excellent resolution of the facility (momentum resolution of $`10^4`$ and angular resolution better than 3 mrad). Liquid hydrogen served as a target, contained in a target cell of 49.5 mm total length with a havar wall thickness of 9 $`\mu `$m. The use of typical electron currents of 30 $`\mu `$A yielded a luminosity of $`=410^{37}`$ cm-2s-1.

In this first VCS experiment below pion threshold, the five-fold differential cross section was measured in a wide angular range at 5 values of the photon momentum $`\mathrm{q}^{}`$: 33.6, 45.0, 67.5, 90.0, and 111.5 MeV/c. Two kinematical variables were kept fixed, namely the virtual photon momentum, q = 600 MeV/c ($`\stackrel{~}{Q}^2`$ = 0.33 GeV2) and the polarization $`ϵ`$ = 0.62. The out-of-plane angle $`\phi `$ range is determined by the acceptance of the two spectrometers around 0 and 180. To ease the presentation (see Fig. 1), the data are plotted with $`\theta `$ ranging from -180 to +180, the negative values correspond to $`\phi =180^{}`$. The wide range of $`\theta `$ from -141 to +6 was covered by only changing the setting of the proton spectrometer. This angular range covers the backward direction relative to the incoming and outgoing electron. Here the VCS contributions are dominant because the BH radiation is mainly emitted in the electrons directions.

The influence of the GPs gives rise to a small deviation of the measured cross section from the BH+Born cross section. Theoretically it is expected to be about 10% at the highest $`\mathrm{q}^{}`$ value. The cross sections at each bin are determined within a statistical accuracy of 3%; however, a careful analysis of possible systematic errors on the above deviation is of particular importance. (1) The proton form factors are not exactly known. Consequently, we also measured the absolute elastic scattering cross section for each kinematical setting of the VCS experiment. These measurements validate the use of the form factor parametrization from Höhler  at a precision better than $`\pm `$ 1%. (2) The luminosity and the detector efficiencies are controlled within the same accuracy. (3) The radiative corrections, which are of the order of 20% of the cross section, have been calculated by Vanderhaeghen et al.  taking into account all the diagrams up to order $`\alpha ^4`$ in the VCS cross section. The corresponding systematic uncertainties are estimated to equal $`\pm 2\%`$. (4) An extensive Monte Carlo simulation  has been performed in order to verify the missing mass spectra and to determine the solid angles within an accuracy of $`\pm `$2%. The code generates events according to the BH+Born cross section and takes into account all resolution deteriorating effects and real photon radiations. (5) While the uncertainties (1)-(4) are constant over the angular range of the real photon, small imperfections in the spectrometer optics calibration produce distortions of the angular distributions: these effects lead to an estimated change in absolute cross section of $`\pm `$2.5%.

Fig. 1 shows the measured differential cross sections (with statistical errors only) as a function of $`\theta `$ for 5 values of the real photon momentum. The cross sections $`d^5\sigma ^{BH+Born}`$, calculated using the Höhler form factors, are presented by the solid lines. At the smallest photon momentum, $`\mathrm{q}^{}`$ = 33.6 MeV/c, the agreement between the data and $`d^5\sigma ^{BH+Born}`$ is excellent, while with increasing $`\mathrm{q}^{}`$ one observes growing deviations from this known cross section. This represents the effect of the polarizabilities.

Fig. 2 shows $`(d^5\sigma d^5\sigma ^{BH+Born})/\varphi \mathrm{q}^{}`$ as a function of the real photon momentum $`\mathrm{q}^{}`$ at the 14 measured angles $`\theta `$, of which the intercept with the ordinate axis directly yields $`\mathrm{\Psi }_0`$ in Eq. (1). Several methods were applied in order to determine the intercept, which will be discussed in a forthcoming paper. The methods that were tried, taking into account a possible $`\mathrm{q}^{}`$ dependence of $`(d^5\sigma d^5\sigma ^{BH+Born})/\varphi \mathrm{q}^{}`$, all show that this dependence is weak. As such, here the most simple approach was used, i.e. no $`\mathrm{q}^{}`$ dependence. $`\mathrm{\Psi }_0`$ is then determined at each angle $`\theta `$ by the weighted mean value of the data at the 5 photon energies. Fig. 3 presents the resulting values of $`\mathrm{\Psi }_0/v_2`$ as a function of $`v_1/v_2`$ (cf. Eq. (2)). The data are well aligned, this allows us to extract the two structure functions $`P_{LL}P_{TT}/ϵ`$ and $`P_{LT}`$ as the slope and intercept, respectively, of a linear fit to the data.

Table 1 shows the results with a statistical error and two systematic errors, the first one corresponds to the normalisation of the angular distributions, the second one stems from the distortion of the distributions. All these effects will be evaluated in detail in a forthcoming paper. The other methods for obtaining the combinations give within the error bars compatible results with the ones that are presented here. The dotted lines in Fig. 1 are the theoretical cross sections calculated using the two structure functions measured in this experiment.

Our results for the two structure functions are also compared in Table 1 with different model predictions: a heavy-baryon chiral perturbation theory calculation (HBChPT) , the linear sigma model (LSM) , an effective lagrangian model (ELM)  and two non relativistic constituent quark models (NRCQM) ). The results seem to favour the HBChPT which at least for $`P_{LT}`$ may not be so surprising since this approach is based on chiral symmetry and therefore should correctly describe the pion cloud which dominates $`P^{(11,00)1}`$ and $`P^{(11,02)1}`$, as it was pointed out in the introduction. Note that the LSM model also respects chiral symmetry but it contains otherwise too restrictive hypotheses. The apparent agreement of the NRCQM models for $`P_{LT}`$ is accidental. It is known that at $`Q^2=0`$ these models overshoot the experimental value of $`P^{(11,11)0}`$ by at least a factor 2.

To summarize, results are reported on the first VCS experiment which allows the determination of two dynamical structure functions which are a linear combination of the GPs. The measurement has been performed at MAMI at $`Q^2`$ = 0.33 GeV2 and gives results compatible with chiral perturbation predictions. Two other experiments, aiming at the determination of the same structure functions are performed at different $`Q^2`$: the Jefferson Lab. experiment  at $`Q^2`$ = 1 and 2 GeV2 and the MIT-Bates experiment  at $`Q^2`$ = 0.05 GeV2. However, to measure independently the 6 GPs and separate non-spin-flip and spin-flip contributions, it will be necessary to perform double polarization experiments.

The authors wish to acknowledge the excellent support of the accelerator group of MAMI. This work was supported in part by the French CEA and CNRS/IN2P3, the Deutsche Forschungsgemeinschaft (SFB 201), the FWO-Flanders (Belgium), the BOF-Gent University, the European Commission ERB FMRX-CT96-0008, the US DOE and the US NSF.


format: text

text:
# Ionic and electronic structure of sodium clusters up to N=59

## I Introduction

Since the pioneering experiments of Knight et al. and their interpretation in terms of the jellium model, sodium clusters have attracted great attention, both experimentally and theoretically. This is due to the fact that sodium is the “simple metal” par excellence, and thus is best suited for the study of fundamental effects. From photoabsorption experiments it is known that small sodium clusters have overall shapes that strongly vary with size and are determined by electronic shell effects. On the other hand, cold clusters with several thousand atoms build icosahedra, i.e., they show ionic shell effects, whereas the bulk material crystallizes in a bcc lattice. Recent experiments indicate that both ionic and electronic degrees of freedom play a role in determining the structural and thermal properties of clusters with several tens of atoms.

On the theoretical side, there have been different approaches to obtain an understanding of the delicate interplay between ions and electrons that is the source of this variety. On one hand, relatively transparent models like the shell and jellium models in several levels of sophistication have been widely used. On the other hand, quantum chemistry and density functional theory offer methods to study clusters on the highest level of sophistication presently possible. But ab initio calculations in the strictest sense, i.e., taking all electrons into account, have only been performed for a few selected cases for the smallest clusters, due to the enormous computational effort that they require. Taking only the valence electrons into account reduces the complexity considerably, but the expense for searching low-energy configurations in three dimensions  still grows rapidly with the system size. The largest ab initio studies of sodium clusters to our knowledge are the recent finite-temperature investigations presented in Refs. .

Several models have been developed to bridge the gap between the ab initio calculations and the shell and jellium models. The ”Spherically Averaged Pseudopotential Scheme” (SAPS) describes the ions by pseudopotentials, in most cases local ones, and the valence electrons are restricted to spherical symmetry. Models based on a volume-averaged or perturbational treatment of ionic effects improve on the treatment of the electrons. Yet further approaches are the Hückel and related models, molecular dynamics based on empirical potentials, and recently, also the extended Thomas-Fermi model combined with a local pseudopotential has been used to study sodium clusters.

From the above examples, two points become clear that considerably facilitate systematical studies of clusters with more than twenty atoms. First, the fact that many calculations make use of phenomenological local pseudopotentials shows the importance of such potentials. This is especially true for sodium, where previous investigations have shown that a local electron-ion interaction can be a good approximation. But especially in cluster physics, where one bridges the region between the atom and the bulk, it is important that a pseudopotential give reliable results for all sizes despite its locality. Second, one needs models that make it computationally manageable to calculate structures and optical properties of clusters systematically for a wide range of considerable sizes, but which on the other hand leave the underlying physics intact and are detailed enough so that relevant information can be drawn from them. Besides extending the computational range, such models will serve the even more important purpose to distill the dominant physical effects from the wealth of details that fully ab initio calculations supply. The results of the model calculations of course must be verified in calculations of higher accuracy and by comparison with experiment. These tasks are addressed in the present article.

In Section II we develop a new phenomenological pseudopotential which meets the just mentioned requirements. Detailed comparisons with ab initio calculations for the smallest clusters in section III verify the validity of the pseudopotential and an axially averaged density functional model for structure calculations. In section IV, a systematic survey over structures and photoabsorption spectra of Na clusters up to $`\mathrm{Na}_{59}^+`$ is given. The results of the survey are summarized and discussed in the concluding section V.

## II A local pseudopotential for structure calculations

Rigorous pseudopotentials in the sense of Phillips and Kleinman  and modern ab initio pseudopotentials are always nonlocal. However, it has been noted early in the development of pseudopotential theory that by relaxing the Phillips-Kleinman condition, one can open up a new class of pseudopotentials. These are also termed “model potentials” because they are constructed by choosing some analytical functions as models for the partial-wave potentials and adjusting their parameters such that some experimentally known quantities, e.g., atomic energy levels, are matched. If several partial-wave potentials can be chosen to be the same, one can construct a local model potential. In fact, such local phenomenological potentials have been and are being used successfully in many branches of physics, see, e.g., Refs. . However, the local approach can only be expected to be good for the so-called “nearly free electron” metals whose atomic structure consists of only $`s`$ or $`p`$ electrons outside of an ionic core with a noble gas configuration. (Lithium, with its lack of $`p`$ electrons in the core, is the prominent exception.) In these metals, the two contributions making up the pseudopotential have a tendency to cancel each other, and the effective potential a valence electron experiences is further diminished by screening effects. Therefore, model potentials have been proposed in the past which tried to exploit this simple electronic structure either by fitting properties of the bulk solid, or of the single atom, or by ad hoc fulfilling desired numerical properties. However, some of these potentials can lead to wrong predictions when they are used in physical surroundings different from the one in which they were set up, or when properties other than the adjusted ones are looked at. With an emphasis on solid state properties, the last point has been discussed in detail previously, see, e.g., Ref.  for an overview. Our aim here is to develop a local model potential that gives a maximum degree of transferability in the sense that potentials constructed according to our scheme should reproduce the important physical properties of a system, irrespectively of its number of atoms or the way in which these are arranged.

### A The Ansatz

The construction of a model potential consists of two steps. The first is to choose the model function. It should allow for a physical interpretation of the final potential, and at the same time should have analytical and numerical properties that allow for an easy application. The parameterization

$$V_{\mathrm{ps}}\left(r\right)=\frac{Ze^2}{r}\left\{c_1\text{erf}\left(\frac{r}{\sqrt{2}\sigma _1}\right)+c_2\text{erf}\left(\frac{r}{\sqrt{2}\sigma _2}\right)\right\},$$
(1)
where

$$\text{erf}(x)=\frac{2}{\sqrt{\pi }}_0^x𝑑y\mathrm{exp}(y^2),$$
(2)
certainly meets the second requirement, since the error function can very efficiently be handled numerically and yields a smooth representation on a grid. That this parameterization also has a transparent physical interpretation will be demonstrated at the end of this section. The second step is the choice of the four parameters $`\sigma _1,\sigma _2,c_1,c_2`$. One necessary condition is that the correct Coulomb limit,

$$\underset{r\mathrm{}}{lim}V_{\mathrm{ps}}\left(r\right)=\frac{Ze^2}{r},$$
(3)
be obtained, which requires

$$c_1+c_2=1.$$
(4)
Thus we are left with three free parameters. Since our aim is to construct a pseudopotential that gives reliable properties for all clusters, i.e., spanning the region from the atom to the bulk, the most important properties of both atom and bulk solid must be reproduced. Therefore, we have chosen to fit the Wigner-Seitz radius $`r_\mathrm{s}`$ and the compressibility $`B`$ of the crystalline metal, together with the energy of the atomic $`3s`$ level $`e_\mathrm{a}`$. These quantities characterize sodium and influence both structure and electronic excitations, thus being of great importance for reliable results. A test for whether we really have captured the relevant physics will be to check if non-fitted quantities (bond-lengths, atomic spectra, bulk binding, dipole resonances) are also reproduced correctly.

### B Determination of basic properties

We calculate the bulk properties in second-order perturbation theory. The unperturbed system is the noninteracting homogeneous gas of valence electrons, the perturbation is given by the potentials arising from the crystal lattice and from the interaction of the electrons with each other. Each crystal ion is described by a pseudopotential $`V_{\mathrm{p}s}`$ centered on a lattice site. Up to second order, the binding energy per valence electron $`e_\mathrm{b}`$ is given by

$`e_\mathrm{b}(r_\mathrm{s})`$ $`=`$ $`e_{\mathrm{k}in}(r_\mathrm{s})+e_{\mathrm{x}c}(r_\mathrm{s})+e_{\mathrm{p}s1}(r_\mathrm{s})+e_\mathrm{h}(r_\mathrm{s})+e_{\mathrm{i}i}`$ (6)
$`+e_{\mathrm{b}s}(r_\mathrm{s}).`$
Here $`r_\mathrm{s}=\left[3/\left(4\pi n\right)\right]^{\frac{1}{3}}`$ is defined in terms of the average valence electron density $`n`$,

$$e_{\mathrm{k}in}(r_\mathrm{s})=\frac{3\mathrm{}^2}{10m}\left(\frac{9\pi }{4}\right)^{\frac{2}{3}}\frac{1}{r_\mathrm{s}^2}$$
(7)
is the noninteracting kinetic energy, and $`e_{\mathrm{x}c}(r_\mathrm{s})`$ is the exchange and correlation energy for which we have employed the LDA functional of Perdew and Wang. The quantity

$$e_{\mathrm{p}s1}(r_\mathrm{s})=\frac{3}{Z4\pi r_{\mathrm{s}}^{}{}_{}{}^{3}}V_{\mathrm{p}s}(r)d^3r$$
(8)
is the first-order contribution of the pseudopotential to the binding energy, $`e_\mathrm{h}(r_\mathrm{s})`$ is the Hartree energy of the valence electron density, and $`e_{\mathrm{i}i}`$ is the electrostatic energy of point ions. Finally, $`e_{\mathrm{b}s}(r_\mathrm{s})`$ is the band structure energy discussed below. Integrations are taken over the (infinite) volume of the crystal. These contributions are rewritten as

$`e_{\mathrm{p}s1}(r_\mathrm{s})+e_\mathrm{h}(r_\mathrm{s})+e_{\mathrm{i}i}(r_\mathrm{s})=`$ (11)
$`{\displaystyle \frac{3}{Z4\pi r_\mathrm{s}^3}}{\displaystyle \left(V_{\mathrm{p}s}+\frac{Ze^2}{r}\right)d^3r}+`$
$`{\displaystyle \frac{3}{Z4\pi r_\mathrm{s}^3}}{\displaystyle \left(\frac{V_\mathrm{h}}{2}\frac{Ze^2}{r}\right)d^3r}+e_{\mathrm{i}i}`$
to obtain the volume-averaged repulsive part of the pseudopotential that defines its strength

$$𝒮_{\mathrm{p}s}=\left(V_{\mathrm{p}s}+\frac{Ze^2}{r}\right)d^3r,$$
(12)
and the total electrostatic energy of point ions in a compensating uniform negative background, called the Madelung energy $`e_\mathrm{m}`$. Separating the Coulomb force into long and short-range components, one can calculate $`e_\mathrm{m}`$, and for the bcc lattice one obtains $`e_\mathrm{m}=0.895929Z^{\frac{2}{3}}e^2/r_\mathrm{s}`$. Thus the only characteristic of the pseudopotential that enters the binding energy in first-order perturbation theory is its strength $`𝒮_{\mathrm{p}𝓈}`$. For fixing the radial dependence of the potential it is therefore essential to take into account the second-order band structure energy

$$e_{\mathrm{b}s}(r_\mathrm{s})=\frac{1}{2}\frac{3}{Z^24\pi r_\mathrm{s}^3}\underset{𝐪\mathrm{𝟎}}{}\left|\stackrel{~}{V}_{\mathrm{p}s}\left(𝐪\right)S\left(𝐪\right)\right|^2\frac{\chi \left(𝐪\right)}{ϵ\left(𝐪\right)}.$$
(13)
In (13),

$$S\left(𝐪\right)=\frac{1}{N}\underset{j=1}{\overset{N}{}}\mathrm{exp}\left(i\mathrm{𝐪𝐑}_j\right)$$
(14)
is the structure factor with the sum running over all ionic positions. For a lattice without basis one obtains $`S\left(𝐪\right)=\delta _{𝐪,𝐊}`$, where $`𝐊`$ denotes reciprocal lattice vectors, i.e., the sum in (13) is running over reciprocal lattice vectors only. Furthermore,

$$\chi \left(𝐪\right)=\frac{mk_F}{2\pi ^2\mathrm{}^2}\left(1+\frac{1\eta ^2}{2\eta }\mathrm{ln}\left|\frac{1+\eta }{1\eta }\right|\right),\eta =\frac{q}{2k_F}$$
(15)
is the Lindhard function with $`k_F=(3\pi ^2n)^{\frac{1}{3}}`$, and

$$ϵ\left(𝐪\right)=1\frac{4\pi e^2}{q^2}\left(1𝒢\left(𝐪\right)\right)\chi \left(𝐪\right)$$
(16)
is the dielectric function including the local field correction

$$𝒢=\frac{q^2}{4\pi e^2}\frac{d^2}{dn^2}\left(ne_{\mathrm{x}c}\left[n\right]\right).$$
(17)
in the LDA. The pseudopotential enters via its Fourier transform, which for our potential is given by

$$\stackrel{~}{V}_{\mathrm{p}s}\left(𝐪\right)=\frac{4\pi e^2Z}{q^2}\underset{i=1}{\overset{2}{}}c_i\mathrm{exp}\left(\frac{q^2\sigma _i^2}{2}\right).$$
(18)
From Eqs. (6) to (16) we calculate $`e_\mathrm{b}`$; the minimum of $`e_\mathrm{b}`$ determines $`r_\mathrm{s}`$, and the bulk modulus is obtained from

$$B=V\frac{P}{V}=\frac{1}{12\pi r_\mathrm{s}}\left(\frac{^2e_\mathrm{b}}{^2r_\mathrm{s}}\frac{2}{r_\mathrm{s}}\frac{e_\mathrm{b}}{r_\mathrm{s}}\right).$$
(19)
The reciprocal lattice vectors are generated numerically from the reciprocal basis, and we carefully checked that in all calculations the sum over reciprocal lattice vectors was numerically converged.

The atomic $`3s`$ level is the lowest eigenvalue $`e_\mathrm{a}`$ of the Schrödinger equation

$$\left(\frac{^2}{r^2}+\frac{l\left(l+1\right)}{r^2}+\frac{2m}{\mathrm{}^2}V_{\mathrm{p}s}\left(r\right)\frac{2m}{\mathrm{}^2}e_\mathrm{a}\right)u\left(r\right)=0$$
(20)
for $`l=0`$, where the pseudo wavefunction has, as usual, been factorized into radial and angular components, $`\psi (r)=\left[u\left(r\right)/r\right]Y_{lm}(\vartheta ,\phi )`$. At this point, a subtlety should be considered. The atomic energy calculated within density functional theory using the LDA will slightly differ from $`e_\mathrm{a}`$ due to the well known fact that in the LDA, the self interactions in the Hartree and exchange energy do not cancel each other exactly and leave a spurious self interaction energy. One therefore might be tempted to “compensate” this self-interaction energy by adjusting the pseudopotential parameters such that the experimental value is matched when the self-interaction energy is included. This, however, would be dangerous for several reasons. First, it must be recalled that for an accurate description of bonding properties, the pseudopotential must lead to an accurate description of the electronic density. If one adjusts the pseudopotential parameters such that a spurious energy is compensated, then the valence electron density resulting from this pseudopotential might correspondingly show spurious, unphysical deformations, leading to wrong bonding properties, as discussed previously. Secondly, the self-interaction energy becomes smaller with increasing delocalization of the electrons. For an electron delocalized over a volume $`\mathrm{\Omega }`$, the self-interaction correction is of order $`\mathrm{\Omega }^{\frac{1}{3}}`$. Since the valence electrons in metal clusters are delocalized to a high degree, and since the cluster volume changes noticeably in the size range from N=8 to N=58 that we are interested in, a compensation of the atomic self interaction energy via the pseudopotential would lead to problems for increasing cluster sizes. And third, building the self-interaction energy into the pseudopotential would be inconsistent with the procedure of fitting to the bulk, because the bulk calculation is less affected by the self-interaction error due to the delocalization of the bulk electrons. Therefore, it is better to exactly solve Eq. (20). This can straightforwardly be done numerically by combining the Runge-Kutta method with adaptive stepsize control with a globally convergent Newton scheme. As a welcomed side effect, our pseudopotential might also be usable in self-interaction free, i.e., beyond-LDA density functional calculations.

### C Results and comparison

The three free parameters of our model potential (1) were now chosen such that the experimental low-temperature value $`r_\mathrm{s}=3.93a_0`$ be reproduced exactly, while at the same time $`e_\mathrm{a}`$ and $`B`$ be as close as possible to the experimental values $`e_\mathrm{a}=0.378Ry`$ and $`B=0.073Mbar`$. These conditions lead to the parameters

$`\sigma _1=0.681\mathrm{a}_0,`$ $`c_1=2.292,`$ (21)
$`\sigma _2=1.163\mathrm{a}_0,`$ $`c_2=3.292.`$ (22)
In Table I we have listed the resulting $`r_\mathrm{s}`$, $`B`$ and $`e_\mathrm{a}`$ for our smooth-core pseudopotential and for other local pseudopotentials that have widely been used in cluster physics.

The empty-core potential with both of the most frequently used choices for its cut-off radius $`r_c`$, and the pseudopotential of Ref.  that was constructed in the same spirit and adjusted in first-order perturbation theory only, lead to considerable deviations from the experimental values for all quantities. The local Heine-Abarenkov potential gives a reasonable $`r_\mathrm{s}`$ and an $`e_\mathrm{a}`$ very close to the experimental one, but 10% error in $`B`$. Our pseudopotential by construction gives no or only very small differences from the experimental values for $`r_\mathrm{s}`$, $`B`$ and $`e_\mathrm{a}`$, showing that it is possible to obtain reasonable results for all these quantities with a local potential. However, a severe test will be whether also non-fitted quantities are reproduced correctly. To this end, we have calculated the bulk binding energy $`e_\mathrm{b}`$, the dimer binding length $`d_{\mathrm{d}imer}`$, the zero of the pseudopotential form factor, given by $`q_0=\sqrt{2[\mathrm{ln}(c_2)\mathrm{ln}(c_1)]/(\sigma _2^2\sigma _1^2)}`$ in our case, and the first seven excited atomic energy levels. Table II shows that also these ten non-fitted quantities come out very close to the measured values, revealing that the model potential defined by (1), (21) indeed incorporates the relevant physical effects.

Fig. 1 shows our pseudopotential in real and in Fourier space, leading us back to the question of the motivation for the ansatz (1). From Eq. (4), one can interpret the coefficients $`c_1`$ and $`c_2`$ as the charges associated with the attractive and repulsive terms, respectively, in the pseudopotential. This interpretation becomes more transparent when one looks at the corresponding pseudocharge density. The pseudocharge $`n_{\mathrm{ps}}`$ is related to the pseudopotential via Poisson’s equation

$$\mathrm{}V_{\mathrm{ps}}\left(r\right)=4\pi e^2n_{\mathrm{ps}}\left(r\right),$$
(23)
and for our pseudopotential is given by

$$n_{\mathrm{p}s}(r)=n_1\mathrm{exp}\left(\frac{r^2}{2\sigma _{1}^{}{}_{}{}^{2}}\right)+n_2\mathrm{exp}\left(\frac{r^2}{2\sigma _{2}^{}{}_{}{}^{2}}\right),$$
(24)
where

$$n_i=c_i\left(2\pi \sigma _{i}^{}{}_{}{}^{2}\right)^{\frac{3}{2}}\mathrm{f}ori=1,2.$$
(25)
Thus, our parameterization describes the ionic core by two overlapping Gaussian charge densities. The different heights and widths of the Gaussians result in a pseudodensity that is negative for small distances, corresponding to a repulsive core, and has a positive tail that compensates the core for larger distances. Equation (24) thus is a generalization of the parameterization that was used in our previous studies, where the pseudodensity has a two step-profile. Besides its transparent physical interpretation, the pseudopotential (1) has excellent numerical properties: in real space the potential is smooth with no extreme peaks and can very efficiently be handled numerically via the pseudodensity, allowing to solve the Coulomb problem for electrons and ions together. At the same time, its rapidly converging Fourier transform makes the potential equally well suited for calculations in reciprocal space.

## III Comparison with ab initio results

In the present work we have employed the “Cylindrically Averaged Pseudopotential Scheme” (CAPS), which has been introduced earlier and which has been modified and improved for the present purposes. Since CAPS treats the electronic degrees of freedom in the Kohn-Sham formalism without limiting them to spherical symmetry, explicitly includes the ionic structure, and furthermore is numerically efficient, it is a very good tool for systematical studies of the interplay between ions and electrons. But before CAPS and the smooth-core pseudopotential are applied on a large scale, we further test their usefulness for cluster structure calculations.

That the smooth core pseudopotential reproduces the experimentally known dimer bond length has already been shown in Table II. For the next larger clusters, no direct experimental information on structures or bond lengths is available, and we therefore resort to comparisons with other theoretical work. Fig. 2 shows the four smallest Na clusters with even electron numbers as they are obtained in CAPS with the smooth core pseudopotential. For these clusters, also three-dimensional geometry optimizations have been performed using ab initio pseudopotentials and Hartree-Fock with configuration interaction (CI) , or DFT with the LDA  for the valence electrons, respectively. Also all-electron Hartree-Fock calculations have been reported. CAPS finds exactly the same structures as the three-dimensional methods, and this in spite of the fact that the ionic configurations of $`\mathrm{Na}_4`$

and $`\mathrm{Na}_5^+`$ might let the cylindrical averaging seem a rather far-fetched approximation. But not only the overall geometry, also the bond lengths are in good agreement with the ones from the three-dimensional calculations. Whereas all-electron Hartree-Fock overestimates bond lengths due to missing correlation effects, the three-dimensional LDA calculations lead to the well-known underestimation. By construction, the phenomenological pseudopotential compensates this underestimation, and it is thus reassuring to see that the resulting bond lengths are close to the ones found in the quantum chemistry calculations. Three-dimensional calculations have been performed for a few other neutral clusters, and CAPS is in agreement with the three-dimensional geometries in these cases, too. A detailed analysis of structures of neutral clusters and their static electric polarizability can be found in Ref. .

A further test for structure calculations is obtained by comparing the photoabsorption spectra corresponding to the theoretically determined structures to the measured ones. To this end, we have calculated the valence electron excitations for our cluster structures in a collective approach. The basic idea of this method has been presented earlier, and a detailed discussion of its extensions and the relation to density functional theory is the subject of a forthcoming publication. Here, we demonstrate the validity of our collective model for the test cases $`\mathrm{Na}_2`$ and $`\mathrm{Na}_8`$, where experimental data and recent ab initio TDLDA calculations are available for comparison. The first row of plots in Fig. 3 shows the excitation spectra obtained with the axial collective model for the CAPS structures

shown in the insets, the second row shows the TDLDA spectra for the corresponding three-dimensional geometries, and the bottom row shows the experiment. The lower two rows were directly adapted from Ref. , and the collective model results were then plotted on the same scale. Besides small differences in the two small peaks seen at high energies for $`\mathrm{Na}_2`$, the agreement between the three sets of data is very good. This demonstrates that our collective model is capable of correctly describing excitations even in small systems. It further is to be noticed that the excitation energies obtained with the smooth-core pseudopotential are $``$ 0.15 eV lower than the TDLDA energies and closer to the experimental ones, which can be attributed to the larger bond lengths that result with the present pseudopotential.

## IV Structures and photoabsorption spectra of singly-charged cationic sodium clusters

Since the comparisons in the previous section have shown that the phenomenological pseudopotential, the collective model and CAPS are well suited for the description of sodium clusters, a systematic study of clusters with even electron numbers between 8 and 58 is presented. For each cluster a great number of simulated annealing runs was started from different random configurations. The search was continued until new runs no longer returned new low-energy structures. Although it must always be kept in mind that no practical algorithm guarantees that all low-energy structures are found, this extensive and unbiased procedure at least gives good hope to do so. The aim of this survey is to investigate how electronic and ionic shell effects play together in determining the cluster structure, and thus gain a better understanding of how matter grows. The geometry optimization was done with CAPS, but the energies of the resulting structures have also been obtained in three-dimensional Kohn-Sham calculations to check the ordering of isomers without axial restriction. Photoabsorption spectra are calculated with the collective model, and comparison with the experimental spectra gives further information on the relevant structures.

The “magic” cluster $`\mathrm{Na}_9^+`$ is spherical according to the jellium model, and the ground state found with CAPS is the $`C_{4v}`$ structure (a) in Fig. 4. Separated by 0.10 eV, CAPS finds the $`D_{4d}`$ isomer (b), and both these geometries were also found in three-dimensional CI and LDA calculations. In addition, CAPS finds the third isomer (c) which is higher by 0.23 eV. All these clusters have nearly spherical valence electron densities, and therefore support the shell model picture. But the ionic configurations are, of course, non-spherical, and this is reflected in the photoabsorption cross sections.

The first thing to be noted is the overall position of the resonance, about which there has been a long-standing debate. In previous work we have shown that when ionic structure is taken into account via a pseudopotential, the detailed form of the potential greatly influences the resonance position. But whereas our previous calculations were mainly compared to high-temperature data, and, due to limited collective model basis sets, fixed the resonance positions only within a few percent (as pointed out in Ref. ), we have now fully converged the basis. The resulting overall resonance positions are very close to the experimental low-temperature positions, showing that consistent inclusion of ionic structure gets the plasmon position in the correct energy interval. The optical response of $`\mathrm{Na}_9^+`$ has also been theoretically investigated previously, and each of these calculations explains the experimental spectrum on the basis of a different isomer. In our calculations, the spectra of isomers (b) and (c) both reproduce the splitting of the main peak seen in the experiment. Isomer (a) does not show this feature, but it leads to the small sub-peak around 3.5 eV that is also seen experimentally. The situation becomes transparent when the results from our three-dimensional calculations, see Table III, are taken into account. They show that isomer (c) is degenerate with isomer (a), and the experimental spectrum can thus be explained even at low temperatures as resulting from a mixture of isomers (a) and (c). (The importance of isomerism on dipole spectra in small sodium clusters has also been pointed out in Ref. .)

The lowest energy structure that we find for $`\mathrm{Na}_{11}^+`$ ($`D_{4d}`$, labeled (a) in Fig. 5) can be understood as resulting from placing one atom in the center of $`\mathrm{Na}_{10}`$, or capping isomer (b) of $`\mathrm{Na}_9^+`$ on the quadratic faces. Separated by 0.11 eV and 0.19 eV, respectively, CAPS finds the $`D_{3h}`$ isomer (b), and its deformed counterpart (c), which is only a shallow local minimum in CAPS and easily transforms into the more stable structure (b). In the three-dimensional calculations, the ordering of isomers is reversed: (b) becomes the ground state,

and our results are thus consistent with the findings of Ref. . We have calculated the photoabsorption spectra with the collective model and obtain two resonances with heights approximately 1:2 for all three geometries. These are centered at 2.41/2.94 eV for structure (a), at 2.39/2.94 eV for structure (b), and at 2.27/2.98 eV for structure (c). The experimental photoabsorption spectrum  of $`\mathrm{Na}_{11}^+`$ at 380 K shows two broad peaks with heights 1:2 around 2.2 eV and 2.85 eV, and a pattern of six peaks when measured at 35 K: three low intensity ones at 1.9 eV, 2.2 eV and 2.4 eV, and three high intensity ones at 2.5 eV, 2.8 eV and 3.0 eV. Obviously, the calculated resonances are blue shifted by a few percent with respect to the hot experiment. This is understandable since the calculations have been done for T=0 K, and thermal expansion of the cluster shifts the plasmon to slightly lower energies. (This effect has recently been put into evidence quantitatively for the static response.) When compared to the cold data, the calculated resonances are in the correct energy range, but the fine structure that the experiment shows is not reproduced. A CI calculation  based on isomer (b) reproduced some of the experimentally observed patterns, but also could not explain all the peaks, and recent three-dimensional TDLDA calculations  have lead to similar results as the collective model. An explanation of all details in the experimental spectrum therefore might require to consider a mixture of low-energy structures, and this will be discussed in detail in a separate publication.

For $`\mathrm{Na}_{15}^+`$ CAPS finds four low-energy structures. The two lowest ones, (a) and (b) in Fig. 6, are very close in energy in both CAPS and the three-dimensional calculations. (a) corresponds to a distorted $`\mathrm{Na}_{11}^+`$ (a) with four atoms added along the $`z`$-axis, and (b) can be understood as a pentagonal bipyramid sharing one ion with $`\mathrm{Na}_9^+`$ (a). Structure (c) results from structure (a) by moving one of the inner single atoms to the bottom face, and the oblate (d) corresponds to a configuration found for $`\mathrm{Na}_{14}`$, but with one atom added along the symmetry axis. In agreement with the experimental photoabsorption spectrum, our results point to prolate clusters as the relevant structures and thus again verify the prediction of the deformed, structure averaged jellium model (SAJM).

The ground state structure found for $`\mathrm{Na}_{17}^+`$, (a) in Fig. 7, is close to the one found in the extended Hückel model for the neutral cluster. It consists of an icosahedral core with a crown of four atoms, and the only difference to the result of Ref.  is that the cylindrical averaging forces the crown atoms into a square. This shows that for such sizes the differences between charged and neutral clusters can already be small. Structure (b), which is slightly higher in energy than (a) in the three-dimensional calculations, has an even more prolate valence electron density, and the collective resonances for these clusters lie at 2.59 eV and 2.86 eV \[structure (a)\] and at 2.34 eV and 2.95 eV \[structure (b)\]. They are thus in contrast to the experimental spectrum that points to an oblate cluster. CAPS finds one more isomer, structure (c), that indeed leads to an oblate valence electron density.

In CAPS, it is 0.16 eV higher than (a), with a half occupied orbital. The collective spectrum for this isomer is very close to the experimental one, as shown in Fig. 7. However, also the three-dimensional calculation gives a non-negligible energy difference between (c) and the prolate isomers, cf. Table III. Therefore, in this case a three-dimensional relaxation of the ions would be necessary to check whether Jahn-Teller distortions lower the energy of isomer (c) so much that it can account for the experimental spectrum.

For $`\mathrm{Na}_{19}^+`$ CAPS finds the four stable geometries shown in Fig. 8. The double-icosahedron (a) is the ground state, and it is separated from structure (b) by 0.23 eV, from (c) by 0.30 eV, and from (d) by 0.35 eV. However, in the calculations without axial restriction, the energetic differences are considerably reduced: isomer (d) becomes equivalent to isomer (b), and both are separated from the ground state by only 0.12 eV. The experimental photoabsorption spectrum sheds further light on the situation. It shows a high peak at about 2.7 eV, followed by a lower one at about 2.9 eV, and is thus very similar to the collective spectrum of isomer (d). The spectra of the other isomers all give the peaks in reversed order and thus do not resemble the experiment. (The collective spectrum of (c), which is not shown in order not to overload Fig. 8, shows a small peak at 2.84 eV and a higher one at 2.93 eV.) Since an internal cluster temperature of about 60 K is sufficient to go from isomer (a) to (d), and since the temperature in the experiment is about 105 K,

it is plausible that isomer (d) will contribute to the spectrum. The experimental spectrum however looks as if it is dominated by isomer (d). One explanation for this can be that structure (d) has a larger “catchment area” than (a) and is entropically favored. The second possibility, which is suggested by the relatively large difference between the energy found in CAPS and the three-dimensional calculation, is that a three-dimensional relaxation of the ions would turn structure (d) into the ground state.

Two degenerate structures are found for $`\mathrm{Na}_{21}^+`$ in CAPS, and also the three-dimensional calculation leads to nearly equal energies. Structure (a) results from the double icosahedron of $`\mathrm{Na}_{19}^+`$ (a) by adding one ion to each of the lower tow pentagons, and (b) results from (a) by moving the lowest ion to the uppermost pentagon. Both geometries lead to collective spectra with two main transitions and some strength above the main resonances. This is in agreement with the experimental data  that also show this high-energy tail. Concerning the main transitions, the two isomers are somewhat different: in structure (a) the main peaks are separated by about 0.1 eV, whereas they nearly fall together for structure (b). Since the structures are extremely close in total energy and since the main peaks of (b) energetically fall exactly between the main peaks of (a), the experimental spectrum can consistently be explained as resulting from a mixture of both isomers.

The simulations performed for $`\mathrm{Na}_{25}^+`$ led to strongly prolate shapes. This is partially understandable on the basis of the SAJM:  the valence electrons force the cluster into the prolate regime. But the inclusion of ionic structure even enhances the prolate tendency, because the lowest configuration, shown in Fig. 10, is a triple icosahedron, and in order to build this configuration, a quadrupole moment larger than the jellium prediction is necessary. From this observation it also becomes clear that the pentagonal bipyramid which has been observed as a building block of the smaller clusters is also important for the larger sizes. For the triple icosahedron, collective resonances are observed at 2.44 eV and 2.96 eV, and also the experiments  indicate peaks at these energies. But the latter shows an additional peak around 2.75 eV which is not found in the calculations. Comparing this to the situation encountered for $`\mathrm{Na}_{27}^+`$ shows that the geometries of the two clusters are closely related: $`\mathrm{Na}_{27}^+`$ is the triple icosahedron with two atoms added to the central pentagons. This shows the consistency of the structure calculations, and for $`\mathrm{Na}_{27}^+`$ the collective spectrum nicely matches with the experimental data. The triple icosahedron is energetically strongly favored over the other structures that CAPS finds. All of them are variations of the ground state geometry with local distortions, and a three-dimensional calculation was performed for the structure that in CAPS is the second lowest. Here, the difference in total energy is a little smaller than in CAPS, but still large, cf. Table III. Since furthermore all low-energy structures give rise to collective spectra that are similar to the one shown in Fig. 10, the explanation for the middle peak observed for $`\mathrm{Na}_{25}^+`$ remains an open question.

Comparing the structures of $`\mathrm{N}a_{25}^+`$ and $`\mathrm{N}a_{27}^+`$ to the ones found for $`\mathrm{N}a_{31}^+`$, $`\mathrm{N}a_{41}^+`$ and $`\mathrm{N}a_{43}^+`$ (Fig. 11, Fig. 12, and Refs. ) allows to identify a growth pattern: the triple icosahedron serves as a building block for the larger clusters. Additional atoms are added on outside faces, as seen in the ground states of $`\mathrm{N}a_{27}^+`$ and $`\mathrm{N}a_{43}^+`$, or can be packed between the pentagonal subunits, as observed for the isomers of $`\mathrm{N}a_{27}^+`$, where the two additional ions are placed at different positions “inside” the $`\mathrm{N}a_{25}^+`$ structure.

The collective spectra of the low-energy structures (a) – (c) of $`\mathrm{Na}_{31}^+`$ reflect the prolate electron densities of these clusters and show one peak around 2.6 eV, carrying about 30% oscillator strength, another one around 2.95 eV carrying about 50% strength, and the rest of the strength scattered at higher energies. Structure (d) is not stable against Jahn Teller distortions in three dimensions. In view of the cluster size, the CAPS differences in total energy for $`\mathrm{Na}_{31}^+`$ are rather small, and they are even smaller in the three-dimensional approach. Since the temperature in the measurement of the photoabsorption cross section was again at least 105 K, it can be concluded from the calculations that a variety of isomers will contribute to the experimental spectrum. It is thus not astonishing that the available experimental data  do not resolve separated peaks but show a rather broad bump.

The pronounced deformation seen for $`\mathrm{N}a_{43}^+`$ in Fig. 12 results from an interplay between ions and valence electrons. The electrons “push” the cluster into the octupole, i.e. pear-shaped form, and the ions arrange under this “constraint”. But the ions favor the icosahedral core, and this increases the deformation. The octupole moment therefore is larger by a factor 1.15 than in the SAJM. The photoabsorption spectrum calculated for $`\mathrm{N}a_{43}^+`$ shows a strong peak at 2.7 eV, followed by a high-energy tail. It is in close agreement with the experimental data, which thus support our structure calculations.

Based on the shell and deformed jellium model a configuration would be expected for $`\mathrm{Na}_{55}^+`$ that gives rise to a prolate valence electron density. The CAPS calculations, however, consistently led to nearly spherical or slightly oblate clusters as low-energy configurations. The lowest energy was found for structure (a) in Fig. 13. It has a fivefold symmetry axis and again shows the close-to-icosahedral core discussed previously. Within 0.4 eV structures (b) and (c) are found with valence electron densities that are more oblate. These structures show half-occupied highest orbitals, indicating that they would Jahn-Teller relax if the axial restriction on the electrons were lifted. The lowest prolate isomer that was found, structure (d), has a quadrupole moment that is close to the one predicted by the SAJM  but in CAPS this configuration is 0.73 eV higher than the ground state. The annealing was also started from an icosahedron with the nearest-neighbor distance of bulk sodium. In this calculation, the strictly geometrical icosahedron transforms into structure (a), which can be regarded as a slightly distorted icosahedron. The experimental photoabsorption spectrum of $`\mathrm{Na}_{55}^+`$ as measured by the Freiburg group is shown in the lowest panel of Fig. 13, and the spectrum has also been measured by Meibom et al. Both spectra have in common that one broad peak is seen, followed by a second, smaller peak or a high-energy shoulder. This is in contrast to the prediction of the jellium model, because the prolate structures found there give the peaks in reversed order, similar to the spectrum of isomer (d). The CAPS low-energy structures, however, lead to collective spectra with a high peak followed by a lower one, and are thus in better agreement with the experiment. The total oscillator strength measured in the experiment  was 70% - 80%. This also agrees with the collective model calculations that find 79% strength within the plotted range for structure (a), 76% for (b) and 80% for (c). As a test for the collective model, a TDLDA calculation with excitation in z-direction was performed for structure (a). In the TDLDA, the high-energy shoulder is more pronounced than in the collective model. To further clear up the situation, the structures (a), (b) and (d) were relaxed in Born-Oppenheimer ab initio molecular dynamics, i.e., fully three-dimensional with the Troullier-Martins pseudopotential. The bond lengths of all three structures shrinked by a few percent due to the different pseudopotential, and isomers (b) and (d) distorted slightly, but otherwise the geometries stayed unaltered. As seen in Table III, the differences in total energy are somewhat smaller in the ab initio calculations, but the ordering of isomers is the same as in CAPS. Finally, as reported in Ref. , the CAPS structure (a) was annealed for 10 ps at about 220 K. In this annealing, the cluster became even more similar to an icosahedron and its overall shape oscillated around the nearly spherical shape of the icosahedron. Thus, the ab initio calculations confirm the finding that in contrast to the deformed jellium model prediction, $`\mathrm{Na}_{55}^+`$ at low to moderate temperatures has a close to spherical valence electron density.

The CAPS results for $`\mathrm{Na}_{57}^+`$, see Fig. 14, are consistent with the results for $`\mathrm{Na}_{55}^+`$: again the low-energy structures are not prolate, but spherical or slightly oblate. The collective spectra somewhat underestimate the high-energy shoulder, and this can be attributed to selective particle-hole effects just as in the case of $`\mathrm{Na}_{55}^+`$. But the overall agreement with experiment  is considerably better than in the deformed-jellium calculation. Structure (a) has the lowest energy, followed by isomers (b) and (c) that are higher by 0.22 eV. For structure (c) we observe a half-occupied orbital. With respect to the ionic geometry, structure (a) is similar to the third isomer of $`\mathrm{Na}_{55}^+`$, and (c) corresponds to the ground state of $`\mathrm{Na}_{55}^+`$ with the two additional ions added on top and bottom faces. Another isomer, not shown in Fig. 14 for the sake of clarity, is found, and there the two additional ions are added in the equatorial plane of $`\mathrm{Na}_{55}^+`$ (a).

In the case of $`\mathrm{Na}_{59}^+`$ the search for low-energy structures has not been as exhaustive as for the other clusters. But several annealing runs were started from random configurations, and the low-energy geometries that were found again showed the ions arranged rather regularly and similar to the structures just discussed. Therefore, further simulations were started from the geometries found for $`\mathrm{Na}_{55}^+`$ and $`\mathrm{Na}_{57}^+`$, plus four or two atoms, respectively, that were added at randomly chosen sites. One of these runs led to a geometry that corresponds to $`\mathrm{Na}_{57}^+`$ (a), capped by two atoms on top and bottom. This structure has the lowest energy of all that were found, and it does not show signs of Jahn-Teller instability. The comparison between the experimental  and theoretical photoabsorption spectra fits into the previous discussion. The strongest collective resonance is seen at 2.8 eV, followed by a smaller one at 2.9 eV and a little strength scattered around 3.25 eV, i.e., the collective model leads to qualitative agreement with the experimental data.

## V Discussion and conclusion

The systematic survey showed that for most of the smaller Na clusters, the overall deformation is determined by electronic shell effects even when the ionic structure is explicitly included. This explains the great success of the deformed jellium model. However, the survey at the same time clearly exhibits the limitations of the jellium picture. Besides the fact that some spectroscopic patterns, e.g., as seen for $`\mathrm{Na}_9^+`$, are directly related to details in the ionic configuration, the fundamental improvement brought about by taking into account the ionic structure is that growth systematics can be identified. In the smallest clusters, the pentagonal bipyramid is a frequently appearing structure and is seen, e.g., in $`\mathrm{Na}_7^+`$, $`\mathrm{Na}_{15}^+`$, $`\mathrm{Na}_{17}^+`$. In $`\mathrm{Na}_{19}^+`$, three of these smallest units are combined to build up the double icosahedron, and from then on, the electronic shell effects and the preference of the ionic structure for icosahedral geometries work hand in hand in determining the cluster structure. This is seen most impressively for $`\mathrm{Na}_{43}^+`$. At $`\mathrm{Na}_{55}^+`$, however, the situation changes. Due to the ionic shell closing, the influence of the ions wins over the electronic shell effects, resulting in nearly spherical structures. This finding is supported by the experimental photoabsorption data and has been verified in ab initio calculations.

In summary, we have presented a local pseudopotential for sodium that accurately models the core-valence interaction and correctly describes the atom, the bulk, and finite clusters. We have demonstrated the influence of pseudopotentials on structural properties and the direct influence of the bond lengths on the resonance positions. This shows that even for the most simple metal, a pseudopotential must be constructed carefully. Cluster structures were calculated in axially averaged and three-dimensional Kohn-Sham calculations. A detailed comparison with ab initio work demonstrated the validity of the CAPS as a tool for the approximate determination of cluster configurations. The systematical survey for clusters with up to 59 ions revealed an icosahedral growth pattern. Collective resonances were calculated for the theoretically determined structures, and comparison with the experimental photoabsorption spectra confirmed the results of our structure calculations. Thus, a distinct step in the growing process of matter has been put into evidence, namely the transition from electronically to ionically determined configurations.

###### Acknowledgements.

We thank J. Akola for the communication of his results prior to publication, and one of us (S. Kümmel) further acknowledges stimulating discussions with A. Aguado and F. Calvo. This work has been partially supported by the Deutsche Forschungsgemeinschaft.


text:
# 1 Introduction

## 1 Introduction

Supersymmetry (SUSY) provides a method of solving the “naturalness” or “hierarchy” problem by introducing a set of new particles which cancel the large radiative corrections to the Higgs mass. The cancellation is achieved by assuming that, for each Standard Model (SM) particle chirality state, there is one additional particle, identical to its SM partner except that its spin differs by 1/2 unit. If SUSY were an exact symmetry, the new SUSY particles would have the same masses as their SM partners. Since this scenario is experimentally excluded, SUSY must be a broken symmetry. It is typically assumed that SUSY is broken in some “hidden” sector of new particles, and is “communicated” (or mediated) to the “visible” sector of SM and SUSY particles by one of the known interactions. Two scenarios for this mediation have been widely investigated: gravity and gauge mediation. In gauge mediated SUSY breaking (GMSB), the hidden111In some GMSB papers this is called the “secluded” sector, to avoid confusion with the hidden sector of gravity mediated SUSY breaking models. sector can lie at energies as low as about $`10^4`$ GeV. In most current GMSB theoretical work , it is assumed that this hidden sector is coupled to a messenger sector, which in turn couples to the visible sector through normal SM gauge interactions. The advantage of GMSB over gravity mediated models is that flavour changing neutral currents cannot be induced by SUSY breaking because the normal gauge interactions are flavour blind. A feature which distinguishes gravity from gauge mediated models is the mass of the gravitino, $`\stackrel{~}{\mathrm{G}}`$. In gravity mediated models, $`\stackrel{~}{\mathrm{G}}`$ is usually too heavy to have a significant effect on SUSY phenomenology, while in GMSB models, the $`\stackrel{~}{\mathrm{G}}`$ is typically light ($`<`$ 1 GeV) and is the lightest SUSY particle, the LSP. While $`\stackrel{~}{\mathrm{G}}`$ is a spin 3/2 particle, only its spin 1/2 component (which has “absorbed” the goldstino associated with spontaneous SUSY breaking via the “superhiggs” mechanism) interacts with weak, rather than gravitational, strength interactions, and contributes to phenomenology.

The next-to-lightest SUSY particle (NLSP) is usually either the lightest neutralino ($`\stackrel{~}{\chi }_1^0`$) or the lightest scalar lepton ($`\stackrel{~}{\mathrm{}}_1^\pm `$), and in a significant fraction of the parameter space the NLSP is the lightest scalar tau lepton ($`\stackrel{~}{\tau }_1`$). The coupling of the SUSY particles to $`\stackrel{~}{\mathrm{G}}`$ is small, and typically SUSY particles will decay to the NLSP, which then decays to the gravitino via $`\stackrel{~}{\chi }_1^0\gamma \stackrel{~}{\mathrm{G}}`$ or $`\stackrel{~}{\mathrm{}}\mathrm{}\stackrel{~}{\mathrm{G}}`$. If the decay to the gravitino occurs with a small lifetime, the distinguishing feature is events with energetic leptons or photons, plus significant missing energy due to the missing gravitinos. OPAL has considered scalar lepton and lightest neutralino pair creation in these scenarios in previous publications . This paper reports the first OPAL results on the systematic search for experimental topologies expected in SUSY models with a light gravitino, assuming prompt decays of the NLSP into the $`\stackrel{~}{\mathrm{G}}`$. It is also possible for the NLSP lifetime to be significant and it may decay near the interaction point, at an observably macroscopic distance, or outside the detector. Many signatures expected with long NLSP lifetimes have been considered in other OPAL publications , and the case of arbitrary lifetime will be considered in a subsequent publication. In addition to the new experimental searches which have sensitivity to general SUSY models with light gravitinos, the results reported in this paper are also combined with those from previous publications to constrain minimal GMSB models. Results from searches for GMSB have also been reported by other collaborations .

## 2 OPAL Detector and Event Simulation

The OPAL detector is described in detail in Reference . The SUSYGEN  event generator was used to simulate most of the signal events. For $`\stackrel{~}{\chi }_1^0`$ pair creation in the $`\stackrel{~}{\chi }_1^0`$ NLSP case, 1000 events were simulated for each $`\stackrel{~}{\chi }_1^0`$ mass, using 6 $`m_{\stackrel{~}{\chi }_1^0}`$ points from 50 GeV to 94 GeV. For $`\stackrel{~}{\mathrm{}}`$ pair production in the $`\stackrel{~}{\chi }_1^0`$ NLSP case, 1000 events were generated for each ($`m_\stackrel{~}{\mathrm{}}`$,$`m_{\stackrel{~}{\chi }_1^0}`$), using a grid of 48 different points for each $`\stackrel{~}{\mathrm{}}`$ flavour. Similar one- or two-dimensional mass grids were used for $`\stackrel{~}{\mathrm{}}`$ and $`\stackrel{~}{\chi }_1^0`$ pair production in the $`\stackrel{~}{\mathrm{}}`$ (and $`\stackrel{~}{\tau }`$) NLSP cases. For chargino pair production, $`\stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^{}`$, and the associated pair production of the lightest and second lightest neutralino, $`\stackrel{~}{\chi }_2^0\stackrel{~}{\chi }_1^0`$, in the $`\stackrel{~}{\chi }_1^0`$ NLSP case, the W and Z boson widths can play an important role and are not fully treated in SUSYGEN. The DFGT generator  is used to simulate these signal events. It includes spin correlations and allows for a proper treatment of the W boson and the Z boson width effects in the chargino and heavy neutralino decays. Both SUSYGEN and DFGT include initial-state radiation. The JETSET 7.4 package  is used for the hadronization. The gravitino mass is set identically to zero in the generation, since a small mass in the range favoured by the models has a negligible effect on the detection efficiencies.

The sources of background include two-photon, lepton-pair, multihadronic, and four-fermion processes. The Monte Carlo generators PHOJET  (for $`Q^2<4.5`$ GeV2) and HERWIG  (for $`Q^24.5`$ GeV2) are used to simulate hadronic events from two-photon processes. The Vermaseren  program is used to simulate leptonic two-photon processes ($`\mathrm{e}^+\mathrm{e}^{}\mathrm{e}^+\mathrm{e}^{}`$, $`\mathrm{e}^+\mathrm{e}^{}\mu ^+\mu ^{}`$ and $`\mathrm{e}^+\mathrm{e}^{}\tau ^+\tau ^{}`$). Four-fermion processes were simulated using KORALW , and with the grc4f  generator, both of which take into account all interfering four-fermion diagrams. The improved simulation of the transverse momentum of photons from initial-state radiation makes the use of KORALW essential for events including photons in the final states. Lepton pairs were simulated using the KORALZ  generator for $`\tau ^+\tau ^{}(\gamma )`$, $`\mu ^+\mu ^{}(\gamma )`$ and $`\nu \overline{\nu }\gamma (\gamma )`$ events, and the BHWIDE  (when both the electron and positron scatter at least 12.5 from the beam axis) and TEEGG  (for the remaining phase space) programs for $`\mathrm{e}^+\mathrm{e}^{}\mathrm{e}^+\mathrm{e}^{}(\gamma )`$ events. Multihadronic, $`\mathrm{q}\overline{\mathrm{q}}(\gamma )`$, events were simulated using PYTHIA .

Generated signal and background events were processed through the full simulation of the OPAL detector  and the same event analysis chain was applied to the simulated events as to the data. A data set of approximately 182 pb-1 at a luminosity weighted centre of mass energy of $`\sqrt{s}=188.7`$ GeV was used for the analysis.

## 3 Analysis

For all the selections, after the event reconstruction, double-counting of energy between tracks and calorimeter clusters is corrected by reducing the calorimeter cluster energy by the expected energy deposition from aligned tracks . The selections for events with leptons or hadronic jets plus photons and missing energy described in Section 3.1, as well as four or six leptons plus missing energy in Section 3.2, use the same lepton identification and isolation requirements and preselection as the OPAL Chargino/Neutralino analysis . The most significant preselection cuts require that there is no significant energy in the OPAL forward calorimeters.

### 3.1 $`\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟏}^\mathrm{𝟎}`$ NLSP

#### 3.1.1 $`\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟏}^\mathrm{𝟎}\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟏}^\mathrm{𝟎}`$ Production with $`\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟏}^\mathrm{𝟎}`$ NLSP

The search for lightest neutralino pair production followed by the decays $`\stackrel{~}{\chi }_1^0\gamma \stackrel{~}{\mathrm{G}}`$ uses the OPAL selection of events with photon pairs and missing energy . The analysis selects events with at least two photon candidates and significant missing energy, along with no other significant energy in the event. A total of 24 events are selected, which is consistent with the expectation of 26.9$`\pm `$1.2 events from Standard Model $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }\gamma \gamma (\gamma )`$ production and 0.11$`\pm `$0.04 from all other sources. The selection efficiency for $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }\gamma \gamma (\gamma )`$ within the kinematic acceptance of the analysis is (66.4$`\pm `$2.9)%. One can calculate  the maximum neutralino mass, $`M_{\stackrel{~}{\chi }_1^0}^{\mathrm{max}}`$, which is consistent with the measured three-momenta of the two photons. A cut on $`M_{\stackrel{~}{\chi }_1^0}^{\mathrm{max}}`$ provides further suppression of the $`\nu \overline{\nu }\gamma \gamma (\gamma )`$ background while retaining high efficiency for the signal hypothesis. We require that the maximum kinematically allowed mass be greater than $`m_{\stackrel{~}{\chi }_1^0}5`$ GeV, which retains $`(95.5_{1.0}^{+2.0})`$% relative efficiency for signal at all values of $`m_{\stackrel{~}{\chi }_1^0}`$ while suppressing much of the remaining $`\nu \overline{\nu }\gamma \gamma (\gamma )`$ background. The number of selected events consistent with a given value of $`m_{\stackrel{~}{\chi }_1^0}`$ varies from 14 for $`m_{\stackrel{~}{\chi }_1^0}`$ 45 GeV to 3 events at the kinematic limit. The expected number of SM background events decreases from $`13.67\pm 0.20`$ at $`m_{\stackrel{~}{\chi }_1^0}45`$ GeV to $`1.34\pm 0.07`$ consistent with $`m_{\stackrel{~}{\chi }_1^0}94`$ GeV.

#### 3.1.2 $`\stackrel{\mathbf{~}}{\mathbf{}}^\mathbf{+}\stackrel{\mathbf{~}}{\mathbf{}}^{\mathbf{}}`$, $`\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟏}^\mathbf{+}\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟏}^{\mathbf{}}`$ and $`\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟐}^\mathrm{𝟎}\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟏}^\mathrm{𝟎}`$ Production with $`\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟏}^\mathrm{𝟎}`$ NLSP

With a $`\stackrel{~}{\chi }_1^0`$ NLSP, scalar lepton, charginos and neutralinos may be observed via
$`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\mathrm{}}^+\stackrel{~}{\mathrm{}}^{}(\mathrm{}^+\stackrel{~}{\chi }_1^0)(\mathrm{}^{}\stackrel{~}{\chi }_1^0)(\mathrm{}^+\gamma \stackrel{~}{\mathrm{G}})(\mathrm{}^{}\gamma \stackrel{~}{\mathrm{G}})`$,
$`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^{}(\mathrm{W}^{()+}\stackrel{~}{\chi }_1^0)(\mathrm{W}^{()}\stackrel{~}{\chi }_1^0)(\mathrm{W}^{()+}\gamma \stackrel{~}{\mathrm{G}})(\mathrm{W}^{()}\gamma \stackrel{~}{\mathrm{G}})`$ and
$`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\chi }_2^0\stackrel{~}{\chi }_1^0\mathrm{Z}^{()}\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0(\mathrm{Z}^{()}\gamma \stackrel{~}{\mathrm{G}})(\gamma \stackrel{~}{\mathrm{G}})`$.
In all cases, the signature is events with two photons plus missing energy, plus other activity in the detector. Scalar lepton production will always lead to a low multiplicity final state, and only events with at most six tracks are considered in the analysis. Chargino pair production and neutralino associated pair production may lead to either low or high multiplicity final states depending on the decays of the $`\mathrm{W}^{()\pm }`$ and $`\mathrm{Z}^{()}`$ bosons. For charginos and neutralinos, the analysis is therefore divided into two categories:

* High-multiplicity topologies, with $`\mathrm{N}_{\mathrm{ch}}\mathrm{N}_{\mathrm{conv}}>4`$, where $`\mathrm{N}_{\mathrm{ch}}`$ is the total number of tracks in the event, and $`\mathrm{N}_{\mathrm{conv}}`$ is the number of tracks originating from identified photon conversions,
* Low-multiplicity topologies, with $`\mathrm{N}_{\mathrm{ch}}\mathrm{N}_{\mathrm{conv}}4`$.

The background composition depends on the event kinematics, which are functions of the mass difference between the produced particles and the lightest neutralino, $`\mathrm{\Delta }m=mm_{\stackrel{~}{\chi }_1^0}`$. The analysis is therefore separately optimized for different $`\mathrm{\Delta }m`$ regions, listed in Table 1.

The analyses select events with significant missing energy and two photons. Cuts are applied on the event acoplanarity222After forcing the event into two jets, the acoplanarity angle is 180 minus the opening angle between the jets in the plane transverse to the beam axis. ($`\varphi _{\mathrm{acop}}`$), polar angle of the missing momentum ($`\mathrm{cos}\theta _{\mathrm{miss}}`$), missing transverse momentum scaled by the beam energy ($`p_T^{\mathrm{miss}}/E_{\mathrm{beam}}`$) and visible energy scaled by the centre-of-mass energy ($`E_{\mathrm{vis}}/\sqrt{s}`$). Additionally, unlike other SUSY searches, the LSP gravitino is essentially massless, and the backgrounds can be further reduced while retaining high efficiency by also imposing a minimum requirement on $`E_{\mathrm{vis}}/\sqrt{s}`$. In the HM analyses, there is significant background from $`\mathrm{e}^+\mathrm{e}^{}\mathrm{W}^+\mathrm{W}^{}`$ in some of the kinematic regions. This background is reduced by removing the most energetic photon and forcing the event into two jets, and then cutting on the two-jet mass ($`M_{2\mathrm{j}\mathrm{e}\mathrm{t}}`$). In the scalar lepton search, it is also required that there be at least one identified, isolated lepton in each event. To maintain a general search, if two leptons are found they are not required to be of the same flavour.

Finally, in all channels, it is required that there be at least two energetic photons in each event. Photons are identified by selecting unassociated clusters in the electromagnetic calorimeter, with the following isolation requirements in a cone centred on the cluster (15 half-angle for the highest energy photon and 10 half-angle for the second photon):

* scalar momentum sum $`<`$ 1 GeV;
* additional electromagnetic calorimeter energy sum $`<`$ 5 GeV;
* hadronic calorimeter energy $`<`$ 5 GeV.

The energies of the most ($`E_{\gamma 1}`$) and second most ($`E_{\gamma 2}`$) energetic photon are also used to select signal events.

In the $`\stackrel{~}{\mathrm{}}^+\stackrel{~}{\mathrm{}}^{}`$ selections, $`\varphi _{\mathrm{acop}}>5^{}`$ is required, while the cut is tightened to $`\varphi _{\mathrm{acop}}>10^{}`$ for the $`\stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^{}`$ and $`\stackrel{~}{\chi }_2^0\stackrel{~}{\chi }_1^0`$ selections. In all channels, $`|\mathrm{cos}\theta _{\mathrm{miss}}|<0.95`$ is required. The values of all other analysis cuts are summarized in Table 1, for the different $`\mathrm{\Delta }m`$ regions.

Examples distributions of the event visible energy and the energy of the most energetic photon are shown in Figure 1. The selection results are summarized in Table 2, and the numbers of selected events are consistent with Standard Model sources.

The detection efficiencies for events from $`\stackrel{~}{\mathrm{}}^+\stackrel{~}{\mathrm{}}^{}`$ production assuming the decays $`\stackrel{~}{\mathrm{}}\mathrm{}\stackrel{~}{\chi }_1^0`$ and $`\stackrel{~}{\chi }_1^0\gamma \stackrel{~}{\mathrm{G}}`$ are typically 30–50% for scalar electrons and muons, and 20–40% for scalar tau leptons. Summing the high- and low-multiplicity selections, the detection efficiencies for events from $`\stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^{}`$ production assuming the decays $`\stackrel{~}{\chi }_1^\pm \mathrm{W}^{()\pm }\stackrel{~}{\chi }_1^0`$ and $`\stackrel{~}{\chi }_1^0\gamma \stackrel{~}{\mathrm{G}}`$ range from 20–50%, depending on the masses of the chargino and lightest neutralino. Summing the high- and low-multiplicity selections, the detection efficiencies for events from $`\stackrel{~}{\chi }_2^0\stackrel{~}{\chi }_1^0`$ production assuming the decays $`\stackrel{~}{\chi }_2^0\mathrm{Z}^{()}\stackrel{~}{\chi }_1^0`$ and $`\stackrel{~}{\chi }_1^0\gamma \stackrel{~}{\mathrm{G}}`$ are typically 20–50%, depending on the mass of the two neutralinos.

### 3.2 $`\stackrel{\mathbf{~}}{\mathbf{}}`$ NLSP

#### 3.2.1 $`\stackrel{\mathbf{~}}{\mathbf{}}^\mathbf{+}\stackrel{\mathbf{~}}{\mathbf{}}^{\mathbf{}}`$ Production with $`\stackrel{\mathbf{~}}{\mathbf{}}`$ NLSP

The search for lightest scalar lepton pair production followed by the decays $`\stackrel{~}{\mathrm{}}^\pm \mathrm{}^\pm \stackrel{~}{\mathrm{G}}`$ uses the OPAL selection of events with lepton pairs and missing energy . The analysis selects events with two lepton candidates and significant missing energy, along with no other significant energy in the event. A total of 301 events were observed in the data, consistent with the 303.3$`\pm `$1.9 events expected from all background sources. A likelihood selection using information about the leptons’ energies, charges and polar angles is used to maximize the sensitivity of the analysis to slepton production for each slepton mass. An optimized cut for each scalar lepton mass on the likelihood is applied, as described in Reference . No evidence for anomalous production of lepton pairs with missing energy is observed.

#### 3.2.2 $`\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟏}^\mathrm{𝟎}\stackrel{\mathbf{~}}{𝝌}_\mathrm{𝟏}^\mathrm{𝟎}`$ Production with $`\stackrel{\mathbf{~}}{\mathbf{}}`$ NLSP and $`\stackrel{\mathbf{~}}{\mathbf{}}^\mathbf{+}\stackrel{\mathbf{~}}{\mathbf{}}^{\mathbf{}}`$ with $`\stackrel{\mathbf{~}}{𝝉}`$ NLSP

With a $`\stackrel{~}{\mathrm{}}`$ NLSP, the large neutralino pair production cross-section may make the 4-lepton plus missing energy signature $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0\stackrel{~}{\mathrm{}}\mathrm{}\stackrel{~}{\mathrm{}}^{}\mathrm{}^{}(\mathrm{}^+\mathrm{}^{}\stackrel{~}{\mathrm{G}})(\mathrm{}^+\mathrm{}^{}\stackrel{~}{\mathrm{G}})`$ the GMSB discovery channel. Because the scalar tau lepton may be the lightest slepton, the signature may predominantly include tau leptons. The selection is sensitive to all 4-lepton plus missing energy final states. Additionally, with a $`\stackrel{~}{\tau }`$ NLSP, if the neutralinos are too heavy to be produced and the scalar tau lepton is significantly lighter than the scalar electron and muon, then the 6-lepton plus missing energy final state may contribute via $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\mathrm{}}^+\stackrel{~}{\mathrm{}}^{}(\mathrm{}^+\stackrel{~}{\tau }\tau )(\mathrm{}^{}\stackrel{~}{\tau }\tau )(\mathrm{}^+\tau ^+\tau ^{}\stackrel{~}{\mathrm{G}})(\mathrm{}^{}\tau ^+\tau ^{}\stackrel{~}{\mathrm{G}})`$.

The analyses select low multiplicity events by allowing at most 10 tracks. The events are required to have significant missing energy by applying cuts on $`\varphi _{\mathrm{acop}}`$, $`\mathrm{cos}\theta _{\mathrm{miss}}`$, $`p_T^{\mathrm{miss}}/E_{\mathrm{beam}}`$ and $`E_{\mathrm{vis}}/\sqrt{s}`$, listed in Table 3. In the 4-lepton analysis, we require two identified, isolated leptons in the event and remove Standard Model $`\tau ^+\tau ^{}\gamma `$ events by vetoing events with photons with more than half the beam energy. In the 6-lepton analysis, only one identified, isolated lepton is required. This is because these events typically have only two high energy leptons, both taus, and they often have nearby tracks from the other decay products and are therefore not isolated. If the event is consistent with one lepton plus two hadronic jets an additional veto on $`\mathrm{e}^+\mathrm{e}^{}\mathrm{W}^+\mathrm{W}^{}\mathrm{q}\overline{\mathrm{q}}\mathrm{}\nu `$ is applied: events are removed if the invariant mass of the most energetic lepton and missing momentum is greater than 60 GeV, or if the mass of the two hadronic jets is greater than 60 GeV. Examples distributions of the event transverse momentum and the number of identified, isolated leptons are shown in Figure 2. The selection results are summarized in Table 4, and the numbers of selected events are consistent with Standard Model sources.

The detection efficiencies for events from $`\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0`$ pair production assuming the decays $`\stackrel{~}{\chi }_1^0\stackrel{~}{\mathrm{}}\mathrm{}\mathrm{}\mathrm{}\stackrel{~}{\mathrm{G}}`$ are fairly uniform for different $`\stackrel{~}{\chi }_1^0`$ and $`\stackrel{~}{\mathrm{}}`$ masses, and about 50% if the neutralino decays into all three slepton generations with equal branching ratios, and 35% if it decays via staus with 100% branching ratio. The detection efficiencies for events from selectron and smuon pair production assuming decays into staus are also fairly uniform for different sparticle masses, and are typically about 50%.

## 4 Systematic Errors and Corrections

Systematic errors on the number of expected signal events arise from the following sources: the measurement of the integrated luminosity (0.5%); Monte Carlo statistics for the signal samples (1–2%), and interpolation errors when determining the efficiencies at arbitrary masses (typically 5%); gaugino field content of the $`\stackrel{~}{\chi }^\pm `$ and $`\stackrel{~}{\chi }^0`$ which can lead to different production and decay angular distributions ($`<5`$%); modelling of the cut variables in the Monte Carlo simulations (5-10%). The cut variable modelling error is determined by shifting each cut by an amount estimated by comparing data and Monte Carlo in high statistics samples.

The systematic errors on the expected number of background events are determined from: Monte Carlo statistics in the simulated background events (typically 5%); modelling of the cut variables (from 10-20%, depending on the analysis and kinematic region).

In the analyses in Sections 3.1.2 and 3.2.2, a common veto on energy deposition in the forward detectors was applied in the preselection. The rate of events in which accidental energy depositions in the forward detectors exceeds the veto thresholds used in the preselection is estimated from luminosity weighted random beam crossing events to be 2.9%. Since this effect is not included in the Monte Carlo simulations, the luminosity is reduced accordingly by this factor when deriving limits using the data.

## 5 Results

No significant excesses are observed in any channels, so limits are derived using the search results. Limits are derived using the method from Reference , including the effects of systematic errors on the signal detection efficiencies and background expectation using the method from Reference .

### 5.1 Model Independent Interpretations

For a $`\stackrel{~}{\chi }_1^0`$ NLSP, assuming the prompt decay $`\stackrel{~}{\chi }_1^0\gamma \stackrel{~}{\mathrm{G}}`$, production cross-section limit contours are calculated in the mass plane of the particle produced vs. the mass of the $`\stackrel{~}{\chi }_1^0`$. In this scenario, limits on the production cross-sections for the processes $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\mu }^+\stackrel{~}{\mu }^{}`$, $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\tau }^+\stackrel{~}{\tau }^{}`$, $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\chi }_1^+\stackrel{~}{\chi }_1^{}`$ and $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\chi }_2^0\stackrel{~}{\chi }_1^0`$ are shown in Figure 3. Typically, production cross-sections in excess of about 0.03 – 0.1 pb are excluded at the 95% confidence level.

For a $`\stackrel{~}{\mathrm{}}`$ NLSP, assuming the prompt decay $`\stackrel{~}{\mathrm{}}\mathrm{}\stackrel{~}{\mathrm{G}}`$, production cross-section limit contours can be calculated in the mass plane of the particle produced vs. the mass of the $`\stackrel{~}{\mathrm{}}`$, shown in Figure 4. Typically, cross-sections for $`\mathrm{e}^+\mathrm{e}^{}\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0`$ larger than 0.05 – 0.06 pb are excluded at the 95% confidence level for the degenerate slepton case, while cross-section larger than 0.07 – 0.15 pb are excluded for the stau NLSP case. For scalar electron or muon pair production with a scalar tau NLSP, the cross-section limits are typically 0.06 – 0.13 pb.

### 5.2 GMSB Model Dependent Interpretations

While there is no single GMSB model, there are typically  six new parameters in addition to those of the SM:

$$F,\mathrm{\Lambda },M,N,\mathrm{tan}\beta \mathrm{and}\mathrm{sign}(\mu ).$$
(1)
The intrinsic SUSY breaking scale is $`\sqrt{F}`$, which also determines the $`\stackrel{~}{\mathrm{G}}`$ mass according to $`m_{\stackrel{~}{\mathrm{G}}}2.5\times F/(100\mathrm{TeV})^2`$ eV. Since $`\sqrt{F}`$ affects primarily the lifetime of the NLSP we do not vary it for this paper, but simply assume that this lifetime is short enough to have no effect on our detection efficiencies. The parameter $`\mathrm{\Lambda }`$ sets the overall mass scale for SUSY particles, $`M`$ is the mass of the messenger particles, $`N`$ is the number of sets333 $`N`$ is technically the Dynkin index of the gauge representation of the messenger fields. To preserve gauge coupling unification, the messengers are assumed to form a GUT representation. In the simplest form, each of the $`N`$ messenger particle sets has the quantum numbers of an $`\mathrm{𝟓}\mathbf{+}\overline{\mathrm{𝟓}}`$ of $`𝑺𝑼\mathbf{(}\mathrm{𝟓}\mathbf{)}`$. The maximum number of messengers can be bounded by requiring the gauge interactions remain perturbative up the GUT scale, although this bound depends on $`M`$. For $`M=`$ 100 TeV, $`N5`$, while for $`M=10^{10}`$ TeV, $`N10`$.  of messenger particles, and $`\mathrm{tan}\beta `$ is the usual ratio of the Higgs vacuum expectation values. The final parameter is just the sign of the Higgs sector mixing parameter, $`\mu `$, which introduces a two-fold ambiguity (the magnitude of $`\mu `$ is calculable from the other parameters in the minimal model by imposing radiative electroweak symmetry breaking). The messenger scale gaugino masses can be calculated using the relation

$$M_i=N\frac{\alpha _i}{4\pi }\mathrm{\Lambda }g(\mathrm{\Lambda }/M),$$
(2)
where the index $`i`$ refers to the $`U(1)`$, $`SU(2)`$ or $`SU(3)`$ gauge group, and the $`\alpha _i`$ are the SM gauge couplings. The function $`g(\mathrm{\Lambda }/M)`$ is always slightly greater than 1, but its effect is only significant when $`\mathrm{\Lambda }M`$. It is also possible to calculate all of the messenger scale scalar SUSY particle masses using $`\mathrm{\Lambda }`$ and $`N`$ via

$$m^2=2N\mathrm{\Lambda }^2f(\mathrm{\Lambda }/M)\underset{i=1}{\overset{3}{}}k_i\left(\frac{\alpha _i}{4\pi }\right)^2,$$
(3)
where $`k_i`$ are multiplicative factors determined by the particle’s SM charge, hypercharge and colour charge. The function $`f(\mathrm{\Lambda }/M)`$ is usually near 1, except when $`\mathrm{\Lambda }M`$. The electroweak scale particle masses may be calculated from the messenger scale masses using the renormalization group equations.

We will work in a theoretical framework based on Reference , extending it by including a full mass treatment for all three generations of sparticles. The theoretical calculations are embedded in the SUSYGEN  generator. The complete interpretation framework is described in Reference . The SUSY breaking scale (or equivalently gravitino mass) is not considered explicitly in this section, although it is assumed that $`M_{\stackrel{~}{\mathrm{G}}}<1`$ GeV for the selection efficiencies to remain valid. The values of the parameters considered in our scan are shown in Table 5.

In Figure 5, 95% C.L. exclusion limits in the $`\mathrm{\Lambda }`$ vs. $`\mathrm{tan}\beta `$ plane are shown for different values of $`N`$. The $`\stackrel{~}{\chi }_1^0`$ NLSP signatures tend to be dominant for small $`N`$ and low $`\mathrm{tan}\beta `$, while the $`\stackrel{~}{\mathrm{}}`$ NLSP signatures are more important for either larger $`N`$ or larger $`\mathrm{tan}\beta `$. Absolute 95% C.L. lower limits on $`\mathrm{\Lambda }`$ of 48, 31, 22 and 19 TeV are established for $`N=`$ 1, 2, 3 and 4, respectively. In Figure 6(a), the excluded region in the $`M_\stackrel{~}{\mathrm{}}`$ vs. $`M_{\stackrel{~}{\chi }_1^0}`$ plane is shown for $`\mathrm{tan}\beta =2`$, corresponding to the case where the three sleptons are degenerate in mass. The dominant exclusion channels in the minimal model are $`\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0\gamma \stackrel{~}{\mathrm{G}}\gamma \stackrel{~}{\mathrm{G}}`$ and $`\stackrel{~}{\mathrm{}}^+\stackrel{~}{\mathrm{}}^{}\mathrm{}^+\stackrel{~}{\mathrm{G}}\mathrm{}^{}\stackrel{~}{\mathrm{G}}`$. In Figure 6(b), the excluded region in the $`M_{\stackrel{~}{\tau }}`$ vs. $`M_{\stackrel{~}{\chi }_1^0}`$ plane is shown for $`\mathrm{tan}\beta =20`$, corresponding to the case where the $`\stackrel{~}{\tau }`$ is significantly lighter than the other sleptons. The $`\stackrel{~}{\chi }_1^0\stackrel{~}{\chi }_1^0\gamma \stackrel{~}{\mathrm{G}}\gamma \stackrel{~}{\mathrm{G}}`$ channel remains powerful in this case, but the light $`\stackrel{~}{\tau }`$ means that only the $`\stackrel{~}{\tau }^+\stackrel{~}{\tau }^{}\tau ^+\stackrel{~}{\mathrm{G}}\tau ^{}\stackrel{~}{\mathrm{G}}`$ channel contributes to the significance from lepton pair final states.

Finally, 95% C.L. limits can be derived on the NLSP mass of $`M_\stackrel{~}{\mathrm{}}>`$ 83 GeV and $`M_{\stackrel{~}{\chi }_1^0}>`$ 85 GeV for $`\mathrm{tan}\beta =2`$, and $`M_{\stackrel{~}{\tau }}>`$ 69 GeV, $`M_{\stackrel{~}{\mathrm{e}},\stackrel{~}{\mu }}>`$ 88 GeV and $`M_{\stackrel{~}{\chi }_1^0}>`$ 76 GeV for $`\mathrm{tan}\beta =20`$.

## 6 Conclusion

We have searched for signatures expected in models with gauge mediated SUSY breaking at a centre-of-mass energy of $`\sqrt{s}=`$ 189 GeV with the OPAL detector at LEP. No evidence in any search channel over the expectations from the Standard Model was observed. Limits are placed on the production cross-sections for a number of processes for the prompt decays of the next-to-lightest SUSY particle to a gravitino. The results are used to constrain minimal models of gauge mediated supersymmetry breaking.

## 7 Acknowledgements

We particularly wish to thank the SL Division for the efficient operation of the LEP accelerator at all energies and for their continuing close cooperation with our experimental group. We thank our colleagues from CEA, DAPNIA/SPP, CE-Saclay for their efforts over the years on the time-of-flight and trigger systems which we continue to use. In addition to the support staff at our own institutions we are pleased to acknowledge the
Department of Energy, USA,
National Science Foundation, USA,
Particle Physics and Astronomy Research Council, UK,
Natural Sciences and Engineering Research Council, Canada,
Israel Science Foundation, administered by the Israel Academy of Science and Humanities,
Minerva Gesellschaft,
Benoziyo Center for High Energy Physics,
Japanese Ministry of Education, Science and Culture (the Monbusho) and a grant under the Monbusho International Science Research Program,
Japanese Society for the Promotion of Science (JSPS),
German Israeli Bi-national Science Foundation (GIF),
Bundesministerium für Bildung und Forschung, Germany,
National Research Council of Canada,
Research Corporation, USA,
Hungarian Foundation for Scientific Research, OTKA T-029328, T023793 and OTKA F-023259.