text:
# Vibrational spectra in glasses

## I Introduction

Propagating density fluctuation of macroscopic size (hydrodynamic limit) are known to exist both in ordered and in disordered materials. Whereas in the ordered ones those excitations (phonons) persist up to momenta of about one tenth of the Debye momentum, the fate of excitations of microscopical size in disordered systems is still a quite puzzling issue, both from the theoretical and the experimental point of view.

Recently, a lot of experimental attention (Buchenau et al. 1986, Benassi et al. 1996, Foret et al. 1996, Masciovecchio et al. 1996, Masciovecchio et al. 1997, Masciovecchio et al. 1998, Monaco et al. 1998, Sette et al. 1998, Fioretto et al. 1999, Ruocco et al. 1999, Sokolov et al. 1999) has been paid to the high frequency dynamics of disordered systems like glasses. As a matter of fact, high-resolution inelastic X-ray scattering (IXS) and neutron scattering techniques have made accessible to experiment the region where the exchanged external momentum $`p`$ is comparable to $`p_0`$, namely the momentum where the static structure factor has its first maximum:

$$0.1p_0<p<p_0$$
(1)
The typical momentum $`p_0`$ can be thought as the generalization of the Debye momentum to the case of disordered systems. With these data available for the study of excitations of microscopic size, one is naturally led to ask which of the well-known features of phonons survive even in disordered systems. A number of facts have emerged from the experiments:

1. The dynamic structure factor $`S(p,\omega )`$ has a Brillouin-like peak for momenta up to $`p/p_00.5`$. This inelastic peak is due to the interaction of the external photon (or neutron) with some excitation of the system. A very controversial issue is the propagating nature of the excitations (see, for example, Caponi et al. and Foret et al. in this same volume) within that range of momenta. As a guideline, one could say that as long as the dispersion relation between the position of the peak (in the frequency domain) and the external momentum $`\omega _P(p)`$ is nearly linear, the excitations are likely propagating. Of course, the study of Anderson localization (see, for instance, Parisi 1994) shows that, strictly speaking, this condition is neither necessary nor sufficient, so care should be taken.
2. The peak has a width $`\mathrm{\Gamma }`$, whose dependence on the momentum $`p`$ within the range (1) has been described by means of the scaling law (for a large variety of materials):

$$\mathrm{\Gamma }p^\alpha $$
(2)
Let us remark that, while there seems to be a quite general agreement about the fact that $`\mathrm{\Gamma }`$ is not affected by changes in the temperature, there is still some debate about the actual value of the exponent $`\alpha `$. More precisely, at low enough momenta (e.g. light scattering) $`\alpha `$ is undoubtedly $`2`$, while at higher momenta (e.g. X-ray scattering) neither the value $`\alpha 2`$, nor the value $`\alpha 4`$, describing instead a Rayleigh scattering regime, can be ruled out on experimental grounds.
3. The density of states (DOS) exhibits an excess respect to the Debye behavior ($`g(\omega )\omega ^2`$), known as the Boson peak. This feature is particularly remarkable for strong glasses. Generally speaking, the typical frequencies of the Boson peak are in the region where the relation between the frequency of the Brillouin-peak $`\omega _p`$ and the momentum $`p`$ is still linear.
4. A secondary peak at frequencies smaller than the Brillouin one develops for large momenta, becoming dominant for $`p/p_00.5`$.

A number of basic insights on the spectral properties of glasses have been obtained by means of molecular dynamics simulations (Mazzacurati et al. 1996, Dell’Anna et al. 1998, Horbach et al. 1998, Ribeiro et al. 1998, Sampoli et al. 1998, Allend et al. 1999, Feldman et al. 1999, Taraskin and Elliot 1999) on systems such as argon, silica and water. Among the others, let us point out the fact that, in the glass phase, the high frequency dynamics is very well described in the framework of the harmonic approximation. Both dynamic quantities like the dynamic structure factor (Ruocco et al. 2000) and thermodynamic observables as the specific heat (Horbach at al 1999) have been shown to be correctly described by models where only vibrations around quenched positions are taken into account. In figure (1) the numerical results obtained by Ruocco et al. (2000) are shown, which illustrate the above-mentioned spectral features.

On the theoretical side, the study of the short time properties of glasses has been attempted from two different (and somewhat complementary) points of view.

* On one hand, it has been shown (Götze and Mayr 2000) that the Mode Coupling Theory (MCT), which usually describes the long time limit of the time correlators, can be modified in the glassy phase so as to describe only the excitations around the quenched structure. Within such an approach, $`S(p,\omega )`$ was computed for a hard-sphere glass by means of a generalized hydrodynamics approximation. In this framework, a Brillouin peak with a linear dispersion relation (propagating excitations) is obtained up to $`p/p_00.5`$, as well as a secondary peak for higher values. Interestingly, the peak width $`\mathrm{\Gamma }`$ was shown to follow the simple scaling law:

$$\mathrm{\Gamma }_{MCT}p^2$$
(3)
only in the very low momenta regime $`p/p_0<0.1`$, whereas a different law applies to the range of momenta comparable with the experimental and simulation values (see the figure (8) in Götze and Mayr 2000) which seems hard to reduce to the simple form (2). Moreover, it has been shown (Ruocco et al. 2000) that, resorting to the harmonic approximation, the memory function can be simply obtained by means of the eigenvalues and eigenvectors of the Hessian matrix.
* A second approach relies on the study of the statistical properties of random matrices (Metha 1991), since within the harmonic approximation the whole dynamical features are encoded in the system’s Hessian matrix. A crucial point is to distinguish two different classes of systems:

1. vibrational systems whose disorder has a topological origin, as in glasses (Elliot 1983)
    2. vibrational systems on a lattice, where random couplings constants (Schirmacher et al. 1998, Montagna et al. 1999, Kantelhardt et al. 2000, Martin-Mayor et al. 2000, Taraskin et al. 2001) are introduced in order to mimic the behavior of real glasses.

The latter class has recently been studied in some detail (Schirmacher et al. 1998, Taraskin et al. 2001) in order to give an insight into the spectral properties, in particular into the Boson Peak. It turns out (Martin-Mayor et al. 2000), however, that these models definitely miss the $`p^2`$ behavior of the peak width. Because of the long range order due to the underlying lattice, the general beahviour is given instead by

$$\mathrm{\Gamma }_{CPA}p^4$$
(4)
even at very low momenta. The above result can be easily obtained within the CPA approximation and can be checked by the direct inspection of the eigenvectors. The discrepancy with the low momenta behaviour of glasses being due precisely to the lack of topological disorder in lattice models, it seems reasonable to address the study of the former class. The problem turns out to reduce to the study of a very special class of random matrices, which have been called Euclidean Random Matrices (ERM) (Mezard et al. 1999). Is is worthwhile to note that this approach would allow to address, besides the study of the high frequency regime of glasses, even the problem of instantaneous normal modes (INM) i.e. the statistical properties of the Hessian matrix of a liquid at equilibrium (Wu and Loring 1992, Wan and Stratt, 1994, Keyes 1997, Biroli and Monasson 1999), the framework being very similar (Cavagna et al. 1999). Nevertheless, in the following we shall focus only on the latter. Recently (Martin-Mayor et al. 2001), the dynamical structure factor $`S(p,\omega )`$ was analytically computed. The computation relies on a pertubative expansion, the expansion parameter being the inverse density particle $`1/\rho `$. Let us recall what those pertubative results undoubtedly have shown:

+ At $`\rho \mathrm{}`$, the dynamical structure factor is composed of a single delta function, whose position changes linearly with the external momentum $`p`$ at low momenta, representing the undamped propagation of a sound wave in an elastic medium. This can be easily understood by considering the infinite number of particles per wavelength in that limit.
    + For high but finite densities, the disorder of the position of the particles involved in the propagation causes a broadening of the Brillouin peak. In other words, the plane waves are no longer eigenstates of the Hessian, hence a finite spreading of the eigenfrequencies involved arises. In principle no analogy with the Rayleigh scattering should be expected. As a matter of fact, the perturbative computation, performed up to $`1/\rho ^2`$ order, shows that at low enough momenta the general behaviour is instead given by:

    $$\mathrm{\Gamma }_{ERM}(1/\rho )^2p^2.$$
    (5)
    Yet at density $`\rho =1`$ a significant deviation from that behaviour is seen as $`p/p_0`$ becomes greater than $`0.1`$, in a good qualitative agreement with the results of the MCT approach.

The perturbative calculation has the problem that both the structure factor and DOS end abruptly at a cut-off frequency. Moreover, the experimental densities are of order $`\rho 1`$ or less (in reduced units), hence the ability of perturbation theory to catch the features of real systems is at least questionable.

In order to overcome those difficulties, we present a resummation scheme for the perturbative results (Grigera et al. 2001) which includes terms at all orders in $`1/\rho `$. Within this scheme, a very good agreement with numerical results, even at large frequencies, is found for values of $`\rho >0.3`$ (reduced units).

## II harmonic approximation

We shall assume that the particles of our system in the glass phase can only oscillate around their equilibrium position, claiming that this is enough to predict the high frequency properties. Here by “equilibrium” we do not mean thermodynamic equilibrium (the thermodynamics of glasses is an entirely different issue than the one addressed here), but rather mechanical equilibrium (i.e. a position where the forces on all particles are zero). For the sake of simplicity, let us consider displacements only along a given direction $`𝒖`$:

$$𝒙_j(t)=𝒙_j^{\mathrm{eq}}+𝒖\phi _j(t)$$
(6)
Since the number of equilibrium positions $`\{𝒙_j^{\mathrm{eq}}\}`$ available to the system is actually infinite, growing exponentially with the number of particles $`N`$, one can ask which one is to take in order to compute the vibrational spectra. The answer lies on the so called self averaging hypothesis, which has proved to be correct for many observables in disordered systems. Hence, we shall assume that in the thermodynamic limit, the spectra obtained for the different realizations of disorder is equal to the one computed by considering the average over the disorder. This statistical approach leads to a tremendous simplification of our task, as shown in the following, because the model is defined simply by the probability distribution of the random variables $`\{𝒙_j^{\mathrm{eq}}\},j=1,\mathrm{},N`$.

Within that harmonic framework, the energy of the system is ($`m`$ is the mass of the particles and $`\mathrm{\Omega }`$ is a frequency scale)

$$V(\{\phi _i\})=\frac{m\mathrm{\Omega }^2}{2}\underset{i,j}{}f(𝒙_i^{\mathrm{eq}}𝒙_j^{\mathrm{eq}})(\phi _i\phi _j)^2.$$
(7)
We choose units such that $`m=1`$ and $`\mathrm{\Omega }=1`$. Thus the Hessian matrix

$$M_{ij}=\delta _{ij}\underset{k=1}{\overset{N}{}}f(𝒙_i^{\mathrm{eq}}𝒙_k^{\mathrm{eq}})f(𝒙_i^{\mathrm{eq}}𝒙_j^{\mathrm{eq}})$$
(8)
is an Euclidean random matrix, whose spectral properties we are studying.

In particular, we are interested both in the dynamic structure factor and in the density of states (DOS).

* The dynamic structure factor $`S(p,\omega )`$, roughly speaking, gives the spectrum of states “excited” by a plane wave with a given momentum $`p`$. It can be obtained from the Hessian matrix of the potential (7) in the classical limit and in the one excitation approximation (see for instance Martin-Mayor et al. 2001 for a detailed discussion). It reads

$$S^{(1)}(p,\omega )=\frac{k_\mathrm{B}Tp^2}{\omega ^2}\overline{\underset{n}{}|n|𝒑|^2\delta (\omega \omega _n)}$$
(9)
where $`|n`$ are the eigenvectors of the Hessian matrix (8), the eigenfrequencies $`\omega _n`$ are the square root of the eigenvalues (which are all positive, see equation(7)), $`T`$ is the temperature, the overline stands for the average on the equilibrium positions $`\{𝒙_j^{\mathrm{eq}}\}`$, while $`|𝒑`$ stands for a momentum plane wave ($`j|𝒑=\mathrm{exp}[\mathrm{i}𝒑𝒙_j^{eq}]/\sqrt{N}`$, where $`|j`$ represents the vector where the $`j`$-th particle has displacement $`𝒖`$ and the other particles do not move).
* The DOS describes instead the density of ’all’ the vibrational states existing within the system. Interestingly, at the level of one-excitation approximation, the DOS can be obtained by:

$$g(\omega )=\frac{\omega ^2}{k_\mathrm{B}Tp^2}S^{(1)}(p\mathrm{},\omega )$$
(10)
The above theoretical result has been pointed out only very recently (Martin-Mayor et al. 2001). Although it can be checked very easily by means of numerical simulations, we believe it would be also very interesting to check it on real systems, using the experimental data at large momenta which are now available. On the other hand, at very high frequencies the one-excitation approximation does not hold, and many-excitations contributions should be taken into account, hence the relation (10) is obviously only an approximation, whose reliability in describing real systems is a very interesting matter. Moreover, we shall see in the following that (10) is crucial in obtaining a model independent derivation of the actual small $`p`$ behaviour of the width of the Brillouin peak, holding only for topologically disordered systems. As a matter of fact, the same argument cannot be applied to lattice models, because the relation (10) is not true in topologically ordered models.

As mentioned, the Hessian is averaged over the disorder, i.e. over all the allowed equilibrium positions. Those are clearly distributed in a highly correlated manner, due to the hard-core repulsion and long range attraction of the potential. It was shown (Martin-Mayor et al. 2001) that, at the level of the superposition approximation, the correlations can be taken into account by doing an uncorrelated average while at the same time suitably renormalizing the interaction $`f`$: if the $`k`$ point static correlation function is written as

$$g(r_1,\mathrm{}r_k)=g(r_1)\mathrm{}g(r_k),$$
(11)
then the inverse-density expansion for the dynamic structure factor of a system with correlation function $`g(r)`$ and force $`F(r)`$ can be obtained from the expansion of a fully uncorrelated system defining a dressed interaction:

$$f(r)g(r)F(r)$$
(12)
where the value of the “spring” constant $`F(r)`$ is weighted by the probability to find two centers of oscillations at the relative distance $`r`$.

## III Non perturbative results

Let $`\widehat{f}(p)`$ be the Fourier transform of $`f(r)`$. The main ’object’ to compute is the resolvent of the Hessian matrix (8), which depends on the complex variable $`z=\omega ^2+\mathrm{i}\eta `$. It can be written in the following way:

$`G(p,z)`$ $`=`$ $`{\displaystyle \frac{1}{zϵ(p)\mathrm{\Sigma }(p,z)}},`$ (13)
$`ϵ(p)`$ $`=`$ $`\rho [\widehat{f}(0)\widehat{f}(p)].`$ (14)
where complex self-energy $`\mathrm{\Sigma }(p,z)`$ has been introduced. Exploiting the well known relation between resolvent and dynamical structure factor (letting aside the prefactor $`k_\mathrm{B}Tp^2`$, unessential for the $`\omega `$ dependence):

$$S^{(1)}(p,\omega )=\frac{1}{\omega \pi }\underset{\eta 0}{lim}\mathrm{Im}G(p,\omega ^2+i\eta ).$$
(15)
the following connections between the main features of the dynamical structure factor and the self-energy are established:

* The ’bare’ dispersion relation $`ϵ(p)`$, which would give the position of the peak in the elastic medium limit, is renormalized by the real part of the self-energy $`\mathrm{\Sigma }^{}(p,z)`$. This gives $`\omega ^{renorm}(p)`$, the position of the maximum of the structure factor in the frequency domain. Let us note that $`\omega ^{renorm}(p)`$ is certainly linear for small $`p`$, as expected.
* The imaginary part $`\mathrm{\Sigma }^{\prime \prime }(p,z)`$ computed at the position of the peak $`\omega =\omega ^{renorm}(p)`$ gives the width, $`\mathrm{\Gamma }(p)`$, of the $`S^{(1)}(p,\omega )`$ by means of:

$$\mathrm{\Sigma }^{\prime \prime }(p,\omega ^{renorm}(p))=\omega ^{renorm}(p)\mathrm{\Gamma }(p)$$
(16)

The self-energy can be obtained as a series in $`1/\rho `$: the $`k^{\mathrm{th}}`$-order corresponds to $`k`$ particle-label repetitions when calculating the moments of the $`S^{(1)}(p,\omega )`$ (the details can be found in Martin-Mayor et al. 2001). It is easy to show that the sum of all the cactus diagrams (see figure (2)) is given by the solution of the following integral equation:

$$\mathrm{\Sigma }(p,z)=\frac{1}{\rho }\frac{d^Dq}{(2\pi )^D}\frac{\left[\rho \left(\widehat{f}(𝒒)\widehat{f}(𝒑𝒒)\right)\right]^2}{zϵ(q)\mathrm{\Sigma }(q,z)}.$$
(17)
Interestingly, Equation(17) provides a model-independent derivation of the $`p^2`$ scaling of width of the peak of the $`S^{(1)}(p,\omega )`$. Indeed, the large $`q`$ contribution to the imaginary part of the integral in equation (17) is, because of (10)

$$\mathrm{\Sigma }_{\mathrm{}}^{\prime \prime }(p,z)=\pi \rho g_\lambda (\lambda )\frac{d^Dq}{(2\pi )^D}\left(\widehat{f}(𝒒)\widehat{f}(𝒑𝒒)\right)^2.$$
(18)
where $`g_\lambda (\lambda )`$ is the density of states in the domain of eigenvalues ($`\lambda =\omega ^2,g_\lambda (\omega ^2)=\frac{g(\omega )}{2\omega }`$). If the spectrum is Debye-like we have $`g_\lambda (\lambda )\lambda ^{0.5}`$, and it is straightforward to show that (18) is proportional to $`\omega ^{renorm}(p)p^2`$. Then the relation (16) implies the scaling:

$$\mathrm{\Gamma }(p)p^2$$
(19)
irrespective of the function $`f(r)`$. Clearly this is only the large $`q`$ contribution to the integral, but it has been shown (Grigera et al. 2001) that it indeed controls the peak width at small $`p`$.

## IV the Gaussian case

In the following we shall study the simple case

$$f(r)=\mathrm{exp}[r^2/(2\sigma ^2)].$$
(20)
This may seem an oversimplification, too distant from any realistic case. However, it is not actually the case, at least for small momenta. For the sake of comparison we can see in figure (3) the Fourier transform of the function $`g(r)F(r)`$ for Argon at very low temperature ($`10K`$) together with that of (20). We thus expect our results to be reasonable for momenta smaller than the first zero of the Fourier transform of $`g(r)F(r)`$, which in general is close to the maximum of the static structure factor, $`p_0`$. Since the Fourier transform of our force decreases an order of magnitude by $`p_0=2/\sigma `$, we shall take this as our $`p_0`$ during the following discussion, and $`\sigma `$ will be our unit of length.

We have numerically solved equation (17) to find the self-energy for several values of $`\rho `$ with the Gaussian choice (20), thus obtaining the structure factor and the density of states, in the eigenvalue domain. In figure (4) we show $`S^{(1)}(p,\lambda )`$ for several values of the momentum as obtained from equation (17) for $`\rho =1`$. We also plot the structure factor obtained numerically by the method of moments (Turchi et al. 1982, Benoit 1989, Benoit et al. 1992) (for lower momenta, the comparison with the method of moments cannot be done due to finite volume effects, Martin-Mayor et al. 2001). A very good agreement with the numerical data is achieved. Note the absence of the secondary peak in this model (see below).

We found that for densities down to $`\rho 0.6`$, the agreement is good for all momenta, in fact good results are also obtained for the density of states (see figure (5)). For lower densities, the DOS starts to deviate from a Debye behavior at small $`\lambda `$, in contradiction with the numerical results. Moreover, below $`\rho 0.31`$ the approximation gives a DOS with an unphysical negative spectrum. The negative spectrum develops continuously as a function of $`\rho `$. Also included in figure (5) is the one-loop result. Notice that even for $`\rho =1`$, the cactus resummation fails to reproduce the exponential decay of the density of states. This is not unexpected, however, due to the non perturbative nature (in $`1/\rho `$) of this tail (Zee and Affleck, 2000).

It is interesting to look at the small $`\lambda `$ limit of $`g_\lambda (\lambda )`$. As seen in figure (5) (top) the behavior is very nearly Debye (i.e. $`\sqrt{\lambda }`$) for small $`\lambda `$. An excess of states relative to the Debye case develops for higher $`\lambda `$, but in a region of eigenvalues well beyond the linear dispersion regime (see bottom of figure (4)). Thus this peak, similar to that found in (Schirmacher et al. 1998), should not be regarded as a Boson peak.

Finally, let us look at the scaling of the position and the width of the peak in the frequency domain. In figure (6) we show the frequency corresponding to the Brillouin peak as a function of the external momentum $`p`$. Let us note that a nearly linear dependance, (pointing likely to a propagating excitation) persists up to $`p/p_00.60.7`$. In figure (7) instead we plot $`\mathrm{\Gamma }(p)`$, obtained by means of (16). As expected, the $`p^2`$ scaling is found for very small momenta, which crosses over to a region where a simple law as (2) is not suitable to describe the real beahviour of the system.

Note that the region where the $`p^2`$ scaling is actually found, i.e. $`p/p_0<0.1`$ is quite different from the region explored by X-rays and neutrons scattering experiments, which rather spans the momenta $`0.1<p/p_0<0.5`$, It is worthwile to note that the same conclusion can be drawn from the results of Götze and Mayr using MCT for hard spheres (Götze & Mayr 2000).

## V Conclusions and outlooks

In summary, we have applied the random-matrix approach to the study of the high-frequency excitations of glassy systems. We have presented a resumation scheme (Grigera et al. 2001) that greatly enhances the predictive power of the $`1/\rho `$ expansion prevoisuly obtained (Martin-Mayor et al. 2001). We have compared our analytical calculation with numerical results, finding that for not too low densities the only failure of the resumation scheme is its inability to repoduce the exponential decay of the density of states.

Our calculations have been performed choosing a gaussian force, in order to modelize the interparticle force dressed with the pair correlation function, and only collinear displacements have been considered. We believe that our results are relevant for realistic glasses, at least in the regime $`p/p_0<1`$. The basic equation (17) shows that, in the general case, the position of the spectral peak is linear at small momenta, and the numerical solution of the integral equation in the Gaussian case shows that the linearity persists up to $`p/p_00.60.7`$.

Moreover, the width of the spectral peak turns out to be proportional to $`p^2`$ in the limit of $`p0`$ (see equation (18)) irrespective of the potential function $`f(x)`$. Interestingly, in the Gaussian case the scaling law holds for a momentum range one order of magnitude smaller than the experimental one. At larger momenta, the law is more complicated (see figure (7)). Being this result quite similar to MCT result for an hard-spheres system (Götze and Mayr 2000), we believe that it could have a certain degree of univerality.

This model lacks an important feature of the experimental spectra, namely the secondary peak of the $`S(q,\omega )`$ (which some authors have interpreted as the Boson peak, see Horbach et al. 1998). We believe this to be related to two important ingredients missing in our model: the vectorial nature of the vibration, and a detailed consideration of the particle correlations. The latter should not affect our results for momenta much smaller than the first maximum of the static structure factor, $`p_0`$. Since the secondary peak appears for $`pp_0`$, we should not really expect to describe it: in this momentum range the dispersion relation is no longer monotonic in real glasses. Furthermore, Dell’Anna et al. 1998 have suggested that transversal excitations may play a prominent role in the raising of the secondary-peak. Our approach can be extended to include transverse displacements, though, and work in this direction is in progress.

VMM was partly supported by CICyT AEN99-0990, AEN97-1693 and M.E.C. TSG was supported in part by CONICET (Argentina).

## VI References

Allen P. B. ,Feldman J. L. Fabian J. and Wooten F. 1999 Phil. Mag. B 79
Angelani L. Di Leonardo R. Ruocco G. Scala A. Sciortino F. 2000 Phys. Rev. Lett. 85, 5356
Benassi P. Krisch M. Masciovecchio C. Mazzacurati V. Monaco G. Ruocco G. Sette F. and Verbeni R. 1996 Phys. Rev. Lett. 77, 3835
Benoit C. 1989 J. Phys.: Condens. Matter 1, 335
Biroli G. and Monasson R. 1999 J. Phys. A: Math. Gen. 32, L255
Buchenau U. Prager M. N cker N. Dianoux A. J. Ahmad N. and Phillips W. A. 1986 Phys. Rev. B 34, 5665
Broderix K. Bhattacharya K. K. Cavagna A. Zippelius A. Giardina I. 2000 Phys. Rev. Lett. 85, 5360
Cavagna A. Giardina I. Parisi G. 1999 Phys. Rev. Lett. 83, 108
Dell’Anna R. Ruocco G. Sampoli M. and Viliani G. 1998 Phys. Rev. Lett. 80, 1236
Ducastelle F. and Treglia G. 1982 J. Phys. C 15, 2891
Elliot S. R. 1983 Physics of amorphous materials (England: Longman)
Feldman J. L. Allen P. B. and Bickham S. R. 1999 Phys. Rev. B 59, 3551
Fioretto D. Buchenau U. Comez L. Sokolov A. Masciovecchio C. Mermet A. Ruocco G. Sette F. Willner L. Frick B. Richter D. and Verdini L. 1999 Phys. Rev. E 59, 4470
Foret M. Courtens E. Vacher R. Suck J. B. 1996 Phys. Rev. Lett. 77, 3831
Götze W. and Mayr M. R. 2000 Phys. Rev. E 61, 587
Grigera T.S., Martín-Mayor V., Parisi G. and Verrocchio P. 2001 cond-mat/0102230
Horbach J. Kob W. Binder K. 1998 J. Non-Cryst solids 235, 320
Horbach J. Kob W. Binder K. 1999 J. Phys. Chem. B 103, 4104
Kantelhardt J. W. Russ S. Bunde A. cond-mat/0012392; Phys. Rev. B (in press)
Keyes T. 1997 J. Chem. Phys. 101, 2921
Martín-Mayor V. Parisi G. Verrocchio P. 2000 Phys. Rev. E 62, 2373
Martín-Mayor V. Mèzard M. Parisi G. Verrocchio P. 2001 J. Chem. Phys. (in press)
Masciovecchio C. Ruocco G. Sette F. Krisch M. Verbeni R. Bergmann U. and Soltwisch M. 1996 Phys. Rev. Lett. 76, 3356
Masciovecchio C. Ruocco G. Sette F. Benassi P. Cunsolo A. Krisch M. Mazzacurati V. Mermet A. Monaco G. and Verbeni R. 1997 Phys. Rev. B 55, 8049
Masciovecchio C. Monaco G. Ruocco G. Sette F. Cunsolo A. Krisch M. Mermet A. Soltwisch M. and Verbeni R. 1998 Phys. Rev. Lett. 80, 544
Mazzacurati V., Ruocco G. and Sampoli M., 1996, Europhys. Lett. 34, 681
Metha M. L. 1991 Random matrices, (Academic Press)
Mèzard M. Parisi G. Zee A. 1999 Nucl. Phys. B559, 689
Monaco G. Masciovecchio C. Ruocco G. Sette F. 1998 Phys. Rev. Lett. 80, 2161
Montagna M. Ruocco G. Viliani G. Dell’Anna R. Di Leonardo R. Dusi R. Monaco G. Sampoli M. and Scopigno T., 1999 Phys. Rev. Lett. 83, 3450
Parisi G. 1994 Field-Theory, Disorder and Simulations, World Scientific.
Ribeiro M. C. C. Wilson M. and Madden P. A. 1998 J. Chem. Phys. 108, 9027
Ruocco G. Sette F. Di Leonardo R. Fioretto D. Krisch M. Lorenzen M. Masciovecchio C. Monaco G. Pignon F. and Scopigno T. 1999 Phys. Rev. Lett.  83, 5583
Ruocco G. Sette F. Di Leonardo R. Monaco G. Sampoli M. Scopigno T. and Viliani G. 2000 Phys. Rev. Lett. 84, 5788
Sampoli M. Benassi P. Dell’Anna R. Mazzacurati V. and Ruocco G. 1998, Phil. Mag. B 77, 473
Schirmacher W. Diezemann G. and Ganter C. 1998 Phys. Rev. Lett.  81, 136
Sette F. Krisch M. Masciovecchio C. Ruocco G. and Monaco G. 1998 Science 280, 1550
Sokolov A. P. Buchenau U. Richter D, Masciovecchio C. Sette F. Mermet A. Fioretto D. Ruocco G. Willner L. and Frick B. 1999 Phys. Rev. E 60, 2464
Taraskin S. N. and Elliot S. R. 1999 Phys. Rev. B 59, 8572
Taraskin S. N. Loh Y. L. Natarajan G. and Elliott S. R. 2001 Phys. Rev. Lett.  86, 1255
Turchi P. Benoit C. Royer E. and Poussigue G. 1992 J. Phys.: Condens. Matter 4, 3125
Wan Y. and Stratt R. 1994 J. Chem. Phys 100, 5123
Wu T. M. and Loring R. F. 1992 J. Chem. Phys. 97, 8368
Zee A. and Affleck I. 2000 J. phys: cond. mat. 12, 8863


format: text

text:
# The Large Area Lyman Alpha Survey

## 1. Introduction

The first major burst of star formation in a young, high redshift galaxy is expected to produce copious numbers of ionizing photons thanks to an abundance of hot, massive main sequence stars. The interstellar medium that gave birth to the stars can then convert 2/3 of these photons into Lyman-$`\alpha `$ photons as neutral hydrogen atoms are ionized and recombine. The strong Lyman-$`\alpha `$ line emission that can result was first proposed as a signpost of young galaxies in formation by Partridge & Peebles (1967; hereafter PP67).

However, searches for Lyman-$`\alpha `$ emission at high redshift failed to discover a field population of such galaxies for thirty years. (See review by Pritchet 1994 and references therein; and more recently, Thompson & Djorgovski 1995; Thommes et al. 1998.) These searches placed stringent upper limits on the abundance and brightness of the Lyman-$`\alpha `$ line at high redshift. More recent searches have finally found Lyman-$`\alpha `$ emitters in the field (Cowie & Hu 1998; Hu, Cowie & McMahon 1998; Pascarelle, Windhorst, & Keel 1998; Dey et al. 1998; Hu, McMahon, & Cowie 1999; Steidel et al. 2000; Kudritzki et al. 2000; Manning et al. 2000), but at line luminosities $`10^2`$ below the early predictions by PP67.

Several factors could be responsible for this paucity of bright Lyman-$`\alpha `$ sources, including dust effects, small star formation units, and long star formation time scales. (PP67 based their estimates on dust free, Milky Way-sized galaxies that form most of their stars in $`3\times 10^7`$ years.)

Dust can quench this line very effectively: Lyman-$`\alpha `$ photons are resonantly scattered by neutral hydrogen atoms, and so traverse a much longer path than continuum photons do in order to escape the galaxy where they are emitted. This long path enhances the optical depth for line photons, reducing the line equivalent width. Quantitative estimates of this effect are however quite sensitive to the spatial distribution of the ISM (e.g. Neufeld 1991). In contrast, star formation in small units would yield small line and continuum luminosities while leaving equivalent widths unaffected. Under this scenario, the expected numbers of Lyman-$`\alpha `$ sources would rise in inverse proportion to their luminosity. Protracted star formation would also reduce both line and continuum luminosities. Additionally, a star formation event lasting longer than the lifetime of a massive star would allow metals and dust produced by the first generation of stars to pollute the galaxy’s ISM, resulting again in dust extinction and potential reductions in Lyman-$`\alpha `$ equivalent width.

Good statistical samples of Lyman-$`\alpha `$ emitters that yield distributions of luminosities and equivalent widths should allow us to distinguish among the above possibilities. We are conducting the Large Area Lyman Alpha (LALA) survey to obtain such a sample.

Whatever the reason for the faintness of Lyman-$`\alpha `$ sources, the effects of dust provide another strong reason for seeking these objects: Truly primordial galaxies, taken here to mean chemically unevolved systems, will have large Lyman-$`\alpha `$ equivalent widths because of their low dust content. Lyman-$`\alpha `$ surveys should therefore preferentially select for this most interesting class of high-redshift galaxy. This is not to say that all Lyman-$`\alpha `$-selected galaxies will be primordial, nor even that a large fraction will be. Dustier galaxies can have substantial Lyman-$`\alpha `$ equivalent widths if their ISM is patchy, or if it is in bulk motion (Kunth et al 1998). However, if chemically primitive galaxies are to be found at all, the Lyman-$`\alpha `$ search is a good tool to discover them.

## 2. The Survey

The Large Area Lyman Alpha (LALA) survey is possible because of new wide field instrumentation. We exploit the 81922 pixel CCD Mosaic camera at the Kitt Peak National Observatory’s 4m Mayall telescope to achieve a survey efficiency (measured in the product $`A\mathrm{\Omega }`$ of collecting area and solid angle) some 6 times greater than comparably deep narrowband surveys at Keck. The “core” survey uses five overlapping $`80`$Å bandpass filters ($`\lambda _c=6559`$, $`6611`$, $`6650`$, $`6692`$, and $`6730`$Å) to look for Lyman-$`\alpha `$ line emission at $`4.37<z<4.57`$. The corresponding survey volume is about $`8.5\times 10^5`$ comoving Mpc3 per 36 field for $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, $`\mathrm{\Omega }=0.2`$, $`\mathrm{\Lambda }=0`$. We integrate for about 6 hours per filter per field, and achieve a limiting sensitivity $`2\times 10^{17}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ in the line. Further details, including a brief description of data reduction, are given in Rhoads et al (2000).

Additionally, we use supporting broad band images from the NOAO Deep Widefield Survey (Bw, R, and I filters; see Schommer, Jannuzi, & Dey 2001, this meeting) and from our own data (V and SDSS z’ filters), and have recently obtained $`8200`$Å narrowband data to extend the LALA survey to $`z5.7`$.

Spectroscopic followup of some LALA candidates has been carried out using the Keck Observatory’s Low-Resolution Imaging Spectrograph. To date, we have confirmed five emission line galaxies at $`z>4`$. Two representative spectra are shown in figures 1 and 2. One of these two is a galaxy with very little continuum that would have gone undetected in a broadband survey of even extreme depth, while the other is a less spectacular emission line source with a clearly detected continuum and Lyman break.

## 3. The Lyman-$`\alpha `$ Source Population

Based on our photometric results and on the statistics of our first spectroscopic followup, we inferred a source density of $`4000`$ Lyman-$`\alpha `$ emitters per square degree per unit redshift at $`z=4.5`$, with an uncertainty dominated by small numbers of spectra (Rhoads et al 2000).

We can estimate the star formation rate density in the Lyman-$`\alpha `$ galaxies in two quasi-independent ways. First, we coadd the observed $`z^{}`$ band flux of all our good candidates and apply starburst models to convert this rest-frame UV continuum into a star formation rate. Correcting this rate for our estimate that $`1/3`$ of the good candidates are actually $`z=4.5`$ Lyman-$`\alpha `$ emitters, we obtain a rate of $`13\times 10^3M_{}\text{yr}^1\text{Mpc}^3`$. Second, we coadd the continuum-corrected line fluxes of all our good candidates and apply a standard conversion factor to go from Lyman-$`\alpha `$ flux to ionizing luminosity and thence star formation rate. Applying again a $`1/3`$ correction factor to account for foreground contamination, we find a lower value of $`2.2\times 10^3M_{}\text{yr}^1\text{Mpc}^3`$. The difference in these two estimates may be due to the effects of dust, which would affect the Lyman-$`\alpha `$ based estimate more strongly than the UV continuum estimate. Alternatively, our simple correction factors of $`1/3`$ for foreground contamination could be wrong if the $`z4.5`$ sources have a systematically stronger line to continuum ratio than the remaining candidates. While these star formation rates are still quite uncertain (pending further spectroscopy), they are comparable to the Lyman break galaxy rate at $`z=4`$, which is $`5\times 10^3M_{}\text{yr}^1\text{Mpc}^3`$ (Madau et al 1996).

The spectra in figures 1 and 2 are representative of our spectroscopically confirmed $`z>4`$ galaxies in the sense that about half have rest frame equivalent widths $`>100`$Å, which is near the upper limit for star forming galaxies ($`100`$ to $`200`$Å; Charlot & Fall 1993). This suggests that some of the LALA sources may indeed be dust-free and chemically primitive. Further spectroscopy to rule out weak narrow-lined AGNs and/or galactic winds (see Kunth et al 1998) will help confirm this possibility.

## 4. Future Directions

Photometric redshift constraints derived from our suite of broadband filters, combined with consistency checks from our growing spectroscopic sample, will enable us to use the entire LALA sample to study population properties of the Lyman-$`\alpha `$ emitters. These include the luminosity function, distribution of equivalent widths, and clustering properties. Planned near-infrared followup of selected sources will allow us to confirm the Lyman-$`\alpha `$ line identification through an \[O II\] $`\lambda `$3727 line search and to look for older stellar populations through a study of rest-frame optical light. Finally, searches in new narrowband filters will allow us to study the evolution of the Lyman-$`\alpha `$ source population. Ultimately, we hope to identify statistical samples of primitive galaxies at a range of redshifts from $`z=4.5`$ to $`z>6`$.

### Acknowledgments.

We thank Andy Bunker and Steve Dawson for help with the spectroscopic observations, and the IRAF team for writing and helping with the MSCRED package. JER’s research is supported by an STScI Institute Fellowship. SM’s research is supported by NASA through Hubble Fellowship grant # HF-01111.01-98A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. The data presented herein were obtained at the Mayall Telescope of the Kitt Peak National Observatory; and at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Keck Observatory was made possible by the generous financial support of the W.M. Keck Foundation.

## References

Charlot, S., & Fall, S. M. 1993, ApJ 415, 580

Cowie, L. L., & Hu, E. M. 1998, AJ 115, 1319

Dey, A., Spinrad, H., Stern, D., Graham, J. R., & Chaffee, F. H. 1998, ApJ 498, 93

Hu, E. M., Cowie, L. L., & McMahon, R. G. 1998, ApJ 502, L99

Hu, E. M., McMahon, R. G., & Cowie, L. L. 1999, ApJ 522, L9

Kudritzki, R.-P., Mendez, R. H., Feldmeier, J. J., Ciardullo, R., Jacoby, G. H., Freeman, K. C., Arnaboldi, M., Capaccioli, M., Gerhard, O., & Ford, H. C. 2000, ApJ 536, 19

Kunth, D., Mas-Hesse, J. M., Terlevich, E., Terlevich, R., Lequeux, J., & Fall, S. M. 1998, A&A 334, 11

Madau, P., Ferguson, H. C., Dickinson, M. E., Giavalisco, M., Steidel, C. C., & Fruchter, A. 1996, MNRAS 283, 1388

Manning, C., Stern, D., Spinrad, H., & Bunker, A. J. 2000, ApJ 537, 65

Neufeld, D. A. 1991, ApJ 370, L85

Partridge, R. B., & Peebles, P. J. E. 1967, ApJ 167, 868

Pascarelle, S. M., Windhorst, R. A., & Keel, W. C. 1998, AJ 116, 2659

Pritchet, C. J. 1994, PASP 106, 1052

Rhoads, J. E., Malhotra, S., Dey, A., Stern, D., Spinrad, H., & Jannuzi, B. T. 2000, ApJ 545, L85

Schommer, R. A., Jannuzi, B. T., & Dey, A. 2001, this meeting.

Steidel, C. C., Adelberger, K. L., Shapley, A. E., Pettini, M., Dickinson, M., & Giavalisco, M. 2000, ApJ 532, 170

Stern, D., & Spinrad, H. 1999, PASP 111, 1475

Thommes, E., Meisenheimer, K., Fockenbrock, R., Hippelein, H., Röser, H.-J., & Beckwith, S. 1998, MNRAS 293, L6

Thompson, D., & Djorgovski, S. 1995, AJ 110, 982


text:
# Influence of statistical fluctuations on 𝐾/𝜋 ratios in relativistic heavy ion collisions

## Abstract

The influence of pure statistical fluctuations on $`K/\pi `$ ratio is investigated in an event-by-event way. Poisson and the modified negative binomial distributions are used as the multiplicity distributions since they both have statistical background. It is shown that the distributions of the ratio in these cases are Gaussian, and the mean and relative variance are given analytically.

Key words: $`K/\pi `$ ratio, statistical fluctuations, relativistic heavy ion collision

The ultimate goal of current and future ultra-relativistic heavy ion program is the production and characterization of an extended volume of deconfined phase of quarks and gluons, the quark-gluon plasma (QGP). The possible existence of QGP as the equilibrium state of strongly interacting matter has been predicted in lattice theories jcc  at sufficiently high temperature and density. So QGP might exist in the early universe and inside the neutron stars. Heavy ion collision at extremely high energies is the only chance to study QGP in a controlled way. In laboratory experiments, however, the QGP can survive only for very short period of time, and it will cool during its expanding and become hadronic matter at freeze-out point. It is very difficult to detect QGP in experiments because most events are, due to collision geometry and quantum fluctuations, without QGP even if the condition for QGP creation can be reached in the experiments. So, a variety of possible signatures for the transient existence of the deconfined state of matter in the collisions have been proposed theoretically and experimentally (see bgsg99 ; pr236 ; na638  for a review). One of the most important and most interesting signatures is the enhanced strangeness production in the process of heavy ion collisions pr142  compared with that in $`pp`$ collisions. The investigation of strangeness production is useful for understanding the mechanism of heavy quark production since no strangeness content exists in the initial state in heavy-ion collisions. Strangeness production can also shed light on the time scale of chemical freeze-out and therefore can carry information about the early stage of a heavy ion collision lpc94 .

Great progress has been made with the proposal and applications of event-by-event analysis zpc54 ; epjc6 ; appb29 ; plb435 ; plb439 ; plb459 ; epjc8  in the study of high energy heavy ion collisions. Event-by-event analysis of heavy ion collisions became possible with the advent of large acceptance detectors la99 . The philosophy of event-by-event physics is based on following speculation: Although the conditions to produce QGP may be reached in every event, the fact that a phase transition is a critical phenomenon implies that it may occur only in a very small sub-sample of events. So the fluctuations accompanying the phase transition will, in effect, be averaged out in the conventional ensemble analyses. The event-by-event analysis searches for fluctuations of observables at the event level, so it can be used to select interesting or anomalous event candidates nato94  with specific dynamical properties. So the new method can provide dynamical information which cannot be obtained from the traditional inclusive spectra.

Recently the event-by-event analysis is used in experiments to analyze strangeness productions in high energy heavy ion collisions. In gr97 ; fs99  the kaon to pion ratio for single event is studied as the event observable with which one can look at fluctuations in the chemical freeze-out stage of the collisions. Of course, the distribution of the $`K/\pi `$ ratio is the result of a combination of statistical and possible dynamical fluctuations in the strangeness production. The method of mixed events, in which there are no momentum correlation among particles by construction, is used in the experimental analyses to estimate the effect of the former on the distribution of $`K/\pi `$ ratio. Both distributions for $`K/\pi `$ ratio from the real data and the mixed events are Gaussian, and the widths are only with a little difference. That means the dominance of statistical fluctuations in strangeness production. Although it is well established that there is no spatial-temporal correlation in the mixed events, we do not know exactly if there exist some other correlations. In the study of $`K/\pi `$ ratio the only quantities concerned are the numbers of produced kaons and pions. The event-mixing technique will not change the total numbers of the particles produced. These numbers carry some important information about the chemical fluctuations at the freeze-out point. So, some information about the strangeness production may still be retained in the mixed events. To investigate the problem theoretically further one should ask: What is the distribution of $`K/\pi `$ ratio if there exist only statistical multiplicity fluctuations in the collisions? Will it be a Gaussian as showed from the mixed events? If yes, what are the mean and the width? Are they in agreement with those obtained from the mixed events? What will be the dependences of them on mean multiplicities of kaons and pions?

In this paper the $`K/\pi `$ ratio for the case with pure statistical fluctuations is studied in an event-by-event way with two assumptions: (1) There is no dynamical fluctuation in the productions of both the kaons and pions, and (2) there is no correlation between the productions of the two kinds of particles. This is, of course, a trivial case, but it is the base for any event-by-event investigation on the influence of dynamical fluctuations on $`K/\pi `$ ratio.

Since we have assumed no dynamical fluctuations in the productions of kaons and pions nor correlation between them, the distributions of the multiplicities of kaons and pions can be specified from some statistical considerations. Two possibilities will be considered in this paper. As a first possibility, the distributions can be Poissonian, i.e. the probabilities with which there are $`k`$ kaons and/or $`\pi `$ pions in an event are, respectively,

$$p_K=\frac{k^k}{k!}\mathrm{e}^k,p_\mathrm{\Pi }=\frac{\pi ^\pi }{\pi !}\mathrm{e}^\pi $$
(1)
with $`k`$ and $`\pi `$ the mean multiplicities of kaons and pions in the collisions. In a global analysis, the global $`K/\pi `$ ratio is then

$$R_g=k/\pi .$$
(2)
From the view point of the theory of coherent state the assumption of Poisson distributions for the multiplicities implies that both kaons and pions are emitted from some coherent sources and that no thermal effect is taken into account. Now we turn to event-by-event analysis of the ratio. For an event with $`k`$ kaons and $`\pi `$ pions the $`K/\pi `$ ratio is

$$r_e\frac{k}{\pi },$$
(3)
and the probability density for $`K/\pi `$ ratio to be in neighborhood of $`r_e`$ can be given by

$$P_{r_e}=\underset{r_e=k/\pi }{}p_K(k)p_\mathrm{\Pi }(\pi )=\underset{\pi =1}{\overset{\mathrm{}}{}}\frac{k^{r_e\pi }}{(r_e\pi )!}\frac{\pi ^\pi }{\pi !}\mathrm{e}^{(1+R_g)\pi },$$
(4)
because the productions of kaons and pions have been assumed to be independent.

Since $`p_K`$ and $`p_\mathrm{\Pi }`$ are peaked at $`k`$ and $`\pi `$ respectively, it is easy to see that $`P_{r_e}`$ will be peaked near $`R_g`$. A natural question is: Can $`P_{r_e}`$ be Gaussian for given $`k`$ and $`\pi `$ ? Since Eq. (4) is too complicated, Monte Carlo simulation can be used to answer this question. To get an event for this case, one needs only to generate two independent random integer numbers (according to two Poisson distributions with means $`k`$ and $`\pi `$ respectively). The random integer number corresponding to $`k`$ is taken to be the kaon multiplicity in the event, the other for pions. Then one gets a corresponding $`K/\pi `$ ratio $`r_e`$ for the event. The distribution of $`r_e`$ from the simulation is given in Fig. 1 with 1 million events used. In the simulation the mean multiplicities for kaons and pions are chosen to be 38 and 200 respectively for illustration. Here we choose $`R_g=0.19`$, the same as shown experimentally in Pb-Pb collisions at 158 $`A`$ GeV, for later comparison with experimental result. One can see that the distribution of the $`K/\pi `$ ratio is indeed approximately a Gaussian. This is consistent with the conclusions given from the mixed events in Refs. gr97 ; fs99 . To get more information about the distribution one may want to know the mean and width of the distribution. From Eq. (4) one can see that the distribution of $`r_e`$ is not symmetric about $`R_g`$, so the mean value of the distribution is, to some extent, different from $`R_g`$. The mean value of the distribution for the trivial case can be calculated analytically as

$$R_e=\frac{k}{\pi }\underset{k=0}{\overset{\mathrm{}}{}}\underset{\pi =1}{\overset{\mathrm{}}{}}\frac{k}{\pi }p_K(k)p_\mathrm{\Pi }(\pi )/\underset{k=0}{\overset{\mathrm{}}{}}\underset{\pi =1}{\overset{\mathrm{}}{}}p_K(k)p_\mathrm{\Pi }(\pi )$$
(5)
$$R_g\underset{\pi =1}{\overset{\mathrm{}}{}}\frac{\pi ^{\pi +1}}{\pi \pi !}\mathrm{e}^\pi .\text{ }$$
(6)
In the transition from Eq. (5) to Eq. (6) condition $`\pi 1`$ has been used. The relative accuracy in Eq. (6) is of the order of $`\mathrm{e}^\pi `$, thus the approximation is quite good for $`\pi `$ a few hundreds which is quite normal in current high energy heavy ion collisions. By using the generating function for the Poisson distribution

$$G_\mathrm{\Pi }(x)=\underset{\pi =0}{\overset{\mathrm{}}{}}\frac{\pi ^\pi x^\pi }{\pi !}\mathrm{e}^\pi =\underset{\pi =1}{\overset{\mathrm{}}{}}x^\pi p_\mathrm{\Pi }(\pi )=\mathrm{e}^{\pi (x1)},$$
(7)
the mean value of the distribution $`R_e`$ can be expressed as

$$R_e=R_g_0^1𝑑x\pi \frac{G_\mathrm{\Pi }(x)G_\mathrm{\Pi }(0)}{x}.$$
(8)
We are particularly interested in the ratio $`R`$ between $`R_e`$ and $`R_g`$, which is a measure of the deviation of the two $`K/\pi `$ ratios obtained from global and event-by-event analyses, respectively. In our trivial case

$$R\frac{R_e}{R_g}=\pi _0^1𝑑x\frac{G_\mathrm{\Pi }(x)G_\mathrm{\Pi }(0)}{x},$$
(9)
which depends only on the mean multiplicity of pions. It can be seen that $`R`$ will approach 1.0 when $`\pi \mathrm{}`$. In this limit the mean values of $`K/\pi `$ ratios from global and event-by-event analyses are equal, as can be naively expected. The change behavior of $`R`$ as a function of $`\pi `$ is shown in Fig. 2. $`R`$ is large for lower mean pion multiplicity and decreases with $`\pi `$. But the values of $`R`$ are always close to the trivial value 1.0 within an accuracy of about 1% in the graphed range of the mean pion multiplicity. So one can hardly find the difference between $`R_e`$ and $`R_g`$. Using the generating function for pion multiplicity distribution, one can also calculate the width of the $`K/\pi `$ ratio distribution. First one can get

$$\left(\frac{k}{\pi }\right)^2=k^2_0^1𝑑x\frac{G_\mathrm{\Pi }(x)G_\mathrm{\Pi }(0)}{x}\mathrm{ln}\frac{1}{x},$$
(10)
where $`k^2=k^2+k`$, since the distribution of the kaon multiplicity is assumed to be a Poissonian. The width of the $`K/\pi `$ ratio distribution is then

$$\sigma =\sqrt{\left(\frac{k}{\pi }\right)^2\frac{k}{\pi }^2}.$$
(11)
The relative variance is then

$$\frac{\sigma }{R_e}=\sqrt{\frac{\left(1+\frac{1}{R_g\pi }\right)_0^1𝑑x\mathrm{ln}\frac{1}{x}\frac{G_\mathrm{\Pi }(x)G_\mathrm{\Pi }(0)}{x}}{\left[_0^1𝑑x\frac{G_\mathrm{\Pi }(x)G_\mathrm{\Pi }(0)}{x}\right]^2}1}.$$
(12)
Different from $`R`$ the relative variance depends on both $`\pi `$ and $`R_g`$ (thus on $`k`$). The behavior of the relative variance as a function of $`\pi `$ for given $`R_g=0.19`$ is shown in Fig. 3. The value of $`R_g`$ is chosen for illustration and is the same as that obtained from current experiments on Pb-Pb collision at 158 $`A`$ GeV. For $`\pi `$ from 100 to 300 the relative variance decreases quite quickly from about 28% to about 16.5%. The relative variance given experimentally from the mixed events is about 23% at $`\pi `$=270.13 croland . Since the relative variance depends on $`R_g`$, different curve for the dependence of $`\sigma /R_e`$ on $`\pi `$ will be obtained if another value of $`R_g`$ is taken as input. For fixed $`\pi `$, the larger $`R_g`$, the smaller $`\sigma /R_e`$. The decrease tendency of $`\sigma /R_e`$ with $`\pi `$ will be similar for all $`R_g`$. The decrease of the relative width $`\sigma /R_e`$ with the increase of $`\pi `$ can be anticipated from the fact that both the distributions for kaons and pions are of Poisson in which the relative widths of the multiplicity distributions are inversely proportional to the root of the mean multiplicities. In Refs. gr97 ; fs99  the mean multiplicity of positively charged pions is 270.13, in the range discussed in this paper. The experimental point is given in Fig. 3 by a solid square. One can see that the theoretically calculated relative width of the $`K/\pi `$ ratio distribution from the Poissonian multiplicity distributions would be too small. The discrepancy of calculated $`\sigma /R_e`$ from the experimental result for the same $`\pi `$ shows the important role played by the decoherent effect in the emission of the pions and kaons. The decoherent effect will broaden the multiplicity distributions of kaons and pions, so that a larger relative width of the $`K/\pi `$ ratio distribution can be expected. Such thermal effect cannot be washed out in the so called mixed events.

As the second possibility for the multiplicity distribution with pure statistical fluctuations, one can work with the modified negative binomial (MNB) distribution mnb1 ; mnb2 . As a generalization of the negative binomial distribution, the thermal effect plays an important role in the MNB distribution. In this case, the generating function is

$$G_{\mathrm{MNB}}(x)=\left(\frac{1\mathrm{\Delta }(x1)}{1r(x1)}\right)^k,$$
(13)
with three parameters $`\mathrm{\Delta },r`$ and $`k`$. As shown in the first paper in mnb1 , $`k`$ can be interpreted as the maximum number of fireballs or clusters in some initial states, $`\mathrm{\Delta }`$ and $`r`$ can be related to the rates of birth and immigration processes. The mean multiplicity in the MNB is $`n=k(r\mathrm{\Delta })`$. The width squared of the multiplicity distribution is $`wnn^2n^2=n(1+r+\mathrm{\Delta })`$. With the MNB distributions for multiplicities of kaons and pions, the distribution of event-by-event $`K/\pi `$ ratio is shown to be also a Gaussian, similar to but broader than that given in Fig. 1. Our calculations show that the mean value and the width of $`K/\pi `$ ratio distribution have weak dependence on the choice of $`k`$. Thus, the parameter $`k`$ is chosen to be 200 in the following. In the calculation of the dependences of the mean value and width of the $`K/\pi `$ ratio distribution on the pion mean multiplicity, there are two free parameters $`w_\pi `$ and $`w_k=(k^2k^2)/k`$ which represents the relative width of the distribution of kaon multiplicity. Then, $`R`$ can be calculated from Eq. (9) with $`G_{\mathrm{MNB}}`$ in place of $`G_\mathrm{\Pi }`$. The width can also be calculated as

$$\frac{\sigma }{R_e}=\sqrt{\frac{\left(1+\frac{w_k}{R_g\pi }\right)_0^1𝑑x\frac{G_{\mathrm{MNB}}(x)G_{\mathrm{MNB}}(0)}{x}\mathrm{ln}\frac{1}{x}}{\left[_0^1𝑑x\frac{G_{\mathrm{MNB}}(x)G_{\mathrm{MNB}}(0)}{x}\right]^2}1},$$
(14)
We do not write $`w_k`$ in Eq. (12) as it is simply 1, because there the kaon multiplicity distribution is assumed to be a Poissonian. Data from NA49 experiments show that $`(N^2N^2)/N`$ is about 2.0 for all charged particles. Considering the fact that most (about 90%) charged particles are pions, one can fix approximately $`w_\pi =(\pi ^2\pi ^2)/\pi =2.0`$ for comparison with the experiment. Numerical calculations show that $`\sigma /R_e`$ has a quite weak dependence on $`w_\pi `$. From the MNB distribution one can get the two parameters in Eq. (13) for pion multiplicity distribution as $`r=(w_\pi 1+\pi /k)/2`$, $`\mathrm{\Delta }=(w_\pi 1\pi /k)/2`$ with $`k`$ preset to be 200. Dependence of $`\sigma /R_e`$ on the global ratio $`R_g`$ can be expected from Eq. (14) in the same way as from Eq. (12), i.e., the larger $`R_g`$, the smaller $`\sigma /R_e`$ for fixed $`\pi `$ and other parameters. It is shown that $`\sigma /R_e`$ increases considerably with $`w_k`$. As argued in Ref. jeon , the main contribution to the fluctuations of the $`K/\pi `$ ratio comes from the fluctuation of multiplicity of kaons which are the particle species with fewer mean multiplicity. Normally, the relative multiplicity fluctuation for such particle species is larger. So in following calculations, three values of $`w_K=2.0,2.5,3.0`$ are chosen to present the dependence of $`\sigma /R_e`$ on the width of kaon multiplicity distribution. Because the mean ratio is independent of $`w_K`$, the calculated $`R`$ are the same and close to 1.0 within 2.5% with the three choices of $`w_K`$, and are given also in Fig. 2. However, $`\sigma /R_e`$ increases considerably with $`w_K`$, as shown in Fig. 3. In those analyses in gr97 ; fs99 ; croland  $`\sigma /R_e`$ from the mixed events is about 23% with $`\pi `$=270.13. Then from Fig. 3 one can see that the calculation from the MNB distribution can give value of $`\sigma /R`$ consistent with the experimental result from the mixed events when $`w_k`$ is suitably chosen (less than 2.5).

However, the value of $`w_K`$ in real experiment is unknown yet and may be different from 2.5 significantly. Thus, the calculated result from pure statistical consideration cannot at present be compared directly with that from the mixed events out of experimental data. The possible difference between them may have deep physical implications. This may show that there exist some other physical effects in the event-by-event $`K/\pi `$ distribution, which might be the correlation between the productions of kaons and pions. Such effects should be studied theoretically and experimentally.

As a summary, the statistical influence on the $`K/\pi `$ ratio is studied in an event-by-event way. It is shown that the distribution of the ratio in the pure statistical case is Gaussian with mean and width depending on the mean multiplicities of kaons and pions. It is interesting to further study the effect on the ratio from correlations between the productions of kaons and pions and the influence of other dynamical contributions.

This work was supported in part by NNSF in China. C.B.Y would like to thank the Alexander von Humboldt Foundation for the Research Fellowship granted to him.


text:
# Towards the fabrication of phosphorus qubits for a silicon quantum computer

\[

## Abstract

The quest to build a quantum computer has been inspired by the recognition of the formidable computational power such a device could offer. In particular silicon-based proposals, using the nuclear or electron spin of dopants as qubits, are attractive due to the long spin relaxation times involved, their scalability, and the ease of integration with existing silicon technology. Fabrication of such devices however requires atomic scale manipulation - an immense technological challenge. We demonstrate that it is possible to fabricate an atomically-precise linear array of single phosphorus bearing molecules on a silicon surface with the required dimensions for the fabrication of a silicon-based quantum computer. We also discuss strategies for the encapsulation of these phosphorus atoms by subsequent silicon crystal growth.
(To appear in Phys. Rev. B Rapid Comm.)

\]

A quantum bit (or qubit) is a two level quantum system that is the building block of a quantum computer. To date the most advanced realisations of a quantum computer are qubit ion trap  and nuclear magnetic resonance  systems. However scaling these systems to large numbers of qubits will be difficult , making solid-state architectures , with their promise of scalability, important. In 1998 Kane proposed a novel solid state quantum computer design  using phosphorus 31P nuclei (nuclear spin I = 1/2) as the qubits in isotopically-pure silicon 28Si (I = 0). The device architecture is shown in Fig. 1a, with phosphorus qubits embedded in silicon approximately 20 nm apart. This separation allows the donor electron wavefunctions to overlap, whilst an insulating barrier isolates them from the surface control A and J gates. These A and J gates control the hyperfine interaction between the nuclear and electron spins and the coupling between adjacent donor electrons respectively. For a detailed description of the computer operation refer to Kane . An alternative strategy using the electron spins of the phosphorus donors as qubits has also been proposed .

One of the major challenges of this design is to reliably fabricate an atomically-precise array of phosphorus nuclei in silicon - a feat that has yet to be achieved in a semiconductor system. Whilst a scanning tunnelling microscope (STM) tip has been used for atomic scale arrangement of metal atoms on metal surfaces , rearrangement of individual atoms in a semiconductor system is not straightforward due to the strong covalent bonds involved. As a result, we have employed a hydrogen resist strategy outlined in Fig. 1b. Here the array is fabricated using a resist technology, much like in conventional lithography, where the resist is a layer of hydrogen atoms that terminate the silicon surface. An STM tip is used to selectively desorb individual hydrogen atoms, exposing the underlying silicon surface in the required array. STM induced hydrogen desorption has been developed and refined over the past ten years  and has been proposed  for the assembly of atomically-ordered device structures. We demonstrate a process to adsorb single phosphine molecules in a predefined array, with atomic resolution, which we have developed specifically for the fabrication of a silicon quantum computer. The incorporation of these arrays in silicon is then discussed.

Whilst hydrogen lithography at the tens of Angstrom linewidth scale has been used to selectively expose the silicon surface to oxygen , ammonia , iron , aluminium , gallium  and cobalt , there has only been one recent report of adsorption, at the atomic scale, in this case individual and clusters of silver atoms . Here we demonstrate the controlled adsorption of a linear array of single phosphine molecules in the extreme case of single hydrogen atom desorption sites for direct application to the fabrication of a scalable silicon quantum computer. This technical achievement has answered the critical questions of whether a hydrogen resist is effective during exposure to phosphine and whether or not phosphine will adsorb to an STM-exposed site sufficiently small to achieve one and only one PH3 molecule at that site.

The requirements for this quantum computer design are very stringent. In order to undertake high resolution lithography the silicon surface must be atomically flat with a low defect density to allow the formation of a near perfect resist layer, where one hydrogen atom bonds to each surface silicon atom. The ability to then desorb individual hydrogen atoms requires a sharp, large cone angle tungsten tip in order to form $`<`$1 nm desorption sites. These sites are subsequently exposed to high purity phosphine gas for the required phosphorus atom placement (Fig. 1b). We demonstrate each of these steps below. Throughout this process particular attention must be made to avoid the introduction of any spin or charge impurities that would be fatal to the operation of the quantum computer .

The final fabrication step outlined in Fig. 1b which is not the subject of this paper, but presents a significant challenge, is to encapsulate the phosphorus qubits in a crystalline lattice of isotopically-pure silicon-28 . The main difficultly in this step is to ensure the phosphorus atoms incorporate into the silicon crystal and remain in their ordered atomic array. The most direct route to achieving this is to desorb the hydrogen resist from the surface by heating to $``$700 K followed by epitaxial silicon growth over the phosphorus array. A possible concern with heating the surface is the potential to induce lateral surface diffusion of the phosphorus atoms in the array. However it is known that the phosphorus atom in the PH2 molecule, with a single bond to the silicon, incorporates into the silicon surface with a more stable threefold coordination geometry over the lower temperature range of 550-650 K . We can avoid heating the surface at all during hydrogen desorption by either direct optical excitation at $`\lambda `$ = 157 nm  or electron bombardment . An alternative room temperature approach that will help to incorporate the phosphorus into the silicon lattice without significant diffusion involves photo-induced excitation of the attached PH2 molecule to increase its chemical activity with the surface . Finally it is also possible to leave the hydrogen resist and grow epitaxial silicon directly on the mono-hydride surface . In all above cases 10-30 Å of epitaxial silicon can be grown at low temperatures, down to room temperature  to encapsulate the phosphorus, followed by elevated temperature growth at $``$500 K to maintain crystallinity . We will return to the issue of phosphorus incorporation in our final discussion.

The first step in fabricating a phosphorus qubit array is to prepare a clean, defect free silicon surface. Fig. 2a demonstrates an optimally-prepared Si(001)2$`\times `$1 surface consisting of rows of $`\sigma `$-bonded silicon-silicon dimers. By heating the surface to 1200 C in a variable-temperature ultra-high vacuum STM system and performing a controlled cool-down  we have achieved large defect free regions. The “bean” shaped protrusions in this filled state image correspond to the charge overlap of the electrons in the dangling bonds on each silicon surface atom giving rise to a weak $`\pi `$-bond. These dangling bonds make the surface reactive, allowing the subsequent adsorption of species such as hydrogen and phosphine.

The adsorption of atomic hydrogen on the Si(001)2$`\times `$1 surface can occur in three ways depending on the surface temperature, forming either a mono-, di- or tri-hydride . An atomic hydrogen source consisting of a heated tungsten filament and water-cooled heat shroud has been used for this work. Fig. 2b shows an image taken at room temperature of a low defect density Si(001)2$`\times `$1 surface after exposure to high purity atomic hydrogen at 600 K. It can clearly be seen that a near uniform coverage of the mono-hydride phase has occurred, where one hydrogen atom bonds to each silicon atom. Sites of di-hydride, where two hydrogen atoms bond to each silicon, and tri-hydride, a mixture of the other two phases, are also indicated. We have found that all three phases passivate the surface and act as an effective resist during subsequent exposure to phosphine. Comparison of current-voltage spectroscopy before and after hydrogen dosing (not shown) confirms the existence of a hydrogen passivated surface, with the silicon $`\pi `$\*-antibonding peak evolving into the silicon-hydrogen antibonding peak after hydrogenation .

An STM tip is then used to desorb hydrogen from the surface by application of a controlled voltage pulse between the tip and sample . Optimisation of the geometry of an oxide-free tungsten tip and controlled voltage pulses (sample bias $``$6 V and tunnelling current $``$1 nA for $``$1 ms) makes atomic resolution desorption possible. In order to allow the adsorption of one phosphine molecule, and therefore only one phosphorus atom, it is necessary to desorb an area that exposes less than or equal to two silicon dimers as shown schematically in Fig. 3a. This is because phosphine bonds to the Si(001)2$`\times `$1 surface with a c(4$`\times `$2) surface periodicity as demonstrated in Fig. 3b where we have dosed a clean Si(001)2$`\times `$1 surface with phosphine . The STM image in Fig. 3c shows three $`<`$ 1 nm diameter hydrogen desorption sites in a row with a pitch of $``$4 nm on a hydrogen terminated Si(001)2$`\times `$1 surface. This image, with such a close spacing between sites, highlights the atomic resolution desorption achieved. The distance between sites can easily be increased to the required qubit spacing of 20 nm, and we have performed controlled lithography of single desorption sites in a line $`>`$ 100 nm in length. The desorption sites in Fig. 3c appear as bright protrusions as a result of the extension of electron density out of the surface due to the silicon dimer surface states of the exposed silicon dangling bonds . Fig. 3 demonstrates that these desorption sites are sufficiently small to allow only one phosphine molecule to bond to the surface at each site.

Following STM lithography to expose small regions of the Si(001)2$`\times `$1 surface it is then necessary to adsorb the individual phosphorus atoms or qubits using phosphine gas. Figs. 4a and 4b show atomic resolution images of two hydrogen desorption sites $``$1.5 nm apart on a dimer row both before and after exposure to phosphine gas. The effectiveness of the hydrogen resist as a barrier to phosphine adsorption is demonstrated by the uniform hydrogen coverage after phosphine dosing except at the previously desorbed hydrogen sites. In order to observe any changes after phosphine exposure we have specifically created single hydrogen desorption sites, rather than larger desorption sites (as in Fig. 3) and present high resolution images where the spacing between sites is very small. Here the bright protrusion at each of the desorption sites in Fig. 4a is the signature of the single silicon dangling bond, after desorption of just one hydrogen atom, in this case from the left side of the dimer. The remaining hydrogen on the silicon dimer is known to be transient  and we have observed it diffusing from one side of the dimer to the other.

In order to obtain high purity phosphine gas delivery, the PH3 micro-dosing system and its connection to the UHV STM employed internally electro-polished gas lines assembled in a clean-room environment. Mass spectra taken in the chamber during the exposure at a pressure of 10-8 mbar reveal no significant increase in the partial pressure of any other species. The sticking coefficient of phosphine on the clean silicon surface is 1 . Fig. 4b shows the same area as Fig. 4a after exposure to phosphine gas at room temperature. Analysis of the line profiles in Fig. 4c and 4d shows a characteristic increase of $``$0.05 nm in the height of the protrusion after phosphine dosing . Figs. 4d and 4e show three desorption sites, in a line perpendicular to the dimer rows, before and after phosphine dosing. The associated line profiles (Figs. 4f and 4g) again show a $``$0.05 nm height increase. This increase, observed at all adsorption sites over several images, was calibrated against an atomic step edge on the same surface (not shown) both before and after phosphine dosing. This reproducible increase confirms the adsorption of a PH3 molecule and corresponds to the difference between the exposed silicon dangling bond and the adsorbed phosphine. The transient nature of the hydrogen atom on the silicon dimer can account for the asymmetry of the image in Fig. 4b where one phosphine molecule has bonded to the left silicon in the dimer (upper) and another phosphine molecule has bonded to the right silicon in the dimer (lower).

Studies of the interaction of phosphine with the clean Si(001)2$`\times `$1 surface  suggest that PH3 molecularly adsorbs to one end of a silicon dimer and can then dissociate to PH2 provided enough silicon dangling bonds are available nearby for the re-adsorption of the dissociated hydrogen. The absence of available dangling bond sites on the hydrogen-terminated surface inhibits this dissociation step. Fig. 5a shows a schematic of a single dangling bond site on a hydrogen terminated surface before phosphine dosing and Fig. 5b shows the proposed structure of Si-PH2 after dissociation of the adsorbed PH3. In this configuration the adsorbed phosphorus atom maintains the tetrahedral geometry and threefold coordination, with the PH2 bonded to one end of the dimer . A similar dissociative process and Si-XH2 configuration is understood to occur in the adsorption of NH3  and AsH3  on the Si(001)2$`\times `$1 surface, where both nitrogen and arsenic are isoelectronic with phosphorus.

In the above discussion we have demonstrated the process of controlled single phosphine molecule adsorption for the extreme case of single hydrogen atom desorption sites. In future work it will be important to maintain the ordered atomic array during the final silicon growth step (final panel of Fig. 1b). Here the optimisation of the size and crystallographic orientation of the desorption sites will be critical. Dissociation of singly-bonded PH2 to threefold coordinate P + 2H can proceed above 500 K , provided there are suitable dangling bond sites that overcome competition with the recombinative desorption of PH3. An alternative option is to dose the hydrogen-patterned surface with phosphine gas at elevated temperatures ($``$550 K) to directly induce phosphorus incorporation into the first atomic layer. We have independently performed elevated temperature studies which indicate that at these temperatures, the silicon surface is stable.

We have demonstrated the effectiveness of a hydrogen resist as a barrier to phosphine adsorption and used STM lithography for the controlled placement of single phosphorus bearing molecules on a silicon surface - central to the construction of a silicon based quantum computer. This process, shown for closely spaced controlled doping, demonstrates the achievability of more widely-spaced ($``$20 nm), precisely positioned phosphorus qubit arrays over large areas. Whilst it is not yet possible to guarantee the subsequent overgrowth steps required for the Kane architecture, this letter demonstrates a significant advance in the solid-state qubit fabrication process and bodes well for the realisation of a scalable silicon quantum computer. Finally, it is important to note that the fabrication strategy demonstrated here is also directly applicable to other silicon based quantum computer architectures .

We would like to thank I. Andrienko, T. Kamins, G. J. Milburn, S. Prawer, C. Thirstrup and S. Williams for useful discussions. This work is funded by the Australian Research Council Special Research Centre scheme, the Australian Government, the US National Security Agency, and the US Advanced Research and Development Activity.