text:
# NEW MECHANISM FOR THE PRODUCTION OF THE EXTREMELY FAST LIGHT PARTICLES IN HEAVY-ION COLLISIONS IN THE FERMI ENERGY DOMAINProceeding of VII International School-Seminar on Heavy-Ion Physics, May 27 – June 1, 2002, Dubna, Russia.

## I Introduction

The production of pre-equilibrium light particles (n, p, t, $`\alpha `$) in nucleus-nucleus collisions depends on the way the nuclear system evolves at the reaction initial stage. For heavy ion collisions, the light particle cross section is known to be a noticeable fraction of the total reaction cross section even at low energies of the order of 10 A MeV, i.e. light particle production is a distinctive feature of all nuclear reactions involving heavy ions. This means that studying the production mechanisms for those particles is likely to provide direct information on both the reaction initial stage dynamics and the potential and dissipative forces of nucleus-nucleus interaction.

Vast bodies of data  demonstrate that in heavy ion reactions at energies per nucleon of the order of the Fermi energy, light particles are produced in a wide angular range, their velocities being two and more times higher than the velocities of the beam particles. Fig. 1 shows the measured differential proton cross section at $`\theta _{lab}=20^{}`$ in the case of the $`{}_{}{}^{16}O`$ (20 A MeV) + $`{}_{}{}^{197}Au`$ collision . What is the mechanism of the production of these extremely fast light particles? The ultimate answer to that question has not been provided yet. Circumstantial evidence suggests that they are of a pre-equilibrium nature. This makes them more difficult to study by a direct experiment since there is currently no way of studying processes taking place in time intervals of the order of $`10^{20}`$ s.

Attempts to interpret the experimental picture have resulted in the creation of a number of theoretical models and approaches. Among these in particular are the models of moving sources , Fermi-jet , dissipative break-up and massive transfer  etc. A comprehensive survey of experimental and theoretical works relating to this problem is given in . Applying these approaches, the authors have succeeded in qualitatively describing the energy and angular dependencies of the spectra of emitted light particles as well as revealing some production mechanisms for them. Among the latter in particular is the mechanism of ”splashing” nucleons out of a retarding projectile nucleus, the high velocities of the light particles being due to the addition of the velocity of Fermi motion within the projectile nucleus and the velocity of its centre of mass. The importance of taking this mechanism into account was shown in refs. , when describing the spectra of the $`\alpha `$-particles produced in projectile break-up and massive transfer at collision energies of the order of several tens of MeV per nucleon.

Nucleon-nucleon collisions are also likely to result in a fast light particle being emitted due to scattering on nucleons of high velocities. The Fermian shape of nucleon momentum distributions explains the existence of such nucleons in a heated nuclear system. The first mechanism, capable of qualitatively describing the experimental results in a transparent way, fails to provide satisfactory quantitative agreement. The nucleon-nucleon collision model is adequate at high energies ($`<`$ 100 A MeV), but in the region of energies $`E_0<`$ 50 A MeV the nucleon-nucleon interaction cross section for heavy ion collisions decreases due to Pauli’s principle. Therefore this mechanism influences the formation of the hard part of the spectrum of light particles to a smaller degree. Recent advances in instrumentation have made possible the precise measurement of the angular and energy spectra of light particles , which is likely to shed light upon yet unresolved problems.

## II Mechanisms of Light-Particles Formation

In ref. , we proposed a four-body classical model of nucleus-nucleus collisions, which permits establishing the role of mean nuclear fields and nucleon-nucleon collisions in the production of light-particles. This model considers the projectile and target nuclei to be bound two-body systems that are composed of a heavy core and a light fragment moving in its field. The interactions of the fragments, which follow classical trajectories, define six pair potentials. The interaction between the light fragments and the cores was described by the Woods-Saxon potential with compilation-based parameters . The interaction between the heavy cores was described by the proximity potential or the Woods-Saxon potential, the parameters of the latter being chosen so as to reproduce the height and position of the Coulomb barrier. The coupling with the reaction channels, which were not taken into account explicitly, was introduced by means of dissipative forces acting between the heavy cores of the target and projectile. The choice of friction coefficients and form-factors for the dissipative forces was based on ref. . To provide a more correct estimate of the channel differential cross sections, the description of the relative motion of the projectile (target) light particle and the target (projectile) core has an absorption probability for them has been introduced as $`P_{ij}^{abs}=1\mathrm{exp}\left(s_{ij}\mathrm{/}\lambda _{ij}\right)`$, where $`s_{ij}`$ is the distance traveled by particle i in the nucleus j, $`\lambda _{ij}`$ the corresponding free path length, which is, as well known, related to the imaginary part of the optical potential, the parameters of which were chosen according to ref. . Specifying the relative distance vector $`𝐫_{ij}`$ for the projectile components, the binding energy $`E_{ij}^{sep}`$ , known by experiment, as well as the orbital momentum $`l_{ij}`$ for their relative motion, completely defines the projectile inner spatial configuration. The components of the vector $`𝐫_{ij}`$ were chosen randomly and equiprobably in the classically allowed region. To have the relative motion momentum $`𝐩_{ij}`$ uniquely determined, one of the components of the vector $`𝐥_{ij}`$ must also be specified (by a random choice) in addition to the $`r_{ij}`$ and $`E_{ij}^{sep}`$ values. Repeating the same operations for the target nucleus and specifying the relative motion of the centers of mass of the nuclei depending on the reaction allow to determine the initial conditions that are required for solving the set of classical equations of motion.

The model outlined above has 15 reaction channels with a different combination of particles in the final fragments. Eight of these channels contribute to the total cross section for pre-equilibrium light particles, which can be divided into two groups: particles emitted from the projectile and from the target. Taking into account all the pair potentials allows us to study the whole range of processes that result in light-particle emission as well as to establish the role played in them by one or another type of interaction. We have investigated heavy ion collisions at energies of several tens of MeV per nucleon and thoroughly studied the production mechanisms for fast light particles. It is remarkable and unexpected that the hardest part of the energy spectrum of emitted light particles corresponds to particles emitted from the target nucleus rather than from the projectile nucleus, as has been assumed up to now. There is a simple enough explanation for this phenomenon, which appears to be unusual at first sight. Let us assume for simplicity that the heavy cores have masses much larger than those of the light particles; the light particle moves in the nucleus at a velocity equal to the Fermi velocity $`v_F`$ ; the nuclei move relative to one another at the velocity $`v_0v_F`$ ; the trajectories of the heavy fragments are taken to be undistorted; the binding energy of the light particle in the nucleus is assumed to be much smaller than the collision energy $`E_{ij}^{sep}E_0=mv_0^2\mathrm{/}2`$, here m is the mass of the light particle. Then, as has been shown in , for example, the neutrons emitted from the projectile (”splashing out”) will have a maximum laboratory system energy expressed as

$$E_{n,\mathrm{max}}\frac{m}{2}\left(v_0^2+2v_0v_F\right)=3\frac{mv_0^2}{2}=3E_0.$$
$`(1)`$
In the framework of this model under the same assumptions, the maximum energy of the neutrons emitted both from the projectile and from the target in elastic nucleon-nucleon collisions is limited by the value

$$E_{n,\mathrm{max}}\frac{m}{2}\left(v_0+v_F\right)^2=4\frac{mv_0^2}{2}=4E_0.$$
$`(2)`$
However, such a mechanism of emitting a fast nucleon suggests that the recoil nucleon is to impart the whole of its kinetic energy to the particle knocked out and, consequently, pass into an occupied lower energy state, which is forbidden by Pauli’s principle. Therefore, this mechanism does not manifest itself in full measure at low and intermediate collision energies and comes into play as the collision energy increases, when states more and more lower in energy become unoccupied in the process of the system getting excited.

The production mechanism for the fast light particles emitted form the target is more complicated. To begin with, let us consider the elastic scattering of a target light particle on a moving infinitely heavy repulsing wall. If before colliding, the light particle and the heavy wall have been moving collinearly in opposite directions with parallel velocities $`v_F`$ and $`v_0`$ respectively, then the centre-of-mass velocity of the light particle has been $`v_{cm}=(v_0+v_F)`$ . After elastically colliding, the light particle will have the velocity $`v_{cm}^{}=v_{cm}=(v_0+v_F)`$ , the lab velocity being $`v^{}=v_{cm}^{}+v_0=2v_0+v_F`$ . Thus the maximum velocity of a neutron emitted from the target will be

$$E_{n,\mathrm{max}}4\frac{m}{2}\left(v_0^2+v_0v_F\right)=8\frac{mv_0^2}{2}=8E_0,$$
$`(3)`$
i.e. the fastest neutrons will have a velocity of more than 2.5 times higher than that of the beam particles. Now, suffice it to say that an elastic collision with a repulsing wall is kinematically equivalent to elastic scattering at the angle $`\theta _{cm}=180^{}`$ (i.e. orbiting) in the attractive field of an incident nucleus. However, orbiting is possible at quite low relative motion energies. At higher energies, a light particle can be deflected only due to the field of the projectile at a certain centre-of-mass angle $`\theta _R`$ called an angle of rainbow scattering. Therefore the maximum energy of a pre-equilibrium light particle emitted will be a function of collision energy, light-particle-incident-nucleus interaction potential, particle binding energy and dissipative forces.

More accurate calculations of the maximum energy of the pre-equilibrium neutrons produced in the reaction $`{}_{}{}^{20}Ne`$ \+ $`{}_{}{}^{165}Ho`$ n + X are presented in Fig. 2. Given here is the $`E_{n,\mathrm{max}}/E_0`$ ratio as a function of the beam energy $`E_0`$ for the above mechanisms: the solid curve is for the energy of the neutrons emitted from the target; the dashed curve for the energy of the neutrons emitted from the projectile; the dash-dotted curve for the energy of neutrons produced in nucleon-nucleon collisions, which was calculated with regard to Pauli’s principle. So, it is seen that at collision energies per nucleon smaller than $`E_F`$ (indicated with an arrow in Fig. 2), the mechanism involving acceleration of neutrons by the moving mean field of the projectile is dominant in the formation of the high-energy portion of neutron spectra. In the high-energy region, the key role is played by nucleon- nucleon collisions. It should be also noted that the curves in Fig. 2 do not represent $`E_{n,\mathrm{max}}`$ in a strict way since our calculations did not take into account the non-zero probability of there existing nucleons of velocities much higher than $`v_F`$. Therefore the boundaries indicated with the curves in Fig. 2 will shift to higher energies.

How should the predictions, made by our model be tested, and how should the yields of pre-equilibrium light particles emitted from a projectile and target be separated? There is currently no direct solution. However the above analysis of various processes of light-particle production in heavy ion collisions shows that the mechanism of emitting fast nucleons from the target is only slightly sensitive to the projectile nucleus mass, and, consequently, the maximum energy of light particles is not likely to change greatly with replacing one projectile by another. On the other hand, the maximum energy of the target nucleons accelerated by the attractive field of the incident nucleus is defined by the corresponding rainbow scattering angle, which depends on the dimensions of the deflecting mean field

$$\theta _R\left(V_C0.56U_0\sqrt{R_U/a_U}\right)/E_{cm},$$
$`(4)`$
where $`V_C`$ is the height of the Coulomb barrier for the emitted light particle; $`U_0`$, $`R_U`$, $`a_U`$ the depth, radius and diffusiveness of the light-particle-projectile interaction potential, and $`R_U1.3A^{1\mathrm{/}3}`$. It should be noted, the coefficient 0.56 in empirical equation (4), obtained in ref. , is not valid in case of nucleon scattering on heavy nucleus. We found that coefficient 0.7 enables to obtain more reasonable value of nucleon rainbow angle $`\theta _R`$.

Now it is evident that the energy distribution of emitted light particles must be more elongated for a reaction with a heavy projectile than for a reaction with a light one. This conclusion is supported by Fig. 3, in which the calculated differential pre-equilibrium proton cross sections are shown for the angle $`\theta _{lab}=51^{}`$ for the $`{}_{}{}^{40}Ar`$ \+ $`{}_{}{}^{51}V`$ and $`{}_{}{}^{132}Xe`$ \+ $`{}_{}{}^{51}V`$ reactions (in normal and inverse kinematics) at the energy $`E_0`$ = 44 A MeV. It is seen, first, that the key role in the formation of the high-energy part of the spectrum is played by the protons emitted from the target nuclei (curves 1 and 2 in Fig 3). Second, the distribution of the protons emitted from projectiles of different masses (curves 3 and 4 in Fig. 3) only changes its shape to a large degree at small energies, i.e. a more massive $`{}_{}{}^{132}Xe`$ projectile loses less energy in a collision with a $`{}_{}{}^{51}V`$ target than does $`{}_{}{}^{40}Ar`$ . Therefore these two reactions have projectile proton distributions with shifted maxima. Third, the $`{}_{}{}^{40}Ar`$ induced reaction in normal kinematics has a pre-equilibrium-proton spectrum, which is several tens of MeV shorter than that for the reaction with the $`{}_{}{}^{132}Xe`$ heavier projectile. So comparing precisely measured data on the distributions of the light particles produced in heavy ion collisions in reactions of normal and inverse kinematics will give us vital information about which mechanism is dominant in emitting a fast light particle.


format: text

text:
# References

Quantum spin liquid in $`\mathrm{𝐒𝐫}_\mathrm{𝟐}\mathrm{𝐂𝐮𝐎}_\mathrm{𝟑}`$ induced by Acoustic Phonons

S.S. Aplesnin

L.V. Kirenskii Institute of Physics, Siberian Branch of the Russian Academy of Sciences, Krasnoyarsk, 660036, Russia , FAX 07 3912 438923, E-mail : apl@iph.krasn.ru

We have carried out theoretical studies of the phonons and spinons excitation spectra of $`Sr_2CuO_3`$ in term of a one-dimensional model with non-adiabatic interaction magnetic and elastic subsystems using continuous time Monte Carlo method. Two gaps in triplet excitation spectrum are found. The bond defects concentration estimated is well agreement with the experimental data. The formation of magnetic moment on site at low temperature $`T_N5K`$ is caused by freezing quantum spin liquid with correlation radius of antiferromagnetic ordering $`\xi /a300500`$.

PACS: 75.10.Jm, 75.30.Fv, 67.20.+k, 63.20.Ry

1. Introduction

Quasi-one-dimensional spin systems with antiferromagnetic interactions have received considerable attention because of a number of basic physical phenomena in low dimensions, such as the existence of a spin gap depending on the spin value or the frustration antiferromagnetic exchange, spin Peierls transitions, accompanied by the structural symmetry breaking. Orthorhombic compound $`Sr_2CuO_3`$ is considered as quasi-1D antiferromagnetic system. However, Raman studies reported that extra peaks were observed along the Cu-O chain directions  which implies that there might be some structural instabilities. This effect is confirmed by polarization-dependent infrared-active phonon spectra which have shown splitting of the stretching mode in the Cu-O chain direction about $`20cm^1`$ and became evident as T decreases  . These unusual behavior are not expected by the previous theoretical calculations.

Another unclear behavior is the uniform spin susceptibility $`\chi ^{}(q=0)`$ have decreased steeply below $`T33K`$ without the signature of three -dimensional short range order approaching the Neel state by NMR Knight shift at the in-chain oxygen. The nuclear spin-lattice relaxation rate offsets the $`1/T`$ increase of the staggered spin susceptibility down to $`10K`$ and depends on frequency. The $`q=0`$ mode of spin transport is ballistic and a diffusion-like contribution at $`T<<J`$  . The bulk susceptibility the measurements have shown decreasing above $`10K`$  . The NMR lineshape have been changed in a very specific way, giving rise to ”features”, which have been explained by the chain edges and the presence of mobile bond-defects. A ”defect” is defined as a local change in the magnetic bond coupling. As the translational symmetry of the spin system is broken, a local alternating magnetization develops around the defects in the direction of the applied field. The density of bond-defects decreases linearly with $`T`$. The source of the bond-defects might be related to the interstitial excess oxygen. However , at the lowest temperature, $`T=20K`$, the density of bond-defects $`\rho 2.810^3`$ remains more than ten times larger than the concentration of the residual spin-1/2 impurities $`(1.310^4)`$ after annealing and the observed $`T`$ dependence of $`\rho `$ is not a constant  .

The concept of mobile ”bond-defects” have been used for explanation of the temperature dependence and the anisotropy of the thermal conductivity below 50 K as the main source of spinon scattering at low temperatures . The calculated mean distance between two neighboring bond defects, consistent with the NMR data between 20 and 60 K. The results of the specific heat measurements have been fit to the data above 4 K using the approximation $`C=\gamma T+\beta T^3+\delta T^5`$. The linear in T contribution is due to spinons excitations of a 1D Heisenberg spin-1/2 antiferromagnetic system with exchange $`J=2620\pm 100`$ which is $`18\%`$ larger value than $`J2200`$ deduced from magnetic susceptibility data  and NMR data . The Neel temperature $`T_N`$ dependents on the measurement time and shifted to lower temperature with increasing time. So Neel temperature is $`T_N=5.4K`$ from neutron scattering  , $`4.15K`$ by $`\mu SR`$  and $`3.5K`$ from the anomalies of the specific heat . Possibly , here it is nonequilibrium process. Physical origin of the bond-defects remains unclear also as abrupt change suseptibility at different temperatures determined by NMR and by external magnetic field.

In the paper we try to answer on these questions basing on non-adiabatic interaction between spins . Numerical methods, such as exact diagonalization and density-matrix renormalization group, face the potential difficulty in dealing with a very large Hilbert space of the phonons. Here , we use Monte Carlo approaches restricted to finite chains $`L=100,200`$ but without any adiabatic approximation and the truncation of the infinite phonon Hilbert space. The method  is based on a path-integral representation for discrete system in which we work directly in the Euclidean time continuous. All the configuration update procedures contain no small parameters . Being based on local updates only, it allows to work with the grand canonical ensemble and to calculate any dynamical correlation function, expectation values.
2. Model and method

The ratio intechain to intrachain exchange in $`Sr_2CuO_4`$ is an extremely small $`10^3`$ and demonstrates a good one-dimensionality. The structural fluctuation were only observed along in chain at low temperature. As a result a spin-phonon interaction were considered along chain. The model Hamiltonian of an spin-phonon system takes the form:

$$H=\underset{i=1}{\overset{L}{}}[J+\alpha (u_iu_{i+1})][S_i^zS_{i+1}^z+(S_i^+S_{i+1}^{}+S_i^{}S_{i+1}^+)/2]+M\dot{u}_i^2/2+K(u_iu_{i+1})^2/2.$$
(1)
Here $`S^{z,\pm }`$ are a spin operator components associated with the site i, $`J>0`$ is the usual antiferromagnetic exchange integral, $`\alpha `$ is the spin-phonon coupling constant, $`u_i`$ is the displacement in the z- direction, $`M`$ is the mass of the ion and $`K`$ the spring constant. Using the quantum representation for phonon operators $`b,b^+`$, the Hamiltonian becomes

$$H=\underset{q}{}\underset{i=1}{\overset{L}{}}[J+\alpha \sqrt{\frac{2\mathrm{}}{M\omega _0}}\sqrt{\mathrm{sin}\frac{q}{2}}\mathrm{cos}(q(i0.5))(b_q+b_q^+)][S_i^zS_{i+1}^z+(S_i^+S_{i+1}^{}+S_i^{}S_{i+1}^+)/2]$$
(2)
$$+\underset{q}{}\mathrm{}\mathrm{\Omega }(q)b_q^+b_q,\mathrm{\Omega }(q)=2\sqrt{\frac{K}{M}}\mathrm{sin}(\frac{q}{2}),\omega _0=2\sqrt{\frac{K}{M}},q=2\pi n/L,n=1,2,\mathrm{}L$$

The phonon frequency $`(\omega _0)`$, spin-phonon coupling $`(\alpha )`$, energy $`E`$, temperature $`T`$ are normalized on the exchange $`J`$ and $`\mathrm{}=1,M=1`$. The temperature used in calculation is $`\beta =J/T=25`$. Our system consists of the two subsystem interacted. The elastic subsystem is described by phonons with the number occupation $`n_{ph}=0,1,2..`$ and magnetic subsystem is characterized by the number occupation $`n_m=0,1`$ and Pauli operators $`a,a^+`$ which coincide with $`S`$ spin operators. We start with the standard Green function of the phonon in the momentum $`q`$\- imarginary-time $`\tau `$ representation:

$$G(q,\tau )=<vac|b_q(\tau )b_q^+(0)|vac>,\tau 0,b_q(\tau )=\mathrm{exp}(H\tau )b_q\mathrm{exp}(H\tau ).$$
(3)
$$G(q,\tau )=\underset{\nu }{}|<\nu |b_q^+|vac>|^2\mathrm{exp}[(E_\nu (q)E_0)\tau ]$$

where $`|\nu >`$ is a complete set of the Hamiltonian $`H`$ in the sector of given $`q`$ , $`H|\nu (q)>=E_\nu (q)|\nu (q)>,H|vac>=E_0|vac>,E_0=0`$. Rewriting Eq.(3) as

$$G(q,\tau )=_0^{\mathrm{}}𝑑\omega A(q,\omega )\mathrm{exp}(\omega \tau ),A(q,\omega )=\underset{\nu }{}\delta (\omega E_\nu (q))|<\nu |b_q^+|vac>|^2$$
(4)
one defines the spectral function $`A(q,\omega )`$ . We calculate a one-particle spin Green function

$$<a(\tau )a^+(0)>=_0^{\mathrm{}}𝑑\omega \rho _t(\omega )\mathrm{exp}(\omega \tau ),$$
(5)
where $`\rho _t(\omega )`$ -spectral density function. Correlation functions are determined on the basis of a complete set of eigenvector of the Hamiltonian:

$$<O>=\frac{_\nu <\nu _i|O|\nu _j>}{_\nu <\nu _i|\nu _j>},$$
(6)
where $`O`$ is the longitudinal spin-spin correlation function $`<S_i^zS_{i+h}^z>`$ and the phonon density- density $`<n_{ph}(q)n_{ph}(q+p)>`$. The distribution of phonons number $`n_{ph}(q)`$ and the magnons number $`n_m(k)`$ as a function of momentum are simulated.

3. Results and discussion

The upper boundary of acoustic phonon band $`\omega _0`$ is in direct proportion to Debye temperature of $`\mathrm{\Theta }_{Sr_2CuO_3}=441\pm 10K`$  in term of Debye approximation. Using known data for $`CuGeO_3`$ , $`\mathrm{\Theta }_{CuGeO_3}=330K`$ , $`\omega _{0,Ge}1440K`$  we estimate $`\omega _02100K`$ for $`Sr_2CuO_3`$. The ratio of $`\omega _0/J0.95`$ will be used for calculation of dynamical characteristics . The spectral density of triplet excitations are shown in Fig.1. At the critical value of spin-phonon coupling $`\alpha _{c1}`$ the spectral density $`\rho _t=0`$ for $`E<\mathrm{\Delta }`$, were $`\mathrm{\Delta }`$ is the gap in the triplet excitation spectrum. The gaps energies presented in inset of Fig.1 fit satisfactory on the straight line $`\mathrm{\Delta }=0.17(2)+1.16(5)\alpha `$ in the range of $`\alpha _{c1}<\alpha <\alpha _{c2}`$ with $`\alpha _{c1}=0.15(1)`$.

For small spin-phonon coupling $`\alpha <\alpha _{c1}`$ the composite quasiparticle consisted spinon and phonon are formed. The inelastic interaction between spinon with momentum $`k`$ and phonon with $`q`$ lead to new quasiparticle with energy $`E=E_q+E_k`$ and momentum $`k+q`$. It follow from the distribution of the phonon number versus momentum as plotted in Fig.2a. The density of spinon excitations exhibits the singularity at the $`E(k=\pi /2)`$ and phonon excitations at the $`E(q=\pi )`$ and summered momentum agrees with MC result $`N_{ph}(q4.5)0`$. Quasiparticles with $`qk`$ are absent . It may be attributed to the repulsion spin-phonon quasiparticle with momentum $`qk`$ from spinon with the same energy $`E=E_kE_q`$. Whereas the density of spinon excitations tends to zero for energies $`E=E_q+E_k`$. The number of spin-phonon quasiparticles is increased against $`\alpha `$ as illustrated in Fig.2 and the bound states of spin-phonon quasiparticls are formed at $`\alpha >\alpha _{c1}`$ that are similar to bipolaron with the symbolic kind $`<b_{\text{q}_i^{n_1}k_i}^+b_{\text{q}+p_i^{n_2}k_i}^+_i^{n_1+n_2}S_{\text{k}_i}^z>`$.

The correlation function of the phonon numbers $`<N_{ph}(k)N_{ph}(k+q)>`$, presented in Fig.3a , is nonzero at the certain wavenumber $`q0.1`$ and $`q\pi /3`$ for $`\alpha \alpha _{c1}`$ and for all $`q`$ at the condition $`\alpha >\alpha _{c2}`$. The elastic strains cause nonuniform distribution of exchange and spin density that confirms spin-spin correlation function. The difference between spin-spin correlation function simulated in chain with spin-phonon coupling and without it is shown in Fig.3b. Small oscillations are observed for $`\alpha <\alpha _{c2}`$. Cycle of the oscillations $`\delta r=9`$ well agrees with cycle of oscillations of the magnons number distribution function $`\delta k=0.36,0.86;\delta r=\pi /\delta k9`$ plotted in Fig.4. The distribution of phonon numbers becomes continuum at the $`\alpha >\alpha _{c2}`$ with one sharp maximum at the wave vector of structure being in the range of $`\pi <Q<2\pi `$ as plotted in Fig.2b. The value $`\alpha _{c2}=0.4\pm 0.02`$ is similar to a critical concentration when the dressed phonons effectively "percolate " along the three low-lying bands. The correlation function $`<N_{ph}(k)N_{ph}(k+q)>`$ is not equal to zero for all momentums as shown in Fig.3.

The transition from one state with localized phonons to delocalized is like to Anderson's transition in disordered systems. The spin density in chain allocates random and nonuniform that follow from the normalized spin-spin correlation function shown in Fig.3b. The estimated values fit by exponential approximation $`<S_0^zS_r^z><S_0^zS_r^z>_{AF}A\mathrm{exp}(r/\xi )`$, where $`\xi `$ is a correlation radius of inhomogeneity. The magnons number distribution function exhibits also the random oscillations (Fig.4c). The transition to delocalized phonons is accompanied by strong increasing of the mean-square amplitude of ions vibration as presented in Fig.5b. At $`\alpha =\alpha _{c1}`$ the function $`<u^2>(\alpha )`$ reveals also a some peculiarity. The average number of phonons , plotted in Fig.5a , increases vs. $`\alpha `$ according to the power law

$$N_{av}=0,0012(\frac{\alpha }{\alpha _{c2}})^{1.75(6)}(7)$$

From the spectral density of phonons shown in Fig.6 and spin excitations it may be stated that the coherent nonlinear bound excitations of phonon and magnon exist at $`\alpha >\alpha _{c1}`$. It may be interpreted as a standing vibration with the spin kink in sites. The disjoint vibrations are performed at the condition that the wave length is changed in two times $`l/a=2;4;8;16;\mathrm{}`$ c $`k=\pi /l,E=w_0sin(k/2)`$. These estimates of energy are in good agreement with Monte Carlo results up to $`l/a=32`$ because the temperature fluctuations $`(T=0.04)`$ cut the kinks interaction radius.
Now we will apply our results to estimate spin-phonon coupling for the one-dimensional spin system $`Sr_2CuO_3`$ with intrachange $`J=2200K,\omega _0=2100K`$. The abrupt change of the bulk susceptibility at $`T10K`$ is associated with gap in the triplet excitation spectrum. Using typical ratio between gap and critical temperature, for example, for $`CuGeO_3,\mathrm{\Delta }/T_c1.8`$ we estimate $`\alpha =0.154`$ from eq. $`0.008=0.17+1.16\alpha `$. The average phonons number is $`N_{ph}=2.210^4`$. These phonons are formed bound states with energy about $`0.55eV,0.76eV,\mathrm{\hspace{0.33em}1}eV`$ as followed from the spectral phonon density plotted in Fig.6 because maximum phonon energy with $`n_q=1,E=0.95J0.18eV`$. This energy of elastic fluctuation is compared to the value of the experimental charge transfer gap $`E_g=1.5\pm 0.3eV`$ derived from He I ultraviolet photoemission combined with bremsstrahlung isochromat spectroscopy  and the on-site Coulomb interaction on oxygen site $`U_p=4.4eV`$  . The elastic fluctuation lead to change of the electron density on oxygen site and discontinuity of antiferromagnetic exchange that causes formation of the paramagnetic spin $`1.310^4`$ observed at low temperature $`4K`$.

There are spin-phonon quasiparticles with energy $`E_{s,ph}=1040cm^1,562cm^1,287cm^1`$,

$`144cm^1,72cm^1,36cm^1,18cm^1`$. The optic phonon mode observed at $`\omega _1=550cm^1`$ corresponds to stretching mode and mode with $`\omega _2=569cm^1`$ is attributed to spin-phonon quasiparticle with $`k=\pi /4,\omega ^{MC}=562cm^1`$. The relaxation time (or lifetime ) of the bound spin-phonon quasiparticle is decreased at heating and at $`TE^{MC}800K`$ it is decomposed that qualitatively agrees with experimental data. The phonon width of optical conductivity spectrum is increased from $`\gamma =17.4cm^1,T=10K`$ to $`\gamma =35.1cm^1,T=300K`$ and the strength of the oscillator is decreased from $`S=0.09cm^1eV,T=10K`$ to $`S=0.06cm^1eV,T=300K`$. The phonon width and the strength of the oscillator of mode $`\omega _1=550cm^1`$ don't change in the range $`10K<T<300K`$ .
The bond defects observed by NMR and heat transport arise from formation of the standing vibration sites. The density of bond defects $`\rho 1/l`$ is in direct proportional to wavevector of the standing vibration $`k=\pi /l`$ and estimated density is satisfactory fit with experimental data at $`T=20K,\rho ^{MC}510^3,\rho ^{ex}2.810^3`$. If discrete spectrum of spin-phonon quasiparticles map to continuum the temperature dependence $`\rho `$ may be presented as $`\rho 2T/\pi \omega _0`$ which well agree with NMR data . It is also explained the linear dependence of the specific heat at low temperatures. Interaction between loosely coupled phonon and spinon modes in the center of Brillouin band may be presented as $`(k\omega /v_{ph})(k\omega /v_{sp})=\stackrel{~}{\alpha }^2`$, where $`v_{ph}=\omega _0/2,v_{sp}=\pi /2J,\stackrel{~}{\alpha }=\alpha /\sqrt{\omega _0}`$. It lead to repulse of these branches and formation of gap in the spinon spectrum with $`\mathrm{\Delta }(k=0)/J=\alpha \sqrt{\pi }/2\sqrt{J}0.019`$. This gap $`\mathrm{\Delta }(k=0)=43K`$ causes abrupt decrease of the dynamical spin susceptibility for $`k=0`$ at $`T33K`$. The $`k=0`$ mode of spin transport has a diffusion like contribution for $`T<<J`$ as a result of the fact that the finite density of bound spin-phonon quasiparticles exists at $`E=0`$ .

At low temperature the spin liquid with non-linear topological excitations having $`S=1/2`$ is freezed that lead to small magnetic moment on site and elastic neutrons scattering. So in compound $`CuGeO_3`$ having transition singlet-triplet at $`T=14K`$ a long range order is observed at low temperatures when ions $`Cu`$ or $`Ge`$ are displaced by $`Mg`$ or $`Sr`$. Breaking a exchange between the nearest neighbors results in a finite spin density near breaking. Between the generated magnetic moments the interaction is existed per phonon field in 3D dimensional space. The critical field may be estimated by mean field (MF) approximation $`T_N=2Sc(J_0+2J_1)`$, where $`c`$ impurities concentration, $`S=1/2`$ , $`J_0`$ and $`J_1`$ is accordingly intra- and interchain exchange in simple cubic lattice. The calculated critical temperatures for $`Cu_{1x}Mg_xGeO_3`$ are in well agreement with experimental data $`T_N^{MF}=2.4K,T_N^{ex}=2.5K,x=0.016;T_N^{MF}=3.2K,T_N^{ex}=3K,x=0.0216`$  . Using similar evaluation for $`Sr_2CuO_3`$ with concentration of bond defect at $`T6K,c210^3`$ we determined the critical temperature $`T_N^{MF}4.5K`$ which is in well agreement with experiment $`3.5K<T_N^{ex}<5.4K`$. Scatter of temperatures may be due to dependence of freezing bond defects from frequency. The average distance between bond defects is $`210^3A`$ and the periodic local alternating magnetization exists which is observed by NMR at low temperatures. Mean square displacement is $`<U^2>10^6`$ that may give rise the structural deformation along chain $`\delta a\sqrt{<U^2>}10^3A`$ for $`a=3.9A,M10^{22}g`$.

In summary, the spectrum of spinon excitations along chain is asymmetric relative wavenumber $`k=\pi /2`$ with two gaps $`\mathrm{\Delta }(k=0)=43K`$ and $`\mathrm{\Delta }(k=\pi )=18K`$ that account for steeply decrease of the bulk susceptibility and the dynamical spin susceptibility for $`k=0`$ at the temperatures $`T10K,T33K`$. The spectrum of bound spin-phonon excitations is due to bond defects the density of which is linearly increased against temperature. At low temperature $`T5K`$ the non-linear spin-phonon excitations are freezed on the breaking exchange and induce magnetic moment on site with correlation radius of antiferromagnetic ordering $`12002000A`$. The magnetic properties of $`Sr_2CuO_3`$ are similar to $`CuGeO_3`$ with impurities having two transitions on temperature. The paramagnetic contribution in the bulk susceptibility at $`T<4K`$ result from a structural defects on oxygen ions due to spin-phonon coupling.

The author are grateful to V.A. Kashurnikov, A.S. Mishchenko, N.B. Prokof'ev, B.S. Svistunov for useful assistance and discussions of Monte Carlo method in the continuous time.

Captions to Aplesnin paper " Quantum spin liquid in $`Sr_2CuO_3`$ induced by Acoustic Phonons "

Fig.1 The spectral density of the one-particle excitations with $`\omega _0=0.95,\alpha =0.05(1),0.15(2)(a)`$; $`\alpha =0.2(1),0.4(2)(b)`$. In insert, the gap energy $`\mathrm{\Delta }`$ in the excitation spectrum as a function of spin-phonon coupling .

Fig.2 Distribution of phonon numbers versus momentum for $`\alpha =0.05(1),0.2(2)(a)`$; $`\alpha =0.4(1),0.6(2)(b).`$

Fig.3 The normalized correlation function of the phonon numbers versus wavenumber for $`\alpha =0.2(1),0.4(2),0.6(3)`$$`(a)`$ and the spin-spin correlation function as a function of distance for $`\alpha =0.2(1),0.6(2)`$ $`(b)`$ . Solid line is shown the exponential dependence $`Aexp(r/\xi ),A=0.008(8),\xi =9.1(9)`$.

Fig.4 Distribution of magnon numbers versus momentum for $`\alpha =0.05,0k\pi (a)`$; $`\alpha =0.2,0<k<\pi (b);\alpha =0.6,0<k<\pi (c).`$

Fig.5 The average phonons number $`<N_{ph,av}>(a)`$ and the mean-square displacement of ions $`<u^2>(b)`$ versus $`\alpha /\alpha _{c2}`$.

Fig.6 The normalized spectral density of phonon excitations $`\rho _{ph}/\rho _{ph,max}`$ versus energy for the spin-phonon coupling $`\alpha =0.25(a);\alpha =0.4(b);\alpha =0.7(c)`$.

$`|`$


text:
# Analysis of directed flow from three-particle correlations

## 1 INTRODUCTION

The measurement of the successive harmonics of azimuthal correlations in a heavy ion collision , $`v_ne^{in(\varphi \mathrm{\Phi }_R)},`$ where $`\mathrm{\Phi }_R`$ is the impact parameter direction, is of paramount importance, since it yields information on the medium created in the collision. In particular, a novel behaviour has been predicted for the first harmonic $`v_1`$, the so-called directed flow, at ultrarelativistic energies .

However, the analysis of directed flow at these energies is highly nontrivial, because $`v_1`$ is very small. Thus, methods based on an analysis of two-particle azimuthal correlations  are likely to be biased by “nonflow” two-particle correlations: quantum (HBT) correlations between identical particles and correlations due to global momentum conservation have been shown to be important at SPS , while correlations due to minijets may be large at RHIC . On the other hand, methods relying on four-particle correlations , which are free from this bias, are plagued by a lack of statistics due to the smallness of $`v_1`$, although they give good results for the analysis of elliptic flow $`v_2`$ (see Sec. 2). To remedy these shortcomings, we proposed a new method of $`v_1`$ analysis , based on the measurement of a mixed three-particle correlation, which involves both $`v_1`$ and $`v_2`$:

$$e^{i(\varphi _a+\varphi _b2\varphi _c)}(v_1)^2v_2,$$
(1)
where $`\varphi _a`$, $`\varphi _b`$, and $`\varphi _c`$ denote the azimuths of three particles belonging to the same event, and the average runs over triplets of particles emitted in the collision, and over events. Once $`v_2`$ has been obtained from a separate analysis, this equation yields $`(v_1)^2v_2`$, thus $`v_1`$.

Here, we apply this method to NA49 data on Pb-Pb collisions at 158 $`A`$GeV. Results obtained using the “standard” flow analysis are given in Ref. . In our method, the first step in the analysis is the measurement of a reference $`v_2`$. Then, an equation analogous to Eq. 1 yields an estimate of the integrated $`v_1`$, i.e., its average value over some phase space region. Finally, restricting $`\varphi _1`$ in Eq. 1 to a small $`(p_T,y)`$ bin allows one to obtain more detailed, differential measurements of $`v_1`$ as a function of transverse momentum or rapidity.

## 2 ELLIPTIC FLOW FROM 2, 4, 6, 8-PARTICLE CORRELATIONS

As stated in the introduction, our method of analysis of directed flow $`v_1`$ requires the preliminary knowledge of an estimate of the elliptic flow $`v_2`$, integrated over some phase space region. Of course, this estimate must be obtained by analysing the same sample of events from which one wants to extract $`v_1`$.

In practice, the average over phase space is a weighted average:

$$w_2v_2w_2e^{2i(\varphi \mathrm{\Phi }_R)},$$
(2)
where $`w_2`$ is the chosen weight. In order to reduce statistical fluctuations, $`w_2`$ must be larger for particles with stronger elliptic flow.

The value of $`w_2v_2`$ is obtained using the cumulant method described in Ref. : one can extract estimates of $`w_2v_2`$ from cumulants of multi- (2-, 4-, 6-…) particle correlations. While two-particle methods are equivalent to the standard flow analysis, higher orders are essentially free from nonflow effects. They were first used in analysing data obtained by the STAR Collaboration at RHIC .

In Fig. 2, we present application of the method to NA49 data, and show the dimensionless quantity

$$v_2\frac{w_2v_2}{\sqrt{(w_2)^2}},$$
(3)
as a function of centrality for charged particles, where we have used $`w_2=p_T`$. We display estimates using cumulants of two-, four-, six-, and even eight-particle correlations , as well as the corresponding quantity for charged pions obtained from the “standard” subevent method .

It is quite remarkable that the four-, six-, and eight-particle estimates all agree: this supports the idea that they are indeed free from nonflow effects, and correspond to a genuine collective motion in the direction of the impact parameter. Moreover, these multiparticle estimates show a slight discrepancy with the two-particle values, as expected if nonflow correlations are sizable . Please note that the statistical uncertainties on high order cumulants remain reasonably small, especially for midcentral collisions. In the following, our reference elliptic flow value will preferably be the estimate from the four-particle cumulant, which is a priori the most reliable since it is free from nonflow correlations and has a smaller statistical error than higher order estimates.

## 3 INTEGRATED DIRECTED FLOW FROM 2, 3, 4-PARTICLE CORRELATIONS

The next step in the analysis is to determine the average value of directed flow over some phase space region. As in the case of elliptic flow, we perform a weighted average

$$w_1v_1w_1e^{i(\varphi \mathrm{\Phi }_R)},$$
(4)
where stronger weight is given to particles with stronger directed flow. In the NA49 analysis, we used a rapidity dependent weight $`w_1=yy_{\mathrm{CM}}`$, where $`y_{\mathrm{CM}}`$ is the centre-of-mass rapidity.

This weighted average is obtained from the following three-particle correlation:

$$w_1(a)w_1(b)w_2(c)e^{i(\varphi _a+\varphi _b2\varphi _c)}=w_1v_1^2w_2v_2,$$
(5)
where $`w_2`$ is the same as in Sec. 2. Using the value of $`w_2v_2`$ obtained in Sec. 2, we thus derive $`w_1v_1`$, up to a global sign.

In practice, the left-hand side of Eq. 5 is constructed using a generating function formalism detailed in Ref. . This procedure is a very efficient way to sum over all possible triplets of particles, and also to remove automatically the effects of slight detector anisotropies.

The method was applied to NA49 data. In Fig. 2, we present as a function of centrality the dimensionless quantity $`v_1w_1v_1/\sqrt{(w_1)^2}`$, analogous to Eq. 3, together with two- and four-particle estimates obtained with the method of Ref. . It is worth noting that the statistical uncertainty on the three-particle estimate is barely larger than that on the two-particle value, while the systematic error due to nonflow correlations (not included in the plot) is a priori much smaller.

## 4 DIFFERENTIAL DIRECTED FLOW

Finally, we want to obtain detailed measurements of $`v_1`$, as a function of $`p_T`$ or $`y`$ and for each particle type. To derive $`v_1(p_T,y)`$, we use a three-particle correlation analogous to Eq. 5, where the particle labeled $`a`$ is restricted to a narrow $`(p_T,y)`$ bin:

$$w_1(b)w_2(c)e^{i(\varphi _a+\varphi _b2\varphi _c)}=w_1v_1w_2v_2v_1(p_{Ta},y_a).$$
(6)
With the previously derived $`w_2v_2`$ and $`w_1v_1`$, this equation yields $`v_1(p_T,y)`$. We illustrate this method on the differential flow of charged pions in midcentral collisions, shown in Figs. 4 and 4 together with two two-particle estimates, either uncorrected or corrected for the effect of momentum conservation .

At high $`p_T`$, $`v_1`$ from three-particle correlations is consistent with the two-particle value corrected for momentum conservation, but significantly lower than the uncorrected one. This shows that correlations from momentum conservation, which are large, are automatically removed in our method. This is also reflected by the behaviour at midrapidity, where $`v_1`$ vanishes, as it should by symmetry, while the uncorrected two-particle estimate does not. The three-particle estimate also vanishes more smoothly at $`p_T=0`$ than two-particle estimates which may be biased by HBT correlations .

## ACKNOWLEDGEMENTS

We thank the NA49 Collaboration for permission to use their data.


text:
# Stability of gap solitons in a Bose-Einstein condensate

## I Introduction

One of the most fundamental facts about Bose-Einstein condensates (BEC) of dilute atomic gases is that they can be very well described by the Gross-Pitaevskii equation (GPE)

$$i\mathrm{}\dot{\psi }(𝐱,t)=\left(\frac{𝐩^2}{2M}+V(𝐱)+\kappa |\psi (𝐱,t)|^2\right)\psi (𝐱,t)$$
(1)
where $`\psi `$ denotes the collective wavefunction of condensed atoms with mass $`M`$ and coupling constant $`\kappa :=4\pi \mathrm{}^2aN_A/M`$, with $`a`$ being the s-wave scattering length and $`N_A`$ denoting the number of atoms in the BEC (see, e.g., Ref. castin01 ). Therefore a BEC provides a physical realization of many nonlinear wave phenomena among which the formation of solitons is particularly interesting.

For our purposes solitons are wavepackets in which the dispersive effect of the kinetic term is exactly cancelled by the nonlinear interaction energy so that their shape does not change. Two fundamental types of solitons have so far been experimentally realized in a BEC: dark solitons sengstock99 ; denschlag00  correspond to a stable density dip in a BEC of repulsive atoms. Bright solitons do exist for atoms with attractive interaction ($`\kappa <0`$) and are described by a one-dimensional solution of the GPE for vanishing potential $`V(𝐱)`$,

$$\psi _{\text{bright}}(z,t)=\frac{1}{\sqrt{2\text{w}}}\text{sech}(z/\text{w})e^{i\omega _st},$$
(2)
where $`\text{w}=2\mathrm{}^2/(M\kappa _{1D})`$ and $`\omega _s=\mathrm{}/(2M\text{w}^2)`$. To relate the one-dimensional GPE to its three-dimensional origin we assumed that the BEC is tightly trapped in the transverse direction so that no transverse excitations can occur. Taking the BEC to be in the transverse ground state amounts to replacing the three-dimensional coupling constant $`\kappa `$ by $`\kappa _{1D}=\kappa /A_{}`$, where $`A_{}=2\pi a_{}^2`$ with $`a_{}:=\sqrt{\mathrm{}/(M\omega _{})}`$ is the harmonic oscillator length of the transverse potential. The attraction between the atoms prevents the dispersion of the sech-shaped wavepacket. Note that for $`\kappa _{1D}>0`$, corresponding to repulsive atom-atom interaction, the kinetic energy and the interaction energy have the same sign and therefore cannot cancel each other, a bright soliton is then not possible. Very recently bright solitons have been created in a quasi one-dimensional setup salomon02 ; hulet02  where the transverse potential $`V(x,y)`$ tightly confines the BEC, thus suppressing transverse excitations and three-dimensional collapse.

In this paper we are concerned with the dynamical stability of gap solitons. This collective state, which has not yet been realized experimentally, exists for repulsive atoms ($`\kappa _{1D}>0`$) in a periodic potential $`V(z)`$ and is related to bright solitons. The basic idea of a gap soliton is the following: as is well known the energy eigenvalues of noninteracting particles in a periodic potential are given by energy bands $`E_n(q)`$, where $`n`$ is the band index and $`q`$ denotes the quasi momentum. If a particle’s state is prepared in the lowest band only the dispersion relation $`p^2/(2M)`$ in free space is replaced by the lowest band energy $`E_0(q)`$. Around the upper band edge, which we take to be at $`q=0`$, this energy can be approximated by $`E_0(q)E_0(0)+q^2/(2M^{})`$ where $`M^{}:=(d^2E_0(q)/dq^2)^1|_{q=0}`$ is the effective mass of the particle in the periodic potential. Since at the upper band edge $`M^{}<0`$ the “kinetic energy” becomes negative and a cancellation with the positive interaction energy becomes possible. The corresponding state is called a gap soliton.

Gap solitons have been realized in nonlinear optics using a periodic modulation of the propagation medium chen87 . In nonlinear atom optics they have first been predicted by Lenz et al. lenz94  in the context of light-induced nonlinearities krutitsky99 ; krutitsky02 . Here we are concerned with the collision-induced nonlinearity appearing in Eq. (1).

We consider a BEC that is placed in a one-dimensional optical lattice, created by far-detuned laser light with wavenumber $`k_L`$, which produces an optical potential of the form $`V(z)=V_0\mathrm{cos}(2k_Lz)`$, and is subject to a tight harmonic transverse confinement of the form $`V_{}(x,y)=M\omega _{}^2(x^2+y^2)/2`$. Here $`\omega _{}`$ is the transverse trap frequency and $`V_0`$ the strength of the optical potential. The derivation of the corresponding gap soliton solution of the one-dimensional GPE is tedious and includes a multiple scales analysis. For nonlinear optics it has been derived by Sipe et al. sipe88 . For nonlinear atom optics a related derivation has been sketched in Refs. lenz94 ; steel98 ; kon01 . One finds that, within the effective mass approximation, the one-dimensional gap soliton is described by

$$\psi _{\text{gap}}(z,t)=\phi _{\text{be}}(z,t)\frac{1}{\sqrt{2\stackrel{~}{\text{w}}}}\text{sech}(z/\stackrel{~}{\text{w}})e^{i\stackrel{~}{\omega }_st},$$
(3)
where $`\phi _{\text{be}}(z,t)=\phi _{\text{be}}(z)e^{iE_0(0)t/\mathrm{}}`$ is the solution of the linear Schrödinger equation that corresponds to the upper band edge. The parameters $`\stackrel{~}{\text{w}}`$ and $`\stackrel{~}{\omega }_s`$ have the same form as w and $`\omega `$ for the bright soliton but with $`M`$ replaced by $`|M^{}|`$ and $`\kappa _{1D}`$ replaced by $`\stackrel{~}{\kappa }_{1D}:=\kappa _{1D}|\phi _{\text{be}}(z)|^4𝑑z`$. Apart from $`\phi _{\text{be}}`$ the solution (3) just corresponds to a bright soliton for a particle of mass $`|M^{}|`$ and coupling constant $`\stackrel{~}{\kappa }`$.

The range of experimentally promising values for the optical potential strength $`V_0`$, the number of condensed atoms $`N_A`$, and the transverse confinement frequency $`\omega _{}`$ has been examined in a study by Brezger et al. brezger02  following standard textbooks agrawal89 . It was found that values around $`N_A=400`$, $`V_0=\mathrm{}^2k_L^2/(2M)`$ and $`\omega _{}=110`$ s<sup>-1</sup> should be optimal for the observation of gap solitons. In the following we will focus on this case. We will also consider the effect of a variation of $`N_A`$ which allows to test the quasi one-dimensionality of the gap soliton.

## II Stability theory of gap solitons and bright solitons

A stationary solution $`\psi _0`$ of the GPE, with chemical potential $`\mu `$, is called dynamically stable if a small deviation $`\delta \psi `$ from $`\psi _0`$ will not grow with time. In this case a small perturbation will not cause the solution to evolve into a completely different wavepacket. To study dynamical instability one can either directly integrate the GPE or solve the associated Bogoliubov-de Gennes equations (BDGE, see, e.g., Ref. lewenstein96 ). The latter arise when one writes the wavefunction $`\psi (𝐱,t)`$ in the form

$$\psi (𝐱,t)=\mathrm{exp}(i\mu t/\mathrm{})(\psi _0(𝐱)+\delta \psi (𝐱,t))$$
(4)
and linearizes the GPE in $`\delta \psi `$. By making the ansatz of a stationary perturbation,

$$\delta \psi (𝐱,t)=u(𝐱)\mathrm{exp}(i\omega t)v^{}(𝐱)\mathrm{exp}(i\omega t),$$
(5)
one arrives at the BDGE

$`\mathrm{}\omega u`$ $`=`$ $`u\kappa \psi _0^2v`$
$`\mathrm{}\omega v`$ $`=`$ $`v\kappa (\psi _0^{})^2u`$ (6)
with $`:=𝐩^2/(2M)+V(𝐱)\mu +2\kappa |\psi _0|^2`$. A solution $`(u,v)`$ with eigenvalue $`\omega `$ corresponds to a quasi-particle mode. The set of all $`\omega `$ forms the quasi-particle spectrum which in general is complex. One can show castin01  that if $`\omega `$ is in the spectrum then so is $`\omega ^{}`$. Using Eqn. (5) it is seen that the existence of a nonzero imaginary part of one quasi-particle frequency implies exponential growth of the mode and hence dynamical instability of the state $`\psi _0`$. To demonstrate that the gap soliton is stable therefore amounts in showing that the associated quasi-particle spectrum is real.

## III Numerical methods

Although the direct numerical integration of the GPE is easy to implement using the split-step method deVries86  it has the disadvantage of not being practical for a 3D study of the gap soliton. The reason is that the state (3) includes two very different spatial scales: the laser wavelength $`2\pi /k_L`$ and the width $`\stackrel{~}{\text{w}}`$ of the soliton’s envelope. To cover both scales simultaneously it was necessary to consider at least 260 periods of the optical potential or 2000 spatial points along the z-axis. Since the number of total points in 3D is restricted by the capacity of the computer a 3D simulation of a gap soliton becomes impractical. We therefore have applied this method only to 1D and 2D simulations.

To derive the spectrum of quasi-particles around the gap soliton we followed the method of Ref. lewenstein96  and expanded the modes $`(u,v)`$ of Eq. (6) as well as the gap soliton wavefunction $`\psi _0`$ in a set of basis functions. For the z-direction we have chosen a number of $`n_z`$ Bloch wavefunctions which are eigenstates of the linear Schrödinger equation with potential $`V(z)=V_0\mathrm{cos}(2k_Lz)`$. As transverse modes we used $`n_xn_y`$ harmonic oscillator eigenstates corresponding to the transverse trapping potential. This turns Eq. (6) into an eigenvalue problem for a $`2n\times 2n`$ matrix, where $`n:=n_xn_yn_z`$ is the total number of mode functions. However, before this eigenvalue problem can be solved, one first has to find the exact wavefunction of the gap soliton in the given set of basis modes (in practice a reasonably large subset is sufficient).

Since the analytical solution (3) is only approximately correct a stability analysis will inevitably lead to a complex quasi-particle spectrum because of the finite difference to the exact solution. To find the exact solution $`\psi _{\text{exact}}`$ we have used a self-consistent field approach (SCF): after expansion of the GPE in the basis modes Eq. (1) is turned into a set of coupled nonlinear algebraic equations for the expansion coefficients of $`\psi _0`$. Using the approximate solution (3) as an ansatz $`\psi _{\text{trial},1}`$ we insert it into the GPE to evaluate the nonlinear terms. The resulting equation,

$$E\psi (𝐱)=\left(\frac{𝐩^2}{2M}+V(𝐱)+\kappa |\psi _{\text{trial}}(𝐱)|^2\right)\psi (𝐱),$$
(7)
represents a linear eigenvalue problem for $`\psi `$ and can easily be solved using standard numerical methods. We pick that solution $`\psi _{\text{trial},2}`$ out of all eigenstates of Eq. (7) which has the least deviation from our previous guess $`\psi _{\text{trial},1}`$ and iterate this procedure until the change in the trial wavefunction is below a given value (we used a relative change of $`10^{14}`$ as accuracy goal). The converged wavefunction may then correspond to the true gap soliton.

## IV One-dimensional results

In one dimension convergence of the SCF algorithm can easily be achieved. We have used up to 180 Bloch modes to expand the exact soliton wavefunction which covered the upper half of the lowest energy band and the lower third of the second Bloch energy band. However about 40 modes covering the effective-mass region around the upper band edge were sufficient to get about the same numerical accuracy. The (real) expansion coefficients for the gap soliton can be seen in Fig. 1. It is interesting to note that the SCF method essentially amounts to removing that part of the approximate solution (3) which corresponds to the second energy band. The exact numerical solution therefore is indeed centered around the upper band edge of the first energy band where the effective mass is approximately constant. We have checked whether the final result of the SCF algorithm indeed describes a soliton by using it as initial condition for the split-step direct integration of the GPE. It was found that the solution does not change its shape for times exceeding 0.2 s.

To verify the stability of the 1D gap soliton we have calculated the quasi-particle spectrum for up to 350 Bloch modes as basis functions. The result for the real part of the spectrum can be seen in Fig. 2 together with the Bloch mode energies. As one can see the two spectra are very similar apart from a constant shift. This shift is given by the chemical potential and arises because of the corresponding phase factor in ansatz (4) The similarity of the curves is due to the small number of condensed atoms ($`N_A=400`$) and the correspondingly small collision effects.

The imaginary part of the quasi-particle spectrum is zero except for a single mode (and the corresponding mode with complex conjugated frequency) for which the frequency is purely imaginary. The value of the imaginary part for this mode depends on the degree of convergence of the SCF wavefunction and on the number of modes included. For the analytical, non-converged solution (3) the decay factor is about $`0.02\times \omega _{}`$, where $`\omega _{}`$ is the transverse trap frequency. For a well-converged solution it is always very small, below $`0.001\times \omega _{}`$.

To understand the unstable mode better a few general facts about the quasi-particle spectrum are helpful: for any potential and any stationary solution of the GPE the mode $`(u,v)=(\psi _0,\psi _0^{})`$ is a quasi-particle mode with frequency $`\omega =0`$. This mode is a Goldstone mode associated with the symmetry of the energy functional

$$E=\left\{\psi _0^{}\left(\frac{𝐩^2}{2M}+V\right)\psi _0+\kappa |\psi _0|^4\right\}d^nx$$
(8)
with respect to a global phase change $`\psi ^{}=\mathrm{exp}(i\alpha )\psi `$. If the external potential $`V`$ is absent then, for any stationary solution $`\psi _0`$ of the GPE, there is a second Goldstone mode $`(u,v)=(\psi _0,\psi _0^{})`$. It arises because of the invariance of $`E`$ against spatial translations. If the effective mass approximation would be exact then the gap soliton would fulfill the same equation as the bright soliton does in free space. Consequently, it would also possess a translational Goldstone mode. However, since the periodic optical potential does explicitly break translational invariance one can expect a shift of the complex Goldstone frequency away from zero. A second effect that explicitly breaks translational invariance is the finite number of Bloch basis functions used in the numerical calculations. This is equivalent to an optical lattice placed in a box whose length is a finite multiple of the lattice period. In our case, the box contained 260 periods.

That the unstable mode indeed corresponds to a modified Goldstone mode can also be seen by looking at its expansion coefficients (Fig. 3). Obviously the shape of $`u(q)`$ is approximately given by $`q\psi _0(q)`$ which would describe the Goldstone mode if the quasi momentum $`q`$ is replaced by the real momentum as it is done in the effective mass approximation. The question remains whether this tiny instability results from numerical aberrations, from the finite number of periods, or whether it is a real physical effect. To shed some light on this question we also have performed a stability analysis of the bright soliton in free space using the same algorithm. We found a similar behaviour: a Goldstone mode develops an imaginary eigenvalue, but this time it is the mode associated with the phase transformation. Since one can prove that this mode has zero frequency we conclude that the tiny imaginary eigenvalue for Goldstone modes is a spurious numerical effect and the gap soliton is dynamically stable in one dimension.

## V Results in two and three dimensions

Having examined the stability of the 1D gap soliton it is of interest whether it will remain stable in a quasi one-dimensional situation. The condition for the latter is usually formulated as follows: the interaction energy, which leads to a coupling between different modes of the corresponding linear Schrödinger equation, should be much smaller than the excitation energy of the transverse trapping potential. If this is fulfilled a BEC in its ground state will effectively behave like a one-dimensional quantum gas since transverse excitations are off-resonant and therefore suppressed.

This is not the case for a gap soliton, however. The reason can be seen by looking at Fig. 4 which displays the mode energies of noninteracting atoms in the optical lattice and with a tight harmonic transverse confinement around the upper band edge. The solid line displays the energy of atoms in the transverse ground state and with longitudinal quasi momentum $`q`$ in the lowest energy band. Each dashed line corresponds to transversally excited atoms with a transverse energy of $`2\mathrm{}\omega _{}`$ to $`6\mathrm{}\omega _{}`$, respectively. It is important to observe that there are resonances between transversally unexcited atoms with $`q=0`$ and transversally excited atoms with $`q0`$. Since the gap soliton is a superposition of Bloch modes around the upper band edge ($`q=0`$), these resonances have the consequence that even for tight transverse confinement a true gap soliton does not exist. This situation is qualitatively the same in two and in three dimensions since in both cases the free energy levels are given by those of Fig. 4. Since in two dimensions the transverse trapping potential is one-dimensional the multiplicity of the energy levels is always one. This is different in three dimensions where a transverse excitation energy of $`n\mathrm{}\omega _{}`$ has an $`n+1`$ fold degeneracy. The number of resonant states is therefore larger than in two dimensions.

Although a true gap soliton does not exist it is of interest to examine a quasi gap soliton of the form

$$\psi _{\text{quasi}}=\psi _{\text{gap}}(z,t)\phi _0(x,y),$$
(9)
with $`\phi _0`$ denoting the transverse ground state. It should be possible to produce a state like $`\psi _{\text{quasi}}`$ using dispersion management treutlein02 . Though the quasi gap soliton is not a true stationary solution of the GPE, it may be sufficiently stable to allow for experimental observation. A signature of a quasi gap soliton would be a strongly suppressed dispersion of the wavepacket along the z axis. To analyze the time scale on which this state decays we first note that the transverse ground state has even parity and because of parity conservation can only couple to even excited levels $`2n\mathrm{}\omega _{}`$. This is the reason why we omitted odd transverse excitations in Fig. 4.

We have used the state (9) as initial condition and studied numerically the time evolution of it in two dimensions. To study the influence of the transverse confinement we have considered three BECs with 400, 1600, and 25600 atoms in Figs. 5 A, B, and C, respectively. The number of atoms $`N_A`$ and the transverse confinement frequency $`\omega _{}`$ have been simultaneously varied keeping the product $`\omega _{}N_A`$ constant. Consequently the interaction energy in Eq. (8) is kept constant, since it is proportional to $`(\omega _{}N_A)^2`$. This also ensures that the first-order soliton condition is fulfilled. The result of the numerical time evolution after 23 ms is shown in Fig. 5. While some excitations are observable, the state still looks very much like a quasi gap soliton. This situation changes after 0.23 s (Fig. 6). Here one can see strong excitations which are growing with decreasing transverse excitation frequency. This is a reasonable result since the ratio of the interaction energy and the transverse excitation energy is greater than one for Fig. 6 B and C hence the nonlinear coupling is strong enough to excite the transverse modes. The figure suggests that in particular the transverse modes with $`2\mathrm{}\omega _{}`$ energy are strongly excited because for a fixed value of $`z`$ there are three density maxima in the x-direction. This is in agreement with the qualitative predictions which we have made above.

To verify this result we also have calculated the quasi-particle spectrum of the state (9) in two and three dimensions. In principle, since the quasi gap soliton is not a stationary state, the spectrum is not enough to predict the evolution of it accurately and a more sophisticated approach is needed castin97 . However, to gain a qualitative understanding of the time scale on which $`\psi _{\text{quasi}}`$ decays the imaginary part of the spectrum is sufficient.

In 2D we used up to 100 Bloch wavefunctions and up to 19 one-dimensional eigenstates of the transverse harmonic trap as basis modes. The quasi gap soliton was expanded using 51 Bloch states. It turned out that there are generally quite many unstable modes, but only few of them do have a considerable overlap with the collective wavefunction. The number of unstable modes depends on the number of basis states used to expand the BDGE since the number of resonant transversally excited states is growing. However, it turned out that this basis dependence does only affect modes with a small instability. Some examples are displayed in Fig. 7. Mode 7A corresponds very well to the anticipated resonant excitation of transverse modes. It has a non-zero overlap with the quasi gap soliton and otherwise only populates even transversely excited basis modes. Correspondingly its instability is rather large; only the instability of mode 7B, which roughly describes the phase Goldstone mode, decays faster. Mode 7B is unstable because, as described above, the quasi gap soliton is not a stationary solution of the GPE. However, due to a coupling to transversally excited states the decay rate is strongly enhanced as compared to the non-converged one-dimensional quasi gap soliton (the modulus of the transversally excited mode coefficients is less than 0.1 and is therefore not visible in Fig. 7). Mode 7C is typical for the many unstable modes which have no overlap with the quasi gap soliton. Therefore, unless transversal excitations are created during the experimental preparation of the quasi gap soliton, they do not contribute to the decay of it. These modes appear because the collective wavefunction provides a linear coupling term between transversally excited states which are in resonance with each other. Such modes can exist for even and odd transverse excitation levels without violating parity. Finally, mode 7D describes a somewhat off-resonant coupling between the transverse ground state and transversally excited states. Because of its off-resonant nature its decay rate Im($`\omega `$) is considerably smaller than for the resonant mode 7A.

To analyse the dynamical instability of a 3D quasi gap soliton we used again Bloch wavefunctions for the longitudinal expansion of the Bogoliubov modes. For the two transversal directions we have chosen a basis of states which are both eigenstates of the Hamiltonian and the angular momentum operator $`L_z`$. These states can be constructed by using the creation operators $`c_\pm :=(a_x^{}\pm ia_y^{})/\sqrt{2}`$ marzlin98 . The basis states are then given by

$$|n,m:=\frac{(c_+^{})^{(nm)/2}(c_{}^{})^{(n+m)/2}}{\sqrt{((n+m)/2)!((nm)/2)!}}|0.$$
(10)
The energy of these states is given by $`n\mathrm{}\omega _{},n=0,1,2,\mathrm{}`$ and their angular momentum by $`m\mathrm{},m=n,(n2),\mathrm{},n`$. We used again 51 Bloch states to expand the Bogoliubov modes along the z-axis and up to 30 transverse modes as given in Eq. (10).

We found again a large number of unstable modes which spuriously depends on the number of transverse basis states. However, similarly to the 2D case, modes with an instability Im$`(\omega )>0.02\omega _{}`$ do not depend significantly on the number of basis states. Some of these unstable Bogoliubov modes in 3D are shown in Fig. 8. Mode 8A is again roughly the phase Goldstone mode plus some small transversal excitations not visible in the figure. Due to the larger number of nearly resonant basis states the decay rate is even larger than in the 2D case. Mode 8B is again a superposition of resonant states and, as in 2D, has the second largest decay rate. These two states as well as most of the other unstable modes are solely composed of states with zero angular momentum, i.e., states which are symmetric in the x-y plane. They are therefore very similar to the 2D case. A new kind of instability in 3D occurs in Mode 8C which has no overlap with the quasi gap soliton. It is composed out of states with $`m=\pm 1`$ and thus describes the decay or increase of rotating states linearly coupled by the collective wavefunction. We found 80 unstable modes of this kind.

## VI Conclusion

We have examined the dynamical instability of gap solitons in a BEC in one, two, and three spatial dimensions under the condition of tight transversal confinement. Using different methods we found that a truly one-dimensional gap soliton is stable. In higher dimensions transversal excitations which are resonant with the upper band edge forbid the existence of a real gap soliton. However, a quasi gap soliton may be experimentally prepared which behaves like the 1D gap soliton for a time smaller than the smallest decay time of one of the unstable Bogoliubov modes. In 3D we numerically found this time to be in the order of $`1/(0.133\times \omega _{})`$ which is the decay time of mode A in Fig. 8. For a transverse trap frequency of $`\omega _{}=110`$ s<sup>-1</sup> we expect the quasi gap soliton to be sufficiently stable for about 70 ms. This should be long enough for experimental observation.

Acknowledgement: We are indebted to Klaus Mølmer and Gora Shlyapnikov for very valuable discussions and Jürgen Audretsch for kind hospitality. This work was supported by the Optikzentrum Konstanz and the Forschergruppe Quantengase. M.K.O. is supported by the Emmy Noether Programm of the Deutsche Forschungsgemeinschaft. K.M.H. would like to thank the following fonds for financial support: Familien Hede Nielsens Fond, Krista og Viggo Petersens Fond, Observator mag. scient Julie Marie Vinter Hansens Fond, Otto Bruuns Fond Nr. 2, Etatsraad C.G. Filtenborg og hustru Marie Filtenborgs studielegat, Rudolph Als Fondet, Frimodt-Heineke Fonden, Sokrates/Erasmus studierejselegat and Professor Jens Lindhards Forskningslegat.


text:
# Hadronic Spectral Functions above the QCD Phase Transition Presented by M. Asakawa

## 1 INTRODUCTION

The spectral functions (SPFs) of hadronic operators play an important role in QCD. The modification of hadrons, which has been suggested, for example, by the dilepton production enhancement in high energy nuclear collisions observed at CERN SPS, can be formulated in terms of SPFs. Although there have been numerous theoretical attempts to understand the modification of SPFs at finite temperature and/or density , the exact nature of the hadronic modes in matter is not understood well.

Recently, the maximum entropy method (MEM) has been used to extract SPFs from lattice QCD data for the first time by the present authors . The Euclidean 2-point functions in the temporal direction and the associated SPFs are related by the Laplace transform. On the lattice, only a finite number of data points with statistical noise are available in the temporal direction. Therefore, direct Laplace inversion from lattice data to SPF is an ill-posed problem. MEM is a method to evade such difficulty on the basis of the Bayes’ theorem in the theory of statistical inference .

## 2 HOW MANY DATA POINTS ARE NECESSARY FOR MEM ?

As we have discussed in detail in , the result of MEM depends strongly on the number of temporal data points. The more data points are used, the closer the MEM result is to the true SPF. Moreover, there is the minimum number of data points $`N`$ to perform the reliable MEM analysis.

In order to find this minimum $`N`$ at $`T=0`$ on the lattice, we have carried out the following analysis . First we calculated hadronic correlators with the quenched approximation on an isotropic $`40^3\times 30`$ lattice at $`\beta =6.47`$ . The number of gauge configurations is 160. Then, MEM analysis has been done for the vector channel by using $`N`$ data points (at $`\tau =2,\mathrm{},N/2+1`$ and $`31N/2,\mathrm{},30`$) out of the total 30 data points. The result is shown in Fig. 1. The figure clearly shows that, SPF changes considerably as $`N`$ and that at least about 30 data points are necessary for the convergence of the result. This minimum $`N`$ would depend on $`\beta `$, and on whether one employs anisotropic lattices, improved actions, and so on. Nevertheless, we use this number as a practical guide in the following.

## 3 ANISOTROPIC LATTICE AT $`T0`$

At finite temperature ($`T`$), because of the relation $`T=1/N_\tau a_\tau `$ ($`N_\tau `$ and $`a_\tau `$ being the temporal lattice size and spacing, respectively), less and less data points are available as $`T`$ increases if $`a_\tau `$ is fixed. Requiring that $`N_\tau >30`$ holds even at the highest $`T`$, e.g. $`2.5T_c`$, we are inevitably led to use an anisotropic lattice.

We have used the bare anisotropy $`\xi _0=3.5`$ and $`\beta =7.0`$ with the naive plaquette action. For the quark part, the standard Wilson action with the quenched approximation is used. The corresponding renormalized anisotropy is $`\xi =4.0`$  and the lattice spacing is $`a_\tau =1/4a_\sigma =9.75\times 10^3`$ fm. Simulations are done on $`32^3\times N_\tau `$ lattices with $`N_\tau =32,`$ 40, 54, 72, 80, and 96. This corresponds to $`T2.3T_c`$, $`1.9T_c`$, $`1.4T_c`$, $`1.04T_c`$, $`0.93T_c`$, and $`0.78T_c`$, respectively. More than 100 configurations are generated for each $`N_\tau `$. The details of the lattice parameters used in the calculation will be given in .

We have studied 2-point correlation functions in the scalar (S), pseudo-scalar (PS), vector (V), and axial-vector (AV) channels. The default models, motivated by perturbative QCD, in each channel are $`m=0.60`$, 1.15, 0.40, and 0.35, respectively.

The lattice data points used in our MEM analysis are $`\tau =5,\mathrm{},21`$ and $`N_\tau 19,\mathrm{},N_\tau 3`$ for $`N_\tau =54`$, 72, 80, and 96. For $`N_\tau =32`$ and 40, $`\tau =5,\mathrm{},N_\tau 3`$ are used. When the source and sink are closer than $`\xi a_\tau `$ in the temporal direction, artifact could arise due to unphysical modes at $`\omega \pi /\xi a_\tau `$. Therefore, we leave out the data points near the edges, i.e., $`\tau =1,\mathrm{},4`$ and $`N_\tau 2,\mathrm{},N_\tau `$.

The number of lattice data points used in the MEM analysis is fixed almost the same at each $`N_\tau `$ except at $`N_\tau =32`$ to keep the resolution in the MEM analysis unchanged. The rest of the MEM procedures such as the averaging over $`\alpha `$ are basically the same as described in  and will be presented in detail in .

Fig. 2 shows the Polyakov loop expectation value $`|L|`$ at each $`N_\tau `$. It starts to deviate from zero around $`N_\tau =80(72)`$, which corresponds to $`T_c=253(281)`$ MeV, and is comparable to the real value 271 $`\pm `$ 2 MeV. We have checked that the Polyakov loop susceptibility also shows a peak in this region. In the following, we use the number $`T_c=271`$ MeV.

## 4 RESULTS OF MEM

In Fig. 3, we show the results of the MEM analysis for SPF $`\rho (\omega )`$ in the S, PS, V, and AV channels on $`32^3\times 54`$ lattice ($`T1.4T_c`$). In this calculation, we have used the hopping parameter corresponding to $`m_\pi /m_\rho 0.7`$ at $`T=0`$.

If the quark-gluon plasma is such an object as described as a free gas of massive quarks, anti-quarks and gluons, $`\rho (\omega )`$ would show a flat behavior with a smooth rise from zero above the $`q\overline{q}`$ threshold as a function of $`\omega `$. To the contrary, Fig. 3 shows that SPFs possess nontrivial structures even above $`T_c`$ in all four channels. There are a sharp peak at about $`\omega =2`$ GeV and two bumps above the peak. Also SPFs in all channels almost degenerate, which shows a strong signature of the restoration of chiral symmetry. The peak and bump structures, and the suppression of SPFs below $`\omega =2`$ GeV are statistically significant according to our error analysis. We have used $`N_\tau `$ twice as large as that employed in . As shown in Sec.2, this is necessary for reliable MEM analysis.

In Fig. 4, we show the results at $`N_\tau =40`$, i.e., $`T1.9T_c`$. There is an apparent peak around $`\omega =0`$. A possible explanation of the peak is the effect of Landau damping . However, this peak is, at the moment, not statistically significant yet. The peaks around $`\omega 4.5`$ GeV are significantly broader than those around 2 GeV at $`T1.9T_c`$.

We have carried out an error analysis for the average of the SPFs in small regions around the peaks and found that the broadening is indeed statistically significant. SPFs in the case of $`N_\tau =32`$ $`(T2.3T_c)`$ have similar structures as those at $`T1.9T_c`$ except the shift of the peak positions. This suggests a possibility of a qualitative change in the state of the deconfined phase between $`1.5T_c`$ and $`2T_c`$.

Acknowledgement This work was supported in part by the Grants-in-Aid by Ministry of Education (No. 12640263, No. 12640296, and No. 14540255). Lattice calculations have been carried out with the CP-PACS computer under the “Large-scale Numerical Simulation Program” of Center for Computational Physics, University of Tsukuba.


text:
# The deconfining phase transition in SU(𝑁_𝑐) gauge theoriesTalk by U. Wenger.

## 1 INTRODUCTION

The large-$`N_c`$ physics of QCD and SU($`N_c`$) gauge theories is of great theoretical and phenomenological interest. Lattice calculations of string tensions, glueball masses etc.  have confirmed that a smooth large-$`N_c`$ limit is achieved by keeping fixed the ’t Hooft coupling $`g^2N_c`$ and that, as expected, the leading corrections are $`O(1/N_c^2)`$.

An interesting question concerns the order of the deconfining phase transition at $`N_c\mathrm{}`$. It has been argued  that the transition in SU(3) may be accidently first order, due to a cubic term in the effective potential, but in general is second order in SU($`N_c`$). Our ongoing lattice study , on which we report here, is addressing this question and thereby improving on earlier results .

We use the standard plaquette action on $`L^3\times L_t`$ lattices with periodic boundary conditions. The deconfining phase transition is studied by keeping $`L_t`$ fixed and varying $`\beta =2N_c/g^2`$ so as to pass the temperature $`T=1/a(\beta )L_t`$ through the transition.

In order to distinguish the confining and the deconfining phases we define the spatial average $`\overline{l}_p`$ of the Polyakov loop and correspondingly the spatial average $`\overline{u}_p`$ of the plaquette variable $`\text{Tr}U_p`$. The two order parameters so defined can be used to characterise the order of the phase transition. For a first order transition, e.g., we should find tunneling between the phases near the critical temperature $`T_c`$ and subsequently a double peak structure in the probability distributions of $`\overline{l}_p`$ and $`\overline{u}_p`$.

In order to determine the critical couplings $`\beta _c(V)`$ at finite volume $`V=L^3`$ we define a normalised Polyakov loop susceptibility

$$\frac{\chi _l}{V}=|\overline{l}_p|^2|\overline{l}_p|^2$$

and, for each volume $`V`$, obtain $`\chi _l`$ as a continuous function of $`\beta `$ using standard reweighting techniques . Then we determine $`\beta _c`$ from the location at which the susceptibility has its maximum.

In complete analogy we can define the specific heat $`C(\beta )`$ from the average plaquette $`\overline{u}_p`$,

$$\frac{1}{\beta ^2}C(\beta )=\frac{}{\beta }\overline{u}_p=N_p\overline{u}_p^2N_p\overline{u}_p^2,$$

where $`N_p=6L^3L_t`$ is the total number of plaquettes. The value of $`\beta `$ where $`C`$ has its maximum provides a different definition for a critical coupling which however should agree with the one from $`\chi _l`$ in the thermodynamic limit. Moreover we can use the finite size scaling behaviour of the peak value $`C(\beta _c,V)`$ as a criterion for the determination of the order of the transition. To be more precise, we expect the leading finite size scaling

$$\frac{C(\beta _c,V)}{V}=\{\begin{array}{cc}a_1L^{\rho d}+a_2L^d,\hfill & \text{2nd order},\hfill \\ a_1+a_2L^d,\hfill & \text{1st order},\hfill \end{array}$$

where $`\rho <3`$ is some combination of critical exponents and $`d=3`$.

An important point we need to address is the influence of a bulk phase transition which can be close to the physical deconfining transition. Our choices of $`L_t`$ for the different gauge groups ensure that the two transitions are well separated and can unambiguously be distinguished.

## 2 $`T_c`$ IN SU(4)

In SU(4) we perform a finite size scaling study at $`L_t=5`$ on $`L=12,14,16,18,20`$ lattices with between 200 and 300 ksweeps at each $`\beta `$-value and lattice size. For $`\beta `$-values close to $`\beta _c`$ we find clear tunneling transitions both in $`\overline{l}_p`$ and $`\overline{u}_p`$. The probability distributions of these quantities exhibit clear double peaking which is characteristic for a first order transition.

In fig. 1 we plot our values of $`\beta _c(V)`$ against the inverse spatial volume expressed in units of the temperature $`T`$ together with the extrapolations to infinite volume.

We note that if the transition is first order, as ours clearly is, then the leading finite-$`V`$ correction should be

$$\beta _c(V)=\beta _c(\mathrm{})\frac{h}{VT^3},$$
(1)
where $`h`$ is independent of the lattice action. We find that all our values of $`\beta _c`$ are consistent with eq. (1) and we obtain $`\beta _c(\mathrm{})=10.63709(72)`$ in the $`V\mathrm{}`$ limit. Moreover the value of $`h=0.09\pm 0.02`$ we thus extract is close to the SU(3) value, $`h0.1`$ . This suggests that $`h`$ depends only weakly on $`N_c`$ and we will therefore use this value to extrapolate to $`V=\mathrm{}`$ in the SU(6) case where we do not perform an explicit finite volume study.

From $`\beta _c(\mathrm{})`$ we obtain the critical temperature $`T_c`$ in terms of the string tension $`\sigma `$ ,

$$\frac{T_c}{\sqrt{\sigma }}=0.6024\pm 0.0045\text{at}a=1/5T_c.$$
(2)
Our preliminary $`L_t=6`$ calculations give a value $`\beta _c(\mathrm{})=10.780(10)`$ which translates to $`T_c/\sqrt{\sigma }=0.597(8)`$. This, taken together with the value in eq. (2), yields the extrapolated continuum value

$$\underset{a0}{lim}\frac{T_c}{\sqrt{\sigma }}=0.584\pm 0.030\text{in SU(4)}.$$
(3)
In fig. 2 we plot $`C(\beta _c)/V`$ against $`1/V`$ and note that the scaling is certainly consistent with a first order one.

The intercept at $`1/V=0`$ provides a measure of the latent heat, $`\mathrm{\Delta }e`$, since

$$\underset{V\mathrm{}}{lim}\frac{C(\beta _c,V)}{\beta _c^2V}(\overline{u}_{p,c}\overline{u}_{p,d})^2=\mathrm{\Delta }e^2,$$

where $`\overline{u}_{p,c}`$ and $`\overline{u}_{p,d}`$ are the average plaquette values at $`\beta _c`$ in the confined and deconfined phases, respectively. We obtain $`\mathrm{\Delta }e(N_c=4)=0.00197(5)`$ and observe that if we take the $`L_t=4`$ SU(3) latent heat in  and naively scale it to $`L_t=5`$, we obtain $`\mathrm{\Delta }e(N_c=3)0.0013`$ which is substantially smaller than our above SU(4) value. This shows explicitly that the SU(4) transition is more strongly first order than the SU(3) one.

## 3 $`T_c`$ IN SU(6)

In SU(6) we perform a study on a $`16^3\times 6`$ lattice at values of $`\beta `$ close to $`\beta =24.845`$. In physical units this is like a $`14^3`$ lattice with $`L_t=5`$ and thus we are confident from our SU(4) calculations to be close to the thermodynamic limit.

As before we observe well defined tunneling between confined and deconfined phases, characteristic of a first order phase transition. Using the reweighting technique we extract a critical value $`\beta _c=24.850(3)`$ which we extrapolate to $`V=\mathrm{}`$ using eq. (1) with the corresponding value of $`h=0.09`$ from SU(4). We obtain $`\beta _c(\mathrm{})=24.855\pm 0.003`$ which translates into

$$\frac{T_c}{\sqrt{\sigma }}=0.588\pm 0.002\text{at}a=1/6T_c$$
(4)
using the values of the SU(6) string tension in . Preliminary $`L_t=5`$ calculations provide an estimate $`\beta _c(\mathrm{})=24.519(33)`$ which translates to $`T_c/\sqrt{\sigma }=0.580(12)`$. Taken together with the value in eq. (4) this gives the continuum value

$$\underset{a0}{lim}\frac{T_c}{\sqrt{\sigma }}=0.605\pm 0.026\text{in SU(6)}.$$
(5)
As for the strength of the first order transition, we note that an estimate of the latent heat, $`\mathrm{\Delta }e`$, in SU(6) is larger than the one on a directly comparable $`16^36`$ SU(4) calculation indicating that the transition is certainly not weakening as $`N_c\mathrm{}`$.

## 4 CONCLUSIONS

We find that SU(4) gauge theory at $`L_t=5`$ reveals a first order deconfining phase transition with a latent heat that is not particularly small. We also find evidence that the SU(6) transition at $`L_t=6`$ is first order as well and certainly not weaker than the one in SU(4). We come to that conclusion by comparing the appropriately rescaled latent heat expressed in terms of $`T_c`$ for the different gauge groups. In doing so we take care to avoid confusion with a nearby bulk phase transition which, on our lattices, can unambiguously be identified.

The $`N_c`$-dependence of the critical temperature when expressed in units of the string tension, $`T_c/\sqrt{\sigma }`$, appears to be weak. From a fit we obtain

$$\frac{T_c}{\sqrt{\sigma }}=0.582(15)+\frac{0.43(13)}{N_c^2}.$$
(6)
Finally, we conclude that the SU($`N_c`$) transition at $`N_c=\mathrm{}`$ is first order and not particularly weak.