# BI-TP 99/31 September 1999 Lattice QCD at Finite Temperature and Density

## 1 Introduction

### 1.1 Preface

There are many reasons to study a complicated interacting field theory like QCD under extreme conditions, eg. at finite temperature and/or non-vanishing baryon number density. Of course, first of all we learn about the collective behaviour of strongly interacting matter, its critical behaviour, equation of state and thermal properties of hadrons as well as their fate at high temperature. However, in doing so we also have a look at the complicated non-perturbative structure of the QCD vacuum from a different perspective and can learn about the mechanisms behind confinement and chiral symmetry breaking. This does allow to check concepts that have been developed to explain these non-perturbative features of QCD.

This review focuses on the former aspect of finite temperature QCD studies; the latter has recently been discussed in . In fact, an important step in the direction of understanding confinement in terms of the dual superconductor picture has been taken recently by the Pisa group . They analysed the scaling behaviour of an order parameter for monopole condensation and could show that it scales in SU(2) and SU(3) gauge theories with the critical exponents expected for the deconfinement transition (Fig. 1). This is a strong hint for the survival of the monopole condensation mechanism as a source for confinement in the continuum limit.

### 1.2 Outline

Significant progress in finite temperature lattice QCD is based on the use of improved actions. In the pure gauge sector we have seen that these actions strongly reduce finite cut-off effects and allow to calculate thermal quantities on rather coarse lattices in good agreement with continuum extrapolated results. A similar approach is now followed in calculations with light quarks. The use of improved staggered and Wilson actions leads to new estimates of the critical temperature and the equation of state which we will discuss in the following sections. Moreover, thermodynamic calculations with domain wall fermions (DWF), which can greatly improve the chiral properties of lattice fermion actions, reached a stage where first quantitative results can be reported on $`T_c`$ and thermal masses. The latter will be discussed in Section 4. In Section 5 we will discuss new developments in the analysis of QCD at finite density. Section 6 contains an outlook.

An important part of finite temperature ($`T`$) and density ($`n_B`$) (chemical potential ($`\mu `$)) calculations consists of trying to understand the QCD phase diagram, i.e. the order of the phase transition as function of the number of flavours ($`n_f`$) and quark masses ($`m_q`$) and its dependence on $`T`$ and $`n_B`$ (or $`\mu `$). We will not discuss these issues here. The finite temperature phase diagram has been discussed extensively in recent reviews  and the possibly rather complex phase structure of QCD at high density and low temperatures has been the topic of E. Shuryak contribution to this conference .

## 2 The Critical Temperature

One of the basic goals of lattice QCD calculations at finite T is to provide quantitative results for the phase transition temperature. In the pure gauge sector this goal has been achieved. The value of the critical temperature for the deconfining phase transition in a $`SU(3)`$ gauge theory is known with small statistical errors. Remaining systematic errors are of the order of 3%. They arise from the calculation of the string tension at $`T=0`$ which is used to set the scale for $`T_c`$. In Table 1 we summarise results for $`T_c/\sqrt{\sigma }`$ obtained with different gauge actions on lattices with varying temporal extent $`N_\tau 1/(aT)`$ and extrapolated to the continuum limit ($`N_\tau \mathrm{}`$).

Using for the string tension the value<sup>1</sup><sup>1</sup>1This may be deduced from quenched spectrum calculations ($`m_\rho /\sqrt{\sigma }=1.81(4)`$ ). Recent estimates from CP-PACS data seem to lead to a somewhat smaller value . $`\sqrt{\sigma }425\mathrm{MeV}`$, we find $`T_c270\mathrm{MeV}`$. This large value can be understood in terms of the particle spectrum of the quenched theory; in the low temperature, confining phase there exist only rather heavy glueballs ($`m_G>1.7`$ GeV). A rather large temperature thus is needed to create a sufficiently dense glueball gas, which can trigger a deconfining transition.

Investigations of the critical temperature in QCD with quarks of finite mass indeed have shown that the transition temperature drops rapidly with decreasing quark masses<sup>2</sup><sup>2</sup>2At intermediate values of $`m_q`$ the transition only reflects a rapid cross-over in thermodynamic observables rather than a true phase transition. The peak in the chiral susceptibility or the Polyakov-Loop susceptibility defines in this case a pseudo-critical temperature. We will continue to call this temperature the transition temperature at finite values of $`m_q`$ For 2-flavour QCD it is expected to be a critical temperature, corresponding to a phase transition, only in the chiral limit. For 3-flavour QCD it is a critical temperature below a certain critical quark mass ..

Unlike in the pure gauge theory the transition temperature for QCD with finite quark masses does seem to be strongly affected by the discretization scheme used for the fermion action. The early calculations with the standard staggered  and Wilson  actions led to widely different estimates for $`T_c`$. Finite-cut off effects are thus expected to be large and improvement of the fermion action should be expected to be important. Indeed, the new calculations performed with improved Wilson (Clover)  and staggered  actions as well as with DWF  yield transition temperatures which are in much better agreement among each other. Unfortunately, this statement is, at present, only partially correct. In fact, for comparable values of $`m_q`$, i.e. fixed ratios of pseudo-scalar (“pion”) and vector meson (”rho”) masses, $`m_{PS}/m_V`$, it only holds when we set the scale for $`T_c`$ using a hadron mass, eg. $`m_V`$. When we follow a similar approach and use the string tension, $`\sqrt{\sigma }`$, to set the scale the agreement is less evident. In this case calculations with the Wilson fermion action typically lead to transition temperatures about 20% below the results obtained with staggered fermion actions. The current status of the determination of $`T_c`$ for 2-flavour QCD is summarised in Fig. 2.

As can be seen the simulations with Clover fermions performed by different groups are consistent with each other. Results on $`T_c/m_V`$ do not seem to depend significantly on the gauge action (one plaquette Wilson , Symanzik improved  or RG-improved ) nor does it seem to be important whether tadpole  or non-perturbative  Clover coefficients are used. Also the DWF results  seem to be insensitive to the gauge action chosen (plaquette or RG-improved) and are consistent with results obtained with the Clover action. Results obtained with an improved staggered fermion action , the p4-action<sup>3</sup><sup>3</sup>3The p4-action improves the rotational symmetry of the quark propagator in $`𝒪(p^4)`$. It strongly reduces the cut-off effects in bulk thermodynamic observables , see Fig. 4., also agree with these data within 10%.

Unfortunately this consistent picture is, at present, not reproduced when calculating $`T_c`$ in units of $`\sqrt{\sigma }`$ (right frame of Fig. 2). A possible source for this discrepancy might be the calculation of the heavy quark potential, which in the case of Wilson fermions so far has only been performed on rather small spatial lattices, e.g. $`8^3\times 16`$. This may lead to an overestimate of the string tension. It, however, also is possible that calculations on $`N_\tau =4`$ lattices are still performed at too strong coupling and do not allow for a unique determination of the scale. A hint in this direction may be the large difference in $`T_c/\sqrt{\sigma }`$ observed in calculations on $`N_\tau =4`$ and 6 lattices with the non-perturbatively improved Clover action . Clearly more work is needed here to establish a unique result for $`T_c`$ using different fermion formulations.

In general, we note that the transition temperature obtained with improved actions tends to be larger than what previously has been quoted on the basis of calculations performed with the standard staggered action. If the quark mass dependence does not change drastically closer to the chiral limit<sup>4</sup><sup>4</sup>4Calculations within the framework of quark-meson models suggest a rapid drop of $`T_c`$ for $`m_{PS}<m_\pi `$ . the current data suggest

$$T_c(170190)\mathrm{MeV}$$
(1)   
for 2-flavour QCD in the chiral limit. In fact, this estimate also holds for 3-flavour QCD. Calculations which are currently performed for $`n_f=2`$ and 3 using the same improved staggered fermion action, suggest that the flavour dependence of $`T_c`$ is rather weak  (see Fig. 3).

It is remarkable that the transition temperature drops significantly already in a region where all hadron masses are quite large. This is apparent from Fig. 3 where we show $`T_c`$ in units of $`\sqrt{\sigma }`$ plotted vs. $`m_{PS}/\sqrt{\sigma }`$ for $`n_f=2`$ and 3. As can be seen the transition temperature starts deviating from the quenched values for $`m_{PS}<(67)\sqrt{\sigma }2.5\mathrm{GeV}`$. We also note that the dependence of $`T_c`$ on $`m_{PS}/\sqrt{\sigma }`$ is almost linear in the entire mass interval. This might be expected for light quarks in the vicinity of a $`2^{nd}`$ order chiral transition where the pseudo critical temperature depends on the mass of the Goldstone-particle like

$$T_c(m_\pi )T_c(0)m_\pi ^{1/\beta \delta }.$$
(2)   
For 2-flavour QCD where the critical indices are expected to belong to the universality class of 3-d, $`O(4)`$ symmetric spin models one would indeed expect $`1/\beta \delta =1.1`$. However, this clearly cannot be the origin of the quasi linear behaviour observed for large hadron masses independent<sup>5</sup><sup>5</sup>5A similar conclusion holds for $`n_f=2`$, when one analyses the unimproved standard staggered fermion data. of $`n_f`$. A resonance gas model would probably be more appropriate to describe the thermodynamics for these heavy quarks.

## 3 The Equation of State

In the pure gauge sector bulk thermodynamic quantities such as the pressure and energy density have been analysed in detail. It has been verified that the most significant cut-off effects result from high momentum modes, which dominate the infinite temperature, ideal gas limit. Improved actions that lead to small cut-off effects in the ideal gas limit<sup>6</sup><sup>6</sup>6In the ideal gas limit the cut-off dependence can be analysed analytically . For a recent analysis of staggered fermion actions see . still do so at finite temperature .

In a recent analysis the CP-PACS collaboration calculated the pressure and energy density using the RG-improved (Iwasaki) action . They confirm that after extrapolation to the continuum limit also this action yields results for the $`SU(3)`$ equation of state which are, within an error of (3-4)%, consistent with the continuum extrapolation obtained from the Wilson , improved  as well as some fixed point actions . In fact, the Wilson and RG-improved actions represent extreme cases for such calculations. Their cut-off corrections have opposite sign which on coarse ($`N_\tau =4`$) lattices leads to deviations of the pressure by about 25% from the continuum extrapolated result. In view of this the agreement reached with different discretization schemes is rather reassuring.

Also in the presence of light quarks the use of an improved action thus seems to be mandatory, if one wants to calculate bulk thermodynamic quantities. The standard staggered and Wilson actions are known to lead to large deviations from the continuum ideal gas behaviour on coarse lattices (small $`N_\tau `$) . For staggered fermions improved actions, leading to smaller cut-off effects in thermodynamic observables, can be constructed by adding suitably chosen three-link terms to the conventional one-link terms . In Fig. 4 (left frame) we show the cut-off dependence of the ideal gas pressure obtained from these minimally improved staggered fermion actions (Naik  and p4  action)<sup>7</sup><sup>7</sup>7Staggered actions with even smaller cut-off effects have been constructed . However, they also require a significantly larger numerical effort..

As expected the systematics seen for the cut-off dependence in the ideal gas limit carries over to finite T. In the right frame of Fig. 4 we show results for the pressure of 2-flavour QCD obtained with the standard staggered action (and Wilson gauge action) on lattices with temporal extent $`N_\tau =4`$ and 6 . These are compared with results obtained with the p4-action (and a Symanzik improved gauge action) . In the latter case also fat links  have been introduced in the one-link terms to improve the flavour symmetry of the staggered action.

While the pressure calculated with the standard staggered action rapidly overshoots the ideal gas limit (as expected from the analysis of the cut-off effects in the ideal gas limit) the results obtained with the p4-action stay below the ideal gas limit and show a temperature dependence very similar to what has been found in the pure gauge sector. We note that the calculations with the p4-action have been performed with rather large quark masses, $`m/T=0.4`$, corresponding to $`m_{PS}/m_V0.7`$, while the calculation with standard staggered action are for $`m_{PS}/m_V0.3`$. In the high temperature phase this does not seem to constitute a major problem<sup>8</sup><sup>8</sup>8Calculations for 4-flavour QCD with quark masses $`m/T=0.2`$ and $`0.4`$ show no significant quark mass dependence . Moreover, also in the continuum the pressure of an ideal gas of fermions with mass $`m/T=0.4`$ deviates by less than 10% from that of a massless gas.. Smaller quark masses are, however, definitely needed close to $`T_c`$ and below in order to become sensitive to the contributions from light pion modes.

In the case of the SU(3) gauge theory the magnitude of the cut-off dependence in the temperature range up to a few times $`T_c`$ has been found to be about half of what has been calculated for the ideal gas, infinite temperature limit. If this continues to hold for the fermion sector one should expect, that the current result, obtained with the p4-action on $`N_\tau =4`$ lattices, underestimates the continuum result by about 15% for $`T>T_c`$. Based on this consideration an estimate for the continuum extrapolated pressure is also shown in Fig. 4. We also note that a calculation performed with identical, improved staggered fermion actions shows that the pressure normalised to the continuum ideal gas value has the same temperature dependence in 2 and 3 flavour QCD .

Unfortunately, Wilson actions with similarly good high temperature behaviour have not been constructed so far. The Clover action does not improve the ideal gas behaviour, i.e. it has the same infinite temperature limit as the Wilson action. Nonetheless, a first attempt to calculate the equation of state, using the tadpole-improved Clover action combined with the RG-improved gauge action, has now been undertaken by the CP-PACS collaboration . A preliminary result for the pressure in 2-flavour QCD is shown in Fig. 5. Similar to simulations with the standard staggered fermion action one observes an overshooting of the ideal gas limit reflecting the cut-off effects in the unimproved fermion sector.

## 4 Thermal Masses

Understanding the temperature dependence of hadron properties, e.g. their masses and widths, is of central importance for the interpretation of heavy ion experiments. Thermal modifications of the heavy quark potential influence the spectrum of heavy quark bound states. Their experimentally observed suppression  thus is expected to be closely linked to the deconfining properties of QCD above $`T_c`$ . Changes in the chiral condensate, on the other hand, influence the light hadron spectrum and may leave experimental signatures, for instance in the enhanced dilepton production observed in heavy ion experiments .

In numerical calculations on Euclidean lattices one has access to thermal Green’s functions $`G_H(\tau ,\stackrel{}{r})`$ in fixed quantum number channels, $`H`$, to which in particular at high temperature many excited states contribute. As the temporal direction of the Euclidean lattice is rather short at finite temperature one usually has not enough information on the correlation function to extract reliably thermal effects on the (pole) masses. A way out may come from the use of anisotropic lattices, which has extensively been explored by the QCDTARO group . In their recent analysis with quenched Wilson fermions they show evidence for a rapid change in the thermal correlation functions across $`T_c`$ which indicates a rapid change in thermal masses above $`T_c`$. At the same time their analysis of pseudo-scalar wavefunctions does, however, show evidence that the correlation among quarks in this quantum number channel remains strong. In how far this hints at the presence of light ”pionic” bound states has to be analysed further.

Anisotropic lattices in combination with NRQCD have also been used to analyse thermal effects on heavy quark bound states . Large mass shifts have been observed for the first excited states.

The missing information on the long distance behaviour of thermal correlation functions may also be overcome by using refined techniques to analyse the numerical data on thermal correlation functions. It recently has been suggested that the maximum entropy method  may help to extract more detailed information on thermal modifications of hadronic spectral functions .

So far the most convincing evidence for modifications of hadron properties at high temperature comes from the analysis of the behaviour of correlation function at large spatial separations, which yields screening masses , as well as thermal susceptibilities,

$$\chi _H=_0^{1/T}d\tau \mathrm{d}^3rG_H(\tau ,\stackrel{}{r}),$$
(3)   
with $`G_H`$ denoting the hadron correlation function in the quantum number channel $`H`$. If there is only a single stable particle of mass $`m_H`$ contributing to $`G_H`$ then $`\chi _Hm_H^2`$. In general the susceptibilities, however, define only effective masses, i.e. they average over contributions from the ground state and all excited states in a particular quantum number channel.

In Fig. 6 we show $`1/\sqrt{\chi _H}`$ for 2-flavour QCD obtained with staggered fermions of mass $`m_q=0.02`$ on lattices of size $`8^3\times 4`$.. This figure is based on data from . As can be seen the $`f_0`$ and $`\pi `$ masses become (almost) degenerate at $`\beta _c`$ while the $`a_0`$ remains heavy at $`\beta _c`$. The former behaviour is expected, $`f_0`$ and $`\pi `$ correlation functions are related through a rotation in $`SU(2)`$ flavour space. The degeneracy thus reflects the restoration of the $`SU(2)`$ flavour symmetry. The difference of $`\chi _{a_0}`$ and $`\chi _\pi `$ on the other hand reflects the persistence of the $`U_A(1)`$ symmetry breaking. A crucial question here is how the gap observed for non-zero $`m_q`$ changes in the chiral limit. At high temperature topologically non-trivial gauge field configurations are expected to be suppressed, which in turn would lead to a strong reduction in the strength of $`U_A(1)`$ symmetry breaking and thus in a strong reduction of the mass splitting between $`a_0`$ and $`\pi `$. For 2-flavour QCD it is expected that the quark mass dependence is quadratic, $`(m_{a_0}m_\pi )(\chi _\pi \chi _{a_0})A+Bm_q^2`$. Previous investigations of this were, however, not conclusive . If a quadratic ansatz is assumed calculations with staggered fermions led to the conclusion that a non-zero mass splitting remains also above $`T_c`$, i.e. the $`U_A(1)`$ remains broken . The problem has now been addressed again by the Columbia group  using DWF. Due to the improved chiral properties of this action one should find a quadratic dependence on the quark mass. For $`T1.2T_c`$ the Columbia group indeed does observe such a quark mass dependence. In the chiral limit they find a non-zero mass splitting from the susceptibilities as well as from the analysis of screening masses ,

$$m_{a_0}m_\pi =0.0606(67)+9.66(58)m_q^2.$$
(4)   
The $`U_A(1)`$ thus remains broken above $`T_c`$, although the mass splitting is strongly reduced (see Fig. 6). This picture is also supported by an analysis of the disconnected parts of flavour singlet correlation functions in quenched QCD using DWF . A vanishing of these would signal a degeneracy between the $`\sigma `$ and $`\eta `$ mesons. The disconnected parts get non-zero contributions only from topologically non-trivial configurations . Indeed, these are strongly suppressed above $`T_c`$. However, at $`T=1.25T_c`$ about 10% of the configurations do still carry a non-zero topological charge indicating the persistence of $`U_A(1)`$ symmetry breaking above $`T_c`$.

## 5 Finite Density QCD

Finite density calculations in QCD are affected by the well known sign problem, i.e. the fermion determinant becomes complex for non-zero values of the chemical potential $`\mu `$ and thus prohibits the use of conventional numerical algorithms. The most detailed studies have been performed so far using the Glasgow algorithm , which is based on a fugacity expansion of the grand canonical partition function at non-zero $`\mu `$,

$$Z_{GC}(\mu /T,T,V)=\underset{B=\alpha V}{\overset{\alpha V}{}}z^BZ_B(T,V),$$
(5)   
where $`z=\mathrm{exp}(\mu /T)`$ is the fugacity and $`Z_B`$ are the canonical partition functions for fixed quark number $`B`$; $`\alpha =3,6`$ for one species of staggered or Wilson fermions, respectively. However, this approach so far did not overcome the severe numerical difficulties.

It thus seems necessary to approach the finite density problems from another point of view. A reformulation of the original ansatz may lead to a representation of the partition function which, in the ideal case, would require the averaging over configurations with strictly positive weights only, or at least would lead to a strong reduction of configurations with negative weights.

A big step away from the original formulation is to start from a Hamiltonian approach . Here it has been shown that problems with a fluctuating integrand can successfully be reformulated in terms of a model where the configurations generated do have strictly positive weights (meron cluster algorithm ). Whether such an approach can be applied to QCD remains to be seen.

An alternative formulation of finite density QCD is given in terms of canonical rather than grand canonical partition functions , i.e. rather than introducing a non-zero chemical potential through which the number density is controlled one introduces directly a non-zero baryon number (or quark number $`B`$) from which the baryon number density on lattices of size $`N_\sigma ^3\times N_\tau `$ is obtained as $`n_B/T^3=\frac{B}{3}(N_\tau /N_\sigma )^3`$.

After introducing a complex chemical potential in $`Z_{GC}`$(Eq. 5) the canonical partition functions can be obtained via a Fourier transform<sup>9</sup><sup>9</sup>9The use of this ansatz for the calculation of canonical partition functions as expansion coefficients for $`Z_{GC}`$ has been discussed in .,

$$Z_B(T,V)=\frac{1}{2\pi }_0^{2\pi }d\varphi \mathrm{e}^{i\varphi B/T}Z_{GC}(i\varphi ,T,V).$$
(6)   
Also this formulation is by no means easy to use in general, i.e. for QCD with light quarks. In particular, it also still suffers from a sign problem. It, however, leads to a quite natural and useful formulation of the quenched limit of QCD at non-zero density .

### 5.1 Quenched limit of finite density QCD

It had been noticed early that the straightforward replacement of the fermion determinant by a constant does not lead to a meaningful static limit of QCD . In fact, this simple replacement corresponds to the static limit of another theory with an equal number of fermion flavours carrying baryon number $`B`$ and $`B`$, respectively . This should not be too surprising. When one starts with and takes the limit of infinitely heavy quarks something should be left over from the determinant that represents the objects that carry the baryon number. In the canonical formulation this becomes obvious. For $`m_q\mathrm{}`$ one ends up with a partition function, which for baryon number $`B/3`$ still includes the sum over products of $`B`$ Polyakov loops, i.e. the static quark propagators which carry the baryon number . This limit also has some analogy in the grand canonical formulation where the coupled limit $`m_q,\mu \mathrm{}`$ with $`\mathrm{exp}(\mu )/2m_q`$ kept fixed has been performed <sup>10</sup><sup>10</sup>10This is a well known limit in statistical physics. When deriving the non-relativistic gas limit from a relativistic gas of particles with mass $`\overline{m}`$, the rest mass is splitted off from the chemical potential, $`\mu \mu _{nr}+\overline{m}`$, in order to cancel the corresponding rest mass term in the particle energies. On the lattice $`\overline{m}=\mathrm{ln}(2m_q)`$ for large bare quark masses..

In the confined phase of QCD the baryon number is carried by the rather heavy nucleons. Approximating them by static objects may thus be quite reasonable and we may expect to get valuable insight into the thermodynamics of QCD at non-zero baryon number density already from quenched QCD.

From a numerical point of view there is hope that we can get control over this limit using different approaches. At non-zero $`\mu `$ a variant of the Glasgow approach  seems to become applicable for large quark masses and also the static limit of Blum et al.  may be explored further .

In the canonical approach simulations at non-zero $`B`$ can be performed on relatively large lattices and the use of baryon number densities up to a few times nuclear matter density may be possible . The simulations performed so far show the basic features expected at non-zero density. As can be seen from the behaviour of the Polyakov loop expectation value shown in Fig. 7 the transition region gets shifted to smaller temperatures (smaller coupling $`\beta `$). The broadening of the transition region may suggest a smooth crossover behaviour at non-zero density. However, in a canonical simulation it also may indicate the presence of a region of coexisting phases and thus would signal the existence of a $`1^{st}`$ order phase transition. This deserves further analysis.

Even more interesting is the behaviour of the heavy quark potential in the low temperature phase. As shown in the right frame of Fig. 7 the potential does get screened at non-vanishing number density. This will have a direct influence on heavy quark bound states at high density.

## 6 Outlook

We have focused in this review on the calculation of basic thermodynamic quantities which are of immediate interest to experimental searches for the Quark Gluon Plasma.

In quenched QCD the critical temperature and the equation of state have been calculated on the lattice and extrapolated to the continuum with an accuracy of a few percent. These calculations set a benchmark for many analytical studies of QCD thermodynamics . The progress made in developing and testing improved fermion actions for thermodynamic calculations shows that a similar accuracy for QCD with light quarks is within reach. The current systematic studies with different improved fermion actions may soon lead to a determination of the transition temperature and the equation of state with similar accuracy.

We do have reached some understanding of thermal effects on hadron properties. In particular, modifications of the light meson spectrum due to flavour and approximate $`U_A(1)`$ symmetry restoration have been established. However, lattice calculations so far did not come up with detailed quantitative results on thermal masses, which could be confronted with experimental data. There are a few promising ansätze which can lead to more detailed information on thermal modifications of hadronic spectral densities.

Of course, there are many more important issues which have to be addressed in the future. Even at vanishing baryon number density we do not yet have a satisfactory understanding of the critical behaviour of 2-flavour QCD in the chiral limit and the physically realized situation of QCD with two light, nearly massless quarks and a heavier strange quark has barely been analysed.

Moreover, the entire phase diagram at non-zero baryon number density is largely unexplored. An interesting phase structure is predicted at high density and low temperatures which currently is not accessible to lattice calculations. This does require new algorithmic developments. 

There are thus many interesting questions waiting to be answered in the next millennium.

Acknowledgements: I would like to thank the CCP at the University of Tsukuba for its kind hospitality during the time this talk has been prepared and written up. I also want to thank N. Christ, S. Ejiri, M. Okamoto, D.K. Sinclair, I.O. Stamatescu, D. Toussaint, and U.-J. Wiese for communication on their results and K. Kanaya and A. Ukawa for comments on the manuscript.