text:
# Low-energy excess of vibrational states in v-SiO2: the role of transverse dynamics

## I Introduction

In this paper we address one of the still open questions in physics of glasses: the origin of the excess of vibrational states in the low energy region with respect to the usual Debye density of states. Buchenau:86 ; Fontana:02  This excess, characteristic of every glassy system, gives rise to the presence, in the region below a few meV, of a more or less pronounced peak, observed both in Inelastic X-ray-Scattering (IXS), Inelastic Neutron-Scattering (INS), and Raman spectra. This peak is universally referred to as Boson Peak (BP). Despite its universality, different interpretations and hypotheses have up to now been proposed both on its origin and on the nature of the excess vibrational states. On the basis of experimental observations and computer simulations, these modes have been catalogued as spatially localized, Buchenau:86 ; Foret:96  spatially delocalized and propagating, Benassi:96 ; DellAnna:98 ; Pilla:00  spatially delocalized but diffusive in character. Feldman:99  Theoretical investigations, on the other hand, have interpreted the BP in glasses as arising from the lowest-energy van Hove singularity of the corresponding crystal; Taraskin:01 ; Taraskin:02  or as the precursor of a dynamical instability expected in a disordered structure; Grigera:01  or as arising from the high frequency dynamics of a model glass within a mode-coupling-like description of Götze Gotze:00  and of Theenhaus. Theenhaus:01

Particularly in the case of v-SiO2 at normal density the most widely accepted explanation is that the BP originates from the piling up of modes near the first van-Hove singularity of the transverse acoustic vibrational branch. Taraskin:01 ; Taraskin:02  These modes should mainly involve relative rotations of almost rigid SiO4 tetrahedra as early pointed out by Buchenau and coworkers. Buchenau:86  In addition to the ambient pressure data, a much useful information is available in the case of v-SiO2 for higher densities, obtained both from in-situ measurements and from permanently densified samples. In general, densification results in a shift of the BP towards higher energies and in a simultaneous decrease of its intensity. Sugai:66 ; Zha:94 ; Inamura:99 ; Inamura:00  In addition to the experimental data, the BP behavior with density has been successfully reproduced by a series of simulations. Lacks:00 ; Jund:00 ; Pilla:03a  Recently, very accurate INS measurements on densified vitreous silica (d-SiO2) have confirmed these effects. Foret:02 ; Courtens:02 ; Ruffle:03  In densified samples the BP lies at higher energies than in vitreous silica at normal density and has a rather weaker intensity. In densified samples the inelastic signal in the BP energy range shows a more marked dependence on the scattering wavevector $`Q`$ as regard its intensity and peak position. We present here simulation results on v-SiO2 at different densities which, when compared with the existing experimental data, can help in clarifying the origin and the nature of the excess modes as well as the intensity and shift effects on the BP as a function of density.

## II Simulation

The simulated system was a glass consisting of 680 SiO2 molecules interacting through the modified BKS two-body interaction potential, vanBeest:90 ; Guissani:96  with periodic boundary conditions. After a long thermalization at a density of $`2.2`$ g cm-3 and at 5000 K the system was slowly cooled down to 300 K at constant volume. At this temperature, an energy minimization was performed to locate the equilibrium configuration at room pressure. The box size was then stepwise scaled by 1.5% amounts up to the final density of $`4.0`$ g cm-3. Once reached this value, the density was decreased in the same way. At each step the energies of the individual realizations were minimized using geometrical methods to obtain the relative equilibrium configurations. This procedure mimics a real hydrostatic pressure experiment at room temperature as described in detail elsewhere. Pilla:03a  The pressure was estimated by using the virial theorem and compared to experimental data. Zha:94 ; Pilla:03b  After the density cycle the zero pressure realization had a density of about $`2.8`$ g cm-3. The dynamical properties of the two zero pressure configurations, and of the higher density ones, were computed in the harmonic approximation by diagonalizing the dynamical matrix. In particular, we have computed the density of states $`g(E)`$ and the longitudinal and transverse currents spectra ($`C(Q,E)`$) which, in the one-excitation approximation, are given by:

$`C_{_{\alpha \beta }}^L(Q,E)=`$ $`{\displaystyle \frac{K_BT}{\sqrt{M_\alpha M_\beta }}}\mathrm{\Sigma }_p|\mathrm{\Sigma }_n(\widehat{Q}\overline{e}_p(n))\mathrm{exp}(i\overline{Q}\overline{R}_n)|^2`$
$`\delta (EE_p)`$
$`C_{_{\alpha \beta }}^T(Q,E)=`$ $`{\displaystyle \frac{K_BT}{\sqrt{M_\alpha M_\beta }}}\mathrm{\Sigma }_p|\mathrm{\Sigma }_n(\widehat{Q}\times \overline{e}_p(n))\mathrm{exp}(i\overline{Q}\overline{R}_n)|^2`$ (1)
$`\delta (EE_p)`$
where $`\alpha ,\beta `$ indicate Si and O, and $`\widehat{Q}`$=$`\overline{Q}/|Q|`$.

The current spectra in transverse polarization are, contrary to the longitudinal currents, not experimentally measurable, and thus they are not directly observable by means of inelastic scattering experiments at least in the low-$`Q`$ and low-$`E`$ ranges. They are nonetheless very useful in understanding the details of the vibrational dynamics, as will be explained in the following section.

## III Discussion

### III.1 Current Spectra

Usually the dynamical data obtained via simulation are reported showing the dynamical structure factor, $`S(Q,E)`$. This is mainly because this procedure has the advantage of being directly comparable to the experimental spectra obtained, for instance, from standard INS and IXS measurements. On the other hand, the current spectra, defined as $`C(Q,E)=E^2S(Q,E)/Q^2`$, are much more useful when one is interested in characterizing the vibrational spectra. In fact, each current spectrum, taken at a fixed $`E`$ and $`C(Q,E)`$, gives the $`Q`$-characterization of the individual vibration with energy $`E`$, being the current strictly related (in harmonic approximation) to the power spectrum of the Fourier transform of the corresponding eigenvector. The computed current spectra, both for longitudinal and transverse dynamics, are shown in Fig. 1(upper and lower panel, respectively) as a contour map plot for v-SiO2 at $`\rho =2.2,\mathrm{\hspace{0.17em}\hspace{0.17em}2.8},\mathrm{\hspace{0.17em}\hspace{0.17em}4.0}`$ g cm-3. The most evident feature is that in the longitudinal currents a periodicity is present, even if much less evident than in a crystal. The same feature is present in the transverse ones, but with a longer quasi period which extends beyond the limit of the figure.

As discussed elsewhere, Taraskin:98 ; Pilla:03a  this difference is due to the fact that the longitudinal dynamics is affected by the average height of the SiO4 tetrahedron, which is the main structural unit of v-SiO2, while the transverse vibrations are sensible to the average height of SiO4 tetrahedron decorated with oxygen atoms. In the low-$`Q`$ longitudinal spectrum, in addition to the quasi-periodic pattern, two well observable features are present in the two lowest-density samples, which consist of protunding peninsulae, nearly parallel to the $`Q`$ axis, centered at about $`20`$ and $`100`$ meV. These features have been interpreted Taraskin:98  as a trace of the vibrational dynamics of the crystobalite structure (the crystalline counterpart of v-SiO2), which shows in that energy range flat dispersion bands. It should be noted that the plots of Fig. 1 are shown as a function of the modulus of the wavevector, which is a relevant quantity in glasses, thus loosing every information on the $`Q`$-direction. If this convention were used also in representing the dispersion curves of a crystal, different $`Q`$-directions in the Brillouin zone would superimpose giving qualitatively the same effect as in Fig. 1. In the $`Q`$-$`E`$ range between $`5`$ and $`15`$ nm-1 and $`10`$ to $`25`$ meV (where the flat band is observed in the longitudinal currents) a much more evident flattening is present in the transverse currents. This coincidence has suggested the interpretation that the flattened band in the longitudinal currents arises from a spilling of transverse dynamics into the longitudinal one, due to disorder. Leadbetter:69 ; Pilla:03a ; Pilla:03b ; Taraskin:98

The maxima of the longitudinal and of the transverse currents in the constant $`Q`$ cuts, are reported in Fig. 2 for three different densities.

In the starting configuration, $`\rho =2.2`$ g cm-3, open and crossed squares of Fig. 2, the main maxima relative to the longitudinal currents have a rather defined dispersive character up to very high energies (about $`100`$ meV), while the transverse dispersion one flattens at small $`Q`$-values. The transverse maxima remain centered, within the uncertainty, at the same energy where a minor peak is observed also in the longitudinal current spectra, whose intensity is roughly 50 to 80% of the main feature. Pilla:03b  By increasing the density up to $`4.0`$ g cm-3, the local structure of the glass, which at normal density consists of fourfold coordinated Si ions, is almost completely substituted by an octahedral one (typical of stishovite, the crystalline form of silica, stable at high densities). As a consequence, also the dynamical properties change and this can be seen in Fig. 2. Both longitudinal and transverse maxima follow a nearly linear dispersion law up to the highest energies and no flattening is observed.

### III.2 Boson Peak

From our simulations we calculated also the density of states as a function of the density. The result are shown in Fig. 3 for selected densities. The low energy part of the total density of states shifts toward higher frequencies and decreases in intensity with increasing density. As a consequence, the BP strongly decreases in intensity and shifts at high energies Pilla:03a .

The reduction of the intensity of the BP with the density can be ascribed to two concurring effects. The first one is an overall reduction of excess states, induced by the pressure, and the second one is the shift towards high energies of the excess which, being weighed by an $`1/E^2`$ factor, appears as having a lower intensity. Actually, both mechanisms are effective in v-SiO2 as we will discuss later. In Fig. 4 we present the maxima of the transverse currents computed at several densities. In the $`\rho =2.2`$ g cm-3 sample, after a linear increase a flattening of the maxima of the transverse currents is observed for $`q`$-values higher than $`6`$ nm-1, at about 15 meV. This plateau progressively shifts, and eventually disappears with increasing density. In the highest density sample it is no more observed and only a linear dispersion relationship is present. By decreasing the density the plateau is at least in part recovered, as shown in the upper panel of Fig. 4. It is worth noting that in the densified sample at $`P=0`$ GPa ($`\rho =2.8`$ g cm-3) the flattening is present, but it saturates at an energy of about 5 meV higher than in the starting realization, indicating an hysteretic behaviour.

The shift and the disapperance of the plateau help in interpreting the intensity and the energy effects on the BP. The excess of states shifts to higher energies and progressively disappears. The fact that these effects take place in an energy region, i.e. 15 meV, higher than the one where the Boson Peak is observed, about 4 meV, is only an apparent contradiction. As a matter of fact, the presence of the $`1/E^2`$ factor severely affects the spectral shape of the excess states by weighting only their low-energy tail. Indeed, the true excess of states should be observed only in the difference between the computed density of states and its Debye approximation.

The determination of the actual shape of the vibrational excess is not a simple task, since it implies the knowledge not only of the density of states, but also of the relative weight of the crystalline one, or at least its Debye approximation. The density of states can be obtained from inelastic neutron scattering data, and also, to a certain extent, from unpolarized Raman data. The Debye density of states can be independently estimated via Brillouin light scattering, while heat capacity experiments give the proper normalization factors. Buchenau:86 ; Carini:95

The result of this procedure is shown in Fig. 5, where INS-deduced $`g(E)`$ is shown with the relative estimate of the Debye contribution. The excess of states (full line in the figure), obtained by subtracting the Debye from the DOS, is a very broad unstructured band centered at energies higher than $`7`$ meV. By dividing this latter by $`E^2`$ its maximum shifts towards lower energies, less than $`4`$ meV, as is shown in Fig. 5 dashed line. This latter peak corresponds to the Boson Peak observed both in INS and IXS measurements of the dynamical structure factor. Actually, the observed BP lies on the low energy tail of the true excess of states, since the contribution of the remaining high-energy vibrations is almost suppressed by the $`1/E^2`$ factor. A similar effect is observed also in the depolarized Raman spectra. Here one does not observe directly the density of states, but this quantity weighed by the Raman coupling function, Fontana:97  as shown in Fig. 5(b). A qualitatively similar behaviour is found even if the scattered intensity is somewhat deformed by the Raman coupling function $`C(E)`$.

The difference between the position of the maxima in the excess of states as found in simulations (about $`17`$ meV as in Fig. 4) and in experiments (about $`7`$ meV as in Fig. 5), can derive from two main causes. From the computational point of view, the discrepancy may derive from the potential used. Moreover, the experimental density of states is obtained through some approximations, like f.i. the incoherent approximation and the estimate of the multiphonon contribution. Also, in the simulations the excess states could be overestimated due to unphysical quenching times currently used in Molecular Dynamics, as pointed out recently by Angell. Angell:03  Finally, the behavior of the transverse current maxima with density shown in Fig. 4, not only accounts for the intensity decrease of the BP, but also explains the coherence effects of the BP measured with in INS and IXS as a function of the density. The BP in normal v-SiO2 has a small coherence effect, i.e. it is rather insensitive to the exchanged wavevector. Benassi:96  On the contrary, densified samples have a marked coherent behavior at least at small $`Q`$-values. Foret:02  From Fig. 4 full circles, one sees that the coherence effects in the BP range, i.e. variation in position and/or intensity, are expected in the linear dispersing range in the case of the $`\rho =2.2`$ g cm-3 sample, being the dispersion curve flat for wavevectors higher than 6 nm-1. In the densified sample at zero pressure, $`\rho =2.8`$ g cm-3, the linear region extends up to 9 nm-1. In the latter sample it is indeed possible to study, using IXS, a more extended linear region and the direct observation of a Brillouin peak directly in the energy cuts is possible in spite of the resolution limitations. Hence, in the densified samples our simulations predict (i) in the low-$`Q`$, low-$`E`$ range a dispersing region where a Brillouin peak is observed extending to wavevectors higher than than in normal silica measurements, (ii) the appearance of the BP (nearly $`Q`$ insensitive) at wavevectors (and corresponding energies) higher than in v-SiO2, (iii) a BP intensity weaker than that of the normal density sample. All these predictions are confirmed by the experiments of Foret and coworkers, Foret:02  while the energy shift and the intensity effect have also been reported by Inamura and coworkers. Inamura:00

## IV Conclusions

In conclusion, the comparison between experiment and simulation allows us to affirm that the BP in v-SiO2 corresponds to the low-energy tail of a flattened transverse branch, whose center lies at relatively high energies. This is also in agreement with earlier studies, Buchenau:86 ; Taraskin:97  which ascribed these vibrations to hindederd librations of coupled SiO4 tethrahedra. The effect responsible for the occurrence of the BP which we have described here is valid for v-SiO2, and it cannot be excluded that other mechanisms are active for different systems, like the idealized ones studied in Ref. Grigera:01, . Nevertheless, it suggests a general approach for the study of the BP in glasses. It is well known that every kind of disorder induces a sort of level repelling and a consequent piling up of modes at the extremes of the vibrational spectrum, as shown by Elliott and coworkers, (see for instance Ref. Simdyankin:02,  and references therein).

## V Acknowledgments

We are grateful to Giancarlo Ruocco, Uli Buchenau and Stephen Elliott for very helpful discussions. This work was financially supported by INFM by PRAGENFDT, Italian Ministero degli Affari Esteri, and MURST Progetto di Ricerca di Interesse Nazionale.


format: text

text:
# Chapter 1 Results from the BAIKAL Neutrino Telescope

## Chapter 1 Results from the BAIKAL Neutrino Telescope

### Abstract

We review the present status of the Baikal Neutrino Project, present updated results on the search for high energy extraterrestrial neutrinos, fast magnetic monopoles and neutrinos induced by WIMP annihilation in the center of the Earth and compare the recorded atmospheric neutrino flux to predictions.

### 1. Introduction

The Baikal Neutrino Telescope NT-200 is operated in Lake Baikal, Siberia, at a depth of 1.1 km. A description of the detector as well as physics results from data collected in 1996 and 1998 (70 and 234 live days, respectively) have been presented elsewhere \[1-3\]. Here we present new limits including data taken in 1999 (268 live days). Data taken in 2000 are presently being analyzed. We also describe NT-200$`+`$ – an upgrade of NT-200 by three sparsely instrumented distant outer strings which increase the fiducial volume for high energy cascades to a few dozen Mtons. A prototype string of 140 m length with 12 optical modules was deployed in March 2003, and electronics, data acquisition and calibration systems for NT-200$`+`$ have been tested.

### 2. Atmospheric Muon Neutrinos as Calibration Tool

The clearest signature of neutrino induced events is a muon crossing the detector from below. Track reconstruction algorithms as well as background rejection have been described elsewhere . Since on the one side the energy threshold for this particular analysis (15-20 GeV) is higher than in underground detectors, and on the other side NT-200 is much smaller than Amanda, event rates of atmospheric neutrinos are small compared to these experiments. Nevertheless, atmospheric neutrinos serve as an important calibration tool and demonstrate the understanding of the detector performance. The data set yields 84 upward going muons. The MC simulation of upward muon tracks due to atmospheric neutrinos gives 80.5 events. The angular distribution for both experiment and simulation as well as the skyplot of upward muons are shown in Fig. 1.

### 3. Search for Neutrinos from WIMP Annihilation

The search for WIMPs with the Baikal neutrino telescope is based on a possible signal of nearly vertically upward going muons, exceeding the flux of atmospheric neutrinos . With no significant excess observed, we derive improved upper limits on the flux of muons from the direction of the center of Earth related to WIMP annihilation. Note that the threshold of 8-10 GeV for this analysis is lower than that for atmospheric neutrinos spread across the full lower hemisphere (see above). Fig.2 compares our new limits to those obtained by other experiments (see  and references given in ).

### 4. Search for Relativistic Magnetic Monopoles

Events due to relativistic monopoles ($`\beta >`$ 0.75) are distinguished by their high light output, allowing identification of events beyond the geometrical boundaries of the detector. The search strategy has been described in . An improved analysis including data from 1996 to 1999 yields a limit about a factor of two below the limit published earlier. This limit is compared to those from other experiments (see  and references in ) in Fig. 3.

### 5. A Search for Extraterrestrial High Energy Neutrinos

The BAIKAL survey for high energy neutrinos searches for bright cascades produced at the neutrino interaction vertex in a large volume around the neutrino telescope . Lack of significant light scattering allows to monitor a volume exceeding the geometrical volume by an order of magnitude. This results in sensitivities of NT-200 comparable to those of the much larger Amanda-B10 detector. The background to this search are bright bremsstrahlung flashes along downward muons passing far outside the array. The method has been described in .

Candidate events do not show a statistically significant excess of hit multiplicity compared to the simulated background from atmospheric muons. Assuming an $`E^2`$ shape of the neutrino spectrum and a flavor ratio $`\nu _e:\nu _\mu :\nu _\tau =1:1:1`$, the new, preliminary 90% C.L. upper limit is $`\mathrm{\Phi }_{\nu _e}E^2=410^7\text{cm}^2\text{s}^1\text{sr}^1\text{GeV}`$, about twice below previous results . The preliminary limit on $`\stackrel{~}{\nu _e}`$ at the W - resonance energy is: $`\mathrm{\Phi }_{\stackrel{~}{\nu _e}}5.4\times 10^{20}\text{cm}^2\text{s}^1\text{sr}^1\text{GeV}^1.`$ These limits do not yet include the effect of systematic uncertainties. Fig.4 (left) shows the experimental upper limits  as well as theoretical limits obtained by Berezinsky (B), by Waxman and Bahcall (WB), by Mannheim et al.(MPR), and predictions for neutrino fluxes from Stecker and Salamon (SS) and Protheroe (P).

Fig.4 (right) shows NT-200+ with its three additional outer (plus one possible central) strings. It will allow a much better vertex identification and hence a significantly more precise measurement of cascade energy in a volume around NT-200. The sensitivity of NT-200+ to high energy cascades will be four times better than that of NT-200, with a moderate 20 percent increase of optical modules only.

This work was supported by the Russian Ministry of Research,the German Ministry of Education and Research and the Russian Fund of Basic Research ( grants  03-02-31004, 02-02-17031, 02-07-90293 and 01-02-17227), Grant of President of Russia NSh-1828.2003.2 and by the Russian Federal Program “Integration” (project no. E0248).

### 6. References

1. Belolaptikov I. et al. 1997, Astr.Ph. 7, 263; and Astr.Ph. 12 (1999) 75.

2. Balkanov V. et al. 2001, Nucl.Ph.(Proc.Sup.) B91, 438.

3. Balkanov V. et al., astro-ph/0211571.

4. Habig A. et al. 2001, Proc.of the XVII ICRC, 1558.

5. Ambrosio M. et al. 2002, Eur.Ph.J. C25, 511.


text:
# Image resonance in the many-body density of states at a metal surface

## I introduction

For the understanding of several fundamental properties of condensed matter surfaces, the knowledge of the electronic density of states and of the parallel-wavevector-resolved density are essential. The dispersion of both occupied and excited surface states, the influence of the changes of the electronic structure in thin film growth and enhanced magnetism are just a few examples of relevant quantities. Experimentally the density of states can be probed by a variety of techniques such as photoemission spectroscopiesMatzdorf (1998) (angle resolved, integrated, inverse and two-photon), and scanning tunneling microscopy.Schintke et al. (2001) To provide an adequate theoretical description of the experimental observables, it is necessary to employ methods which retain the continuous character of the spectrum. In other words one has to take into account that certain quantities may be defined for any energy in a given energy interval. Furthermore, the formalism must be capable of describing excited-state properties.

Regarding excited states, considerable experimental interest has been devoted to image-potential induced (IPI) statesDose et al. (1984); Straub and Himpsel (1984); Giesen et al. (1985); Hertel et al. (1996) and resonances.Straub and Himpsel (1986); Quiniou et al. (1993); Yang et al. (1993); Petaccia et al. (1999) IPI states are present in systems where a bulk band gap provides a barrier, trapping electrons in the image tail of the surface potential. If no gap is present at the IPI energies, an electron is not reflected completely at the bulk barrier and hybridization with surface truncated bulk states becomes possible. This results in the formation of resonances for some materials. A comprehensive theoretical description of IPI states has been given by Echenique and coworkers,Echenique and Pendry (1978); Echenique et al. (1985); Bausells and Echenique (1986); Chulkov et al. (1998) whereas the situation for IPI resonances is much less satisfactory. In principle, IPI resonances are intrinsically contained in the many-body framework, already at the level of the $`GW`$ approximation,Hedin and Lundqvist (1969) but to our knowledge surfaces have only been investigated in this context using a repeated slab geometry.Eguiluz et al. (1992); Deisz et al. (1993); Heinrichsmeier et al. (1998) Such a simplified treatment cannot capture the continuous spectrum of a real surface, because the spectral function will inevitably be composed of a limited number of sharp, discrete peaks in place of the resonance.

In this Article we present ab initio many-body calculations of the local density of states (LDOS) of a semi-infinite jellium surface and demonstrate the presence of a broad IPI resonance. We calculate the LDOS – in its many-body generalization, the spectral functionHedin and Lundqvist (1969) – decomposed according to the surface parallel wave-vector $`𝐤_{}`$ within the surface $`xy`$-plane as

$$\sigma (𝐫,\omega )=\frac{\mathrm{d}^2𝐤_{}}{(2\pi )^2}A(z,𝐤_{},\omega ),$$
(1)
where $`A`$ is the $`𝐤_{}`$-resolved LDOS or spectral weight function, defined by

$$A(z,𝐤_{},\omega )=\frac{1}{\pi }\mathrm{}G(z,z,𝐤_{},\omega )\text{sgn}(\omega \mu ).$$
(2)
Here $`G`$ is the one-particle Green’s function in the representation indicated, and Hartree atomic units ($`a_0=0.529`$ Å, 1 hartree$`=27.2`$ eV) are used.

$`G`$ is obtained from Dyson’s equation

$$G=G^{\text{DFT}}+G^{\text{DFT}}\left[\mathrm{\Sigma }_{\text{XC}}v_{\text{XC}}\right]G,$$
(3)
where $`v_{\text{XC}}`$ is the exchange and correlation DFT potential and $`\mathrm{\Sigma }_{\text{XC}}`$ is the $`GW`$ electron self-energy.

Equation (3) is solved using a recent method developed to perform $`GW`$ calculations in infinite, non-periodic geometriesFratesi et al. (unpublished) based on the embedding method.Inglesfield (1981) The advantage of this approach is that the semi-infinite substrate, surface and vacuum regions are treated equally without the need for any fitting parameters or a repeated cell geometry.

The main steps of the computation are as follows: (i) The Kohn-Sham equationKohn and Sham (1965) is solved self-consistently within LDA.Perdew and Zunger (1981) In this framework, the embedding method permits an exact treatment of the semi-infinite substrate. $`G^{\text{DFT}}`$, obtained by numerical inversion of the Kohn-Sham Hamiltonian, describes a continuous spectrum of fictitious non-interacting electrons. (ii) From $`G^{\text{DFT}}`$ we compute the polarization $`P`$ in the random phase approximation. (iii) The inverse dielectric response, $`(1vP)^1`$, yields the effective interaction $`W^{\text{RPA}}`$. (iv) The self-energy $`\mathrm{\Sigma }_{\text{XC}}=iG^{\text{DFT}}W^{\text{RPA}}`$ is calculated using a real-space, imaginary-frequency representation and one obtains the self-energy on the real frequency axis by means of analytic continuation. (v) Equation (3) can now be solved to update $`G`$.

In Sec. II we present and discuss the results. Section III is devoted to conclusions.

## II Results and discussion

We examine first the $`z`$-dependence of the spectral weight function $`A(z,𝐤_{},\omega )`$. This quantity is proportional to the probability amplitude for a particular wave-vector $`𝐤_{}`$ and energy $`\omega `$. In Fig. 1 we report $`A(z,𝐤_{},\omega )`$ for $`k_{}`$ equal to zero and $`\omega `$ equal to the chemical potential $`\mu `$.

All results shown in this Article are obtained for a jellium substrate of aluminum density ($`r_s=2.07`$ $`a_0`$). The spectral weight outside the surface is enhanced by the improved description of exchange-correlation effects in $`GW`$, in common with the states of Al(111) studied in Ref. White et al., 1998. By varying the energy $`\omega `$ we have verified that this feature is common to all bound states. In the bulk, the $`GW`$ spectral weight is lower than that calculated by the LDA because some weight is transferred to lower energies through electron-plasmon coupling. The larger amplitude of $`A`$ in the surface layer, together with the absence of decay into the bulk, identifies the states in this part of the surface band structure as forming a surface resonance.

In Fig. 2 we plot $`\mathrm{\Delta }\sigma =\sigma ^{GW}\sigma ^{\text{LDA}}`$, the difference between our many-body LDOS and that in the LDA.

The energy-dependence of $`\mathrm{\Delta }\sigma `$ as $`z`$ moves from the bulk towards the vacuum is indicated by contour levels. We recapitulate that, for a jellium surface, the LDA effective potential $`v_{\text{eff}}(z)`$ approaches the constant limits $`V_\text{B}`$ and $`V_\text{V}`$ for $`z\mathrm{}`$ (bulk) and $`z+\mathrm{}`$ (vacuum) respectively. As a consequence, the bulk and vacuum limits of $`\sigma ^{\text{LDA}}`$ are simply proportional to the well-known expressions $`\sqrt{\omega V_\text{B}}`$ and $`\sqrt{\omega V_\text{V}}`$. In the bulk, the LDA therefore predicts no states below $`V_\text{B}`$. The electron-plasmon coupling, automatically included in $`GW`$, is responsible for moving states down in energy from above to below $`V_\text{B}`$, yielding the sub-band shown in Fig. 3 for energies around 1 hartree below the chemical potential, as already demonstrated for the homogeneous electron gas.Hedin and Lundqvist (1969)

The presence of the surface is still noticeable even some atomic units into the substrate through Friedel oscillations, which are visible for bound energies from the contour levels of Fig. 2 and in the LDOS at a distance from the surface $`z=10`$ $`a_0`$ in Fig. 3.

When $`z`$ approaches the surface, coupling with plasmons considerably reduces, eventually becoming negligible a few atomic units outside the surface. As $`z`$ moves into the vacuum, the electron density decays to zero. Self-energy effects are no longer present and the LDA description becomes exact. We then obtain $`\mathrm{\Delta }\sigma =0`$.

Apart from surface-truncated bulk structures, Fig. 2 also highlights a completely new feature in the form of a teardrop-shaped enhancement of the LDOS, localized in the near-surface region at energies between the chemical potential and the vacuum level. To highlight this feature we plot the energy dependence of $`\mathrm{\Delta }\sigma `$ for $`z=3.5`$ $`a_0`$ in Fig. 4 (note that this corresponds to a cross section plot along a horizontal line through Fig. 2).

The peak in $`\mathrm{\Delta }\sigma `$, located at about $`1/32`$ hartree below the vacuum level (the position of the first Rydberg stateEchenique and Pendry (1978)), clearly identifies an IPI surface resonance.

The presence of IPI resonances, even in the limiting case of a substrate without bulk reflectivity (the bulk reflectivity is associated with scattering by the atomic-coresPapadia et al. (1990)), has been unclear. Simple model potentials in independent-particle approximations produced IPI resonances for jellium substratesLindgren and Walldén (1989) in some cases, but the presence of a clear resonance depended sensitively on the precise details of the chosen model. We remind that IPI states and resonances are a many-body phenomenon and can thus be included only in an ad hoc way in single-particle dynamic approximations. In contrast to these earlier, parameter-dependent results, our method is fully ab initio and includes many-body correlations.

In the case of semi-infinite jellium, the zero bulk reflectivity yields a large linewidth of the resonance, since hybridization with bulk states is not prevented at all. We estimate of the full width at half maximum (FWHM) from $`\mathrm{\Delta }\sigma `$ as $`0.07`$ hartree (about $`2`$ eV). 111Experimentally, the FWHMs of surface resonances at metal surfaces are usually of the order of tenths of eV.Quiniou et al. (1993); Bulović et al. (1994) The difference may be explained by the non-negligible reflectivity of the atomic cores, reducing the hybridization with the underlying bulk continuum of states and thus the FWHM. We remark that in systems without a bulk band gap the hybridization with the continuous bulk band is the major contributor to the linewidth, as the self-energy alone accounts only for about one tenth of the total FWHM.Chulkov et al. (1998)

Owing to its large linewidth, the resonance shown in Fig. 4 encompasses the whole series of Rydberg image states,Quiniou et al. (1993) which in the infinite bulk barrier model extends from $`V_\text{V}1/32`$ hartree to $`V_\text{V}`$. The resonance is clearly visible in the difference plot between the many-body LDOS and the LDOS in the LDA. In the many-body LDOS itself, on the other hand, it resolves as a shoulder feature just below the vacuum energy (Fig. 5) rather than a distinct peak, because the LDA LDOS is already a strongly increasing function of energy in this range.

The quantity reported in Fig. 4 is the difference between two theoretical quantities and is therefore not directly accessible to experiment. However, it allows us to isolate an important contribution to the phenomenon – the many-body interaction – which could improve the description of the electronic structure of real semi-infinite surfaces beyond the usual single-particle approaches such as Ref. Ishida, 2001.

We will now examine the origin of the resonance emerging in Fig. 4. For this reason it is important to look at the spectral weight function $`A`$, whose integral over $`𝐤_{}`$ is the LDOS \[Eq. (1)\]. In Fig. 6 we present the energy dependence of $`A(z,𝐤_{},\omega )`$ for $`𝐤_{}=0`$ at the same position outside the jellium edge as Figs. 4 and 5 ($`z=3.5`$ $`a_0`$).

We first consider the LDA result. For non-interacting particles in a constant potential $`V`$ the $`z`$-resolved spectral weight function is essentially one-dimensional and hence proportional to $`(\omega k_{}^2/2V)^{1/2}`$. Thus in the vacuum a singularity is present at $`V_\text{V}`$ for $`k_{}=0`$. The spectral weight then diminishes with increasing energy. But near the surface weight is transferred to bulk states decaying out into the vacuum: the spectral function has a peak very close to $`v_{\text{eff}}(z)`$, but singularities are no longer present. We emphasize that in a slab geometry Fig. 6 would appear as a collection of delta functions, each delta corresponding to one of the discrete eigenstates.

When we describe the system in the interacting picture ($`GW`$), we still observe bulk states that spill out into the vacuum. But we also have a new class of states, constituting the IPI resonance. They have substantial weight at and outside the surface. As a consequence spectral weight from the LDA states is transferred into these quasiparticle states. Spectral weight hence moves down from higher energies to the energies of the IPI resonance. This produces the shoulder in $`A`$ at $`V_\text{V}`$, and the displacement of the peak energy to lower energies.

If $`k_{}`$ is increased from its zero value, the profile of $`A`$ shown in Fig. 6 is shifted towards higher energies by $`k_{}^2/2`$, without being distorted too much (i.e., the dispersion of the states is nearly parabolic with effective mass equal to $`1`$). When we integrate over $`k_{}`$ to get the LDOS, one notices that for energies roughly below $`V_\text{V}`$ the $`GW`$ spectral weight $`A^{GW}`$ is greater than $`A^{\text{LDA}}`$ for any value of $`k_{}`$, leading to the positive $`\mathrm{\Delta }\sigma `$ shown in Fig. 4. For higher energies, $`A^{GW}A^{\text{LDA}}`$ changes from negative to positive as $`k_{}`$ increases. Negative contributions dominate and the corresponding $`\mathrm{\Delta }\sigma `$ is negative.

## III Conclusions

In summary, we have evaluated from first principles the local density of states of a semi-infinite simple metal surface. Many-body correlations are included via the $`GW`$ self-energy. The substrate is described as truly semi-infinite, thus enabling a calculation of the continuous spectrum necessary to properly account for the hybridization between surface electronic states and bulk states. We demonstrate the presence of an IPI resonance just below the vacuum energy, encompassing the full series of image-states that would be present in a system with a surface bandgap. The origin of the resonance is explained in terms of the spectral weight tranfer to lower energies due to the inclusion of electron-electron correlation.

###### Acknowledgements.

We are grateful to S. Modesti, L. G. Molinari, G. Onida, and M. I. Trioni for useful discussions. This work was supported by the Italian MIUR through Grant No. 2001021128. Patrick Rinke acknowledges the support of EPSRC and the DAAD.


text:
# Fiber-assisted detection with photon number resolution

## Abstract

We report the development of a photon-number resolving detector based on a fiber-optical setup and a pair of standard avalanche photodiodes. The detector is capable of resolving individual photon numbers, and operates on the well-known principle by which a single mode input state is split into a large number (eight) of output modes. We reconstruct the photon statistics of weak coherent input light from experimental data, and show that there is a high probability of inferring the input photon number from a measurement of the number of detection events on a single run.

Many quantum information strategies require the preparation of nonclassical states. For example, the method of linear optical quantum computing proposed by Knill, Laflamme, and Milburn KLM  demands the preparation of single photon states as well as maximally entangled photon multiplets. A number of schemes have been proposed for the preparation of such states, including single photon emitters singleph ; YamamotoSP ; Rempe  and conditionally prepared photon pairs from parametric downconversion Konrad ; Pittman ; Loock . Conditional state preparation requires the ability to distinguish states of different photon number, which is not possible using conventional photodetectors. Photon number resolution is also desirable to enhance the security of quantum cryptographic schemes Hwang ; Christine . In this case it is important to measure the photon statistics of the source at the sending and receiving stations. For implementations that use weak coherent states this means distinguishing between the detection of, say, one or two photons.

According to the quantum theory of photodetection, the signal obtained from an ideal noise-free detector has a discrete form corresponding to the absorption of an integer number of quanta from the incident radiation. In practice, however, the granularity of the output signal is concealed by the noise of the detection mechanism. When very low light levels are detected using devices with single-photon sensitivity such as photomultipliers or Geiger-mode operated avalanche photodiodes (APDs), the electronic signal can be reliably converted into a binary message telling us with high efficiency whether an absorption event has occurred or not. However, the intrinsic noise of the gain mechanism necessary to bring the initial energy of absorbed radiation to the macroscopic level completely masks the information on exactly how many photons have triggered that event.

There are several methods for constructing photon number resolving detectors. Among those demonstrated to date are the segmented photomultiplier yamamoto , the superconducting bolometer NIST , and the superconducting transimpedance amplifier sobolewski . These detectors operate at cryogenic temperatures and have single photon quantum efficiencies ranging from about 20% for the superconducting devices to approximately 70% in the case of the segmented photomultiplier. On the other hand, conventional room temperature APDs have intrinsic quantum efficiencies up to 80%, though they respond only to the presence or absence of radiation. The ease of use of these devices suggest that it is worth exploring ways to develop photon resolving capability.

In this paper, we present an experimental detection scheme that implements a full photon-number resolved measurement of light intensity for optical pulses. Our experimental setup follows the main idea outlined in BanaszekWalmsley . The detection is done by splitting the input pulse into separate parts that are expected to contain no more than one photon and then detecting them with conventional avalanche photodiodes. The setup we describe here is constructed only from standard passive fiber optic components, but nevertheless permits the splitting into arbitrary many time-delayed pulses, while retaining a fixed number of output spatial modes. The design is depicted schematically in Fig. 1. The basic elements are 50/50 couplers and fibers of variable lengths. In the first step the input signal is divided into two and launched into two fibers of unequal lengths, such that the partite pulses are delayed with respect to each other. Their subsequent combination at another 50/50 coupler then yields two pulses in each of the two channels. A second stage in which the relative delay is twice as long as the previous stage provides four pulses per output channel. Further doubling of the number of temporal output modes can be achieved by adding more stages. Thus the photon number of the incident pulse can then be detected using only two APDs, if the time separation is ensured to exceed the dead times of the APDs. A less efficient scheme using a single coupler and a fiber loop has also been investigated czech .

For our implementation we used single mode fibers at 780nm (Lucent SMC-AO780B) to construct a detector with two stages, i.e. eight temporal modes. After the fiber setup the pulses were separated in time by 108 ns to 164 ns. Without insertion losses the transmission through the complete fiber system was measured to be higher than 56%. A standard laser diode (Thorlabs V3-780-TO-DA) driven in pulsed operation at 777 nm served as a source of coherent states. These pulses had a width of approximately 14 ns and a repetition rate of 10 kHz. The detection of the attenuated signals at the fiber outputs was carried out with two standard APDs (Perkin Elmer SPCM-AQR-13-FC), which are specified to have an efficiency of 66% and typical deadtimes of 50 ns. Hence the time bins with non-zero photons could be identified using a digital oscilloscope.

A careful analysis of the measured data is crucial since the statistics can emulate “non-classical” properties such as sub-Poissonian photon number distributions. This is mainly due to the fact that high photon numbers result in a non-negligible probability that two photons remain together in one pulse and are counted as one. Therefore the number of stages of the fiber configuration essentially limits the incident photon numbers that can be reliably distinguished by the detector. The probability $`p(k)`$ of the detected counts is linked to the signal photon-number distribution $`\rho (n)`$ by

$$p(k)=\underset{n}{}p_{kn}(k|n)\rho (n).$$
(1)
The photon number distribution can be reconstructed, if the conditional probabilities $`p_{kn}(k|n)`$ of obtaining $`k`$ counts for $`n`$ incident photons are known. The exact values of $`p_{kn}(k|n)`$ depend on the detailed experimental setup including non-ideal 50/50 splitting and unbalanced losses. For known parameters of the fiber system they can be calculated by a basic stochastic model that takes into account the different possibilities for spreading the incident photons into the output modes. The inversion of the associated matrix allows the identification of the incident photon number distribution $`\rho (n)`$ from the measured distribution $`p(k)`$. For our detector with eight output pulses we only expected a good performance for distributions with mean photon numbers smaller than three. Hence we could restrict our data processing to $`p_{kn}(k|n)`$ with $`n`$ smaller than eight, which gives an estimated error being smaller than one percent.

To simplify the inversion the influence of the convolution and the losses can be considered separately. This is equivalent to assuming that the probability of sending an input photon to all the output channels sums to one and introducing an overall loss factor. Our reconstructions of the photon number from the experimental data do not correct for losses, since these leave the form of the Poissonian distributions unchanged and reduces the means only. In addition, by not including the losses we demonstrate the potential capability of the detector to distinguish between different photon numbers if the transmission through the fibers is optimized and highly efficient detectors are available. Fig. 2 shows our experimental results for coherent states with mean photon number of 0.79 and 2.00 respectively. The insets depict the actual count statistics we directly obtained from the two APDs for samples of $`10^4`$ single measurements. To eliminate dark counts and afterpulsing we applied a temporal gating such that only counts within the expected time windows of $`45`$ ns were accepted. The bars in the main graphs represent the normalized photon number distributions we acquired from the detected data using the inversion of $`p_{kn}(k|n)`$. We determined the given error bars by running 1000 Monte Carlo simulations for our measurement samples. The dots indicate the theoretical fits for the respective distributions with the same average photon numbers as the experimental data. The analysis reveals excellent agreement between experiment and theory for distributions with average photon numbers smaller or equal to one (see Fig. 2a). Moreover, as seen in Fig. 2b, the direct inversion still works reasonably well for experimental data corresponding to a coherent state with mean photon number of $`2.00`$, though in this case negative probabilities may appear. Note, however, that we did not optimize the inversion in any way nor did we put any constraints that would exclude negative probabilities.

We tested our inversion for increasing mean photon numbers with simulated Monte Carlo data for input Poissonian statistics, which we truncated at a photon number of eight. In this way we were able to confirm that direct inversion is only possible up to a mean photon number two. The error in these reconstructions, based on the mean square differences between the distributions, is $`4.17\times 10^4`$. The comparison of the experimental data with simulated data revealed that the inversion is very sensitive to events that measure counts with $`k>6`$. These events cause negative probabilities to emerge in the photon distribution due to the inversion. Using advanced estimation techniques that exploit known probability properties should allow accurate reconstruction of distributions with higher mean photon number KonradEstimation . Following that line we investigated coherent states with averages up to four. If we assume Poissonian distributions, we can estimate the mean of the detected data by a simple maximum likelihood estimation, taking into account the conditional matrix $`p_{kn}(k|n)`$. Using this approach we found good consistency between experimental data and the corresponding theoretical distribution for a mean photon number of $`3.78`$, where the deviations for all possible photon numbers summed up to be smaller than $`0.012`$. Thus the photon resolving detector proves to be a valuable tool to explore quantum states in the photon number basis.

Generally speaking there exist two main tasks where photon number resolving detectors are needed. So, far we have discussed one case when the initial photon distribution is either unknown or—like in quantum cryptography—is to be confirmed by ensemble measurements. In the second case the incident state is well known from the beginning, and one wants to perform the detection of the photon number on a single shot basis to address states individually. In such applications, e.g. conditional state preparation Konrad , the performance of the detector has to be evaluated on the confidence that $`l`$ counts have actually been triggered by $`l`$ photons. In our experiment this confidence can be characterized by the conditional probabilities $`\stackrel{~}{p}_{nk}(n=l|k=l)`$, which describes the probability that $`l`$ photons cause $`l`$ detection events. The $`\stackrel{~}{p}_{nk}(n=l|k=l)`$ always depend on the losses of the complete system as well as on the photon statistics, i.e. on the input state itself. To illustrate the reliability of our detector in such settings we have calculated the relevant conditional probabilities for input Poissonian photon statistics and different losses (see Table 1). We find that a lossless detector with eight temporal modes allows appropriate discriminations of low photon numbers for coherent states with mean photon numbers up to $`\overline{n}=1.5`$. The distribution of the input light into a limited number of modes can therefore be tolerated if the ratio of maximum photon number to outgoing modes is sufficiently small. As expected losses decrease the confidence for measurements on single quantum systems, which restricts the usage of the detector to photon number distribution with low enough mean values.

In summary we have demonstrated a fiber based photon detector, which is capable of resolving multiple photons. In our particular setup with eight temporal modes we could measure and recover the photon number statistics by direct inversion for coherent states with mean photon numbers up to two. The presented design for the photon-counting detector can easily be extended to further split the incident pulse in time without the need to increase the number of spatial modes, which are monitored by the two APDs. Adding on stages allows both an increase in the reliability of the detector and permits the resolution of higher photon numbers. Because of the comparative simplicity of the setup it is feasible to establish this photon-counting detector as a standard device in quantum optics. In the context of quantum information processing this is a step forward in implementing more sophisticated schemes based on conditional state preparation.

We are grateful for conversations with Czesław Radzewicz. This research was supported by the US Department of Defense through the Army Research Office via grant number DAAD19-02-1-0163. C. Śliwa’s current address is: Centre for Theoretical Physics, Al. Lotnikow 32/46, 02-668 Warszawa, Poland. After completion of this experiment, we became aware of very similar work performed independently Fitch .


text:
# A comparison of statistical hadronization models

## 1 Statistical hadronization models

There is increasing interest to apply the statistical model  in a description of both particle yields  and spectra  produced in heavy ion collisions in terms of a statistically hadronizing transversely expanding fireball. The fitted parameters, and in particular the temperature, however, have varied considerably, ranging from as low as 110 MeV  to 140 MeV  to as high as 160 and 170 MeV . This report examines these differences in terms of the different possible approaches to statistical hadronization, and studies data sensitivity to model choice.

The most commonly used prescription for statistical emission of hadrons is the truncated Cooper-Frye formula .

$$\left(E\frac{dN}{d^3p}\right)_{cf}=p^\mu d^3\mathrm{\Sigma }_\mu f(p^\mu u_\mu ,T,\lambda )\theta (p^\mu d^3\mathrm{\Sigma }_\mu )+\left(E\frac{dN}{d^3p}\right)_{Res}$$
(1)
Where $`p^\mu `$ is the particle’s four-momentum, $`u^\mu `$ is the system’s transverse flow, T is the temperature, $`\lambda `$ is the normalization and/or chemical potential, $`f(E,T,\lambda )`$ is the statistical distribution of the emitted particles in terms of energy and conserved quantum numbers, $`\mathrm{\Sigma }^\mu `$ describes the hadronization geometry and $`\theta (p^\mu d^3\mathrm{\Sigma }_\mu )`$ is a step function, which corrects for inward emission .

The Cooper-Frye formula is appealing because it’s the most general covariant way of describing statistical particle emission. However, its form is compatible with a variety of freeze-out scenarios.

### 1.1 Resonances and normalization

As Fig. 1 shows, the resonance contribution $`(E\frac{dN}{d^3p})_{Res}`$ is a non-negligible and strongly temperature dependent admixture to the total multiplicity. Its role within particle spectra, however, is a matter of controversy. If an interacting hadron gas phase re-equilibrates resonance decay products , it might be more realistic to not include resonances when analysing spectra . Direct detection of resonances  through invariant mass reconstruction, however, suggests that at least in some cases their feed-down has to be taken into account. These resonances will have a different distribution than those produced directly, and their inclusion can shift the temperature by as much as 20 MeV .

If an incoherent system is considered, one can factor out the decay matrix elements (the angular distribution of decay products and lifetime), and arrive at the distribution of “daughter” particles ($`N`$ with momentum $`p`$) from “parent particles” ($`N_R`$ with momentum $`p_R`$) by simply integrating over the available phase space

$$\left(E\frac{dN}{d^3p}\right)_{res}=\underset{j=2}{\overset{N}{}}\frac{d^3p_j}{E_j}\delta (p_R^\mu p^\mu \underset{2}{\overset{N}{}}p_j^\mu )\left(E_R\frac{dN_R}{d^3p_R}\right)$$
(2)
In general, this expression can get very complicated, and Montecarlo integration  becomes necessary. For most cases considered here, where there is one feed-down and two or three body decays, eq. 2 can be integrated semi-analitically .

While including resonances in this way is computationally demanding , it can be implemented with no extra degrees of freedom. It is therefore straightforward to check the effect of resonance inclusion on $`\chi ^2`$. A large increase in statistical significance would be a strong indication that resonance decay products do in fact emerge from the system with little rescattering. Otherwise, rescattering has to be taken into account.

The handling of resonance decays is connected with the way particle spectra are normalized. Several ways to normalize spectra are possible, and the fitted temperature could change significantly depending on the chosen method, since temperature becomes correlated to a different set of parameters.

The simplest approach is to treat normalization of each particle as a free parameter . If the post-hadronization interacting hadron gas phase is long enough to allow inelastic interactions to alter the chemical composition of the fireball, this may be the correct approach. In this case, temperature will be correlated to flow but not to normalization.
The arbitrary normalization means there is no way to include feed-down from resonances without introducing many more parameters (the resonance normalizations) into the fit. (If rescattering is significant, however, these parameters would be physically justified).

Alternatively, a common normalization volume, together with the introduction of flavor chemical potentials, can be used . This approach, justified by the success in handling short-lived resonances through thermal models , is consistent with a scenario in which post-hadronization dynamics does not change particle abundances or distributions significantly. It does not require extra degrees of freedom for resonances to be included.

However, feed-down from resonances means the fitted temperature will be strongly correlated with chemical potentials and reaction volume, in addition to transverse flow: Most short-lived resonances have the same quark content as the lighter daughter particles, making their yield relative to the lighter particle dependent on temperature only. Hence, temperature will influence both absolute normalization (correlating it with volume, chemical potentials) and the slope (correlating it with flow).

If chemical potentials are used, the eventuality of chemical non-equilibrium raises the possibility that, additionally to chemical potentials, the flavor phase space occupancy parameter $`\gamma `$ needs to be used for either strange  or light  flavors. If these parameters are used, they will correlate temperature with chemical potential even more strongly, since, just like temperature, $`\gamma `$ affects particles and antiparticles abundance in the same way.

### 1.2 Freeze-out geometry and flow profile

Several choices of the freeze-out geometry $`\mathrm{\Sigma }^\mu `$ are possible. This choice is important both for getting a full understanding of how freeze-out happens  and because it correlates with the shape of particle spectra, and hence on parameters such as transverse flow and temperature.

Measured RHIC rapidity distributions  indicate that around midrapidity Bjorken boost-invariance  holds. This, together with cylindrical symmetry (appropriate for the most central collisions) constrains $`d^3\mathrm{\Sigma }^\mu `$ to be of the form

$`\mathrm{\Sigma }^\mu =[t_f\mathrm{cosh}(y_L),r\mathrm{cos}(\theta ),r\mathrm{sin}(\theta ),t_f\mathrm{sinh}(y_L)]`$ (3)
$`p^\mu d\mathrm{\Sigma }_\mu =(m_T\mathrm{cosh}(yy_L){\displaystyle \frac{t_f}{r}}p_T\mathrm{cos}(\theta \varphi ))rdrdy`$ (4)
where $`y,\varphi `$ are the particle’s rapidity and momentum directions, while $`y_L,\theta `$ are the longitudinal rapidity and emission angle of the freeze-out hypersurface. $`t_f`$ is the freeze-out time in a longitudinally co-moving frame.

It is clear that several freeze-out models can be produced with different choices of $`t_f/r`$, giving quantitatively different results. It is also clear that the choice of $`t_f/r`$ might influence fitted parameters such as temperature, transverse velocity and normalization. Recently published fits use a variety of freeze-out geometries:

*  assume $`t_f/r=0`$, i.e. emission is independent of radius. Emission happens in the same lab momentum time.
*  fits $`t_f/r`$ as a free parameter, with the result being $`1`$.
This indicates explosive freeze-out, driven by negative vacuum pressure
* In  freeze-out happens at the same proper time $`\tau `$ throughout the fireball, ($`d^3\mathrm{\Sigma }^\mu u^\mu `$). Freeze-out might have this form (“Hubble” freeze-out) if dynamical effects are negligible during hadronization .

An additional physical uncertainty, which is found to be correlated to flow and $`t_f/r`$ in fits, is the flow profile. Hydrodynamics requires that each part of the fireball volume will in general have a different density and transverse expansion rate. For this reason, the integral over $`d^3\mathrm{\Sigma }`$ will in general span a range of flows, weighted by density function $`\rho (r)`$.

$$E\frac{dN}{d^3p}=r𝑑r(Ep_T\frac{dt_f}{dr})f(T,y_T(r),\lambda )\rho (r)$$
(5)
Just like the form of $`d^3\mathrm{\Sigma }`$, the flow profile depends on the hadronization conditions. For instance, assuming that freeze-out happens at a critical energy density yields $`y_Tr`$ , while assuming freeze-out happening at approximately the same time in the lab frame leads to a power distribution for the transverse expansion velocity $`v_Tr^\alpha `$
One can approximate the flow profile with one “average” flow and assume

$$E\frac{dN}{d^3p}(Ep_T\frac{dt_f}{dr})f(T,<y_T>,\lambda )$$
(6)
However, this approach might result in a systematic shift of fitted parameters such as flow and $`t_f/r`$

## 2 Sensitivity to model choice

The ambiguities presented in sections 1.1 and 1.2 mean that it is important to study their effect on the statistical model’s fitted parameters. One way to do this is to use a Montecarlo to generate data according to a particular freeze-out model, and to see what happens if the “wrong” model is used to perform the fit. We have written a MonteCarlo program which can be used for this purpose. While a detailed description of the algorithm is outside the scope of this write-up , we shall provide a short summary of its features here. An acceptance/rejection algorithm is used to generate particles in a statistical distribution in the volume element’s rest frame. The accepted particles are then Lorentz transformed to the lab frame. (Any flow and density profile, as well as any freeze-out surface can be accomodated). Resonance decays are handled through Eq. 2, using the MAMBO algorithm  to generate points in phase space. Output can be used to generate spectra or fed into a microscopic model such as uRQMD .

The Montecarlo output was used to produce the data points in Fig. 1. It can be seen that a boost-invariant statistical hadronization can explain the global properties of the system such as $`dN_{tot}/d\eta `$. It can also be verified that the role of resonances is absolutely crucial.

We proceeded to generate three datasets of particles. Each data set had a temperature of 140 MeV, a maximum transverse flow of 0.55, and out of equilibrium chemistry ($`\gamma _q=1.4,\gamma _s/\gamma _q=0.8`$). Generated particles include $`\pi ,K^+,K^{},p,\overline{p},\mathrm{\Lambda },\mathrm{\Xi },\mathrm{\Omega }`$ and their resonances. The three samples differ in their choice of freeze-out geometry (specifically $`t_f/r`$) and flow profile:

$`t_f/r=0`$ and no flow profile, as fitted in

$`t_f/r=0`$ and a quadratic flow profile, as fitted in

$`t_f/r=1`$, the boost-invariant analogue of .

We have fitted the three samples to a variety of models, producing $`\chi ^2/DoF`$ profiles for freeze-out temperature and the $`t_f/r`$ parameter. Fig. 3 shows the profiles resulting in the fit to sample III while samples I and II are shown in Fig. 2. A full chemistry model with resonances (solid black line) seems to be equivalent (as far as the position of the $`\chi ^2`$ minimum and the value of $`\chi ^2/DoF`$) to a fit in which normalization is particle-specific (dashed red line). However, the chemical potential fit has greater statistical significance since considerably more degrees of freedom are required for arbitrary normalization.

If chemical potentials are included resonances become essential since a fit with chemical potentials but no resonances (dotted blue line) loses all statistical significance. Similarly, the physical presence of non-equilibrium ($`\gamma _{q,s}1`$) means chemical potentials have to include the non-equilibrium parameter for the fit to be meaningful (black dot-dashed line). The freeze-out geometry does not seem to impact the temperature minimum that much. However, the correct freeze-out geometry can be picked out by a comparison of fits to different models by choosing the model with the lower $`\chi ^2/DoF`$. Moreover, data sensitivity to temperature is strongly affected by freeze-out geometry: Comparing the temperature profiles for different choices of $`t_f/r`$ (Figs. 2 and 3) it is apparent that the temperature $`\chi ^2`$ minimum is more definite in the $`t_f/r=0`$ case. In the case of explosive freeze-out, the correlation between temperature and other parameters in the fit (notably flow) increases, resulting in a shallow $`\chi ^2`$ increase at larger than minimum temperatures.

Finally it should be noted that flow profile, freeze-out geometry and temperature appear to be strongly correlated. If data sample II is fitted with a distribution with no flow profile (such as Sample I) there is a non-negligible shift in both the fitted temperature and $`t_f/r`$, and a small rise in $`\chi ^2/DoF`$ (Fig. 2, brown dot-dashed line).

## 3 Fit to RHIC data ($`\sqrt{s_{NN}}=130GeV`$)

Finally, we have perfomed a fit to the available RHIC data. The data sample we used is the same as the one used for the Montecarlo, but, since the STAR and PHENIX spectra had different trigger requirements (notably centrality) and acceptance regions, we have used two different system volumes, one for STAR particles and another one for PHENIX.

The results are similar to the Monte Carlo data in many ways. The $`\chi ^2/DoF`$ was only slightly larger.
A fit with particle-specific normalization gives a very similar $`\chi ^2`$ and fitted temperature to a fit including chemistry and resonances. If chemistry is included in particle spectra analysis, resonances and non-equilibrium are essential. The fit to freeze-out geometry weakly points to $`t_f/r=1`$, a picture that is supported by the temperature $`\chi ^2`$ profile, virtually identical to data sample III. However, we can not claim our study to be complete in this respect, since we have not yet investigated the effect of including flow profile in the models. As Montecarlo simulations have shown, the conclusion can differ once these are taken into account.

Fig. 5 shows the hyperon and pion spectra from the global fit of Fig. 4. A comparison of the fits on the left panel confirms that a model with no chemistry is about as good at fitting particle spectra as a model with resonances and chemical potentials. However, the second fit also has predictive power: We have used the fitted parameters to predict the $`m_T`$ spectrum for the $`\mathrm{\Sigma }^{}`$. Unsurprisingly, we found that the $`\mathrm{\Sigma }^{}`$ should have roughly the same slope as the $`\mathrm{\Xi }`$s, but its total multiplicity should be about three times as big. We therefore suggest that a greater sample of spectra, in particular more spectra of heavy resonances taken within the same centrality bin as light particles (to make sure both flow and emission volume are the same for each particle) would help in establishing whether chemical potentials are a good way to normalize hadron spectra or not.

The only spectrum which presents a significant systematic deviation for most models is the $`\pi ^{}`$. As Fig. 5 (right panel) shows, most models fail to catch the upward dip of the low momentum pions, and indeed simple transverse expansion predicts the reverse trend . Including resonances, and allowing for $`\gamma _q>1`$ helps (the latter is equivalent to postulating a pion “chemical potential” ) However, to fully account for the lowest momentum pions, even addition of resonances and $`\gamma _q>1`$ are not enough. One has to add a power-law component to the pion spectrum

$$E\frac{dN}{d^3p}=\left(E\frac{dN}{d^3p}\right)_{cf}+\frac{A}{(p_T+p_{T0})^\alpha }$$
(7)
This contribution (roughly $`6\%`$ of the total pion yield in the best fit) also accounts for the highest $`p_T`$ pions. Such a parametrization has been justified  and successfully used  before.

In conclusion, we have presented a comparison of different statistical models currently used to fit spectra in heavy ion collisions. We have described how these different models arise from different freeze-out scenarios of a system hadronizing from a thermalized quark gluon plasma. We have used MonteCarlo simulated data to study the sensitivity to model choice of presently available experimental data, and have evaluated different models ability to fit presently existing RHIC data. While data slightly favors a chemical non-equilibrium explosive freeze-out, there is not enough evidence to make a definitive conclusion about this issue. We hope the forthcoming 200 GeV results will clarify this further.

Supported by a grant from the U.S. Department of Energy, DE-FG03-95ER40937 .
GT acknowledges partial conference support provided by NSF grant PHY-03-11859
We thank Magno Machado (IF-UFRGS) for helpful revision suggestions.

## References


text:
# Field-doping of C60 crystals: A view from theory

## 1 Effect of a strong electric field

Reaching substantial charging (of the order of $`n`$ electrons per C60 molecule) in a field-effect device requires enormous electric fields. As the induced charge is basically restricted to one monolayer of C60 wehrli , a rough estimate can be obtained from simple electrostatics: For neutrality the charge on the gate must equal that on the monolayer, hence the field originating from the gate electrode is given by $`E_{\mathrm{gate}}=2\pi n/A_{\mathrm{mol}}`$, where $`A_{\mathrm{mol}}`$ is the area per C60 molecule in the monolayer and $`n`$ is the number of induced electron charges per molecule. Thus the external field is about 1 V/Å per induced charge, corresponding to a voltage drop of about 10 V across the molecule. In such a strong external field the C60 molecules are strongly polarized. Nevertheless, we find that their response is still in the linear regime. Furthermore, in the charged monolayer, the field experienced by a molecule is screened by the polarization of the neighboring molecules. Taking this into account, we find that the field is reduced by about a factor of two. Calculating the splitting of the molecular levels in this screened homogeneous field, we find a quadratic Stark effect, with the splitting of both, the $`t_{1u}`$ and the $`h_u`$ level, becoming of the order of the band-width for a field corresponding to an induced charge of three to four charges per molecule. This seems to be consistent with the typical doping levels that had been reported.

For a more realistic description of the electrostatics in the field-effect device, we have, however, to go beyond considering only a homogeneous field. Given the spherical shape of C60, the natural approach is via a multipole expansion multipole : We choose an external field and the corresponding induced charge per molecule. We determine the multipole expansion of the field generated by all other molecules about the molecule centered at the origin. Using the linear response of a C60 molecule to multipole fields (calculated ab initio), we determine the new charge distribution on the molecules and repeat the procedure until self-consistency is reached.

Figure 1 shows the splitting of the molecular levels in the self-consistent multipole field for different doping levels. While for an external homogeneous field ($`l=1`$ multipole) the splitting is independent of the direction of the field, including the effect of the induced charge on the neighboring molecules breaks this symmetry. Surprisingly, the asymmetry in the splitting is quite strong, even though the fields that break the symmetry are fairly weak compared to the homogeneous field. This is because the multipole potentials with even $`l`$ give rise to a linear Stark effect, which changes sign with the external potential and which gives rise to a strong splitting even for weak fields. In addition it turns out multipole  that the splitting due to the $`l=1`$ and $`l=2`$ potentials add or subtract, depending on the sign of the external field: When inducing electrons they add for the $`t_{1u}`$ level and almost cancel for the $`h_u`$, while when inducing holes the situation is reversed. I.e., when a molecular level is filled, the splitting is substantially enhanced. It reaches the order of the band width when inducing about two electrons, or somewhat more than two holes per molecule. Beyond these fillings the effect of the Stark splitting on the electronic structure of the C60 monolayer that carries the induced charge will clearly be very large, and one can definitely no longer speak of doping.

## 2 HOMO-antiscreening

As we have seen, the charge density of the C60 molecule is strongly polarized in an electric field, and one might expect that the main contribution comes from the polarization of the highest molecular orbital (HOMO). Calculating the change in the HOMO charge density, we find, however, that the dipole moment of the HOMO charge density is of the same order of magnitude as that of the total charge – but of the opposite sign (see Fig. 2). This surprising result can be understood, e.g., in terms of perturbation theory: Expanding the wave function to first order in the external field $`V=E_zz`$ and calculating the dipole moment $`p=ez`$, we find that the leading term is given by a sum over the matrix element, squared, with all unperturbed molecular orbitals of different parity divided by the characteristic energy denominator. Hence the main contribution comes from energetically close-by levels, and the sign of their contribution is determined by whether they are energetically above or below the level under consideration. For the HOMO in a molecule with a large HOMO-LUMO gap this means that the contributions mainly come from the molecular levels below – implying antiscreening. We thus see that HOMO-antiscreening should be quite general for molecules with large HOMO-LUMO gap, and, in fact, it can also be found, e.g., in the series of polyacenes: benzene, naphthalene, anthracene, tetracene, and pentacene.

## 3 Mott transition

Since the bands in the fullerenes are narrow, while the Coulomb repulsion between two electrons on the same molecule is sizable, the doped fullerenes show effects of strong correlation. It is, e.g., only due to orbital degeneracy that A3C60 is metallic and not a Mott insulator c60mott . In field-doped fullerenes the electrons are restricted to a monolayer wehrli . Hence the number of nearest neighbors to which an electron can hop is reduced and the bands are even more narrow. It is therefore expected that the Mott transition occurs at critical values of the Coulomb interaction $`U`$ below those found in the bulk. To determine the transition point, we have performed quantum Monte Carlo calculations filling  for a doped (111)-layer (without Stark splitting) and find that the Mott transition occurs between $`U_c=0.81.2`$ eV for doping with three electrons, and $`U_c=1.21.6`$ eV for the half-filled $`h_u`$ band (see Fig. 3). For integer dopings other than half-filling the transition is expected to occur for even smaller values of $`U`$ filling . Furthermore, the splitting of the molecular levels in the electric field should weaken the effect of the degeneracy on the Mott transition and lead to still smaller values of $`U_c`$ manini . In particular for the $`h_u`$ orbital, the strong electron-phonon coupling might lead to a further reduction of $`U_c`$ JThund . One has, however, to keep in mind that the Coulomb interaction $`U`$ depends on the environment of the molecule. Screening due to the polarization of the neighboring molecules is, e.g., responsible for a large reduction of $`U`$ in the crystal as compared to the value for an isolated molecule c60U . Likewise, it is to be expected that for a molecule in the monolayer next to the gate dielectric, $`U`$ might be substantially changed from the bulk value. This effect is, however, hard to quantify without knowing the microscopic structure of the oxide-C60 interface. Nevertheless it seems likely that field-doped C60 should be insulating at integer fillings.

## 4 Enhancement of transition temperature

In A3C60 the superconducting transition temperature $`T_c`$ increases with increasing lattice constant, i.e., with increasing density of states at the Fermi level rmp . It is therefore natural to try the same for field-doped C60. The simplest way to increase the distance of the molecules in the conducting monolayer is to apply uniaxial stress uniax . An alternative approach is the intercalation of the crystal with inert molecules. In fact, for field-doped C60 intercalated with haloform molecules spectacularly increased transition temperatures have been reported latticeexp ; retract . A subsequent analysis of the lattice structure of these crystals revealed, however, that the lattice is mainly expanded perpendicular to the conducting layer, and that the density of states in the doped layer shows no correlation with the reported $`T_c`$ jansen . Therefore, the additional coupling to the vibrations of the haloform molecules has been proposed as an alternative explanation of the enhancement of the transition temperature bill . It has, however, turned out that such a coupling is very small and, for the two-fold degenerate modes, is even excluded by symmetry CHX3 . The reported enhancement of $`T_c`$ in haloform intercalated C60 is therefore not understood.


text:
# Two–dimensional clusters of liquid 4He

## I Introduction

In recent years, a great deal of work has been devoted to study quantum liquids in restricted geometries.kro02  One important feature of these systems is that their internal structure becomes more easily observable than in bulk liquids due to the restricted motion of the particles in the confining potential. Among these systems the study of quantum films has received particular attention. They consist of liquid helium adsorbed to a more-or-less attractive flat surface. In 1973 M. Bretz et al.bre73  observed for the first time the adsorption of 4He onto the basal plane of graphite. In the last few years, adsorption properties of helium on different substrates such as carbon, alkali and alkaline-earth flat surfaces, carbon nanotubes and aerogels have become a fertile topic of research.

The structure and growth of thin films of 4He adsorbed to a substrate was studied by Clements et al.cle93  employing the optimized hypernetted-chain Euler-Lagrange theory with realistic atom-atom interactions. It turns out that films with low surface coverages, where all atoms cover the surface with a thickness corresponding to a single atom, can be approximated reasonably well by a 2D model. In connection with these systems, an interesting question naturally arises as how physics depends on the dimensionality of the space.

The homogeneous 2D liquid has been studied using different theoretical methods, such as molecular dynamicscam71  and quantum Monte Carlo simulations either Green’s Functionwhi88  or diffusiongio96  techniques. The inhomogeneous case was studied by Krishnamachari and Chester,kri99  who used a shadow variational wave function to describe 2D puddles of liquid 4He. In this work we report energies and density profiles of puddles calculated within the diffusion Monte Carlo (DMC) method. Our main objective is to give an accurate estimate of the line energy or the line tension of the 2D liquid 4He. As atom-atom interaction we have used the revised version of the Aziz potential dubbed as HFD-B(HE).azi87  This potential has been used to study ground–state properties of 3D bulk 4He and 3He,bor94 ; cas00  within the DMC framework, and it has proved to accurately reproduce the ground–state properties of both liquids at zero temperature.

The trial wave function used for the importance sampling in the DMC calculation is introduced in Section II, where the variational Monte Carlo results for this wave function are also reported. A brief explanation of the DMC techniques used in the present paper is presented in Section III. Section IV contains the DMC results and their analysis in terms of a mass formula in 2D. The line tension is extracted from this mass formula. Properties characterizing the puddles, such as the density profiles, are discussed in Section V. Finally, the main conclusions are summarized in Section VI.

## II Variational ground–state energies

To study a system of $`N`$ 4He atoms in two dimensions we start from the following trial wave function

$$\mathrm{\Phi }_T(𝐑)=\underset{i<j}{}\mathrm{exp}\left[\frac{1}{2}\left(\frac{b}{r_{ij}}\right)^\nu \frac{\alpha ^2}{2N}r_{ij}^2\right],$$
(1)
written in the same way as in the 3D case.panda86  The coordinate $`𝐑`$ indicates the set of coordinates of all the particles $`\{𝐫_1,𝐫_2,\mathrm{},𝐫_N\}`$, while $`r_{ij}`$ stands for the interparticle distance, $`r_{ij}=|𝐫_j𝐫_i|`$. The trial wave function contains the simple McMillan formmcm65  to deal with the very short-range part of the interaction, and the translationally invariant part of a harmonic oscillator (HO) wave function with parameter $`\alpha `$, to roughly confine the system.

In our calculations the value $`\mathrm{}^2/m_4=12.1194`$ K Å2 has been employed for the atom mass and the parameters $`b`$ and $`\nu `$ have been fixed to the values 3.00 Å and 5, respectively, the same values as in 3D calculations. The variational search has thus been restricted to the HO parameter $`\alpha `$, whose optimal value is given in Table 1. The expectation value of the Hamiltonian, as well as the separate contributions of kinetic and potential energies are given in the same Table for puddles with $`N`$ atoms. It can be seen that the total energy results from an important cancellation between kinetic and potential energies, which is in fact larger than in the 3D case. Let us recall that in 3D bulk, the energy per particle results from adding $``$ 14 K of kinetic energy with $`21`$ K of potential energy. In 2D, both kinetic and potential contributions are very close to each other, which makes the calculation very delicate.

In the last column of Table 1 the VMC results of Krishnamachari and Chesterkri99  are reported. As compared with their results, our calculations provide smaller binding energies, in spite of the fact that the interaction used in Ref. kri99,  is an older version of the Aziz potential, which tends to underbind the systems. In fact, the shadow wave function used in Ref. kri99,  contains more elaborated correlations not present in our simple trial wave function. The VMC energy for the bulk system corresponds to the saturation density, $`\rho _0=0.04344`$ Å-2, taken from the DMC calculation of Ref. gio96, .

We have also performed calculations using a different trial wave function, replacing the translationally invariant gaussian part by an exponential one, i.e.

$$\mathrm{\Phi }_T(𝐑)=\underset{i<j}{}\mathrm{exp}\left[\frac{1}{2}\left(\frac{b}{r_{ij}}\right)^\nu \frac{\alpha }{2}r_{ij}\right],$$
(2)
expecting that this larger tail in the wave function will result in more binding. Actually, we do not find significant differences for small values of $`N`$. For instance, in the case $`N=8`$, using the same values for $`b`$ and $`\nu `$ as before, we get $`E/N=0.2178(5)`$ K, $`T/N=1.266(2)`$ K, $`V/N=1.484(2)`$ K for $`\alpha =0.035`$ Å-1. When the values of $`b`$, $`\nu `$ and $`\alpha `$ are optimized, we obtain a slightly larger binding energy, $`E/N=0.2267(8)`$ K for $`b=3.04`$ Å, $`\nu =5.0`$ and $`\alpha =0.035`$ Å-1. For greater values of $`N`$, the gaussian ansatz tends to provide more binding than the exponential. For example, with the exponential ansatz, for $`N=16`$ we get $`E/N=0.1816(7)`$ K for $`\alpha =0.023`$ Å-1, and optimizing the different parameters one gets $`E/N=0.2514(6)`$ K, with $`b=3.04`$ Å. In conclusion, the gaussian wave function seems appropriate to be used as importance sampling in the DMC calculations.

## III Diffusion Monte Carlo ground–state energies

Quantum Monte Carlo (QMC) methods provide the exact ground–state energy of a boson system, except for statistical errors. These techniques solve numerically the Schrödinger equation by means of a statistical simulation. They have been widely described in the literature; hence we briefly recall here the main ideas, referring the reader to, for example, Ref. hlr94,  for a more detailed description on QMC techniques. In this work we use the diffusion Monte Carlo (DMC) method to solve the Schrödinger equation in imaginary time ($`\tau =it`$) for the function

$$f(𝐑,\tau )=\mathrm{\Phi }_T(𝐑)\mathrm{\Psi }(𝐑,\tau )$$
(3)
where $`𝐑`$ represents all the particle coordinates and is usually called “the walker”, $`\mathrm{\Psi }(𝐑,\tau )`$ is the wave function of the system, and $`\mathrm{\Phi }_T(𝐑)`$ is the previously determined trial wave function (Sect. II), used here as importance sampling. It is convenient to write the solution of the time-dependent Schrödinger equation in the following form

$$f(𝐑,\tau +\mathrm{\Delta }\tau )=𝑑𝐑^{}G(𝐑,𝐑^{},\mathrm{\Delta }\tau )f(𝐑^{},\tau )$$
(4)
where $`G`$ is the time-dependent Green’s function and is formally written as

$$G(𝐑,𝐑^{};\tau )=𝐑e^{H\tau }𝐑^{}$$
(5)
where $`H`$ is the Hamiltonian of the system. The function $`G(𝐑,𝐑^{};\tau )`$ represents the amplitude probability for the transition from an initial state $`𝐑^{}`$ to a final one $`𝐑`$ in a time $`\tau `$. In the limit $`\tau \mathrm{}`$, Eq. (4) gives the exact ground–state wave function. Thus, knowing $`G`$ for infinitesimal time steps $`\mathrm{\Delta }\tau `$, the asymptotic solution for large times $`f(𝐑,\tau \mathrm{})`$ can be obtained by solving iteratively the above equation. To this end, the exponential entering Eq. (5) is approximated to some fixed order in $`\mathrm{\Delta }\tau `$. Both first and second orderbor94  propagators have been implemented in the present work and both of them provide the same extrapolated energy, within statistical errors, using the trial function $`\mathrm{\Phi }_T`$ of Eq. (1) as guiding function. Our simulations have been carried out with a population of typically 400 walkers. As usual, some runs are first done to establish the asymptotic region of the short time propagator, then several values of the time step have been used, and finally a fit, either linear or quadratic, has been carried out to obtain the extrapolated energy. For example, for $`N=16`$ the time steps 0.0001, 0.0002, 0.0003, and 0.0004 have been used to perform the extrapolation. In general, the statistical error is of the order of the time step error in our calculations.

In Table 2 we present the results of our linear DMC calculations of the total energy per particle for puddles containing $`N`$ atoms. We have also reported and reproduced the results of the binding energy per particle of homogeneous 2D liquid 4He at the equilibrium density $`\rho _0^{DMC}=0.04344(2)`$ Å-2; obtained in Ref. gio96, , where the same version of the Aziz potential was used. For this case, the simulations have been carried out for a system of 64 atoms with periodic boundary conditions, for which the errors due to finite size effects are smaller than the statistical oneswhi88 . We have also performed quadratic DMC calculations for some puddles, and found results which are compatible with the linear DMC ones within their error bars. For example, the quadratic algorithm provides $`E_{quad}=0.2612(2)`$ K for $`N=8`$ and $`0.652(4)`$ for $`N=64`$. As expected, the DMC results lower the corresponding energies obtained by VMC either with our simple variational wave function or with a shadow wave function,kri99  by up to $`25\%`$ in the case of the bulk system. It is worth mentioning that the final DMC result for the energy does not depend on the trial wave function for a boson system like the studied here, a fact that in the present case has been numerically checked for the Gaussian and the exponential ansatzs, Eqs. (1) and (2). Indeed, for boson systems the DMC method provides exact ground–state energies, within statistical errors.

## IV Energy and line tension

For a saturating self-bound system, the ground–state energy per particle can be expanded in a series of powers of the variable $`N^{1/D}`$, where $`N`$ is the number of constituents and $`D`$ is the dimensionality of the space. This is the well-known mass formula, which in the present case writes

$$E(N)/N=\epsilon _b+\epsilon _lx+\epsilon _cx^2+\mathrm{}$$
(6)
with $`x=N^{1/2}`$. The two first coefficients of this expansion are the bulk energy $`\epsilon _b`$ and the line energy $`\epsilon _l`$, out of which the line tension $`\lambda `$ is defined by $`2\pi r_0\lambda =\epsilon _l`$. Here $`r_0`$ is the unit radius, defined as the radius of a disk whose surface is equal to the inverse of the equilibrium density of the 2D bulk liquid, i.e. $`\rho _0\pi r_0^2=1`$. Finally, $`\epsilon _c`$ is the so-called curvature energy.

Our calculated ground–state energies (Tables 1 and 2) are plotted in Fig. 1 as a function of $`N^{1/2}`$. One can see that the differences between our VMC and DMC energies increase with the number of atoms in the puddle. This clearly mirrors the simplicity of the trial wave function, which could be improved by including, for example, three-body correlations. Nevertheless this trial function is adequate for the importance sampling in the DMC calculation.

We have fitted these energies to a parabolic mass formula like Eq. (6). The coefficients of the fit are given in Table 3, together with the deduced line tension. Notice that the coefficient $`\epsilon _b`$ is identical, within statistical errors, to the bulk energy per particle of Table 2. In fact, the $`\chi ^2`$ of the fit is very small, $`\chi ^2=5.7\times 10^6`$. Regarding the line tension, and despite using a different version of the Aziz potential and a different trial function, we notice that our VMC estimate is rather close to the one reported in Ref. kri99, , $`\lambda =0.07`$ K/Å. However, both VMC results are remarkably different from the DMC line tension, $`\lambda =0.121`$ K/Å.

To stress the curvature effect we have also plotted in the figure a straight line between the $`N=8`$ and bulk DMC values. A linear fit of the DMC energies gives coefficients $`\epsilon _b=0.885`$ K and $`\epsilon _l=1.80`$ K, which are appreciably different from the previous ones. The bulk energy extrapolated from this linear fit differs from the directly calculated value, and the corresponding line energy is closer to the variational one, thus giving a bad estimation for the linear tension. In all cases, the line energy coefficient is approximately minus twice the volume energy term, similarly to the 3D case,panda86  and therefore one expects curvature effects to be important. In both VMC and DMC cases the extracted $`\epsilon _c`$ is negative, i.e. the binding energy is a convex function of $`x`$ as it also happens for the 3D clusterspanda86 . This is in contrast with the value of $`\epsilon _c`$ reported in Ref. kri99,  which was positive but rather smaller in absolute value and with larger error bars. Actually, as argued in Ref. panda86,  for 3D clusters, one would expect the curvature correction to the energy of a circular 2D cluster to be positive. Therefore, one should take the extracted value for $`\epsilon _c`$ with certain caution and not to emphasize its physical significance. However, it turns out that the value and sign of $`\epsilon _c`$ are stable against different possible fits, e.g., changing the number of points to build the fit, or using a cubic mass formula. In any case the two first coefficients $`\epsilon _b,\epsilon _l`$ are quite robust against all performed fits. As an illustration, if one takes out the bulk point, the predicted bulk energy per particle and surface tension are equal to the reported ones within 1% and 5% respectively. Therefore, the extracted line tension should be reliable, as it also happens for 3D clusterspanda86 .

## V Density profiles

The calculation of observables given by operators that do not commute with the Hamiltonian poses new problem to the DMC method. After convergence, the walkers are distributed according to the so-called mixed probability distribution given by the product of the exact and the trial wave functions. Therefore averaging the local values of the operator does not give the exact expectation value unless the operator commutes with the propagator. The result obtained by straightforward averaging is the mixed estimator which is of first order error in the trial wave function. Several options have been proposed in the literature in order to obtain unbiased (trial function independent and exact) values. In this work we have used the so called forward or future walking techniquehlr94  to calculate unbiased, also called pure, density profiles. The key ingredient to correct the mixed estimator is to include as a weight in the sampling the quotient $`\mathrm{\Phi }_{\mathrm{exact}}(𝐑)/\mathrm{\Phi }_{\mathrm{trial}}(𝐑)`$ for each walker, given by the asymptotic number of walkers. Several algorithms have been proposed in order to compute this quantity. In this work we use the algorithm developed in Ref. cas95,  that constitutes a simple and efficient implementation of the future walking method.

The pure DMC estimates of the density profiles for several puddles are plotted in Fig. 2. The figure also contains an horizontal line which indicates the saturation density ($`\rho _0^{DMC}`$) of the homogeneuos system. For the puddle containing 36 atoms, the VMC profile obtained from a Gaussian ansatz (Eq. (1)) is also shown for comparison as a dotted line. As one can appreciate in the figure, the process of optimization implied by the DMC method changes the profile reducing its thickness, i.e., producing a sharper surface. It can be seen that for the smaller clusters the central density is below $`\rho _0`$, while for the larger values of $`N`$ shown in the figure the central density is above $`\rho _0`$, indicating a leptodermous behaviour. One expects that, increasing the number of particles, the central density will approach $`\rho _0`$ from above, as in the 3D case.str87 ; chi92  It is worth noticing the oscillating behaviour in the interior part of the density profile for $`N=121`$. It is difficult however to decide whether these oscillations are genuine or are simply due to a poor statistics in evaluating the pure estimator. Unfortunately, to discard this last option would require an exceedingly long computing time, within the scheme of this work.

The solid lines plotted in Fig. 2 are fits to our DMC densities provided by a generalized Fermi profile of the form:

$$\rho (r)=\frac{\rho _f}{\left(1+\mathrm{exp}\left(\frac{rR}{c}\right)\right)^\nu }$$
(7)
The parameters defining the Fermi profile are given in Table 4 together with the thickness $`t`$ and the root mean square (rms) radius $`r^2^{1/2}`$. We have checked that the rms radius calculated within the DMC code and the one derived form the fit agree to better than 0.5%, except for the $`N=121`$ case, where the difference is 1%. The rms radius grows with the number of particles as $`N^{1/2}`$, as expected. Therefore it grows faster than in 3D, in which case grows as $`N^{1/3}`$. This behaviour allows for an alternative determination of the saturation density by performing a linear fit to the relation

$$\sqrt{r^2^{1/2}}=\sqrt{\frac{N}{2\pi \rho _0}}$$
(8)
The extracted value of $`\rho _0`$ from the mean square radius reported in Table 4 is 0.043 Å-2, in good agreement with the determination from the calculation for the homogeneous system.

In the interval of $`N`$ considered, the thickness $`t`$, defined as the distance over which the density falls from 0.9 of its value at origin to 0.1, is continuously increasing. However, as the finite value of the thickness for the semiinfinite system should define the asymptotic value of $`t`$, one expects that for larger puddles the thickness will probably have a maximum and smoothly approach this asymptotic value, as happens in the 3D case.str87  Finally, one also observes the asymmetric character of the density profiles with respect to the point at which the density falls at half its value at the origin. This can be appreciated by looking at the value of $`\nu `$, which grows with $`N`$, and also to the increasing difference between the quantities $`R`$ and $`r^2^{1/2}`$.

## VI Summary and conclusions

In this work we have considered strictly two–dimensional systems of liquid 4He, which are of course an idealization of a real quantum film. They are nevertheless interesting because their study can enlighten the underlying structure of real quasi-2D systems. Of course, in the latter case, one has to take also into account the interaction with the substrate, which basically provides a global attractive potential. In the ideal 2D case, the suppression of the wave function component in the third dimension, produces an increment of the global repulsion between atoms, resulting in a smaller binding energy per particle, and a decrease of the equilibrium density.apa

The binding energies of two–dimensional 4He clusters, calculated by means of a diffusion Monte Carlo method, are well fitted by a mass formula in powers of $`x=N^{1/2}`$. The analysis of the mass formula provides the main result of this paper, namely the value of the line tension $`\lambda =0.121`$ K/Å, which significantly differs from the one obtained from a similar analysis of VMC data and the one previously reported in the literaturekri99 . The quadratic term of the mass formula cannot be neglected and results in a negative value of the curvature energy as in the 3D case.panda86 ; chi92  However, the studied clusters may be too small to give physical significance to this result.

The density profiles obtained with the pure estimator have been fitted by a generalized Fermi function, and the behaviour of the rms radius and the thickness as well as the asymmetry character of the profile as a function of $`N`$ have been discussed. However, calculations for larger puddles, which are out of the scope of the present paper due to limitations in computing time, would be necessary to describe the complete $`N`$-dependence of the density profiles.

###### Acknowledgements.

Fruitful discussions with X. Vi as are gratefully acknowledged. This work has been supported by DGICYT (Spain) project PB98-1247, DGI (Spain) grants BFM2001-0262 and BFM2002-00200, Generalitat de València grant GV01-216, Generalitat de Catalunya project 2001SGR00064, and MURST (Italy) grant MIUR-2001025498. J. M. P. acknowledges a fellowship from the Generalitat de Catalunya.