text:
# 1 Introduction: why neutrinos are interesting?

## 1 Introduction: why neutrinos are interesting?

At the present time neutrinos are the only elementary particles that indicate new physics beyond the Minimal Standard Model (MSM). Otherwise all experimental data are in perfect agreement with MSM. It is well established that all constituents of matter are three families (or flavors) of quarks and leptons. Each family of quarks may have one of three colors, while leptons are white or colorless. Both quarks and leptons participate in electroweak interactions which are realized by the exchange of massless photons ($`\gamma `$), coupled to electric charge of particles, and heavy $`W^\pm `$ and $`Z^0`$ bosons. Since neutrinos do not have any electric charge, they interact only with W and Z. That’s why their interactions are so weak at low energies ($`E<M_{W,Z}`$). Quarks, due to their color charges, interact in addition with 8 colored gluons. These interactions are described by quantum chromodynamics (QCD). This theory well describes all observed phenomena in particle physics and only neutrinos disturb the peace. On the other hand, astronomy and cosmology strongly demand new physics beyond MSM.

Neutrinos play a very important role in astrophysics. Stellar evolution and e.g. supernova explosions would be drastically different without neutrinos. Since heavy elements, necessary for life, are produced in SN explosions, one may say that life itself strongly depends upon existence of neutrinos. Due to their large penetrating ability neutrinos comes out from the deep stellar interior and could give information about processes in stellar cores. Neutrino telescopes will open (in fact have already opened) a new window in astronomy. Observation of solar neutrinos created serious puzzles that most possibly could be resolved by new phenomena in particle physics (neutrino oscillations) or by something unexpected (really drastically) in solar astrophysics.

In cosmology, neutrinos have strong impact on Big Bang Nucleosynthesis (BBN) and thus determine chemistry of the universe. They could even create chemically inhomogeneous, though energetically smooth, universe. Neutrinos may participate in formation of the large scale structure of the universe, making hot, warm and possibly self-interacting dark matter. Detailed mapping of the large scale structure will permit to measure neutrino mass with the accuracy better than 1 eV. Massive or massless, neutrinos have a noticeable impact on the angular spectrum of cosmic microwave background radiation (CMB) and the future missions, MAP and especially Planck, will determine the number of light species in the universe and measure their masses.

## 2 What do we know about neutrinos.

In the standard model neutrinos possess only weak (and of course gravitational) interactions. Anomalous, even stronger than weak, interactions are not excluded but only between neutrinos and/or some other new not yet discovered particles. Neutrinos are coupled to the intermediate bosons according to: $`W^+\nu _ll^+`$ ($`l=e,\mu ,\tau `$) and $`Z^0\nu _l\overline{\nu }_l`$. Measurement of the decay width of $`Z^0`$ permitted to conclude that the total number of different neutrino species is $`N_\nu =3.07\pm 0.12`$, while the combined fit to LEP data gives the more accurate result, $`N_\nu =2.994\pm 0.01`$ (for this and other limits below see Particle Data Group, 1998).

Electric charge of neutrinos are believed to be exactly zero, while magnetic moment should be non-vanishing. Direct measurements give the upper bounds: $`\mu _{\nu _e}<1.810^{10}\mu _B`$, $`\mu _{\nu _\mu }<710^{10}\mu _B`$, $`\mu _{\nu _e}<610^7\mu _B`$, where $`\mu _B`$ is the Bohr magneton. Studies of the stellar cooling rate permit to obtain much stronger limit: $`\mu _\nu <310^{12}\mu _B`$ (see e.g. Raffelt, 1996). These results are are about 10 orders of magnitude above the standard theoretical expectations, but in some extensions of MSM $`\mu _\nu `$ might be large.

Neutrinos may be massless, but there is no theoretical principle that requests the vanishing of $`m_\nu `$, so it is natural to expect that $`m_\nu 0`$. All direct experiments are compatible with $`m_\nu =0`$ but with a quite different level of precision: $`m_{\nu _e}<3`$ eV, $`m_{\nu _\mu }<160`$ keV, and $`m_{\nu _\tau }<18`$ MeV. Indication for $`m_\nu 0`$ comes from possible observations of neutrino oscillations (see below). Two different types of mass are possible for neutrinos: the Dirac mass, that distinguishes particles from anti-particles (such a mass have e.g. electrons), and Majorana mass, that mixes neutrinos and antineutrinos. Majorana mass would induce non-conservation of leptonic charge. From the absence of the double beta-decay of heavy nuclei one can conclude that $`m_M<0.1`$ eV.

Leptonic charge is conserved in all observed up to now reactions, however it is not excluded that it may be strongly non-conserved in neutrino oscillations. The hypothesis that neutrinos might oscillate was proposed in 1957 by B. Pontecorvo and its verification now presents a major challenge in experimental particle physics. Oscillation is a generic phenomenon if neutrinos are massive. Indeed, mass eigenstates normally do not coincide with the interaction eigenstates because the mechanism of mass generation does not know anything about gauge interactions with $`W`$ and $`Z`$. Thus for example the electronic neutrino, i.e. the one that coupled to electrons through $`W\nu _ee`$ could be a superposition of two mass eigenstates:

$`|\nu _e=\mathrm{cos}\theta |\nu _1+\mathrm{sin}\theta |\nu _2`$ (1)
and e.g. $`\nu _\mu `$ is the orthogonal combination of $`\nu _{1,2}`$. The mixing may be more complicated and include all three active neutrinos, $`\nu _e,\nu _\mu `$, and $`\nu _\tau `$. Moreover new sterile neutrinos, $`\nu _s`$, may be involved. If all existing data on neutrino anomalies are correct (see below) and the anomalies are explained by the oscillations, sterile neutrinos seem to be necessary.

A very interesting effect may take place in neutrino oscillations in inhomogeneous or non-stationary matter, e.g. in stellar interior or in primeval cosmic plasma. Though neutrinos very weakly interact with matter, their refraction index is still different from unity. And if the mass difference is sufficiently small, the matter effects could be significant and in particular a resonance transition between neutrino flavors could be possible (Miheev and Smirnov, 1985; Wolfenstein, 1978). In this case, even for a very small mixing angle in vacuum, the mixing in matter could be of order unity.

## 3 Neutrino anomalies

The following phenomena in neutrino physics do not fit the MSM expectations:

1. Deficit of solar neutrinos: the measured flux of solar neutrinos is approximately twice smaller than theoretical predictions. If this deficit is explained by the oscillations then the following values of oscillation parameters are possible: 1) $`\delta m^2510^6`$ eV2, $`\mathrm{sin}^22\theta 610^3`$; this is the MSW resonance solution; the mixing could be either with an active neutrino or a sterile one; 2)$`\delta m^2510^5`$ eV2, $`\mathrm{sin}^22\theta 0.75`$; 3) “just so” vacuum solution: $`\delta m^2810^{11}`$ eV2, $`\mathrm{sin}^22\theta 0.75`$ (for the review see e.g. Bahcall et al, 1998). Long base-line experiments, in particular, CERN neutrinos registered in Gran Sasso, could fix the parameters.

Deficit of atmospheric $`\nu _\mu `$:  the flux of $`\nu _\mu `$ observed in cosmic rays is twice smaller than predicted. The recently measured angular dependence (SuperKamiokande, 1998) gives a strong argument in favor of the mixing of $`\nu _\mu `$ with $`\nu _\tau `$ or $`\nu _s`$ with $`\delta m^2=2.210^3`$ eV2 and $`\mathrm{sin}2\theta 1`$.

Direct observation of neutrino oscillations (LSND, 1998): $`\nu _\mu `$ produced by the decay $`\pi \mu +\nu _\mu `$ induce the reactions $`\nu +Ce^{}+X`$, i.e. electronic neutrinos were generated by muonic ones: $`\nu _\mu \nu _e`$.

KARMEN anomaly (1995, 1999): neutrinos produced by pion decay at rest were registered at a certain distance behind a shield. The moment of the pion production was well fixed. The time distribution of the events should follow the pion decay exponent. However a considerable excess of the events was observed. This excess can be explained by a production of a heavy ($`m=33.9`$ MeV) sterile neutrino that slowly propagated to the detector and decayed there (Barger et al, 1995). This interpretation however meets serious problems with cosmology and astrophysics (Dolgov et al, 2000a).

## 4 Cosmological limits on neutrino mass

Cosmology permits to put strong limits on neutrino mass. The first limit was derived by Gerstein and Zeldovich in 1966. The number density of relic neutrinos at the present time can be calculated if they were in thermal equilibrium in the early universe: $`n_{\nu _a}=n_{\overline{\nu }_a}=3n_\gamma /22=56/\mathrm{cm}^3.`$ This result is true for any stable neutrino with normal weak interactions and the mass below a few MeV. From this expression follows

$`{\displaystyle \underset{a}{}}m_{\nu _a}=94\mathrm{eV}\mathrm{\Omega }_\nu h_{100}^2`$ (2)
where $`h_{100}=H/100`$ km/sec/Mpc $`0.65`$ and $`\mathrm{\Omega }_\nu =\rho _\nu /\rho _c`$. Definitely $`\mathrm{\Omega }_\nu <1`$ and thus $`m_\nu <40`$ eV. However, neutrinos cannot dominate the energy density of the universe, otherwise the structure formation at large $`z`$ (at early time) would be suppressed. For a flat $`\mathrm{\Lambda }`$ \+ CDM model with $`\mathrm{\Omega }_m=0.3`$ the sum of neutrino masses is bounded by $`m_\nu <2`$ eV (Fukugita et al, 2000). The future sensitivity of detailed galaxy surveys and of CMB angular spectrum is expected to be sensitive to $`m_\nu 1`$ eV.

Very heavy neutrinos with mass above 45 GeV are not excluded by LEP measurements and in a large mass interval they are permitted by cosmology. Heavy neutrinos more efficiently annihilate in the early universe. Their annihilation cross-section is proportional to $`m_\nu ^2`$ for $`m_\nu <m_{W,Z}`$, while at higher energies it is inversely proportional to $`m_\nu ^2`$. Correspondingly neutrinos with mass in the interval from 2 GeV up to approximately 100 TeV are allowed by cosmology. There are no reliable calculations of the cosmological mass limits at the heavy end because such a large mass of neutrinos corresponds to the strong coupling regime of the electroweak theory.

If neutrinos constitute all dark matter in galaxies then their mass should be larger than 100-200 eV (Tremain and Gunn, 1979). This limit is an impressive manifestation of quantum mechanics at astronomical scales. Since neutrinos are fermions, they obey the Fermi exclusion principle and hence one cannot pack too many neutrinos into a galaxy. Correspondingly their mass should be bounded from below to account for all dark matter.

## 5 Neutrinos and big bang nucleosynthesis

There are three physical effects through which neutrinos can effect the creation of light elements:

Cooling rate of the universe. Comparing the cosmological critical energy density, $`\rho _c=3m_{pl}^2/32\pi t^2`$, with the equilibrium energy density of thermal plasma, $`\rho _T=\pi ^2g_{}T^4/30`$, one can find the law of the temperature evolution, $`tg_{}^{1/2}T^2`$. Here $`g_{}`$ is the number of particle degrees of freedom in the plasma. During primordial nucleosynthesis $`g_{}=10.75`$. Each neutrino plus antineutrino species contributes into this number 7/4. A change in the number of particle species would change the temperature of neutron-proton freezing, $`T_fg_{}^{1/6}`$, and also change the moment of the onset of nucleosynthesis. The latter takes place when $`T65`$ keV but the time when such temperature is reached depends upon $`g_{}`$. Correspondingly the number of neutrons that survived the decay till beginning of light element formation would be different. Combination of these two effects results in a change of mass fraction of $`{}_{}{}^{4}He`$ by 4% if the number of neutrino species is changed by 1. This phenomenon was first noticed by Hoyle and Taylor in 1964 and later by Peebles (1966). Detailed calculations were done by Shvartsman (1969) and Steigman et al., (1977).

According to the recent analysis (Tytler et al 2000) the data permits not more than 0.2 extra neutrino species at BBN. However the conclusions of the earlier paper (Lisi et al, 1999) is that one should be rather cautious in interpretation of the data and a safer limit is $`\mathrm{\Delta }N_\nu <1`$. Moreover there is a conflicting evidence of possible high and low primordial deuterium (see the next section).

Lepton asymmetry. The standard calculations of primordial abundances of light elements are done under assumption of vanishingly small lepton asymmetry. However the cosmological asymmetry in lepton sector is unknown and comparison of the BBN calculations with observations permits to put some limits on the value of the asymmetry. If one parameterizes neutrino distribution by the equilibrium function with a non-zero chemical potential, $`f_\nu =[1+\mathrm{exp}(E/T\xi )]^1`$ the following upper limits on values of $`\xi `$ can be established: $`\xi _{\nu _e}<0.03`$ and $`\xi _{\nu _\mu ,\nu _\tau }<3`$. The results are more sensitive to the asymmetry in the sector of electronic neutrinos because the latter directly influence the $`n/p`$-transformation through the reactions $`n+\nu _ep+e^{}`$ and $`n+e^+\overline{\nu }_e+p`$, while $`\nu _\mu `$ and $`\nu _\tau `$ effect the cooling rate only.

Neutrino spectrum. It is assumed normally that $`\nu _e`$ participating in $`np`$-reactions have thermal spectrum. However even in the standard model their spectrum differs from the thermal one by roughly 1% (Dolgov and Fukugita, 1992a; Dodelson and Turner, 1992; Hannestad and Madsen, 1995; Dolgov et al, 1997). The impact of this spectral distortion on the light element abundances is very small, at the level of $`10^4`$, but in some modification of the standard model (in particular if there are oscillations between $`\nu _e`$ and $`\nu _s`$) the effect may be significant.

## 6 Possible spatial variation of primordial abundances

The standard cosmological model satisfies cosmological principle stating that the universe is the same everywhere. A strong argument in favor of this hypothesis is that the universe is very smooth energetically as is shown by the almost perfect isotropy of of CMB. However it is possible to formulate a model that gives a strong variation of the primordial chemical content of the universe at different space points and very low density variation (Dolgov and Pagel, 1999). This work was stimulated by the observation of primordial deuterium at large red-shifts, $`z=0.53`$. Different groups observed three different values of the number density of primordial deuterium Most observations give “normal” abundance, $`D/H=(34)10^5`$ (see Tytler et al, 2000). This value is observed also in our neighborhood. There are a few observations of high deuterium, $`D/H=(1020)10^5`$ (Webb et al, 1997, and references therein) and there are even recent observation of the region with a low deuterium, $`D/H=(12)10^5`$ (Molaro et al, 1999).

If the effect is real its importance for the cosmology is difficult to overestimate. Since the characteristic scale of variation is very large (more than several hundred Mpc), the proper conditions should be arranged during inflationary stage. A model of leptogenesis, that can explain the observed variation (Dolgov and Kirilova, 1991; Dolgov, 1992b), is based on the Affleck and Dine scenario (1985) that gives a large lepton asymmetry together with a small baryon asymmetry. Moreover, the cosmological lepton asymmetry in this model could be non-uniform with a characteristic scale of variation bigger than several 100 Mpc. To escape large angular fluctuation of CMB temperature one has to assume that there exists “lepton conspiracy” (Dolgov and Pagel, 1999), ensuring a symmetry with respect to permutation of different chemical potentials (electronic, muonic, and tauonic).

The model predicts that there should be three regions on the sky: 2/3 are normal with the standard abundances; 1/6 has low deuterium and the other 1/6 has high deuterium. In the regions with normal deuterium the mass fraction of primordial $`{}_{}{}^{4}He`$ is also normal, $`Y_p25\%`$. In deuterium poor regions $`Y_p12\%`$, and in deuterium rich regions there is a lot of helium, $`R50\%`$. Surprisingly, no data at the present time exclude such unusual regions situated at the distances of several hundred Mpc (or larger).

Independently on the model and the data on primordial deuterium, it is an important question: what are the limits on the mass fraction of the second most abundant element in the universe, $`{}_{}{}^{4}He`$, at large distances? A possible indication for a high or low fraction of $`{}_{}{}^{4}He`$ could be a different stellar evolution. However stars are not resolved at large distances. One can observe only galaxies and it is difficult to say, if one sees an unusually blue or red galaxy, whether the effect is attributed to an abnormal helium content or to galaxy evolution. Detailed calculations of stellar evolution with high and low initial helium-4 are necessary. For more discussion see the papers by Dolgov and Pagel (1999) and by Dolgov (1999).

Possibly more promising is a search for the variation of helium through the angular spectrum of CMB. The model predicts some peculiar features in CMB angular fluctuations both at high and low angular scales. Especially interesting is the decay of the fluctuations related to diffusion (or Silk) damping (Silk, 1968). The diffusion rate depends upon the number density of electrons prior to hydrogen recombination. In the region with high $`{}_{}{}^{4}He`$ the number of electrons is smaller because helium recombined earlier and took some electrons from the plasma. Thus, the damping in such regions would be stronger. In the regions with low helium the effect is the opposite (Hu et al, 1995). According to the calculations made by P. Naselsky, RATAN-600 is sensitive enough to observe the effect at $`l\mathrm{a}\mathrm{few}10^3`$. Planck could observe the effect if additional information about the position of poor or rich regions is known. The effect is rather striking: the slope of the angular spectrum of CMB would be very different at different directions on the sky.

## 7 Neutrino oscillations in the early universe and BBN

Neutrino oscillations in the early universe differ from those in stellar interior by two important features: 1) neutrinos may change the medium by back reaction from the oscillations; 2) the loss of coherence is essential and the density matrix formalism should be employed (Dolgov, 1981; Sigl and Raffelt, 1993).

Neutrino refraction index $`n`$ in the cosmological plasma was calculated by Nötzold and Raffelt in 1988. There are two types of contribution into $`n`$. The first one comes from the averaging of the neutrino interaction current $`J_\alpha `$. Due to homogeneity of the universe only the time component of the current gives non-zero contribution proportional to the charge asymmetry of the plasma. The second term comes from non-locality of neutrino interactions due to exchange of intermediate bosons. If charge asymmetry has the normal value $`10^{10}`$, the first term is subdominant. However in the case of resonance transition the effect of asymmetry can be strongly amplified. Indeed, the resonance may first take place for neutrinos so that more neutrinos would be transformed into sterile ones by oscillations. This would give rise to an increase of the asymmetry in the active neutrino sector and in turn would enhance the oscillations. In this way an exponential instability with respect to generation of lepton asymmetry is developed. The effect was found by Barbieri and Dolgov in 1991 but the authors concluded that the back reaction from the plasma quickly turn off the instability. This conclusion was reconsidered by Foot et al (1996) who argued that the back reaction is not so important and that the asymmetry could rise almost to unity. However recent semi-analytical calculations by Dolgov et al (1999c) found a much milder enhancement, by 5 orders of magnitude only. At the moment different groups are in disagreement with each other and the issue remains unresolved. The discussion and the list of relevant papers can be found in Dolgov et al (1999c), Dolgov (2000b), and Di Bari and Foot (2000).

The situation is somewhat simpler in the non-resonance case but still there is a disagreement between the first paper by Barbieri and Dolgov (1990) and all the subsequent ones. While in the first paper it was assumed that the probability of sterile neutrino production is proportional to the inverse annihilation rate of active neutrinos, in the other ones it was argued that the breaking of coherence was determined by the total reaction rate (which is approximately 10 time larger) and this rate must determine the production of $`\nu _s`$. The problem was reconsidered recently by Dolgov (2000) and it was shown that the annihilation rate plays a decisive role in production of $`\nu _s`$ in accordance with conjecture of Dolgov and Barbieri (1990).

A resolution of this issues is essential for derivation of BBN bounds on oscillation parameters. According to Dolgov (2000c) the limits are:

$`\delta m^2\mathrm{sin}^42\theta |_{\nu _\mu ,\nu _\tau }<3.310^4\mathrm{\Delta }N_\nu ^2,`$ (3)
$`\delta m^2\mathrm{sin}^42\theta |_{\nu _e}<510^4\mathrm{\Delta }N_\nu ^2.`$ (4)
The last result should be somewhat stronger because electronic neutrinos effect the BBN not only producing a sterile partner but also through their spectrum. Numerical calculations of the last effect were done by Kirilova and Chizhov (1997, 1998a,b) for small $`\delta m^210^7`$ eV2. As is argued by Dolgov (2000c) for a larger $`\delta m^2`$ the $`\nu _e`$ spectrum remains of thermal form but with a non-zero chemical potential. Its impact on $`n/p`$-ratio permits to exclude mixing of $`\nu _e`$ with $`\nu _s`$ to the region:

$`\mathrm{sin}^22\theta <0.32\mathrm{\Delta }N_\nu `$ (5)
for all $`\delta m^2>510^6`$ eV2. BBN limits may essentially reduce the permitted parameter space for the oscillations and more work is necessary to resolve the existing controversies.


format: text

text:
# Collapse of a randomly forced particle: a comment

## Abstract

We refute the arguments by Anton in cond-mat/0004390, which set out to disprove the existence of a collapse transition for a randomly forced inelastic particle.

An eprint by Anton  has recently appeared on this archive, criticizing our proposal of a collapse transition for an inelastic, randomly-forced particle . Our conclusions were that a particle forced by Gaussian white noise that rebounds from a wall with coefficient of restitution $`r`$ will, with probability 1, dissipate all its energy and come to rest at the wall after an infinite number of collisions in a finite time, provided $`r<r_c=e^{\pi /\sqrt{3}}`$. It has been pointed out that the transition is not present for certain discretizations of the Langevin equation , and it is also not clear how it manifests itself in experimentally realizable systems . However, ref.  goes further, and argues that the transition is also absent for the ideal case of a pure white-noise force in continuous time. The purpose of this Comment is to point out substantial errors in Anton’s arguments, which, we believe, invalidate his conclusions.

The existence of the collapse transition is supported by analytical calculations  by Burkhardt, Franklin, and Gawronski (henceforth BFG), who constructed a steady-state solution of the Fokker-Planck equation for a particle confined in a finite spatial interval. They observed that the solution is well-behaved when the coefficient of restitution $`r`$ is in the range $`r_c<r<1,`$ where $`r_c`$ is the critical value for collapse proposed in , whereas the rate of collisions of the particle with the boundary diverges as $`r`$ approaches $`r_c`$ from above. Anton’s criticisms of this work, together with our responses, are as follows:

* *Burkhardt et al.’s probability density function is not normalizable*. This is incorrect. BFG state explicitly that their solution is normalizable for $`r_c<r<1`$, and it is in fact the onset of non-normalizability as $`rr_c^+`$ that lead them to support our conclusions.
* *The system does not approach a steady state, since the moments of a time-dependent solution are divergent in time*. Anton bases this on two calculations. The first shows that the energy of an *elastic* particle ($`r=1`$) diverges linearly in time, a result which is also stated in ref. . The second is a calculation (incorrect, as we later shall argue) of spatial moments for an inelastic particle ($`r<1`$) in a *semi-infinite* region, which diverge in time. Neither of these calculations are relevant to the case of an *inelastic* particle in a *finite* region, as studied by BFG. On the other hand, physical considerations show that these moments cannot diverge in that case. Firstly, the spatial moments are restricted due to the spatial confinement of the particle. Secondly, the energy is prevented from diverging by the fact that, at high speed $`v`$, energy is dissipated at a rate $`(1r^2)v^3`$ (collision rate $`v`$, energy dissipation per collision $`=(1r^2)v^2`$), whereas energy is being supplied at a constant rate .
* *The boundary conditions used by BFG lead to a contradiction*. The argument is as follows. The boundary condition for the probability density function $`G(x,v)`$ used by BFG, which corresponds to the stipulation that particles incident on the boundary at velocity $`v`$ rebound at $`rv`$, is

$$vG(0,v)dv=(rv)G(0,rv)d(rv).$$
(1)
Setting $`v=1`$, and then replacing $`rw`$ leads to

$$G(0,w)=\frac{1}{w^2}G(0,1).$$
(2)
However, contrary to Anton’s claims, this does not imply that the flux $`vG(0,v)`$ has a non-integrable singularity at $`v=0`$. The reason is that the solution must depend upon the coefficient of restitution—i.e., $`G`$ has a parametric dependence upon $`r`$, which has been suppressed in the above notation. This means that, for a given value of $`r`$, eqn. (2) is only valid for a single value ($`=r`$) of $`w`$. If we include explicitly the $`r`$-dependence in $`G(x,v;r)`$, eqn (2) becomes

$$G(0,w;w)=\frac{1}{w^2}G(0,1;w),$$

which does not imply that there is a non-integrable singularity in $`vG(0,v,r)`$ at constant $`r`$.

The final point leads Anton to propose alternative boundary conditions, which contradict eqn. (1) above. However, we have shown that there is no reason for believing eqn. (1) to be wrong, whereas Anton’s alternative is not justified on physical grounds. We therefore believe that he has not constructed a valid solution for this system, and the ensuing conclusions in sections III and IV of ref.  are erroneous.

Ref.  contains a further, independent argument against the existence of collapse. Anton states that, when collisions with the boundary occur, the Langevin equation for a damped particle near an elastic wall reduces to the equation of motion for the inelastic, undamped particle. Comparing the discretized versions of these equations at the instant of collision with the boundary , we find

$`\text{Damped,\hspace{0.17em} elastic:}v(t+\mathrm{\Delta }t)`$ $`=`$ $`v(t)\mathrm{\Delta }t\left\{\eta _1(t)(1r)v(t)\right\}`$
$`\text{Undamped,\hspace{0.17em} inelastic:}v(t+\mathrm{\Delta }t)`$ $`=`$ $`rv(t)\mathrm{\Delta }tr\eta (t)`$
These are far from equivalent—the inelastic particle loses a finite fraction of its energy at the instant of collision, whereas the damped particle loses an infinitesimal fraction of its energy. The equations become the same if one sets $`\mathrm{\Delta }t=1`$, but this does not allow for a systematic approach to the continuum limit since the Langevin equation has already been adimensionalized (i.e., the timescale has already been set in eqn. (6) of ref. ). In any case, it is already known  that inelastic collapse is destroyed by such crude discretization schemes.

To summarize, we do not believe that Anton has provided a substantive argument against the existence of the collapse transition. However, we agree that further work is needed to investigate whether, in the ideal system, inelastic collapse is definitive (i.e., is a particle that comes to rest at the wall capable of escaping again in a finite time?), and also what remnant of the collapse transition occurs in real systems.

We would like to thank Lucian Anton for interesting discussions.


text:
# Quantum Suppression of the Rayleigh Instability in Nanowires

## I Introduction

Plateau’s celebrated study  of the stability of bodies under the influence of surface tension established a fundamental result of classical continuum mechanics: a cylinder longer than its circumference is energetically unstable to breakup. Here we consider a quantum mechanical generalization of Plateau’s study—the stability of a metallic nanowire. Methods from the study of quantum chaos  are used to examine the energetics of a free electron gas in a cylindrical jellium filament, and show that an infinite filament can be stable against all axisymmetric perturbations if its electrical conductance $`G`$ \[in units of $`G_0=2e^2/h=(12,906.4\text{Ohm})^1`$\] belongs to a set of “magic numbers” $`n=1,3,5,6,\mathrm{}`$ and is otherwise unstable. Here $`e`$ is the charge of an electron and $`h`$ is Planck’s constant. Such magic numbers are analogous to those found in the shell model of atomic nuclei  and in metal clusters . Our stability analysis elucidates and confirms the stability properties of simple metal nanowires conjectured by Yanson, Yanson, and van Ruitenbeek  based on their observation of shell structure in the conductance statistics of sodium nanowires.

Structural multistability of metallic nanowires was previously postulated based on classical molecular dynamics simulations . However, the magic numbers observed in conductance histograms in alkali metals  are clearly an electronic shell effect, as shown below, and cannot be explained with classical molecular dynamics. The common occurence of multistability in two such radically different models likely stems from the fact that both models introduce an additional length scale (the Fermi wavelength $`\lambda _F`$ in the free electron model, the atomic diameter in molecular dynamics simulations), which leads to commensurability effects.

The properties of nanowires formed from simple monovalent metals, like the bulk properties of these materials , are determined largely by the delocalized conduction electrons. A free electron jellium model, treating the electrons as a noninteracting Fermi gas confined within the wire by hard-wall boundary conditions, provides an intuitive and even quantitative explanation of observed quantities like conductance , force , and shot noise . To examine energetic stability, we consider a weakly deformed cylinder, and find its thermodynamic potential to quadratic order in the amplitude of the deformation. A contribution proportional to the filament area, naturally identified as the surface tension, is only one of several competing terms in the free energy, which consists of terms which vary smoothly with the geometry of the filament and an oscillatory contribution directly connected with the discrete energy levels that are solutions to Schrödinger’s equation. The smooth terms appear in descending powers of length (proportional to volume, surface area, mean curvature, etc.), and are quite analogous to those found in the study of classical wave equations in curved domains  and the related problem of classically screened Coulomb (Debye-Hückel) interactions of curved surfaces in an electrolyte . The oscillatory part of the free energy, in particular, alters Plateau’s classical stability analysis in an essential way. The present work is thus fundamentally distinct from that which has dealt with the quantum mechanical origin of surface tension itself in metallic fluids , as well as those which consider quantum mechanical corrections to classical surface tension due to the quantization of capillary waves .

Long gold nanowires suspended between gold electrodes have been produced and imaged with a transmission electron microscope by Kondo and Takayanagi ; in particular, they have observed wires which are stable and almost perfectly cylindrical. The wires of Ref.  are of finite length, and attached to electrodes at either end. It seems reasonable to consider an infinite wire as a theoretical starting point, and to assume that the length of the wire acts as a cut-off: deformations with a wavelength longer than the wire length cannot occur. Thus one does not have to model the boundary conditions at the ends of the wire, which would introduce additional parameters.

To be specific, we examine an infinite cylindrical wire and its sole classically unstable deformation—an axially symmetric one (see Fig. 1). Any such deformation can be written as a Fourier series:

$$R(z)=R_0\left(1+_0^{\mathrm{}}𝑑qb(q)\mathrm{cos}(qz+\varphi (q))\right),$$
(1)
where $`R(z)`$ is the radius at position $`z`$, $`R_0`$ is the unperturbed radius, $`b(q)`$ is the (infinitesimal) perturbation coefficient, and $`\varphi (q)`$ is an arbitrary phase shift. The coefficients $`b(q)`$ are chosen such that the total volume of the wire is unchanged by the deformation. Other physically reasonable constraints  are also possible, but lead to similar results.

The metallic nanowire is an open system, connected to macroscopic metallic electrodes at each end . Therefore the change of the grand canonical potential $`\mathrm{\Omega }`$ under the perturbation determines its stability. $`\mathrm{\Omega }`$ is related to the electronic density of states $`g(E)`$ by

$$\mathrm{\Omega }=k_BT_0^{\mathrm{}}𝑑Eg(E)\mathrm{ln}\left(1+\mathrm{e}^{(E\mu )/k_BT}\right),$$
(2)
where $`k_B`$ is Boltzmann’s constant, $`T`$ is the temperature, and $`\mu `$ is the electronic chemical potential specified by the macroscopic electrodes. Our aim is to expand $`\mathrm{\Omega }`$ up to second order in the coefficients $`b(q)`$ characterizing the deformation. As we will show, this yields

$$\mathrm{\Omega }[b]=\mathrm{\Omega }[0]+_0^{\mathrm{}}𝑑q\alpha (q)[b(q)]^2+𝒪(b^3).$$
(3)
The change in the grand canonical potential is of second order in $`b`$ and contributions from deformations with different $`q`$ decouple. If the prefactor $`\alpha (q)`$ is negative for any value of $`q`$, then $`\mathrm{\Omega }`$ decreases under the deformation and the wire is unstable.

For an infinite cylindrical wire, the transverse motion is quantized, leading to the formation of discrete electronic subbands. The total density of states is the sum of the contributions from each subband (see Fig. 2): every subband begins to contribute at a threshold energy equal to the energy of its quantized transverse motion with a sharp spike, falling off smoothly for increasing energy. If the Fermi energy $`E_F`$ lies near one of these sharp peaks, certain small deformations of the wire can dramatically increase the density of states. According to (2), this lowers the grand canonical potential, leading to an instability. On the other hand, if there is no subband threshold sufficiently close to $`E_F`$, the density of states will instead decrease with any deformation, implying the existence of stable regions intervening between the instabilities associated with the opening of each subband.

## II Quantum Chaos Approach

In order to examine this picture quantitatively, we use a semiclassical approach , which enables $`g(E)`$ to be split into a smooth average contribution $`\overline{g}(E)`$, referred to as the Weyl contribution, and an oscillatory part $`\delta g(E)`$, whose average value is zero:

$$g(E)=\overline{g}(E)+\delta g(E).$$
(4)
The Weyl term $`\overline{g}`$ contains terms proportional to the volume of the nanowire, and to the area and curvature of its surface: $`\overline{g}(E)=(1/E_F)\stackrel{~}{g}(E/E_F)`$, where

$$\stackrel{~}{g}(x)=\frac{x^{1/2}}{2\pi ^2}k_F^3V\frac{1}{8\pi }k_F^2S+\frac{x^{1/2}}{6\pi ^2}k_FK,$$
(5)
and $`k_F=2\pi /\lambda _F`$ is the Fermi wavevector. The volume $`V`$, surface area $`S`$, and integrated mean curvature $`K`$ can be calculated for arbitrary perturbations by simple geometric considerations. The oscillatory part $`\delta g`$ of the density of states is a quantum correction, and may be calculated in the semiclassical approximation as a sum over all periodic classical orbits of the system .

Let us first consider an undeformed cylindrical wire. Each periodic orbit lies in a plane transverse to the wire’s axis, and there is a correspondence to the orbits  in a circular billiard (see Fig. 3). Note that the action $`S_{vw}`$ of each orbit is invariant under both translations of the orbit parallel to the $`z`$-axis and rotations about the $`z`$-axis. In a system with continuous symmetries, when taking the trace of the electron Green’s function to calculate the density of states, one must first integrate exactly over these symmetries  before employing the semiclassical (stationary phase) approximation. For electrons in a cylinder of length $`L`$, one obtains the following trace formula

$$\delta g(E)=\frac{mL}{\pi \mathrm{}^2}\underset{w=1}{\overset{\mathrm{}}{}}\underset{v=2w}{\overset{\mathrm{}}{}}\frac{f_{vw}L_{vw}}{v^2}\mathrm{cos}\left(\frac{S_{vw}(E)}{\mathrm{}}\frac{3v\pi }{2}\right),$$
(6)
where $`m`$ is the electron mass, $`v`$ and $`w`$ are defined in Fig. 3, $`f_{vw}=1+\theta (v2w)`$ counts the discrete degeneracy of the orbit under time reversal, $`L_{vw}=2vR_0\mathrm{sin}(\pi w/v)`$ is the length of a periodic orbit, and the action $`S_{vw}/\mathrm{}=k_FL_{vw}\sqrt{E/E_F}`$.

When the wire is deformed, the translational symmetry is broken, and we use a semiclassical perturbation theory  to calculate $`\delta g`$. This theory assumes that for small perturbations, the amplitudes and topology of the unperturbed orbits can still be used in (6), but that the lengths of the orbits—and hence their actions—change. The trace formula becomes

$$\delta g(E)=\frac{m}{\pi \mathrm{}^2}\underset{w=1}{\overset{\mathrm{}}{}}\underset{v=2w}{\overset{\mathrm{}}{}}\frac{f_{vw}}{v^2}\text{Re}\left\{e^{i(S_{vw}/\mathrm{}3v\pi /2)}_0^L𝑑z(L_{vw}+\mathrm{\Delta }L_{vw})e^{i\mathrm{\Delta }S_{vw}/\mathrm{}}\right\},$$
(7)
where $`\mathrm{\Delta }S_{vw}/\mathrm{}=k_F\mathrm{\Delta }L_{vw}\sqrt{E/E_F}`$, and

$$\mathrm{\Delta }L_{vw}=2v\mathrm{sin}(\pi w/v)(R(z)R_0),$$
(8)
which may be expressed in terms of the perturbation coefficients using Eq. (1). Combining this result for $`\delta g`$ with (5), it is straightforward to calculate the density of states up to second order in the coefficients $`b(q)`$ for a deformed wire. The result can be integrated to obtain the grand canonical potential, which indeed has the expansion (3).

## III Stability Analysis

Let us first discuss the stability of a nanowire at zero temperature. Fig. 2 shows the stability coefficient $`\alpha `$ (lower diagram) and the density of states at the Fermi energy (upper diagram) as functions of the radius of the unperturbed wire. The wavelength of the perturbation was taken to be $`qR_0=1`$, the critical wavelength for stability in Plateau’s classical analysis of a body under the influence of surface tension. In addition to surface tension, the present model for $`\mathrm{\Omega }`$ has a curvature energy, which enhances the instability for small $`R_0`$, and an oscillatory component associated with the opening of successive subbands as $`R_0`$ increases. As discussed above, $`\alpha `$ has sharp negative peaks—indicating strong instabilities—at the subband thresholds, where the density of states is sharply peaked. Under surface tension and curvature energy alone (dashed curve in Fig. 2), the wire would be slightly unstable at the critical wavevector $`qR_0=1`$, since the curvature term is negative. However, the quantum correction is positive in the regions between the thresholds to open new subbands, thus stabilizing the wire. Since the oscillatory contribution to $`\alpha `$ is independent of $`q`$, we find that regions of stability persist for arbitrarilly long wavelength perturbations, indicating that an infinitely long cylindrical wire is a true metastable state if the radius lies in one of the windows of stability.

With these results, we can construct the zero temperature stability diagram for the wire \[see Fig. 4(a)\]. In contrast to Plateau’s classical stability analysis, an additional quantum length scale arises here, namely the Fermi wavelength $`\lambda _F`$. The stability problem is now determined by two dimensionless parameters: $`qR_0`$ and $`k_FR_0`$. In Fig. 4, regions of instability, where the coefficient $`\alpha (q)`$ is negative, are shaded grey, while the stable regions are shown in white. Note that many of the white regions of stability persist all the way down to $`q=0`$. The multistability of the system, indicated by the alternating stable and unstable stripes, reflects commensurability effects between $`\lambda _F`$ and $`R_0`$.

We remark that in addition to the axially symmetric modes considered here, which are the sole unstable modes in the classical limit, perturbations which break axial symmetry may also become unstable near the subband thresholds. However, these Jahn-Teller like modes  are likely less unstable than the axisymmetric modes, and will not destroy the regions of stability shown in Fig. 4(a).

## IV Conductance Magic Numbers

The electrical conductance $`G`$ of a perfect cylindrical nanowire is quantized  in units of $`G_0=2e^2/h`$. The quantized conductance values of the stable cylindrical configurations are indicated by bold numerals in Fig. 4(a). For comparison, the right vertical axis of the figure is labeled with the corrected Sharvin conductance

$$G_S=\left(\frac{k_FR_0}{2}\right)^2\left(1\frac{2}{k_FR_0}\right),$$
(9)
which gives a smooth approximation to $`G/G_0`$. The conductance values of the stable configurations are somewhat analogous to the magic numbers of enhanced stability in atomic nuclei  and metal clusters . An important distinction is that the magic numbers in nuclei and clusters refer to the number of fermions in a finite system, while we consider an infinite, open system, with magic numbers describing the number of conducting transverse modes  which hold the wire together. The number of conducting modes is approximately equal to the number of atoms which fit within a cross section of the wire.

The integer magic numbers identified here must be distinguished from the concept of “magic wire configurations” proposed by van Ruitenbeek and others . The latter represented discrete minima of the energy per unit volume of a jellium wire as a function of its cross-sectional area, but the stability with respect to deformations with $`q>0`$ was not considered. A similar scenario was advanced in Ref.  based on a partial summation of the periodic orbits in the expression for $`\mathrm{\Omega }`$ itself. While the theoretical results of Refs.  seem to imply that stability occurs only for a discrete set of radii, of measure zero, our analysis finds stability with respect to small perturbations over broad intervals of radius. Furthermore, the very existence of discrete energetic minima of the type discussed in Refs.  depends sensitively on the numerical value of the surface tension in the model; for example, they do not occur whatsoever in the free-electron model at constant volume for $`G>G_0`$. In contrast, the finite windows of stability in our analysis are robust with respect to variations of the surface tension.

## V Comparison to Experiments in Alkali Metals

The sequence of magic numbers $`G/G_0=1,3,5,6,\mathrm{}`$ is consistent with the observation of conductance quantization in alkali metal point contacts . Recently, conductance histograms for sodium nanowires with pronounced peaks up to $`G/G_0100`$ were obtained by Yanson et al. . They argued that these peaks could not be understood based on conductance quantization alone, but rather reflected energetically preferred wire configurations. In order to construct a theoretical conductance histogram, we need to know the a priori probability of occurence of a nanowire of a given cross-sectional area. Here, we make the simplest hypothesis: that nanowires of different cross-sectional areas occur with equal a priori probability if stable, and with zero probability otherwise. On this hypothesis, the relative probability of observing a contact with a given quantized conductance value is proportional to the width $`\mathrm{\Delta }G_S`$ of the corresponding stable region, $`G_S`$ being a dimensionless measure of the contact area.

The conductance histogram from Ref.  taken at a temperature of $`T=80\text{K}`$ is reproduced in our Fig. 5(a). For comparison, the theoretical magic numbers at the same temperature are plotted as vertical bars, with height equal to the width $`\mathrm{\Delta }G_S`$ of the corresponding stable region. Unlike the idealized wires in our analysis, the experimental wires are of finite length, and contain imperfections. Thus the peaks in the experimental histogram are shifted  below the theoretical integer values due to backscattering, and are broadened  due to tunneling, disorder-induced conductance fluctuations, and inelastic processes. Furthermore, the relative heights of the peaks may be influenced by dynamical as well as energetic effects; in particular, the peak at $`G=G_0`$, which is quite pronounced at the lowest temperatures , is absent from the experimental histogram at $`T=80\text{K}`$—presumably indicating that thermally activated processes lead to the rupture of this metastable configuration.

While the first four theoretical peaks at $`G/G_0=1,3,5,6`$ can be identified unambiguously with the narrower peaks in the low temperature data of Refs. , it is not entirely clear whether one can match up the broad peaks in the 80K experimental histogram with individual theoretical magic numbers. In particular, there appears to be fine structure in the theoretical histogram which is not present in the 80K histogram, possibly because it is obscured due to broadening. In order to see whether the theoretical histogram nonetheless correctly describes the gross features of the experimental histogram, it is useful to take the Fourier transforms of the two histograms, singling out the contributions of the shortest periodic orbits.

The actions of the periodic orbits are proportional to the radius of the wire, and hence are approximately proportional to the square root of the conductance, by Eq. (9). The contributions of the various periodic orbits to the conductance histogram can thus be extracted  by taking a Fourier transform with respect to $`\sqrt{G/G_0}`$. The Fourier power spectrum of the theoretical conductance histogram is shown as a solid curve in Fig. 5(b), where all magic numbers up to $`G/G_0=207`$ were included. It shows clear peaks corresponding to the three shortest periodic orbits. For comparison, the experimental Fourier spectrum for sodium nanowires at 90K obtained by Yanson et al.  is shown as a dashed curve. The experimental spectrum contains two broad peaks, one centered at the frequency of the diameter orbit, and a second which spans the frequencies of the triangle and square orbits. The experimental peaks in Fig. 5(b) are broader than the theoretical peaks since the oscillatory structure in the experimental histogram is damped for $`G/G_0>50`$. The widths of the theoretical peaks, in contrast, are determined by the maximum conductance included, $`\text{max}\{G/G_0\}=207`$; increasing this cutoff decreases the widths of the peaks. We have simply taken a cutoff larger than the maximum measured conductance, so that the cutoff does not introduce any artifacts in the theoretical Fourier spectrum. The overall vertical scale of the power spectra shown is arbitrary, so a comparison between theory and experiment must be based on the relative spectral weights of the various peaks. In both spectra, the weight of the diameter orbit is roughly half the combined spectral weights of the triangle and square orbits—a rather good agreement between theory and experiment.

## VI Classical Limit

At zero temperature, the pattern of stable regions separated by unstable stripes shown in Fig. 4(a) continues up to arbitrarilly large radii. However, at any finite temperature $`T`$, the quantum oscillations in $`\alpha `$ are smoothed out, and the classical stability criterion $`qR_0>1`$ is recovered asymptotically for sufficiently large radii. The crossover from the $`T=0`$ result to the classical limit occurs when $`k_BTE_F(G_0/G)`$, i.e., when the thermal energy $`k_BT`$ is comparable to the average transverse level spacing. Fig. 4(b) shows the stability diagram for $`T/T_F=.05`$, where $`T_F=E_F/k_B`$ is the Fermi temperature. One sees that the stability boundary indeed begins to cross over to the classical line $`qR_0=1`$ for $`G_S>20`$.

In Fig. 4(b), there are no true metastable configurations with $`G_S>25T_F/T`$, indicating that all thicker wires would be dynamically unstable (like a column of fluid) at this temperature, once the electronic shell effects have been smoothed out. However, $`T_F=3.75\times 10^4\text{K}`$ in sodium, so multistability from electronic shell effects can be expected to occur in sodium contacts with $`G/G_0125`$ up to at least 300K.

## VII Discussion

It should be pointed out that thermal averaging is not the only mechanism which can suppress electronic shell effects. Disorder also tends to smooth out the sharply peaked structure in the density of states, so that one can expect a reduction of shell effects when the diameter of the wire exceeds the mean free path. Furthermore, the tendency of the positive ions to order themselves into regular arrays  will certainly affect the stability of metallic nanowires. Indeed, pioneering theoretical investigations  of the dynamics of nanowires focused exclusively on the arrangement of the ions. Based on the relative importance of electronic shells and crystal structure in metal clusters , one would expect electronic shell effects to dominate the energetics of very thin wires, particularly in the alkali metals, with crystal structure becoming increasingly important for thicker wires, and for metals where the bonding is more directional.

Although conductance histograms of gold nanowires  do not exhibit the sequence of magic numbers predicted here and observed in alkali metals, there is some evidence that gold nanowires can otherwise be adequately described using the free electron model . It is thus worthwhile to speculate about a possible electronic origin for the remarkable stability of wires of individual gold atoms. Linear chains composed of four to seven gold atoms suspended between two gold electrodes, with a conductance $`G=G_0`$, were found to be stable in the laboratory for hours at a time . Given that such a configuration has an enormous surface energy, its stability is at first sight surprising. However, in our free electron model, we find that an infinitely long wire with a conductance of $`G_0`$ is indeed stable with respect to small perturbations.

Our results may also be relevant to explain the observed stability of gold wires with larger radii . The atomic-scale structure of these wires exhibits a remarkable diversity, ranging from helical, multishell structures  to crystalline structures with surface reconstruction ; but a common factor in the observed stable structures is that they are almost perfectly cylindrical. A direct comparison to our theoretical stability analysis would be facilitated by measurements of the conductance of the stable wires, which have not yet been carried out. However, a rule of thumb is that the conductance $`G/G_0`$ of a metal wire made of monovalent atoms is roughly equal to the number of atoms in the cross section. Thus the hexagonal prism structure with a cross section of 30 atoms determined in Ref.  does correspond to a stable configuration in our analysis \[vertical bar at $`G/G_0=30`$ in Fig. 5(a)\]. Similarly, the thinnest of the helical wire structures determined in Ref. , with a cross section of 8 atoms, also corresponds to a stable structure in our analysis, while the two thickest wire structures determined in Ref. , with cross sections of 22 and 24 atoms, respectively, appear to straddle the predicted island of stability at $`G/G_0=23`$. Further work is needed to elucidate the interplay of atomistic stacking effects and electronic shell effects in these structures.

Finally, let us comment on the dynamical evolution of a nanowire under elongation or compression. Consider stretching a nanowire that is initially in a stable configuration (white areas in Fig. 4). Under elongation, the radius of the wire decreases, so that one moves downward on the stability diagram. When a stability boundary is encountered, it becomes energetically (and dynamically) favorable for the wire to deform spontaneously, until another stable configuration of smaller radius is reached, thus causing the conductance to jump abruptly from one magic number to a smaller one, and conversely under compression. This scenario is consistent with the claim  that the structure of a metallic nanowire undergoes a sequence of abrupt changes as a function of elongation or compression. The finite widths of the unstable tongues in Fig. 4 also provides a possible explanation for the hysteresis  observed in the conductance as a function of elongation: the critical radius at which the wire’s conductance jumps between neighboring magic numbers is different, depending on whether the tongue is approached from above or below, i.e., depending on whether the wire is stretched or compressed.

## Acknowledgments

After completion of this work, we became aware of additional relevant experimental results and a novel data analysis by A. I. Yanson et al. . We thank J. M. van Ruitenbeek and A. I. Yanson for extensive discussions, and for permission to reprint the experimental data shown in Fig. 5. F.K. and H.G. were supported by Grant SFB 276 of the Deutsche Forschungsgemeinschaft, C.A.S. by NSF Grant DMR0072703, and R.E.G. by NSF Grant DMR9812526. This research was supported by an award from Research Corporation.


text:
# Are there universal freeze-out parameters for central 158𝐴GeV Pb+Pb collisions?

## 1 Introduction

Experiments with ultrarelativistic heavy ions at the CERN SPS which are on the quest for a new state of strongly interacting matter have obtained a wealth of information during the last years (see e.g. ). While there are several interesting non-trivial observations in these experiments, our understanding of the evolution of the reactions is far from complete. However, for an interpretation of the promising possible signatures of the quark-gluon plasma a better understanding of the reaction dynamics is mandatory.

Hydrodynamical models or simple parameterizations have been used as tools to help understanding some of the issues in heavy ion reactions. While even the more complex versions of such models are believed to be crude approximations, it is worthwhile to use them and also the simpler parameterizations to test strongly different scenarios and provide first guidelines on the possible dynamics, if the reaction mechanisms are not too far away from the concept of local thermal equilibrium. One of the key questions in this context is the partition of the excitation energy into collective and thermal (i.e. random) motion of the produced particles. This is mostly characterized by the average values of the temperature and the radial flow velocity of the finally observed hadrons, the so-called kinetic freeze-out parameters.

Still, in such comparisons one has to be cautious to

1. use only those experimental results which are most likely too share the same equilibrium conditions and
2. avoid approximations that are wrong even in the most optimistic cases.

In the following I will study some of the experimental results on momentum spectra and interferometry that have been used in applications of these hydrodynamical models and check their implications for the estimates of freeze-out parameters. I will first use estimates based on effective slopes of transverse momentum spectra for different particle species and compare them with the results of hydrodynamical parameterizations of neutral pion spectra. I will then turn to the interpretation of pion interferometry results. There the momentum dependence of effective radii extracted from fits to two-particle correlation functions is used to extract information about temperature and flow. As the method to obtain this information is much more indirect, it is naturally more complicated to interpret the corresponding results. I will try to demonstrate that some of the recipes established should be critically revisited.

## 2 Analysis

### 2.1 Inverse Slopes

The transverse momentum spectra of hadrons are believed to be frozen in a relatively late stage of the nuclear collision. If there is local equilibration, the elementary fireball volumes can be associated with local temperatures. In a finite system there will necessarily be hydrodynamic motion because of the pressure gradient. Therefore, at freeze-out the system may be characterized by local temperatures and transverse flow velocities which will determine the shape of the transverse momentum spectra.111There is an additional dependence of the shape on particle production by resonance decay which I will not discuss in detail here. If the composition of the system is relatively homogeneous, the freeze-out may occur at a relatively similar temperature throughout the whole volume. It is therefore reasonable to speak of one general freeze-out temperature. One should keep in mind, however, that this is a crude approximation: It is more a question of how different all the local temperatures are and for which particle species they are relevant. This simplification may be used for the discussion of momentum spectra, where a superposition of distributions with different temperatures may well be approximated by one distributions with an average temperature. This simplification is much more crucial in the case of pion interferometry, as I will discuss below.

The general idea in using different particle species to estimate these parameters relies on the fact that in a system with a temperature and a collective flow velocity, heavier particles will acquire higher momenta from the shared velocity because of their larger mass and will exhibit flatter momentum spectra, i.e. larger inverse slopes. In  the authors have presented inverse slopes $`T_{\mathrm{eff}}`$ of transverse mass spectra of different particle species and used those to estimate a freeze-out temperature and flow velocity from the relation:

$$T_{\mathrm{eff}}=T+m\beta _T^2.$$
(1)
They obtain a freeze-out temperature of $`T140\mathrm{MeV}`$. This analysis has been criticised as neglecting relativity , the authors of the latter publications present a better estimate of the slopes:

$$T_{\mathrm{eff}}=\sqrt{1\beta _T^2}\frac{T}{1\beta _T\sqrt{1+m^2/p_T^2}}.$$
(2)
Then the limit of zero particle mass:

$$T_{\mathrm{eff}}=T\sqrt{\frac{1+\beta _T}{1\beta _T}}$$
(3)
is used to reestimate the freeze-out temperature. They obtain $`T90\mathrm{MeV}`$ using the data of . Both groups use also somewhat more sophisticated models to obtain freeze-out parameters which support their simple temperature estimates.

I would like to demonstrate the ambiguities in such estimates by using the correct relativistic formula 2 to calculate freeze-out temperatures from the same slope parameters. For this purpose I have to assume $`p_T=p_T`$. This is of course an approximation which should however not be as severe as those used in  and . The slopes and the values of $`p_T`$ used are given in table 1. It turns out that not all three particle species investigated can be described with the same set of parameters. Comparing positive kaons and protons yields $`T=139\mathrm{MeV}`$ and $`\beta _T=0.33`$, while for negative kaons and antiprotons one gets $`T=152\mathrm{MeV}`$ and $`\beta _T=0.28`$. However, estimates using pions and (anti)protons yield $`T100\mathrm{MeV}`$ and $`\beta _T0.41`$ and those for pions and kaons yield $`T84\mathrm{MeV}`$ and $`\beta _T0.51`$.

Although the authors of  exclude the region of $`m_Tm<0.25\mathrm{GeV}/c`$ from the fit for pions, it is still very likely that the slopes for pions are influenced by resonance decays which would result in a lowering of the inverse slope and thus implicitly lower the extracted freeze-out temperature. Kaons and protons should be relatively free from such influence, so that the estimate using these particle species should be more reliable. The temperature $`T=140150\mathrm{MeV}`$ is also in line with the hydrodynamical fit analysis in .

The detailed analysis in  which obtains lower freeze-out temperatures uses additional information, namely the normalization of the spectra and two-particle correlations. It is known from other analyses (see e.g. ) that two-particle correlation which can provide constraints on $`\beta _T^2/T`$ within hydrodynamical models effectively gives very little information on the temperature but constrains mostly the flow velocity (see light grey area in figure 1). On the other hand the normalization of spectra of different particle species is related to *chemical freeze-out*, while the spectral shape is influenced by *kinetic freeze-out*.222There is an influence of the chemical temperature on the spectral shape via the resonance decays – this should however be relevant only for pions. To assume that these two corresponding temperatures are the same, as was apparently done in , is very likely wrong.333A similar assumption is used in the analysis in  which will not be discussed here.

From the large pion-nucleon cross section one would assume that these two particle species are most likely to share the same freeze-out conditions. However, the influence of resonance decays in the pion spectrum forbids the simple use of slope parameters to estimate freeze-out conditions. To perform such a comparison the slope parameters for heavier particles will be used together with a full hydrodynamical parameterization of neutral pion spectra as described in . The hydrodynamical computer program described in  has been modified to describe neutral pions and accommodate a Woods-Saxon density profile in addition to the default Gaussian profile and was used to obtain best parameter sets for $`T`$ and $`\beta _T`$. The $`2\sigma `$-allowed values obtained for different profiles (Gaussian and Woods-Saxon of different diffuseness) are shown in figure 1 as a grey band. Among these fits, those obtained with large diffuseness (as e.g. the Gaussian) tend to yield smaller flow velocities for the same temperature but also constrain the allowed parameters to the high $`T`$ region . In addition, the slope parameters for kaons, protons, antiprotons and deuterons from NA44 and NA49 as summarized in table 1 have been used with the help of equation 2 to obtain similar relations between $`T`$ and $`\beta _T`$. All these estimates are also displayed in figure 1.

The NA49 collaboration has also published freeze-out constraints obtained from negative hadron and deuteron transverse mass spectra . The negative hadrons are dominated by pions which should allow to use them in hydrodynamical descriptions assuming the pion mass, however, the analysis was performed using a non-relativistic approximation and ignoring resonance production. The deuteron spectra are not affected by resonances but still suffer from the non-relativistic approximation. The corresponding parameters are therefore not reliable. In fact, the estimate from the slope included in figure 1 differs considerably from the constraints in , esp. at large velocities. While the curves for kaons, protons and deuterons seem to converge around $`T=140\mathrm{MeV}`$ and $`\beta _T=0.33`$, the neutral pion results have a further bias to somewhat higher temperatures. One should note that the analysis of the pion spectra differs from the others, as

1. a full parameterization including the effects of resonance decays has been performed compared to using an effective slope together with an approximate relation connecting $`T`$ and $`\beta _T`$ and
2. they cover a different $`m_T`$ range – the high $`m_T`$ tail of the neutral pion spectrum has also been used as an upper limit as described in .

Still this figure suggests that a common parameter set for all the momentum spectra considered here can only be obtained for $`T135\mathrm{MeV}`$ and $`\beta _T0.35`$. This appears to contradict constraints obtained from pion interferometry measurements  (see light grey area in figure 1), which favor larger transverse flow velocities.

### 2.2 Interferometry

It is well known that correlations between position and momentum space modify the radii measured in pion interferometry. One of the most interesting examples of such a correlation is collective flow , both longitudinal and transverse, although there are other examples, like a source with local temperatures $`T=T(x^\mu )`$, as I will discuss below. In the idealized case of pure collective motion (no thermal contribution) only the source element moving towards the detector will be measured, so the radius extracted from interferometry will correspond to this small elementary volume. If there is additional thermal motion, the probability of a particle to still be measured while being radiated from another source element at another space point and directed slightly away from the detector will increase. The effective source size will therefore decrease, if the amount of collective motion relative to thermal motion is large, and vice versa. In addition, this effect is a function of the momentum of the particle, and in most scenarios the effective radii will decrease with increasing transverse mass of the particle.

Most quantitative studies of these effects start from the Cooper-Frye formula  for the momentum distribution of the source:

$$E\frac{d^3N}{dp^3}=\frac{2J+1}{(2\pi )^3}_\mathrm{\Sigma }d^3\sigma _\mu \frac{p^\mu }{\mathrm{exp}\left\{\left[p_\mu v^\mu \mu \right]/T\right\}1}.$$
(4)
In this formula the factor $`p^\mu `$ in the numerator and the exponent $`p_\mu v^\mu \mu `$ originate from the Lorentz transformation with $`v^\mu =v^\mu (x^\mu )`$ being the local four-velocity of the source elements. In general also $`T`$ is again a local variable. $`\mu `$ is the chemical potential and the $``$ signs are relevant for bosons and fermions, resp.

There has been substantial work on extracting expressions for observables in interferometry, the most systematic approach using a parameterization of the source distribution and a saddle point approximation to extract approximate analytic relations for the effective radii . A relatively general result is obtained using a source distribution expressed as:

$$S(x,K)=\frac{m_t\mathrm{ch}(\eta Y)}{(2\pi )^3\sqrt{2\pi (\mathrm{\Delta }\tau )^2}}\mathrm{exp}\left[\frac{Ku}{T}\frac{x^2+y^2}{2R_G^2}\frac{(\tau \tau _0)^2}{2(\mathrm{\Delta }\tau )^2}\frac{\eta ^2}{2(\mathrm{\Delta }\eta )^2}\right].$$
(5)
Gaussians are used for the space-time and rapidity distributions, the momentum spectrum is approximated as a Boltzmann distribution:

$$E\frac{d^3N}{dp^3}\mathrm{exp}\left(\frac{Ku}{T}\right).$$
(6)
For a model with transverse expansion, the authors obtain approximate relations for the transverse and longitudinal radii:

$`R_T`$ $``$ $`R_G\left[1+{\displaystyle \frac{m_T\beta _T^2}{T}}\mathrm{ch}(\overline{\eta }Y)\right]^{1/2}`$ (7)
$`R_L`$ $``$ $`\tau _0\left[{\displaystyle \frac{m_T}{T}}\mathrm{ch}(\overline{\eta }Y){\displaystyle \frac{1}{\mathrm{ch}^2(\overline{\eta }Y)}}+{\displaystyle \frac{1}{(\mathrm{\Delta }\eta )^2}}\right]^{1/2}.`$ (8)
For a boost invariant longitudinally expanding source the expressions are more simple, e.g. for the longitudinal radius one obtains:

$$R_L\tau _0\sqrt{\frac{T}{m_T}}.$$
(9)
As it turns out  that the source rapidity and the pair rapidity are similar $`\overline{\eta }Y`$, also the expression for the transverse radius may be simplified:

$$R_TR_G\left[1+\frac{\beta _T^2}{T}m_T\right]^{1/2}.$$
(10)
These two latter equations are mostly used to extract information from experimental two-particle correlation functions. There are a number of assumptions which have to be made to obtain equations 7 and 8 and even more for equations 9 and 10. Some of these are discussed in , however there are further assumptions which I will try to concentrate on here.

To be able to use equation 10 to extract $`\beta _T^2/T`$ as in , one of the major assumptions is that the parameter $`R_G`$, which contains information on the geometrical size of the source, does not depend on the transverse momentum of the emitted particles. This is not necessarily correct. In general, there will not be a universal value for the temperature, it will be a local parameter depending on the position of the source element with respect to the center or the surface as discussed above.

This induces an additional correlation between momentum and coordinate space. If the variations in $`T`$ are very strong, this would of course invalidate all the attempts to extract *the* freeze-out temperature. If the variations are relatively weak, it may still be reasonable to use one effective temperature. It may e.g. be used to characterize momentum spectra, as for those the spatial distribution and therefore also these correlations are not directly relevant.

However, they are crucial to the analysis of pion interferometry, especially to the above mentioned method of extracting $`\beta _T^2/T`$. Consider e.g. that a certain fraction of pions is emitted before freeze-out, from the surface of the source, and that this fraction increases with increasing momentum. One might try to describe this with a superposition of two sources. Just as one example one may consider the following – arbitrary – scenario: One source of effective size (Gaussian) $`R_1=10\mathrm{fm}`$ and temperature $`T_1=200MeV`$ and a second with $`R_2=2\mathrm{fm}`$ and $`T_2=400MeV`$.444For simplicity exponential distributions in transverse mass are assumed. If one investigates produced particles in different intervals of transverse mass, the relative contributions of these two sources will change, and implicitly the effective source size as seen in interferometry. This is shown in figure 2, where the effective source size – again approximated by a Gaussian – is shown as a function of the transverse mass for the case where the two sources discussed above both contribute the same fraction to the total multiplicity of particles. It can be seen easily that the effective source size decreases significantly with increasing transverse mass. Of course, this scenario has not been derived from physical arguments here, so it should only serve as an illustration of possible effects of such correlations between momentum and coordinate space. How large such effects might be would require a detailed model calculation which is not the scope of this paper. One should only note that in general $`R_G`$ in equation 7 is expected to be a function of $`m_T`$ within the model used. Values of $`\beta _T^2/T`$ should therefore be considered as upper limits unless more detailed model comparisons are performed. Pion interferometry is much more sensitive to such effects than momentum spectra.

A more fundamental caveat to the use of the $`m_T`$ dependence of HBT to estimate freeze-out parameters relates to another approximation used in extracting the analytical formula which are applied in most of the analyses. Specifically, the authors of  use the Boltzmann approximation to the distribution function. While this may be a good approximation for most of the pions used in experimental HBT analyses, it is well known that the deviations to the true Bose-Einstein distributions are crucial at low momenta. This is illustrated in figure 3, where the inverse slopes

$$T_{eff}=\left[\frac{d}{dE}\mathrm{log}f(E)\right]^1$$
(11)
of Bose-Einstein and Fermi-Dirac distributions

$$E\frac{d^3N}{dp^3}\frac{1}{\mathrm{exp}\left(\frac{E}{T}\right)1}.$$
(12)
are shown. Strong deviations from Boltzmann behaviour are seen for both distributions for $`E/T<1`$.

As in the analyses in question most of the $`m_T`$ dependence is seen at relatively low momenta, one should carefully reinvestigate the interpretation of these data with the correct statistical distributions. The trend of such a corrected analysis may be guessed from the behaviour of the slopes. Using the crude assumption that the BE-distributions may be locally approximated by a Boltzmann distribution, one might then use the formalism developed in  to estimate the $`m_T`$ dependence. However, now the effective slope for a Bose-Einstein distribution ($`T_{BE}`$) as shown in figure 3 should be used instead of the Boltzmann temperature $`T`$. In figure 4 the transverse radii as a function of transverse mass are shown for three different estimates:

1. Following the formalism as in , i.e. using the Boltzmann-like temperature in equation 10 and the results obtained in the article, i.e. $`T=120\mathrm{MeV}`$, $`\beta _T=0.55`$ and $`R_G=6.5\mathrm{fm}`$,
2. substituting the temperature by the corresponding effective slope of a Bose-Einstein distribution, but using the same parameters as above and
3. substituting the temperature by the corresponding effective slope of a Bose-Einstein distribution, but now using $`T=85\mathrm{MeV}`$, $`\beta _T=0.55`$ and $`R_G=7.1\mathrm{fm}`$.

It can be seen from this example that cases 1 and 2 show significant differences, and that for a rough agreement with 1 the temperature in the corrected version has to be modified from 120 MeV to 85 MeV as in case 3.

In this context it is interesting to investigate also results for proton source sizes. While the mechanism generating the correlation of particle pairs is different compared to pions – a strong interaction quasi-resonance instead of Bose-Einstein statistics, the picture of effective source sizes modified by the interplay of random thermal and collective hydrodynamical motion remains valid. People have therefore applied similar ideas to pp-correlation radii. The experiments NA49  and NA44  have measured proton correlations at midrapidity in Pb+Pb collisions - they obtain Gaussian radii in the range of $`R_{pp}=3.54.0\mathrm{fm}`$ which are considerably smaller that those extracted from pion interferometry ($`R_{\pi \pi }=67\mathrm{fm}`$). Similar values have been estimated by applying coalescence models to measured proton and deuteron multiplicities . The smaller effective radii have been attributed to the much larger transverse mass of the protons, if a scaling with $`1/\sqrt{m_T}`$ which may be derived as an approximation from equations 9 and 10 is used for both protons and pions.

From the above argument on BE-distributions it is obvious that equivalent modifications are important here. However, the effective slope for a FD-distribution is larger at low momenta. In addition, the exponent of the FD-distribution contains a non-vanishing chemical potential which makes it much more likely that the exponent is small and that modifications are important.

Actually, another important argument is relevant here. If, as commonly believed, chemical freeze-out occurs relatively early before kinetic freeze-out, the relevant picture for the latter is the canonical ensemble, as the particle number for each species is fixed, so that there is no chemical potential as a parameter. The relevant energy in the statistical distributions is the kinetic energy and not the total energy, as no new particles are produced and the ground state energy of each particle is the rest mass. The Cooper-Frye formula should then be changed to:

$$E\frac{d^3N}{dp^3}=\frac{2J+1}{(2\pi )^3}_\mathrm{\Sigma }d^3\sigma _\mu \frac{p^\mu }{\mathrm{exp}\left\{\left[p_\mu v^\mu m_0\right]/T\right\}1}.$$
(13)
The major change for the transverse degrees of freedom would be to use $`m_Tm_0`$ instead of $`m_T`$. This should hold for both bosons and fermions, but while it might have only small effects on the pions, it will be dramatic for protons because of their large mass. It is beyond the scope of this paper to estimate the effect of such modifications on the $`m_T`$ spectra in general or, more specific, on the $`m_T`$ dependence of proton source sizes. However, it is also premature to conclude that proton and pion source sizes are compatible without performing detailed calculations. As the values of $`m_Tm_0`$ for protons and pions are roughly similar in the experiment, no difference in the radii would be expected between protons and pions from their transverse masses. This would however contradict the measured source sizes.

## 3 Conclusion

I have tried to illustrate above that the freeze-out parameters in ultrarelativistic heavy ion collision are not well established quantitatively. One should be careful to use only data sets which very likely share the same conditions. Also, the used approximations should be thoroughly investigated. The simple re-analysis of some of the available data performed here indicates the following:

* Transverse momentum spectra of deuterons, protons, kaons and pions seem to indicate that compatible parameter sets can be obtained only for $`T135\mathrm{MeV}`$ and $`\beta _T0.35`$.
* The $`m_T`$ dependence of interferometry radii has to be reinvestigated using true statistical distributions instead of Boltzmann approximations at low momenta. This would also require to study the role of modifications related to the change at chemical freeze-out to a system with fixed particle number.
* Other possible sources of correlations between momentum and coordinate space may cause further alterations of the $`m_T`$ dependence of interferometry radii.
* The difference in pion and proton source sizes is very likely not explained by the difference in transverse mass of the two particles.

The three last observations shed doubt on the use of pion interferometry as the main ingredient to obtain constraints on the amount of transverse flow at freeze-out. A more thorough analysis of the pion interferometry results is necessary to clarify the situation.

Two possible scenarios may be envisaged which explain the data discussed here:

1. Pion HBT does not yield reliable constraints on kinetic (thermal) freeze-out.555It may be that a treatment of collective flow influence on interferometry using correct statistical distributions can yield a consistent picture for pions and protons and can serve to obtain freeze-out constraints – these will likely be considerably different from the ones presently used. It may actually be that pion interferometry measures a different, later time than proton correlations and momentum spectra. This late “wave function freeze-out” may be related e.g. to very soft or mean field interactions of pions which do not alter the momentum distributions as mentioned in . Momentum spectra, however, favour $`T135\mathrm{MeV}`$ and $`\beta _T0.35`$.
2. Freeze-out does not occur with one even approximately universal temperature for different species and also different parts of the source.