format: text

text:
# 1 INTRODUCTION

## 1 INTRODUCTION

The field of relativistic heavy-ion collisions is at a unique time in its development. The Relativistic Heavy-Ion Collider (RHIC) at BNL started regular operation in the summer of year 2000. In a very short, but extremely successful, run with an integrated luminosity of only a few $`\mu b^1`$, the four RHIC experiments delivered an impressive amount of results offering a first glimpse at the physics of Au-Au collisions at the energy of $`\sqrt{s_{NN}}=130GeV`$, almost one order of magnitude larger than the highest energy of $`\sqrt{s_{NN}}=17.2GeV`$ available at CERN. Some of these results are new and exciting, some are puzzling and some follow the pattern already established at lower energies. A lot more is expected to come from the second run at the design energy of $`\sqrt{s_{NN}}=200GeV`$, completed by the time of this conference with almost two orders of magnitude larger integrated luminosity. In addition to that, the 15 year programme at CERN SPS and BNL AGS has produced a wealth of very interesting and intriguing results and more is still to be expected from the SPS in the next few years. This combination results in very exciting opportunities for the study of matter under extreme conditions and in particular for the quest of the phase transition associated with quark-gluon plasma formation and chiral symmetry restoration.

About a dozen of different experiments at the AGS and SPS and four experiments at RHIC are or have been involved in this endeavour covering a very broad range of observables. In the limited space of this paper it is not possible to do justice to the vast amount of available results . A selection is therefore unavoidable and this review is restricted to selected topics on recent results and highlights of the field on: (i) global observables where systematic comparisons are made from AGS up to RHIC energies (Section 2), (ii) hadron spectra and yields (Section 3), (iii) excess emission of low-mass lepton pairs (Section 4) and $`J/\psi `$ suppression (Section 5) which are among the most notable results hinting at new physics from the SPS programme (results on these two topics are not yet available at RHIC as they require a much higher luminosity than achieved in the first year), and (iv) suppression of large $`p_T`$ hadrons (Section 6) which is one of the highlights of the first RHIC run pointing to new phenomena opening up at RHIC energies.

## 2 GLOBAL OBSERVABLES

### 2.1 $`N_{ch}`$ and $`E_T`$ distributions

Global observables, like multiplicity and transverse energy, provide very valuable information. Besides defining the collision geometry, they can be related to the initial energy density, e.g. using the well-known Bjorken relation , shed light into the mechanisms of particle production and provide constraints to the many models aiming at describing these collisions .

The charged particle rapidity density at mid-rapidity has been measured by the four RHIC experiments -BRAHMS, PHENIX, PHOBOS and STAR- with very good agreement among their results . For central Au-Au collisions at $`\sqrt{s_{NN}}=130GeV`$, the global average is $`dN_{ch}/d\eta =580\pm 18`$. For the transverse energy rapidity density there is only one measurement, by PHENIX, with a value $`dE_T/d\eta =578_{39}^{+26}GeV`$ for the most central 2% of the inelastic cross section . Using the Bjorken formula  this translates into an initial energy density of $`ϵ=5.0GeV/fm^3`$ 222using the canonical formation time $`\tau =1fm/c`$. Much shorter formation times, as advocated in saturation models , would result in energy densities as high as $`ϵ20GeV/fm^3`$. i.e. 60-70% larger than the corresponding values at the full SPS energy.

The energy dependence of the charged particle density $`dN_{ch}/d\eta `$, normalized to the number of participating nucleon pairs ($`N_{part}/2`$), exhibits an almost logarithmic rise with $`\sqrt{s_{NN}}`$ from AGS up to RHIC energies, as shown in Fig. 1. However, the dependence is very different from the $`p\overline{p}`$ systematics also shown in Fig. 1. At $`\sqrt{s_{NN}}=200GeV`$, $``$65% more particles per pair of participants are produced in central Au-Au collisions than in $`p\overline{p}`$. Note also that the large increase predicted by the HIJING model with jet quenching  (upper curve in Fig. 1) is not observed, and from $`\sqrt{s_{NN}}=130`$ to $`200GeV`$ the particle density increases by only 15% . (See also Sections 2.2 and 6 for further discussion on jet quenching).

The centrality dependence of the charged particle density has been proposed as a sensitive tool to shed light on the particle production mechanism: soft processes are believed to scale with the number of participants $`N_{part}`$, whereas hard processes, expected to become more significant as the energy increases, lead to a scaling with the number of binary collisions $`N_{bin}`$. PHENIX  and PHOBOS  have reported an increase of $`dN_{ch}/d\eta `$ and $`dE_T/d\eta `$ strongly than linear with $`N_{part}`$ (the PHENIX results are shown in Fig. 2 left panel) and in the framework of models with these two components, like HIJING , such an increase is interpreted as evidence of the contribution of hard processes to particle production 333This is, however, not a unique interpretation. The centrality dependence results of PHENIX shown in Fig. 2 have also been explained with a model based on purely soft processes .. This is to be contrasted to SPS results where WA98 reports a much weaker increase (also shown in Fig. 2)  and WA97 an even weaker increase consistent within errors with proportionality with $`N_{part}`$ . Both the PHENIX and WA98 data, extrapolated to peripheral collisions, are in very good agreement with the $`p\overline{p}`$ result at the same $`\sqrt{s}`$ derived from the UA5 data . The increased role of hard processes at higher energies and their scaling with $`N_{bin}`$ could then explain the stronger increase of particle production in $`AA`$ collisions compared to $`p\overline{p}`$ previously discussed in the context of Fig. 1. The importance of hard processes at RHIC energies is a very interesting issue which will be further discussed in this review.

More surprising is the behavior of the ratio ($`dE_T/d\eta `$) / ($`dN_{ch}/d\eta `$), the average transverse energy per charged particle, shown in Fig. 2 (right panel). This ratio is found to be independent of centrality and approximately equal to $`0.8GeV`$. Within errors of the order of 10-20% this ratio is also independent of $`\sqrt{s_{NN}}`$ from AGS up to RHIC, implying a constant or a very moderate increase of the average $`p_T`$ per particle.

It is also interesting to note that UA1 quotes a very similar ratio at $`\sqrt{s}=200GeV`$  but with a rather large error of $`20\%`$. This seemingly universal behaviour of constant energy/particle is one of the most puzzling results. The increase in $`dE_T/d\eta `$ with $`\sqrt{s_{NN}}`$ translates into an increase in the number of produced particles rather than in the production of particles with larger $`E_T`$.

### 2.2 Elliptic flow

Elliptic flow originates from the spatial anisotropy of the diamond shaped system formed in non-central collisions, later converted into momentum anisotropy if enough rescatterings among particles occur in the further evolution of the system. Elliptic flow provides therefore an opportunity to learn about the degree of thermal equilibrium achieved with emphasis on the early collision dynamics.

Elliptic flow has been measured at RHIC by STAR , PHOBOS  and PHENIX  using different analysis methods and with good agreement between their results. The magnitude of elliptic flow, quantified e.g. by the second Fourier coefficient $`v_2`$ of the azimuthal particle distribution with respect to the reaction plane, is shown in Fig. 3 (left panel) as function of beam energy. The strength of $`v_2`$ at RHIC follows the systematic trend established over a very broad range of lower energies from SIS to SPS. The STAR value of 6% is approximately 50% larger than at the SPS, suggesting a much stronger early approach to local thermal equilibrium at RHIC.

Other interesting features appear in the $`p_T`$ dependence of $`v_2`$ shown in Fig. 3 (right panel) . Up to $`p_T`$ of $`1GeV/c`$, hydrodynamic calculations, which implicitly assume local thermal equilibrium, are in remarkably good agreement with the data, contrary to the observations at the SPS where the experimental results are significantly below the hydrodynamic limit . At higher $`p_T`$, when pQCD starts to become significant, $`v_2`$ appears to saturate. This behaviour can be reproduced by incorporating in the hydrodynamic model the energy loss of hard partons (jet quenching) in a dense medium with an initial gluon density of $`dN^g/dy5001000`$  444This might not be a unique interpretation. In a very recent paper, the $`v_2`$ saturation at high $`p_T`$ is explained by the two-particle correlation of produced minijets .. This novel feature of jet quenching is further discussed in Section 6.

### 2.3 Particle interferometry

Two-particle correlation with small relative momentum (HBT) is a useful tool to learn about the lifetime and size of the source at the time of particle emission. Since pions and kaons are emitted rather late, these correlations provide information about the space-time extent of the hadron system in the later stages of the collision.

The HBT technique has been widely used in ion experiments over the past decade at the AGS and SPS  and first results from RHIC on pion interferometry are already available . The systematic behaviour of the source parameters –the radii $`R_o,R_s,R_l`$ and the coherence parameter $`\lambda `$–, as function of energy is displayed in Fig.4. The RHIC results are a big puzzle. All HBT parameters have almost the same values as at the SPS. In particular, the larger source size $`R_s`$ and longer source lifetime reflecting itself in $`R_o/R_s2`$, expected if a QGP is formed , are not observed at RHIC. Instead, $`R_s`$ is very close to the SPS value and $`R_o/R_s1`$ again as at the SPS. These results are surprising and raise questions about the present understanding of the space-time evolution of relativistic heavy-ion collisions.

## 3 HADRONIC OBSERVABLES

Hadronic observables (spectra of identified particles, ratios of particle/antiparticle production, absolute particle yields) contribute, together with the global observables discussed in the previous section, to complete our view of the space-time evolution of RHI collisions. From their systematic studies, the picture has emerged of an initial high density state, which quickly reaches equilibrium and cools down while undergoing collective expansion. Particles decouple when the temperature and thus the particle density are low enough that interactions stop. We distinguish the chemical freeze-out when inelastic collisions cease to occur and particle abundances are frozen and the kinetic freeze-out, occuring later and at a lower temperature, when elastic processes stop and transverse momenta are frozen.

### 3.1 Chemical freeze-out

If we assume that an ideal hadron gas in thermal equilibrium is formed, then, in a statistical model using the grand canonical ensemble, the final particle production ratios are governed by only two independent variables, the temperature T and the baryochemical potential $`\mu _B`$ characterizing the system at chemical freeze-out . This simple model works extremely well and reproduces all data available from SIS, AGS, SPS and also the new data from RHIC. It is interesting to note that the enhanced production of strange and multi-strange hadrons, one of the first proposed signatures of QGP formation, is also described by the statistical model nearly as well as all other species.

The parameters T and $`\mu _B`$ derived from SPS and RHIC data are listed in Table 1 .

The temperatures are practically the same at RHIC and SPS, whereas $`\mu _B`$ is considerably lower at RHIC reflecting a lower net baryon density (see below). However, at both SPS and RHIC, the freeze-out parameters in the T vs $`\mu _B`$ plane are close to the expected boundaries of the QGP phase transition. So presumably, prior to freeze-out, the system was in a deconfined or at least in a mixed phase.

### 3.2 Kinetic freeze-out

The conditions at kinetic freeze-out are derived from the transverse momentum spectra of identified particles. As an example, Fig. 5 shows the positive and negative $`p,K,\pi `$ spectra measured by PHENIX in central Au-Au collisions at $`\sqrt{s_{NN}}=130GeV`$ .

The spectra exhibit an almost exponential shape with a temperature T (the inverse slope parameter) increasing with particle mass. This qualitative behaviour is independent of $`\sqrt{s_{NN}}`$. It has been observed all the way from BEVALAC up to RHIC energies and can be understood if one postulates that there is collective radial flow on top of the thermal chaotic contribution , so that that the measured inverse slope $`T=T_{fo}+m<\beta _T>^2`$. Particle spectra are then reproduced with only two parameters, characterizing the system at kinetic freeze-out, the freeze-out temperature $`T_{fo}`$ and the collective radial flow velocity $`<\beta _T>`$. $`T_{fo}`$ is basically constant at $`130MeV`$ from AGS to RHIC as expected if particles decouple at a given density. A constant radial flow velocity of $`<\beta _T>0.4`$ has been established at AGS and SPS energies, whereas the first RHIC results point to a higher value of $`0.6`$ indicating higher initial pressure in the system. The resulting larger inverse slopes lead to an interesting feature not observed at lower energies. The $`\overline{p}`$ and $`\pi ^{}`$ (as well as the $`p`$ and $`\pi ^+`$) yields become comparable at $`p_T2GeV/c`$.

At the SPS, the $`\mathrm{\Omega }`$ and $`J/\psi `$ (and probably also the $`\mathrm{\Xi }`$) particles follow a different pattern . They exhibit stronger slopes than expected from the previously quoted formula. This is taken as evidence that these particles decouple earlier from the system due to their relatively large mean free path. It will be very interesting to see whether similar data from RHIC will follow the same pattern.

### 3.3 Baryon yields

Baryon yields have consistently attracted interest: (i) the net baryon yield at mid-rapidity is an important observable allowing to study the amount of baryon stopping (the transport of participating nucleons to mid-rapidity) that occurs early in the collision. (ii) the mechanism of baryon stopping is a an unsolved question and a topic of recent debate . (iii) the total baryon density plays an important role in connection with in-medium modification of light vector mesons and low-mass dilepton production as emphasized recently  (see Section 4).

The $`\overline{p}/p`$ and $`\overline{\mathrm{\Lambda }}/\mathrm{\Lambda }`$ ratios increase dramatically from values of $`0.1`$ at SPS to values of $`0.7`$ at RHIC . Clearly the ratios do not yet reach unity. The system approaches net baryon free conditions, but remarkably, still $`5\%`$ of the participating nucleons are transported over 5 units of rapidity as shown in Table 2.

From the inclusive (i.e. non-corrected for feed-down) measurement of $`p`$ and $`\overline{p}`$ one can get a crude estimate of the total baryon density (under the assumption that $`p`$ and $`n`$ have the same stopping probability, that $`p\overline{p}`$ and $`n\overline{n}`$ pairs are produced in equal quantities and that higher mass baryons (mainly lambdas and sigmas) finally end up as $`p,\overline{p},n`$ or $`\overline{n}`$). The relevant quantities for RHIC are listed in Table 2. A similar estimate at

the maximum SPS energy of $`\sqrt{s_{NN}}=17GeV`$, derived from the feed-down corrected $`p`$, $`\overline{p}`$ yields and the $`\mathrm{\Lambda }`$, $`\overline{\mathrm{\Lambda }}`$ yields  gives a total baryon density of $`dN_B/dy=110`$. It appears therefore that, contrary to previous expectations, the total baryon density is approximately the same at SPS and RHIC, the strong decrease in nuclear stopping being compensated by copious production of baryon-antibaryon pairs.

## 4 LOW-MASS $`e^+e^{}`$ PAIRS and PHOTONS

Dileptons and photons are unique probes to study the dynamics of relativistic heavy-ion collisions . The interest stems from their relatively large mean free path. As a consequence, they can leave the interaction region without final state interaction, carrying information about the conditions and properties of the matter at the time of their production and in particular at the early stages when deconfinement and chiral symmetry restoration have the best chances to occur.

The prominent topic of interest, both in dileptons and photons, is the identification of thermal radiation emitted from the collision system. This radiation should tell us the nature of the matter formed, a quark-gluon plasma (QGP) or a high-density hadron gas (HG). The physics potential of low-mass dileptons is further augmented by their sensitivity to chiral symmetry restoration. The $`\rho `$-meson is the prime agent here. Due to its very short lifetime ($`\tau =1.3fm/c`$) compared to the typical fireball lifetime of $`10fm/c`$ at SPS energies, most of the $`\rho `$ mesons decay inside the interaction region providing a unique opportunity to observe in-medium modifications of particle properties (mass and/or width) which might be linked to chiral symmetry restoration. The situation is different for the $`\omega `$ and $`\varphi `$ mesons. Because of their much longer lifetimes they predominantly decay outside the interaction region after having regained their vacuum properties. The $`\omega `$ and $`\varphi `$ mesons remain nevertheless important messengers: the undisturbed $`\omega `$ can serve as a reference and the $`\varphi `$ with its $`s\overline{s}`$ content is a probe of strangeness production.

For the moment, the stage for low-mass dileptons belongs to CERN and more specifically to the CERES experiment. CERES has systematically studied the production of $`e^+e^{}`$ pairs in the mass region $`m1.0GeV/c^2`$, with measurements of p-Be, p-Au, , S-Au  and Pb-Au . Whereas the $`p`$ data are well reproduced by the known hadronic sources, a strong enhancement is observed in the mass region $`200<m<600MeV/c^2`$ both in the S and Pb data with respect to those sources scaled to the nuclear case with the event multiplicity. Fig. 6 (left panel) shows the results from Pb-Au collisions at $`40AGeV`$ with an enhancement factor of $`5.1\pm 1.3(stat.)`$ relative to the hadronic cocktail.

Within uncertainties, this is consistent with, or even larger than, the enhancement factor of $`2.9\pm 0.3(stat.)\pm 0.6(syst.)`$ obtained from the combined 95-96 Pb-Au data at $`158AGeV`$. Further studies on the latter demonstrate that the enhancement is more pronounced at low pair $`p_T`$ and increases faster than linearly with the event multiplicity. In all cases, with the S and Pb beam, the enhancement sets at $`m2m_\pi `$.

The enhancement of low-mass dileptons has triggered a wealth of theoretical activity (for a comprehensive review see ). There is consensus that a simple superposition of $`pp`$ collisions cannot explain the data and that an additional source is needed. The pion annihilation channel ($`\pi ^+\pi ^{}\rho l^+l^{}`$), obviously irrelevant in $`pp`$ collisions, has to be added in the nuclear case. This channel accounts for a large fraction of the observed enhancement (see line cocktail + free $`\pi \pi `$ in Fig. 6 right panel) and provides first evidence of thermal radiation from a dense hadron gas. However, in order to quantitatively reproduce the data in the mass region $`0.2<m_{e^+e^{}}<0.6GeV/c^2`$, it was found necessary to include in-medium modifications of the $`\rho `$ meson. Li, Ko and Brown , following the original Brown-Rho scaling , proposed a decrease of the $`\rho `$-mass in the high baryon density of the fireball, as a precursor of chiral symmetry restoration, and achieved excellent agreement with the CERES data (see the line cocktail + dropping mass in Fig. 6 right panel).

Another avenue, based on effective Lagrangians, uses the broadening of the $`\rho `$-meson spectral function resulting from its propagation in the medium, mainly from scattering off baryons , and achieves also an excellent reproduction of the CERES data (see line cocktail + in-medium $`\pi \pi `$ in Fig. 6 right panel). The success of these two different approaches, one relying on quark degrees of freedom and the other on a pure hadronic model, has attracted much interest. Rapp raised, and provided support for, the hypothesis of quark-hadron duality down to low-masses by showing that in a high density state the dilepton production rates calculated with the in-medium $`\rho `$-meson spectral function are very similar to the $`q\overline{q}`$ annihilation rates calculated in pQCD . It is interesting to note that the $`40AGeV`$ data is also equally well reproduced by the two approaches . The accuracy of the data does not allow to discriminate among the two. The CERES run of year 2000 taken with improved mass resolution and good statistics may allow to do that.

The extension of these studies to RHIC energies should be very interesting. The total baryon density, which is the key factor responsible for in-medium modifications of the $`\rho `$ meson both in the dropping mass and the collision broadening scenarios, is almost as high at RHIC as at SPS (cf. Section 3.3), contrary to previous expectations. Updated RHIC predictions that incorporate this and other acquired knowledge on global and hadronic observables, show indeed that the enhancement of low-mass $`e^+e^{}`$ pairs persist at the collider with at least comparable strength . The calculations further predict in-medium modifications of the $`\omega `$ and $`\varphi `$ mesons. These are much less dramatic than in the case of the $`\rho `$ meson but should nevertheless be readily observable with the excellent mass resolution of the PHENIX detector. First results on single electrons from PHENIX are already available  and results on the $`\varphi `$ meson from the second RHIC run are expected soon. The combined analysis of $`\varphi K^+K^{}`$ and $`\varphi e^+e^{}`$ within the same apparatus should provide a very powerful diagnostic tool to make evident in-medium effects.

Direct photons are expected to provide analogous information to thermal dileptons. However, the physics background for real photons is larger by orders of magnitude compared to dileptons making the measurement of photons much tougher and less sensitive to a new source. And indeed, only a small direct photon signal of $`20\%`$ has been observed in central Pb-Au collisions by the WA98 experiment . The better conditions at RHIC, larger particle densities and hence larger initial temperatures should provide a more compelling evidence of this key signal.

## 5 $`J/\psi `$ SUPPRESSION

The melting of charmonium states in a deconfined state of matter is one of the earliest signatures of QGP formation  and has provided one of the most exciting sagas of the SPS programme. Already in the first runs with O and S beams, the NA38/NA50 experiment observed a suppression of $`J/\psi `$ which immediately captured the interest. Intensive theoretical efforts were devoted to explain the effect within conventional physics

and it soon appeared that all experimental data, including systematic studies of $`pp`$, $`pA`$ and $`SU`$ collisions, can be reproduced by invoking the absorption in nuclear medium of the $`c\overline{c}`$ pair before it forms a $`J/\psi `$, with a cross section $`\sigma _{abs}=6.4\pm 0.8mb`$ .

However, a different behavior is observed with the $`158AGeV`$ Pb beam. Whereas peripheral collisions seem to follow the regular absorption pattern, an anomalously larger suppression occurs at central collisions characterized by impact parameters $`b<8fm`$ or $`E_T>40GeV`$ which can be translated into a Bjorken energy density $`ϵ>2.3GeV/fm^3`$ (see Fig. 7) . The data in this figure exhibit a two-step suppression pattern which NA50 attributes to the successive melting of the $`\chi _c`$ and directly produced $`J/\psi `$ mesons in a quark-gluon plasma scenario.

The same data normalized to the Drell-Yan cross section are shown in Fig. 8 as function of $`E_T`$. The two-step pattern can be discerned here at $`E_T`$ values of $`30`$ and $`100GeV`$. Most published calculations, based on conventional hadronic models including the effect of absorption by comovers, fail to reproduce the results as illustrated in the left panel of Fig. 8. The recent calculations of Capella et al.  are in much better agreement with the data and fail only at the most central collisions (middle panel).

On the other hand, quite good agreement over the whole $`E_T`$ range is achieved by models assuming QGP formation (right panel) . The model of Blaizot et al. reproduce the data remarkably well by invoking $`J/\psi `$ suppression whenever the energy density exceeds the critical value for deconfinement (first step) together with fluctuations of the transverse energy for the most central collisions (second step) .

The extrapolation of these results to RHIC energies is not at all clear and ranges from total suppression up to enhancement !  thus ensuring another exciting chapter on the $`J/\psi `$ saga. First glimpse on $`J/\psi `$ production at RHIC is expected from the PHENIX results of year 2 run.

## 6 SUPPRESSION OF LARGE $`p_T`$ HADRONS

With the energies available at RHIC, one order of magnitude higher than at the SPS, new channels and probes become accessible for matter diagnostics. In particular, energy loss through gluon radiation of high $`p_T`$ partons resulting from initial hard scattering and jet production has been predicted as a possible signature of deconfined matter. The phenomenon, commonly referred to as jet quenching, has been extensively studied over the last decade  and should manifest itself as a suppression of high $`p_T`$ particles. Such a suppression has been observed already in the first low-luminosity RHIC run by the two large experiments STAR  and PHENIX  and certainly constitutes the most interesting result from RHIC so far. The suppression is evidenced by plotting the so-called nuclear modification factor defined as the ratio of $`AA`$ to $`pp`$ $`p_T`$ spectra, scaled by the number of binary collisions:

$$R_{AA}(p_T)=\frac{d^2\sigma _{AA}/dp_Td\eta }{<N_{bin}>d^2\sigma _{pp}/dp_Td\eta }$$

In the absence of any new physics this ratio should be equal to 1 at the high $`p_T`$ characteristic of hard processes. At low $`p_T`$, dominated by soft processes which scale with the number of participants, the ratio is expected to be lower, e.g. for central collisions $`R_{AA}(0)0.5N_{part}/N_{bin}0.2`$. The STAR and PHENIX results for central Au-Au collisions at $`\sqrt{s_{NN}}=130GeV`$ are shown in Fig. 9. In both cases the $`pp`$ data at $`130GeV`$ were derived from interpolation of $`pp`$ data  at lower and higher energies.

STAR shows the ratio for negative hadrons and PHENIX for negative hadrons as well as identified $`\pi ^o`$. As expected, for low $`p_T`$ the ratio is small, $`0.2`$. It raises with $`p_T`$, getting close to unity and then falls off at higher $`p_T`$ in the STAR data whereas it seems to flatten in the PHENIX data. The suppression appears to be stronger for $`\pi ^o`$ than for charged particles. This implies that the ratio $`\pi ^o/h`$ is considerably smaller than unity and reflects the large $`p,\overline{p}`$ contribution to the charged particle spectra at high momentum (cf. Section 3.2). A similar suppression pattern appears when the ratio of central to peripheral collisions, each divided by its corresponding value of $`N_{bin}`$, is plotted as function of $`p_T`$ . This behaviour is in marked contrast to the observations in Pb-Pb and Pb-Au collisions at the SPS (the solid lines in Fig. 9 right panel represent the band of uncertainty of the CERN results) where the ratio overshoots unity and saturates at high $`p_T`$, like in pA collisions  due to the well known Cronin effect . The high $`p_T`$ suppression observed at RHIC can be quantitatively reproduced in terms of jet quenching if an average energy loss of $`dE/dx0.25GeV/fm`$ is assumed within a parton model.

In the context of the new RHIC data, jet quenching has been discussed in three different topics (cf. Sections 2.1, 2.2 and this Section). Whereas the phenomenon helps describing the observed effects in elliptic flow and high $`p_T`$ suppression, the expected increase in particle production has not been observed (cf. Section 2.1). Fine tuning of the theoretical treatment to achieve consistent description of all these observables together with new RHIC data allowing to reach much higher $`p_T`$ values and reference data on $`pp`$ and $`pA`$ measured within the same apparatus will be very valuable to consolidate these intriguing results.

## 7 CONCLUSION

In a relatively short run, the RHIC experiments have produced an impressive amount of results and much more is expected over the next years. The second RHIC run has just been completed with a recorded luminosity almost two orders of magnitude larger than in year-1. This, together with the still ongoing yield of interesting results from the SPS programme, places the field of relativistic heavy-ion collisions at a very unique phase with exciting physics prospects.


text:
# Gluon Correlators in the Kogan–Kovner Model

## 1 Introduction

Strongly interacting non-Abelian gauge theories require non-perturbative methods. One such, the application of the variational approximation to the Schrödinger formulation, goes back to the days of Feynman . Though such methods are approximate, they incorporate many of the desirable features of non-Abelian gauge theories, in particular the gluon condensate, asymptotic freedom, instantons, dynamical mass generation and confinement . Variational methods have also met success applied to Abelian theories .

In the Schrödinger formulation of pure $`SU(N)`$ Yang–Mills theory, the wavefunctional $`\mathrm{\Psi }[A_i^a(x)]`$ is an eigenstate of the Hamiltonian

$$H=d^3x\frac{1}{2}[E_i^{a2}(x)+B_i^{a2}(x)],$$
(1)
where $`E_i^a(x)=i\frac{\delta }{\delta A_i^a(x)}`$, $`B_i^a(x)=\frac{1}{2}ϵ_{ijk}(_jA_k^a(x)_kA_j^a(x)+gf^{abc}A_j^bA_k^c)`$ and physical states satisfy Gauss’ law as a constraint:

$$[_iE_i^a(x)gf^{abc}A_i^bE_i^c]\mathrm{\Psi }[A]=0.$$
(2)
Kogan and Kovner  proposed the gauge invariant Gaussian ansatz

$$\mathrm{\Psi }[A_i^a(x)]=DU\mathrm{exp}\left[\frac{1}{2}d^3xd^3yA_i^{Ua}(x)\delta ^{ab}\delta _{ij}G^1(xy)A_j^{Ub}(y)\right],$$
(3)
where under a local gauge transformation $`U(x),`$

$$A_i^a(x)A_i^{Ua}(x)=S^{ab}(x)A_i^b(x)+\lambda _i^a(x),$$
(4)
and $`DU`$ is the $`SU(N)`$ group invariant measure. Here $`S^{ab}=\frac{1}{2}\mathrm{tr}(\tau ^aU^{}\tau ^bU)`$, $`\lambda _i^a=\frac{i}{g}\mathrm{tr}(\tau ^aU^{}_iU)`$ and $`\frac{\tau ^a}{2}`$ form an $`N\times N`$ Hermitian representation of $`SU(N)`$: $`[\frac{\tau ^a}{2},\frac{\tau ^b}{2}]=if^{abc}\frac{\tau ^c}{2}`$ with normalisation $`\mathrm{tr}(\tau ^a\tau ^b)=2\delta ^{ab}.`$ The $`A`$-field propagator has a mass parameter $`M`$, with respect to which the Hamiltonian is minimised. The authors found that $`M`$ is of the correct size to generate the observed gluon condensate .

In sections 3 and 4 of this paper, the Lorentz invariant gluon correlators $`F^2(x)F^2(x^{})`$ and $`F\stackrel{~}{F}(x)F\stackrel{~}{F}(x^{})`$ are calculated using the same ansatz, and their scaling properties at large separations are analysed. Naïvely, one expects that the correlators will decay as $`\mathrm{exp}M|xx^{}|`$, with mass scale $`M`$. Firstly though, let us recall the details of the ansatz and calculation of the correlation functions therein.

## 2 Correlator Rubric

To calculate correlation functions, consider the generating functional

$$Z=DA\mathrm{\Psi }[A]\mathrm{\Psi }[A];$$
(5)
the correlator of operator $`O`$ is then

$$O=\frac{1}{Z}DA\mathrm{\Psi }[A]O\mathrm{\Psi }[A].$$
(6)
For gauge invariant $`O`$, one of the gauge integrations is redundant. Since the gauge transformation is linear on the gauge fields, the gauge field path integral is Gaussian. Using matrix notation for integration, *viz.*

$$(AB)_{ik}^{ac}(x,z)=\mathrm{\Sigma }_{b,j}d^3yA_{ij}^{ab}(x,y)B_{jk}^{bc}(y,z),$$
(7)
the generating functional is proportional to

$$DADU\mathrm{exp}[\frac{1}{2}(A+a)^T(A+a)]\mathrm{exp}[\frac{1}{2}\lambda ^T(G+SGS^T)^1)\lambda ].$$
(8)
Here, matrices are defined as

$`S_{ij}^{ab}(x,y)`$ $`=S^{ab}(x)\delta (xy)\delta _{ij},`$
$`G_{ij}^{ab}`$ $`=G(xy)\delta ^{ab}\delta _{ij},`$
$``$ $`=S^TG^1S+G^1\mathrm{and}`$
$`a`$ $`=^1S^TG^1\lambda .`$ (9)
Note that $`S`$ is orthogonal whereas $``$ and $`G`$ are symmetric.

To calculate a correlator, substitute the operator in (6) to obtain the correlator of some string of $`A`$ fields. Then do the $`A`$-field path integral; this is a free field theory. Thus

$$A_i^a(x)A_j^b(y)=a_i^a(x)a_j^b(y)+_{ij}^{1ab}(x,y),$$
(10)
and so on via Wick’s theorem.

This done, one is left with an integral over gauge transformations $`U`$. The stratagem here is to factor the gauge transformation $`U`$ into a product $`U_LU_H`$ with parts dependent on high or low momentum modes of magnitude $`k`$: $`U_H`$ with $`k>M`$ and $`U_L`$ with $`k<M`$ respectively. Both the measure and action factorise similarly, since the group measure is invariant and the high and low modes are orthogonal. The $`A`$-field propagator takes the form $`G^1(k)=\sqrt{k^2+M^2}`$. For $`k>M`$, take $`G^1(k)k`$ and treat the theory perturbatively. With $`U=\mathrm{exp}[ig\varphi ^a\tau ^a]`$, to $`O(g^2)`$ the high mode action takes the form

$$\mathrm{\Gamma }_H[\varphi ]=\frac{1}{4}d^3xd^3y_i\varphi ^a(x)G^1(xy)|_{k>M}_i\varphi ^a(y).$$
(11)
One again obtains a free field, with propagator

$$\varphi ^a(x)\varphi ^b(y)=2\delta ^{ab}[_i^x_i^yG^1(xy)|_{k>M}]^12\delta ^{ab}H^1(xy).$$
(12)
For $`k<M`$, take $`G^1(k)M`$ such that $`G^1(xy)|_{k<M}`$ is local; the matrices $`U_L`$ vary only on scales greater than $`1/M`$ and can be considered constant in the action.111In the action, $`U_L`$ are multiplied by $`G^1(x)0`$ for $`|x|>1/M`$ The resulting leading order low mode action is

$$\mathrm{\Gamma }_L[U]=\frac{M}{2g^2}\mathrm{tr}d^3x_iU_L^{}(x)_iU_L(x).$$
(13)
The integration over high modes amounts to a renormalization group transformation. The effect on the low modes is to replace $`g`$ in (13) by the coupling at the scale $`M`$. This sigma model is calculated in the mean field approximation: the action is replaced by

$$\frac{M}{2g^2(M)}\mathrm{tr}d^3x\left[_iU_L^{}(x)_iU_L(x)+\sigma ^2(U_L^{}U_L1)\right],$$
(14)
subject to the constraints

$$U_L^{}U_L=1,$$
(15)
$$\sigma U_L=0.$$
(16)
So the $`U_{mn}(x)`$ act as $`N^2`$ free charged scalar fields with propagator

$$U_{Lkl}^{}(x)U_{Lmn}(y)=\frac{2g^2}{M}\delta _{kn}\delta _{lm}(_i_i+\sigma ^2)^1\delta (xy).$$
(17)
The momentum split allows the further simplifications

$`S^{ab}(x)`$ $`=S_L^{ab}(x),`$
$`^{ab}(x,y)`$ $`=2\delta ^{ab}G^1(xy),`$
$`\lambda _i^a(x)`$ $`=\lambda _{iL}^a(x)+\lambda _{iH}^a(x),`$
$`a_i^a(x)`$ $`={\displaystyle \frac{1}{2}}\lambda _{iL}^a(x)+{\displaystyle \frac{1}{2}}\lambda _{iH}^bS_L^{ba}(x).`$ (18)
Before tackling specific correlators, let us note also the following points:

* The non-perturbative physics is captured in the low mode action; elsewhere we work to $`O(g)`$.
* A factor $`S_L`$ in a correlator contributes $`U_L^{}U_L`$ in the Wick expansion, and therefore a factor of $`g^2`$, from (17). The only leading order terms are those in which the orthogonality of the $`S_L`$ can be used, such that they cancel.
* Orthogonality of the $`S_L`$ can be used whenever they are separated by a scale less than $`1/M`$. Thus, wherever, $`S_L`$ appears, it is sufficient to assume orthogonality, and integrate up to distance scale $`1/M`$. Contributions on larger scales are higher order in $`g`$.
* The action is quadratic in the decoupled $`\lambda _L,\lambda _H`$; correlators with an odd number of $`\lambda _L`$ or $`\lambda _H`$ vanish.
* A correlator with an odd number of $`a`$ and *ergo* $`A`$ will therefore vanish to leading order.
* Since $`\lambda _{iH}^a`$ is just $`_i\varphi ^a`$, $`ϵ_{ijk}_j\lambda _{kH}^a=0`$.
* $`\lambda _L`$ is pure gauge and Abelian magnetic terms like $`(ϵ_{ijk}_j\lambda _{kL}^a)^2`$ vanish to $`O(g)`$.
* Thus, any correlator with two or more $`ϵa`$ terms must vanish to $`O(g)`$.
* In seeking the scaling properties of correlators, constants and contact terms may be discarded.

## 3 Scalar gluon correlator

The Lorentz invariant gluon correlators are $`F^2(x)F^2(x^{})`$ and $`F\stackrel{~}{F}(x)F\stackrel{~}{F}(x^{})`$, where

$$F^2=F_{\mu \nu }^aF^{a\mu \nu }=2(B_i^{a2}E_i^{a2}),$$
(19)
$$F\stackrel{~}{F}=\frac{1}{2!}ϵ_{\mu \nu \sigma \rho }F_{\mu \nu }^aF_{\sigma \rho }^a=4B_i^aE_i^a.$$
(20)
Discarding a contact term, the scalar correlator is thus

$$F^2(x)F^2(x^{})=4B^2(x)B^2(x^{})+E^2(x)E^2(x^{})2B^2(x)E^2(x^{}).$$
(21)
To $`O(g)`$, the chromo-magnetic term is

$`B_i^{a2}(x)B_i^{{}_{}{}^{}a2}(x)`$ $`=ϵ_{ijk}ϵ_{ilm}ϵ_{ijk}^{^{}}ϵ{}_{}{}^{}{}_{ilm}{}^{}[_jA_k^a_lA_m^a_j^{^{}}A_k^{}_{}{}^{}a_l^{^{}}A_m^{}_{}{}^{}a`$
$`+gf^{abc}_jA_k^aA_l^bA_m^c_j^{^{}}A_k^{}_{}{}^{}a_l^{^{}}A_m^{}_{}{}^{}a+gf^{{}_{}{}^{}abc}{}_{}{}^{}{}_{j}{}^{}A_{k}^{{}_{}{}^{}a}A_l^{}_{}{}^{}bA_m^{}_{}{}^{}c_jA_k^a_lA_m^a],`$ (22)
where a prime on a field indicates that its parameters carry primes. The terms of $`O(g)`$ have five $`A`$-fields and these vanish. Using (10), the only terms with less than two $`ϵa`$s contain $``$s only. To wit,

$`B_i^{a2}(x)B_i^{{}_{}{}^{}a2}(x)=ϵ_{ijk}ϵ_{ilm}ϵ_{ijk}^{^{}}ϵ_{ilm}^{^{}}[`$ $`_j^x_l^y_{km}^{1aa}(x,y)_j^{}_{}{}^{}x_l^{}_{}{}^{}y_{km}^{{}_{}{}^{}1aa}(x,y)|_{x=y}`$
$`+`$ $`_j^x_j^{}_{}{}^{}x_{kk^{}}^{1aa^{}}(x,x^{})_l^x_l^{}_{}{}^{}x_{mm^{}}^{1aa^{}}(x,x^{})`$
$`+`$ $`_j^x_l^{}_{}{}^{}x_{km^{}}^{1aa^{}}(x,x^{})_l^x_j^{}_{}{}^{}x_{mk^{}}^{1aa^{}}(x,x^{})].`$ (23)
Using $`ϵ_{ijk}ϵ_{ilm}=\delta _{jl}\delta _{km}\delta _{jm}\delta _{kl}`$ and $`_{ij}^{ab}(x,y)=2\delta ^{ab}\delta _{ij}G^1(xy)`$ yields (restoring primes and discarding a constant)

$`+{\displaystyle \frac{1}{2}}`$ $`(N^21)_j^x_j^{}^x^{}G(xx^{})_j^x_j^{}^x^{}G(xx^{})`$
$`+{\displaystyle \frac{1}{2}}`$ $`(N^21)_j^x_j^x^{}G(xx^{})_j^{}^x_j^{}^x^{}G(xx^{}).`$ (24)
How do these terms scale? Recall

$$G(x)=_{\mathrm{}}^{\mathrm{}}\frac{d^3k}{(2\pi )^3}G(k)e^{ik.x},$$
(25)
where $`G(k)^1=\sqrt{k^2+M^2}`$. Choose spherical polars with $`k_3`$-axis along $`x`$. Then

$$G(x)=\frac{1}{2\pi ^2x}_0^{\mathrm{}}𝑑k\frac{k}{\sqrt{k^2+M^2}}\mathrm{sin}kx=\frac{1}{4\pi ^2ix}_{\mathrm{}}^{\mathrm{}}𝑑k\frac{k}{\sqrt{k^2+M^2}}e^{ikx}.$$
(26)
The contour can be deformed in $``$ to surround the branch cut going from $`iM`$ to $`\mathrm{}`$ (*cf.* \[12, pp. 27-28\]) and for large $`x`$ approaches

$$\frac{1}{2\pi ^2x}_M^{\mathrm{}}𝑑y\frac{ye^{xy}}{\sqrt{y^2M^2}}=\frac{M}{2\pi ^2x}K_1(Mx)(\frac{M}{8\pi ^3x^3})^{1/2}e^{Mx}.$$
(27)
Similarly,

$$G^1(x)=\frac{1}{4\pi ^2ix}_{\mathrm{}}^{\mathrm{}}𝑑kk\sqrt{k^2+M^2}e^{ikx}=\frac{M^2}{2\pi ^2x^2}K_2(Mx)(\frac{M^3}{8\pi ^3x^5})^{1/2}e^{Mx},$$
(28)
where $`K_{1,2}(Mx)`$ are modified Bessel functions .

Carrying out the differentiation in (3), one finds that

$$4B^2(x)B^2(x^{})\frac{4(N^21)M^5e^{2M|xx^{}|}}{(2\pi )^3|xx^{}|^3},$$
(29)
where the expected mass scaling is indeed present.

Consider now the chromo-electric correlator. Here things become more complicated. Doing the functional differentiation in (6), one obtains

$`E_i^{a2}(x)E_i^{{}_{}{}^{}a2}(x)`$ $`=6(N^21)[G^1(xx^{})]^2`$
$`4G^1(xx^{}){\displaystyle d^3yd^3y^{}G^1(y^{}x^{})G^1(yx)A_i^a(y)A_i^a(y^{})}`$
$`+{\displaystyle d^3yd^3y^{}d^3zd^3z^{}G^1(y^{}x^{})G^1(z^{}x^{})G^1(yx)G^1(zx)}`$
$`A_i^{}^a^{}(y^{})A_i^{}^a^{}(z^{})A_i^a(y)A_i^a(z),`$ (30)
plus constants. Using (10) and (2), and the orthogonality properties of the $`S_L`$, this becomes

$`{\displaystyle \frac{3}{2}}(N^21)[G^1(xx^{})]^2`$
$``$ $`{\displaystyle \frac{1}{2}}G^1(xx^{}){\displaystyle d^3yd^3y^{}G^1(y^{}x^{})G^1(yx)[\lambda _{iL}^a(y^{})\lambda _{iL}^a(y)+\lambda _{iH}^a(y^{})\lambda _{iH}^a(y)]}`$
$`+`$ $`{\displaystyle \frac{1}{16}}{\displaystyle d^3yd^3y^{}d^3zd^3z^{}G^1(y^{}x^{})G^1(z^{}x^{})G^1(yx)G^1(zx)}`$
$`[\lambda _{iL}^{}_{}{}^{}a(y)\lambda _{iL}^{}_{}{}^{}a(z)\lambda _{iL}^a(y)\lambda _{iL}^a(z)+\lambda _{iH}^{}_{}{}^{}a(y)\lambda _{iH}^{}_{}{}^{}a(z)\lambda _{iH}^a(y)\lambda _{iH}^a(z)+4\lambda _{iH}^{}_{}{}^{}b(y)S_{iL}^{{}_{}{}^{}ba}(y)\lambda _{iL}^{}_{}{}^{}a(z)\lambda _{iH}^b(y)S_{iL}^{ba}(y)\lambda _{iL}^a(z)].`$ (31)
The inverse relation $`d^3yG^1(xy)G(yz)=\delta ^3(xz)`$ has been used.

The averaging over the gauge transformations $`U`$ remains. The orthogonality of the $`k`$-modes ensures the following property: in any integral of the form $`d^3xf(x)g(x)`$, where $`g`$ contains only high (low) mode contributions, only high (low) mode parts of $`f`$ contribute to the integral. Then, since $`G^1(xy)|_{k<M}M\delta (xy)`$ is local, the integrations in (3) involving $`\lambda _L`$ are trivial.

The integrations involving $`\lambda _H`$ are a little more tricky. Consider, for example, the simplest one:

$$\frac{1}{2}G^1(xx^{})d^3yd^3y^{}G^1(y^{}x^{})G^1(yx)\lambda _{iH}^a(y^{})\lambda _{iH}^a(y).$$
(32)
Using (12), this is

$$(N^21)G^1(xx^{})d^3yd^3y^{}G^1(y^{}x^{})G^1(yx)_i^y_i^y^{}H^1(yy^{}).$$
(33)
Fourier decomposition of the propagators yields the simple result

$$(N^21)G^1(xx^{})G^1(xx^{})|_{k>M}.$$
(34)
The quartic term is similar. Since

$`\varphi ^{}_{}{}^{}a(y)\varphi ^{}_{}{}^{}a(z)\varphi ^a(y)\varphi ^a(z)=4[(N^21)^2H^1(y^{}z^{})H^1(yz)`$
$`+(N^21)H^1(y^{}y)H^1(z^{}z)+(N^21)H^1(y^{}z)H^1(z^{}y)],`$ (35)
the terms dependent on $`xx^{}`$ become (Fourier decomposing again)

$$\frac{(N^21)}{4}_{k>M}\frac{d^3k}{(2\pi )^3}\frac{d^3k^{}}{(2\pi )^3}\frac{G^1(k)G^1(k^{})}{k^2k^{}_{}{}^{}2}(k.k^{})^2e^{i(k+k^{}).(xx^{})}.$$
(36)
These integrals can be evaluated, but are not very illuminating, amounting to Fourier transforms of various simple powers of $`k`$. If the theory is regularised with a cut-off, these are just spurious oscillatory functions of $`M|xx^{}|`$. One obtains the same generic oscillatory result if one tries to calculate $`G^1(xx^{})|_{kM}`$ explicitly; what is important for our purposes is that the only mass scale present is $`M`$.

The remaining term in $`\lambda _H`$ in (3) is quadratic. The high mode propagator contributes a factor $`\delta ^{bb^{}}`$ and the orthogonality property of the $`S_L`$ can then be used.

The low mode expectation values remain in (3). The simplest terms are quadratic:

$$\lambda _{iL}^a(x^{})\lambda _{iL}^a(x).$$
(37)
The Fierz identity for $`SU(N)`$: $`\tau _{ij}^a\tau _{kl}^a=2(\delta _{il}\delta _{jk}\frac{1}{N}\delta _{ij}\delta _{kl})`$ implies

$$\lambda _{iL}^a(x^{})\lambda _{iL}^a(x)=\frac{2}{g^2}\mathrm{tr}(U_L^{}(x^{})_iU_L(x^{})U_L^{}(x)_iU_L(x)).$$
(38)
The $`U_L`$ contain modes with $`k<M`$; for $`xx^{}`$, the unitarity of the $`U_L`$ cannot be invoked. Doing the Wick expansion, one sees that these term are $`O(g^2)`$.

The quartic low mode term is (using (17))

$`{\displaystyle \frac{M^4}{16}}\lambda _{iL}^{}_{}{}^{}a(x)\lambda _{iL}^{}_{}{}^{}a(x)\lambda _{iL}^a(x)\lambda _{iL}^a(x)`$
$`=`$ $`{\displaystyle \frac{M^4}{4g^4}}\mathrm{tr}(_i^{^{}}U(x^{})_i^{^{}}U^{}(x^{}))\mathrm{tr}(_iU(x)_iU^{}(x))`$
$`=`$ $`2M^2N^2\left[{\displaystyle _{k<M}}{\displaystyle \frac{d^3k}{(2\pi )^3}}{\displaystyle \frac{k_ik_i^{}}{k^2+\sigma ^2}}e^{ik.(xx^{})}\right]^2.`$ (39)
This low mode term is very interesting. Consider the integral

$$_{k<M}\frac{d^3k}{(2\pi )^3}\frac{e^{ik.x}}{k^2+\sigma ^2}=\frac{1}{4\pi ^2iy}_M^M𝑑k\frac{ke^{iky}}{k^2+\sigma ^2},$$
(40)
where ($`y=|xx^{}|`$). In the limit $`M\mathrm{}`$ this is a simple contour integral. Closing the contour in the upper half plane, and enclosing the pole at $`k=i\sigma `$ one obtains

$$\frac{1}{4\pi y}e^{\sigma y}.$$
(41)
Moreover, one can evaluate the integral for finite $`M`$ using the exponential integral function (*cf.* \[14, p. 228\]); one finds precisely this term, plus complicated terms involving $`M`$. But this term takes the form of a mass term: $`\sigma `$ appears as an *additional* mass scale in gluon correlation functions.

So now there are two mass scales, $`M`$ and $`\sigma `$, in marked contrast to what one expects. How do we get ourselves out of this quandary? One way is to choose $`\sigma M`$. Then the mass scale $`\sigma `$ is present and yet impossible to detect. The other way out is to force $`\sigma `$ to be zero.

It seems then, that if one makes a somewhat arbitrary choice of $`\sigma =0`$ or $`M`$ to fit the expected behaviour, then things work out satisfyingly. But these choices look distinctly *less* arbitrary in the context of Kogan’s and Kovner’s model. They consider the action (13) as a statistical mechanical system with $`SU(N)\times SU(N)`$ symmetry at temperature $`g^2(M)`$. Since QCD is asymptotically free, large $`M`$ corresponds to low temperatures, and conversely. One expects that the system will undergo a phase transition from an ordered state with spontaneous symmetry breaking to a disordered state with the full symmetry restored at some critical temperature. One is then left to speculate as to the critical temperature (coupling) and the nature of the phase transition. In the context of mean field theory, $`\sigma `$ acts as a mass gap; for a second order transition it vanishes continuously at the critical point, and is otherwise non-zero.

Kogan and Kovner argue that the system’s internal energy will be minimised at, or close to, the critical point, with $`\sigma `$ taking its critical value, corresponding to the mass gap. We have argued that $`\sigma `$ can only take the values zero or $`M`$. Thus, the case $`\sigma =0`$ corresponds to no mass gap, and a second order phase transition, as predicted by the mean field theory. The case $`\sigma M`$ corresponds to a mass gap as big as the ultraviolet cut-off in the low mode model: clearly this corresponds to a strongly first order phase transition. Intermediate behaviour (*e.g.* a weakly first order transition) has observable consequences for the gluon correlators.

Collecting up the terms generated by (3), that is (34), (36) and (3), the chromo-electric correlator becomes, for $`\sigma =0`$

$`4E_i^{a2}(x)E_i^{{}_{}{}^{}a2}(x)`$ $`=6(N^21)[G^1(xx^{})]^2`$
$`4(N^21)G^1(xx^{})G^1(xx^{})|_{k>M}`$
$`+\left[8N^2{\displaystyle _{k<M}}+(N^21){\displaystyle _{k>M}}\right]{\displaystyle \frac{d^3k}{(2\pi )^3}}{\displaystyle \frac{d^3k^{}}{(2\pi )^3}}G^1(k)G^1(k^{}){\displaystyle \frac{(k.k^{})^2}{k^2k^{}_{}{}^{}2}}e^{i(k+k^{}).(xx^{})}.`$ (42)
The final part of $`F^2(x)F^2(x^{})`$ is (up to a contact term)

$$8B^2(x)E^2(x^{})=8ϵ_{ijk}ϵ_{ilm}_j^xA_k^a_l^xA_m^ai\frac{\delta }{\delta A_i^{}_{}{}^{}a(x)}i\frac{\delta }{\delta A_i^{}_{}{}^{}a(x)},$$
(43)
where terms of $`O(g^2)`$ or with an odd number of $`A`$-fields have been discarded. Calculating as before, we find that this is a contact term.

After a good deal of effort, we have arrived at an expression for the correlator $`F^2(x)F^2(x^{})`$ for $`\sigma =0`$, namely the sum of (29) and (3). The explicit form of the correlator is not especially illuminating. What *is* significant is that unless $`\sigma =0`$ or $`M`$, there is an additional mass scale in the correlator.

## 4 Pseudo-scalar gluon correlator

Let us proceed to calculate the scaling properties of the pseudo-scalar glueball correlator. In calculating the scalar correlator, the $`\sigma `$ mass scale came from the quartic term in $`\lambda _L`$. Such a term is not present in the pseudo-scalar correlator, suggesting that only the expected mass scale $`M`$ will be present.

Discarding contact terms, the correlator can be written unambiguously as $`16BB^{}EE^{}`$. Doing the functional differentiation in (6) and retaining terms even in $`A`$ of $`O(g)`$ gives

$$16ϵ_{ijk}_j^xA_k^aϵ_{ijk}^{^{}}_j^{}_{}{}^{}xA_k^{}_{}{}^{}a[\delta ^{aa^{}}\delta _{ii^{}}G^1(xx^{})d^3y^{}d^3zG^1(y^{}x^{})G^1(zx)A_i^{}_{}{}^{}a(y)A_i^a(z)].$$
(44)
The first term contains two $`ϵ`$s only, so its only non-zero part has $`AA^1`$, *i.e.*

$$16(N^21)G^1(xx^{})_j^x_j^x^{}G(xx^{}),$$
(45)
with the expected mass scale $`M`$. The second term is more complicated. The possible non-zero parts in the four-point $`A`$ correlator in (44) are

$`_j^xA_k^a_j^{}_{}{}^{}xA_k^{}_{}{}^{}aA_i^{}_{}{}^{}a(y)A_i^a(z)`$ $`=_j^x_j^{}_{}{}^{}x_{kk^{}}^{1aa^{}}(x,x^{})_{i^{}i}^{1a^{}a}(y^{},z)`$
$`+_j^x_{ki^{}}^{1aa^{}}(x,y^{})_j^{}_{}{}^{}x_{k^{}i}^{1a^{}a}(x^{},z)`$
$`+2_j^x_{ki^{}}^{1aa^{}}(x,y^{})_j^{}_{}{}^{}xa_k^{}_{}{}^{}a(x)a_i^a(z)`$
$`+_j^x_j^{}_{}{}^{}x_{kk^{}}^{1aa^{}}(x,x^{})a_i^{}_{}{}^{}a(y)a_i^a(z).`$ (46)
The first term yields $`8(N^21)G^1(xx^{})_j^x_j^x^{}G(xx^{})`$ and the second term turns out to be a contact term.

The third and fourth terms on the RHS of (4) have both high and low mode contributions. For the low mode parts, the integrations in (44) are trivial. they result in a factor $`\lambda _L(x)\lambda _L(x^{})`$, which is $`O(g^2)`$.

The only remaining contribution comes from the high modes in the last term. This is

$$16(N^21)[_j^x_i^x^{}G(xx^{})_{k>M}\frac{d^3k}{(2\pi )^3}\frac{k_ik_j}{k^2}G^1(k)e^{ik.(xx^{})}_j^x_j^x^{}G(xx^{})G^1(xx^{})|_{k>M}].$$
(47)
Thus, the only mass scale in the pseudo-scalar glueball correlator is $`M`$.

## 5 Discussion

In the above, an extensive analysis of the gluon correlators in the Kogan-Kovner ansatz has been made. It has been shown that the gauge and Lorentz invariant gluon correlators contain the expected dynamically generated mass scale $`M`$. Moreover, it has been shown that an additional mass scale $`\sigma `$ is present in the scalar glueball, and should be observable, *unless* it vanishes, or lies close to $`M`$. This constrains the phase transition in the $`\sigma `$-model part of the Kogan–Kovner ansatz to be of second order, or strongly first order, respectively. If this is the case, then the scalar and pseudo-scalar glueball masses are degenerate.

## 6 Acknowledgements

I would like to thank I. I. Kogan and A. Kovner for their comments.


text:
# How the quasispecies evolution depends on the topology of the genome space

## 1 Introduction

The conceptual simplicity and readiness for mathematical description, along with obvious practical relevance makes the Darwinian evolution a prominent subject of theoretical biology. Moreover, the formalisation in terms of stochastic processes reveals close relations to many problems studied before in theoretical physics. Therefore, it is natural to look for physical tools which may answer biologists’ questions about natural evolution. A lot of effort was devoted to this area recently .

The process of the biological evolution consists in three steps, though two of them are simultaneous: $`reproductionmutationselection.`$ The important thing to note is that the *biological fitness* function which denotes individuals’ ability to produce viable offspring is dependent on their phenotype. On the other hand, the mutations occur in the genotype — information stored in a sequence of DNA. The overwhelming complexity of the fitness function, which assigns to each microscopic genotype corresponding macroscopic phenotype, makes the theoretical treatment of evolution an extremely complicated task.

Simplifications of the problem are necessary. One of the most prominent schemes was developed within the *quasispecies* model of the biological evolution, introduced by Eigen . It showed to be plausible and very useful when investigating the mechanisms of microevolution. In the language of physics we can roughly speak of diffusion in a complicated potential. Originally it was introduced as a model of chemical prebiotic evolution, but it is also reasonable to use it for viruses and bacteria which do not haven too complex genome. Among the host of other approaches, let us mention at least the model of coevolution described by Kauffman’s NKCS model .

Very clear and visual presentation of the evolutionary process is to consider it as a random walk in the *fitness landscape*. This is a surface above the genome space, where the altitude of each point is determined by the fitness attributed to the corresponding genotype. In fact, the fitness landscape is not static, as it strongly depends on the ever changing environment, which includes interactions with (co)evolving species as well as abiotic influences . The movement in the fitness landscape is driven by the mutations and optimal genotypes are situated on tops of peaks in the landscape.

A broad set of different fitness landscapes was used recently to study the behaviour of the evolutionary system. These models included both static  and dynamic landscapes . They include the *sharply-peaked landscape* (SPL) , the *Fujiyama landscape* or the *holey landscape* (HL). Several recent reviews summarise various approaches explored .

In the Eigen’s quasispecies model an individual is characterised by a sequence of $`d`$ nucleotides (or loci) $`𝐬=\{s_1,s_2,s_3,\mathrm{},s_d\}`$. In real situation we should consider that there are 4 main nucleotides A, G, C, T(U), so the alphabet contains four letters. Two of them are pyrimidines and the other two are purines. To make the treatment mathematically simpler, we will use only two-letter alphabet $`\{0,1\}`$. We may interpret this simplification in two alternative ways. We can consider putting the digit $`0`$ in the site of the sequence when in the same position in the nucleotide acid is a pyrimidine and $`1`$ if there is a purine. Or, we can consider a set of loci and for each locus suggest exactly two possible alleles, denoted by $`0`$ and $`1`$.

In every position of the sequence a mutation can occur with a probability $`\mu `$. In our work we consider only point mutations (transversions and transitions or an exchange of alleles in a certain locus) which are only a part of all processes which change genotype and consequently the phenotype. For example we neglect duplications which,it seems, are important for our understanding of protein interaction networks, see . The mutation is then an exchange of $`1`$ and $`0`$ in a particular site in the sequence or in a certain locus. The on site mutation rate $`\mu `$ is considered to be uniform for all sites of the sequence, although this is not the case which we see in fact. It is considered that this simplification does not change the results of the model.

The properties of the evolutionary process, which we are going to compare in the article, are the critical coefficients of the error threshold transition. The *error threshold* is the phase transition found in SPL which separates two regimes of the quasispecies evolution — the *adaptive regime* and *wandering regime*. In the adaptive regime the cloud of quasispecies is formed around the *wildtype*, whereas in the wandering regime no wildtype is formed due to a mutation load.

The error threshold phenomenon comes out of the competition between the selective advantage of the wildtype $`\alpha `$ and the mutation load presented by the rate $`\mu `$ and the genome length $`d`$. The threshold is characterised by the specific value of the selective advantage $`\alpha _c`$. From the physical point of view, the error threshold may be thought of as a dynamic phase transition from localised to extended ground-state eigenvector. Recently it was shown  that coevolution of several species may drive the system close to the error threshold. It is therefore highly biologically relevant to study the properties of the error threshold transition in detail.

It is evident that the native geometry of the genome space is the hypercube $`𝒮=\{0,1\}^d`$. So, this is the natural starting point when building a model for a fitness landscape. The main focus will be aimed at the presence or absence of adaptive regime induced by a single maximum in the fitness landscape. Therefore, the first thing to try is the sharply–peaked landscape on the hypercube. This model was solved exactly by Galluccio et al. .

However, the real landscapes are far more complicated. Indeed, it is expected that rugged landscapes with many competing maxima represent a realistic picture . Fortunately enough, such landscapes are well–known in the theory of spin glasses  and a “spin–glass” theory of evolution was investigated .

The main result coming from the spin–glass theory is the ultrametric structure of the minima in the free–energy landscape of spin glasses, which can be translated into ultrametric or hierarchical structure of maxima in the corresponding fitness landscape.

Somewhat more coarse–grained description of the evolution may arise, if we assume the ultrametric structure of fitness maxima. Indeed, we may think of the evolution as a hopping process between various peaks, where the system spends most of its time, while the detailed course of the hops, composed of many individual point mutations, can be accounted for phenomenologically by tuning appropriately the hopping probabilities.

Here, we will be interested again in the error–threshold transition, now interpreted in such a way that all but one of the peaks have the same weight, while the singled–out peak is higher and corresponds to the “wildtype” of the sharply–peaked landscape. Therefore, the treatment can be very similar to SPL from the mathematical point of view, but differs significantly in that now we take into account, at least approximately, the ruggedness of the real fitness landscape.

Yet another approach to the modelling of the fitness landscape is possible, namely considering the holey landscape. Indeed, most of the point mutations which may occur at the basic level of the evolutionary picture are lethal for the individual. Therefore, the hypercube does not represent a good approximation to the evolutionary dynamics, because only few of the hypercube’s edges represent paths to possible new genomes. A sparsely connected set of points selected at some of the hypercube corners is perhaps better choice.

The problem resembles the study of topologically disordered solids, where the random network of bonds is often well modelled by the Bethe lattice. Therefore, third geometry we investigate in our work is that of the Bethe lattice. Again, as in the previous cases, one of the points will be considered as the “wildtype” with larger fitness, while all other points will have uniform fitness (smaller than the wildtype).

## 2 General scheme for sharply–peaked and related landscapes

In the quasispecies model, evolution of the size $`p_i`$ of an asexual population for the genome represented by the lattice site $`i`$ is governed by the equation

$$\dot{p}_i(t)=\underset{j}{}T_{ij}p_j(t).$$
(1)
describing a Markov process in the genome space.

The quantity $`p_i`$ will represent the relative population size, or it will be proportional to the probability of a certain individual to belong to the population in the site $`i`$. In so doing, we tacitly assume to work within the *infinite population* limit.

The time evolution may take place in discrete time; in this case the derivative corresponds to $`\dot{p}(t)p(t+1)p(t)`$. However, in our work the continuous–time description will be used, $`\dot{p}(t)\frac{p(t)}{t}`$, unless explicitly stated the opposite.

The dynamical matrix can be written in the form $`T_{ij}=Q_{ij}+\delta _{ij}(A_iD)`$, where $`A_j`$ is proportional to the fitness of the genome $`i`$ and so it represents a reproductive potential of the quasispecies. The mutation matrix $`Q_{ij}`$ represents the rate of mutation from the point $`i`$ in the genome space to another point $`j`$. The diagonal elements of the matrix $`Q_{ii}`$ are the probability rates of the perfect replication of the individual $`i`$. We may also introduce the homogeneous external stress which conserves the population size, through the decay rate $`D`$.

The dynamical matrix $`T_{ij}`$ is supposed to be symmetric. This allows us to apply the common machinery of Green functions. Indeed, the main object of interest will be the resolvent (Green function) of the dynamical matrix

$$𝒢(\zeta )=\frac{1}{\zeta T}$$
(2)
This formalism is particularly suitable for the case of the sharply–peaked landscape. Indeed, we suppose that the dynamical matrix can be decomposed as $`T=T^0+V`$, where $`T^0`$ corresponds to the flat landscape of neutral evolution and the perturbation $`V`$ describes the presence of the peak localised at point $`i=1`$,

$$V_{ij}=\alpha \delta _{ij}\delta _{i1}$$
(3)
Indeed, supposing that we are able to fully solve the neutral evolution, i. e. to find the “free” resolvent $`𝒢^0=(\zeta T^0)^1`$, the error threshold is closely related to the existence of the split-off pole in the resolvent $`𝒢(\zeta )`$. For the position of the pole $`z`$ we have the equation (Koster–Slater condition)

$$\frac{1}{\alpha }=G(z)$$
(4)
where $`G(\zeta )=𝒢_{11}^0(\zeta )`$ is the diagonal matrix element of the free resolvent at the site of the peak. If the equation (4) has a real solution, there is a split-off pole, the stationary state is localised and we are in the adaptive phase of the evolution. On the contrary, the absence of a real solution indicates that the stationary state is extended and we are in the wandering regime.

Our strategy in the following will consist in calculating the free resolvent for various models of the genome space and then analysing the solution of the equation (4). To start with, we reexamine the previously obtained results on the hypercube.

## 3 Exact solution on the hypercube

In the article  Galluccio et al. investigated the evolution of the quasispecies in the sharply–peaked landscape on the hypercubic lattice. The model is exactly soluble. Let us start by re-derivation of their main results.

The following assumptions are made in :

1. The topology of the genome space is the hypercube.
2. The evolution takes place in discrete time, i.e. $`\dot{p}p(t+1)p(t)`$.
3. At one time step, only single point mutation occurring with probability $`\mu `$ is allowed. Therefore, the mutation matrix connects only the nearest neighbours along the edges of the hypercube.
4. The fitness landscape is flat with only a single point which takes the additive selective advantage $`\alpha >0`$ (sharply–peaked landscape).

The consequence of allowing single point mutations in discrete time evolution, is that they had to restrict a mean number of mutations per genome — genomic mutation rate $`\overline{\mu }=\mu d`$ — to be less than 1. This need not to be valid in general, for example in natural viral populations, as is shown in paper  where the genomic mutation rates are estimated for measles virus, poliovirus, VSV and rhinovirus. The values are in the range from 0.13 to 1.15 mutations per replication. In other work  the genomic mutation rates were estimated in the range from 0.475 to 4.28 mutations per replication. This is the reason why we decided to use a continuous time approximation. The value $`\mu `$ can be easily found if we know the length of the life cycle and the mutation rate per site per replication. We took the rate from mentioned articles. The “lifetime” of measles virus or more precisely the time required for production of a new burst of viral particles is in the order of hours or days. So, for example, the mutation rate in RNA virus is around $`\mu =10^4`$ per site per generation . With the considered generation time $`20`$ hours we get approximately $`\mu =1.5\times 10^9`$ per site per second. The presented models of quasispecies have to respect the values of $`\mu `$ and the length of the sequence which is in order of kilobases. For instance, $`d=15,894`$ bases for the measles virus . Therefore, we choose the continuous-time description of the evolution process instead. The rest of the assumptions listed above remain valid. Let us proceed with formalisation of the above assumptions.

The genome space is the $`d`$-dimensional hypercube, $`𝒮=\{0,1\}^d`$. The points in the genome space are the sequences of base pairs $`i=[i_1,i_2,\mathrm{},i_d]𝒮`$. Let us denote $`\overline{i}(k)𝒮`$ the sequence which differs from $`i`$ by the point mutation at $`k`$-th base pair, $`(\overline{i}(k))_l=i_l+(12i_l)\delta _{kl}`$. Then the dynamic matrix of the neutral evolution corresponds to the diffusion in the hypercube

$$T_{ij}^0=\mu \underset{k=1}{\overset{d}{}}\delta _{j\overline{i}(k)}\mu d\delta _{ij}$$
(5)
The sharp peak located at the wildtype sequence $`0[0,\mathrm{}0]𝒮`$ contributes by the term expressed by Eq. (3).

The general condition (4) assumes very simple form, when we use the Fourier transform in the space $`𝒮`$ and divide the space with respect to the Hamming distance $`h`$ from the wildtype. The number of sequences in each Hamming distance class is given by the binomial coefficient $`\left(\genfrac{}{}{0pt}{}{d}{h}\right)`$. We obtain

$$\frac{1}{\alpha }=\frac{1}{2^d}\underset{h=0}{\overset{d}{}}\left(\genfrac{}{}{0pt}{}{d}{h}\right)\frac{1}{z+2\mu h}.$$
(6)
In the limit $`2^d\mathrm{}`$ the calculation further simplifies. The point is that the binomial distribution of classes reduces in that limit to a very thin and high band which we approximate by a $`\delta `$-function in the mean value, $`\left(\genfrac{}{}{0pt}{}{d}{h}\right)2^d\delta (hd/2)`$.

However, we should be aware that this approximation holds only for $`z`$ not too close to $`0`$. On the other hand, the error threshold itself is characterised by the fact that the solution $`z`$ of the equation (6) approaches to $`0`$. Therefore, we have to consider in (6) separately the term with $`h=0`$, while the rest of the sum, for $`h>0`$, is approximated by the $`\delta `$-function as hinted above. Hence, the approximate form of the Eq. (6) is111we have done numerical solutions of (6) and (7) for $`d=44`$ and even for such a short sequence the approximative formula fits very well. So, we expect very good agreement for the genomes with the length comparable to e. g. the measles virus’s genome.

$$\frac{1}{\alpha }=\frac{1}{2^d}\frac{1}{z}+\frac{2^d1}{2^d}\frac{1}{z+\mu d}.$$
(7)
The approximation we use here goes along different line than the one used in  but it gives the same results in lowest order.

The solution $`z`$ of Eq. (7) is

$$z=\frac{1}{2}\left(\alpha \mu d+\sqrt{(\alpha \mu d)^2+\frac{4\alpha \mu d}{2^d}}\right).$$
(8)
Now we perform the limit $`d\mathrm{}`$. It requires some care, as the product of $`\mu d`$ must stay finite, otherwise diagonal elements of the transfer matrix diverge. So we fix $`\overline{\mu }=\mu d`$ constant. Performing the limit we have to distinguish two different cases,

$`z=`$ $`{\displaystyle \frac{1}{2^d}}{\displaystyle \frac{2\alpha \overline{\mu }}{\overline{\mu }\alpha }}`$ $`\mathrm{for}(\alpha \overline{\mu })<0`$
$`z=`$ $`\alpha \overline{\mu }`$ $`\mathrm{for}(\alpha \overline{\mu })>0.`$ (9)
As we see behaviour of the largest eigenvalue $`z`$ changes sharply in the point $`\alpha _c=\overline{\mu }`$. Here we recognise the *error threshold*.

Our next aim is to calculate the *susceptibility* of the quasispecies to a change of the sequence length. Elongation of the sequence increases a number of possible sequences $`N=2^d`$. The susceptibility, then, shows us how the largest eigenvalue and thus a relative fitness of the wildtype changes with the “volume” of the sequence space.

We differentiate (8) with the respect to $`N`$. Then the limit $`N\mathrm{}`$ is done. Defining

$$\chi =\underset{N\mathrm{}}{lim}N^2\frac{z}{N}$$
(10)
we get in the vicinity of the error threshold

$$\chi =\frac{\alpha \alpha _c}{|\alpha \alpha _c|}$$
(11)
Now we introduce critical exponents $`\beta ^\pm `$, $`\gamma ^\pm `$. They characterise behaviour of the largest eigenvalue $`z`$ and the susceptibility, respectively, close to the error threshold, from the left $`()`$ and right $`(+)`$ side. In the vicinity of the error threshold we assume the critical behaviour

$`\alpha <\alpha _c`$ $`zN^1|\alpha \alpha _c|^\beta ^{}`$ $`\chi |\alpha \alpha _c|^\gamma ^{},`$
$`\alpha >\alpha _c`$ $`z|\alpha \alpha _c|^{\beta ^+}`$ $`\chi |\alpha \alpha _c|^{\gamma ^+}.`$ (12)
From (9) we see that $`z`$ behaves close to the error threshold as a power function with the exponents $`\beta ^{}=1`$ and $`\beta ^+=1`$, and the Eq. (11) implies that $`\gamma ^{}=\gamma ^+=1.`$

## 4 The Quasispecies Evolution in an Ultrametric Space

The sharply–peaked landscape is certainly an oversimplified approximation. It describes the adaptation close to a single selected peak within the true fitness landscape. Indeed, in reality the landscape is rugged and consists of very many concurrent peaks. In order to investigate their combined effect, we must turn to some model scheme.

It was realised quite a long ago that complicated free–energy landscapes with many minima are characteristic for disordered spin systems, namely spin glasses . This analogy was brought to biology within the spin–glass model of evolution . The main result of the spin–glass theory, which will be relevant for us here, is the tree–like (or ultrametric) structure of minima in the free–energy landscape . We therefore assume that the same ultrametric structure of peaks will be characteristic for the rugged fitness landscape. Then, we will observe the evolution on a “mesoscopic” time–scale, where a single event will be the jump from one peak to the other. Clearly, it consists of many “microscopic” events, i. e. point mutations in the DNA sequence.

The evolution takes place in the space $`𝒮`$ constructed as the set of endpoints of a regular three with $`K`$ levels, where at each level every branch splits into $`p`$ new branches. The ultrametric tree is illustrated in the left panel of Fig. 2.

From the mathematical point of view, the ultrametric structure is described through $`p`$-adic numbers and we will use the $`p`$-adic Fourier transform as the main tool for the solution. The ultrametric space is represented by the set $`𝒮=\{0,1,2,\mathrm{},p^K1\}`$ endowed with ultrametric distance between each pair of points, denoted $`|ij|_p`$ for $`i,j𝒮`$. We refer the reader to the Appendix for necessary definitions and formulae. The diffusion in the ultrametric space was already thoroughly studied before in the context of spin-glass dynamics  and the concept of $`p`$-adic Fourier transform was developed quite long ago in the mathematic literature , and was used implicitly (see e. g. ) but its use in physics was made explicit only recently .

The ultrametric structure will be reflected by the fact that the probability of jumps (mutations) from one point to the other will depend only on the ultrametric distance between the two points. It is analogous to the hypercube, where the mutation probability depends only on the Hamming distance.

Movement of the individual on the ultrametric lattice is driven by the dynamic matrix $`T`$ which is similar to the one used on the hypercubic lattice. There is one important difference — in the previous case jumps only to the nearest neighbours were considered, but now we must allow a possibility to jump further, otherwise we get evolution only in one isolated branch of the tree.

The evolution in ultrametric analog of the sharply–peaked landscape is described again by the equation (1) where now the “unperturbed” matrix $`T^0`$ reflects the symmetry of the ultrametric space. This means that the value of the matrix element $`T_{ij}^0`$ depends only on the ultrametric distance between the points $`i`$ and $`j`$. The “perturbation” $`V`$ corresponds to the single site with selective advantage. We can write explicitly

$$\begin{array}{ccc}T_{ij}^0=\hfill & \stackrel{~}{T}\hfill & \text{if}i=j\hfill \\ T_{ij}^0=\hfill & T_m\hfill & \text{if}|ij|_p=p^m,\hfill \\ V_{ij}=\hfill & \alpha \delta _{ij}\delta _{i0}\hfill & .\hfill \end{array}$$
(13)
Moreover, we will require that the bare transition matrix conserves probability. This fixes the value of the diagonal element,

$$\stackrel{~}{T}=\underset{b(a)}{}T_{ab}=\underset{m=0}{\overset{K1}{}}(p1)p^mT_m$$
(14)
This requirement is not substantial for our treatment, but simplifies somewhat the formulae. Especially it has the consequence that the matrix $`T^0`$ has an uniform eigenvector $`(1,1,\mathrm{},1)`$ with corresponding eigenvalue equal to 0.

Thus, the fitness landscape is fully described by the sequence of $`K`$ numbers $`T_0,T_1,\mathrm{},T_{K1}`$ and the parameter $`\alpha `$. As $`T_m^0`$ is the probability of hopping to the distance $`p^m`$, it is natural to require that

$$T_{m+1}^0>T_m^0$$
(15)
To follow the general scheme of section 2 we should first find the unperturbed resolvent $`𝒢^0(\zeta )=(\zeta T^0)^1`$. So, we need to diagonalise the matrix $`T^0`$.

To this end we use the $`p`$-adic Fourier transform. For details we refer the reader to the Appendix. The basic result is, that the eigenvalues of the matrix $`T^0`$ are $`\widehat{T}_m`$, $`m=0,1,\mathrm{},K`$, where

$$\begin{array}{cc}\hfill \widehat{T}_0=& \stackrel{~}{T}+\underset{l=0}{\overset{K1}{}}T_lp^{Kl1}(p1)\hfill \\ \hfill \widehat{T}_m=& \stackrel{~}{T}+\underset{l=m}{\overset{K1}{}}T_lp^{Kl1}(p1)T_{m1}p^{Km},m=1,\mathrm{},K.\hfill \end{array}$$
(16)
From the conditions (14) and (15) we can immediately deduce that $`\widehat{T}_0=0`$ is the maximum eigenvalue, all other eigenvalues are negative and $`\widehat{T}_{m+1}^0<\widehat{T}_m^0`$. Using the inverse $`p`$-adic Fourier transform we get for the unperturbed resolvent the expression

$$G(z)=\frac{1}{p^K}\left[\frac{1}{z}+\underset{l=1}{\overset{K}{}}\frac{p^{l1}(p1)}{z\widehat{T}_l}\right].$$
(17)
Introducing the normalised density of states (DOS) as

$$𝒟(\widehat{T})=\frac{1}{N}(\delta (\widehat{T})+\underset{j=1}{\overset{K}{}}p^{j1}(p1)\delta (\widehat{T}\widehat{T}_j))$$
(18)
we can write the Koster-Slater condition (4) in general form

$$\frac{1}{\alpha }=\frac{𝒟(\widehat{T})}{z\widehat{T}}.$$
(19)
Clearly, the factor $`p^{j1}(p1)`$ is the multiplicity of the eigenvalue $`T_j`$.

Let us come back to biological aspects. Each possible quasispecies has got a certain number of alleles, (or sets of subsequences) but we can not consider that every combination of those alleles forms a quasispecies which is able to survive. In fact, the ultrametric tree need not be regular but bears some level of randomness. The ultrametric space considered so far was regular, characterised by fixed branching ration $`p`$. We shall take into account the randomness by the following approximate procedure.

Let us take an arbitrary tree, maybe random. With the number of distinct end-points $`N=p^K`$ kept fixed we enlarge the number of levels, $`K\mathrm{}`$. We can freely do so, as the definition of levels is arbitrary: if no branches appear at certain level, we simply have $`p=1`$ locally. With this definition, $`p`$ may assume different values at different levels as well as at different branches within the same level. The randomness of the tree is then translated into the set of random $`p`$-s. We illustrate the situation in the Fig. 2. The approximation consists in characterising the tree by an average value of $`p`$, which now may be any real number larger than $`1`$. Furthermore, we will assume that the properties of trees with non-integer values of $`p`$ can be approximately calculated by analytical continuation from the results obtained for regular, non-random trees and arbitrary integer $`p`$. Moreover, we should keep in mind that with $`N`$ fixed and diverging number of levels $`K\mathrm{}`$ we should make the limit $`p1^+`$. At the end of the calculations, we will make also the thermodynamic limit $`N\mathrm{}`$.

Therefore, we take the formulae obtained before for general $`p`$ and suppose the variable $`p`$ is now a real number. Then, we will perform the limits

$$\begin{array}{cc}\hfill \text{first:}& K\mathrm{},p1^+\text{with}N=p^K\text{fixed}\hfill \\ \hfill \text{next:}& N\mathrm{}\hfill \end{array}$$
(20)
Now we should specify the form of the jump rates $`T_m`$. The rate $`T_m`$ of a certain individual to mutate to some other with the sequence which is in the ultrametric distance $`p^m`$ from the original one will be considered in the form

$$T_m=\mu p^{(mK)\tau }.$$
(21)
We solve this case for the value of $`2>\tau >1`$. We will see later that this is indeed the interval relevant for the error-threshold phenomenon.

The choice of an arbitrary power $`\tau `$ allows simulation of a quite general situation. In the limit $`p1`$ possible mutation rates to distinct sites are in a very wide range, starting in the value $`\mu p^{(K1)\tau }`$ and ending with the highest possible jump rate $`\mu `$. And this is exactly, in our opinion, what we find in reality — a probability rate of the jump from one codon sequence (or a certain set of alleles) to other is not uniform.

In order to solve the equation (19) we need to compute the density of states (18). The eigenvalues of $`T`$ can be easily found using the Fourier transform. A value of the diagonal element $`\stackrel{~}{T}`$ is again established with the help of the condition that the sum of elements of the matrix’s column equals zero, i. e. $`\widehat{T}_0=0`$.

$$\widehat{T}_m=\frac{\mu }{p}\frac{p1}{p^{\tau 1}1}p^{(\tau 1)K}\left(p^{(\tau 1)m}\left(1+\frac{p^{\tau 1}1}{p^{\tau 1}(p1)}\right)1\right)$$
(22)
The density of states $`𝒟(\widehat{T})`$ will be calculated in the limit (20). Its support can be found easily by computing the lowest eigenvalue $`\widehat{T}_K`$, which in the limit (20) is equal to $`\widehat{T}_{\mathrm{}}=\frac{\mu \tau }{\tau 1}`$. Substituting $`x=p^m`$ into (22) and inverting the dependence we obtain the function $`x(\widehat{T})`$, hence

$$𝒟(\widehat{T})=\frac{1}{N}\left|\frac{\mathrm{d}x}{\mathrm{d}\widehat{T}}\right|=\frac{1}{\tau 1}\left(\frac{\tau 1}{\mu \tau }\right)^{\frac{1}{\tau 1}}\left(\widehat{T}\right)^{\frac{2\tau }{\tau 1}}$$
(23)
for $`\widehat{T}(\frac{\mu \tau }{\tau 1},0)`$ and $`𝒟(\widehat{T})=0`$ outside this interval.

The equation for the split-off eigenvalue $`z`$ now becomes

$$\frac{(\tau 1)\overline{\mu }^{\frac{1}{\tau 1}}}{\alpha }=_0^{\overline{\mu }}\frac{y^\epsilon }{z+y}dy$$
(24)
where we denoted

$$\epsilon =\frac{2\tau }{\tau 1}.$$
(25)
and

$$\overline{\mu }=\frac{\mu \tau }{\tau 1}$$
(26)
The integral on the RHS of (24) is, for non-integer $`\epsilon `$,

$$_0^{\overline{\mu }}\frac{y^\epsilon }{z+y}dy=\frac{1}{\epsilon }\overline{\mu }^\epsilon \frac{\pi }{\mathrm{sin}\epsilon \pi }z^\epsilon +z\overline{\mu }^{\epsilon 1}\underset{k=0}{\overset{\mathrm{}}{}}\frac{1}{k+1\epsilon }\left(\frac{z}{\overline{\mu }}\right)^k.$$
(27)
(For integer values of $`\epsilon `$ the infinite series becomes a finite sum accompanied by a logarithmic correction). When analysing the equation (24) we must distinguish two regimes. For $`0<\epsilon <1`$ the leading contribution in (27) is of order $`z^\epsilon `$, while for $`\epsilon 1`$ the leading contribution is linear in $`z`$. The former regime applies for $`2>\tau >\frac{3}{2}`$ while the latter holds for $`1<\tau <\frac{3}{2}`$. This gives for the critical exponent at the error threshold

$$\beta ^+=\mathrm{max}(1,\frac{\tau 1}{2\tau }).$$
(28)
The location of the error threshold is determined by the first term on the RHS of the Eq. 27

$$\alpha _c=\mu \tau \frac{2\tau }{\tau 1}$$
(29)
We can see that the error threshold vanishes in the limit $`\tau 2`$ and diverges for $`\tau 1`$ which justifies the choice of $`\tau `$ within the interval $`(1,2)`$. We investigated also the case $`\tau =1`$, where we send $`\mu 0`$ simultaneously with $`N\mathrm{}`$ with $`\stackrel{~}{\mu }=\mu \mathrm{ln}N`$ kept fixed. In this case we also observe the error threshold for certain value of $`\stackrel{~}{\mu }`$, with exponent $`\beta ^+=1`$.

## 5 Evolution on Bethe lattice

Another approach which takes into account the complicated structure of real fitness landscapes uses the so-called holey landscape. This concept reflects the fact that most of the random point mutations which may occur are lethal and the evolution ends in a dead end there. So, we may start with the hypercubic lattice $`\{0,1\}^d`$ as in section 3 but leave only viable sites, while all lethal genomes are deleted. The remaining lattice has highly random geometry and the symmetry properties of the hypercube, which enabled us to use the Fourier transform, are no more applicable here. Thus, another scheme of the solution should be looked for.

In fact, topologically disordered lattices are widely studied in solid-state theory. The Bethe lattice is a well-known model of random lattices. It has got that large advantage that many problems can be solved analytically. It is widely used in statistical physics, see .

We take the infinite Bethe lattice with the connectivity $`k`$. In each site there is a population of individuals. With the uniform probability rate $`\mu `$ an individual may jump (mutate) to one of $`k`$ nearest neighbours’ site. Again one site with a selective advantage $`\alpha `$ is present and we will study creation of the quasispecies around that site. The unperturbed dynamical matrix $`T^0`$ corresponds to the topology of the Bethe lattice. The off-diagonal terms will be $`T_{ij}^0=\mu `$ if $`j`$ and $`i`$ are neighbours on the lattice and $`T_{ij}^0=0`$ elsewhere. The diagonal terms on a lattice with connectivity $`k`$ are obviously $`T_{ii}^0=k\mu `$ for all $`i`$.

Following the general scheme of Sec. 2 we calculate the diagonal element of the resolvent $`𝒢^0=(\zeta T^0)^1`$ at the selected site. As in the infinite Bethe lattice all sites are equivalent, we can calculate the diagonal element at any site inside the Bethe lattice. This will be done using the partitioning (projector) method.

The procedure is based in splitting the Bethe lattice into separate parts. This will be done in two steps, as illustrated in Fig. 3. We calculate the diagonal element $`G(\zeta )=𝒢_{11}^0(\zeta )`$ at site $`1`$. First, we separate the central site $`1`$ from the rest of the lattice. Because there are no loops, the rest of the Bethe lattice now consists of $`k`$ disconnected and identical infinite branches. The branches start at *terminal* sites denoted by $`2`$.

Let $`P`$ be the projector to the site $`1`$ and $`Q=1P`$ the complementary projector, which separates the rest of the lattice. Then, we get for the diagonal element of the resolvent

$$G(\zeta )=P𝒢(\zeta )P=\frac{P}{\zeta PT^0P\mathrm{\Sigma }(\zeta )},$$
(30)
where the self-energy has the form

$$\mathrm{\Sigma }(\zeta )=PT^0Q\frac{Q}{\zeta T^0}QT^0P.$$
(31)
By extracting the site $`1`$ the lattice consists of $`k`$ identical parts. So, $`\mathrm{\Sigma }(\zeta )`$ is a sum of $`k`$ identical contribution, each of them coming from one isolated branch. We can write

$$\mathrm{\Sigma }(\zeta )=\mu ^2k\mathrm{\Gamma }(\zeta ),$$
(32)
where $`\mathrm{\Gamma }(\zeta )`$ is the matrix element of the projected resolvent $`(\zeta QT^0Q)^1`$ at any of the terminal sites denoted by $`2`$ in Fig. 3.

The next step of the partitioning takes one branch and separates the terminal point $`2`$ from it. As a result, $`k1`$ disconnected branches arise. Their terminal points are denoted as $`3`$ in Fig. 3. An equation analogous to (30) is found also for $`\mathrm{\Gamma }(\zeta )`$.

We may iterate the same procedure arbitrary number of times. But if we consider infinite Bethe lattice, the branch starting with site $`2`$ must be identical to any of the branches starting at site $`3`$. This leads to closure in the relation for the quantity $`\mathrm{\Gamma }(\zeta )`$.

$$\mathrm{\Gamma }(\zeta )=\frac{1}{\zeta +\mu k\mu ^2(k1)\mathrm{\Gamma }(\zeta )}.$$
(33)
Solving the equations (33) and (30) we find the function $`G(\zeta )`$. From the two roots of the quadratic equation the proper one is chosen by the requirement that the resolvent must have correct asymptotic behaviour $`G(\zeta )1/\zeta `$ for $`|\zeta |\mathrm{}`$.

However, one complication arises here. If we calculate the density of states $`𝒟(\zeta )=\frac{1}{\pi }lim_{\epsilon 0}\text{Im}G(\zeta \text{i}\epsilon )`$, we observe for $`k>2`$ a gap arising between the upper edge of the density of states and the point $`\zeta =0`$. This reflects the well-known pathology of the Bethe lattice, which stems from the fact that the surface points of the Bethe lattice constitute finite fraction of the whole even in the thermodynamic limit, unlike lattices embedded in $`d`$-dimensional Euclidean space, where the fraction of the surface vanishes as $`N^{1/d}`$ when the number of sites $`N`$ goes to infinity. When investigating the diffusion on the Bethe lattice, it has the unnatural consequence that the probability density flows out constantly toward the surface and tends to zero in any finite portion of the lattice. In order to compensate for this unnatural effect, we put uniform probability flux into the lattice. In the context of biological evolution, this amounts simply to addition of a constant to the fitness of all genomes. Such an uniform change dos not affect the natural selection and evolutionary competition, but simplifies the mathematical treatment.

Indeed, introducing the uniform probability influx amounts simply to shifting the variable $`\zeta \zeta +\mu k2\mu \sqrt{k1}`$ so that the upper edge of the density of states is located at $`\zeta =0`$. The shifted resolvent is then

$$G(\zeta )=\frac{2(k1)}{(k2)(\zeta +2\mu \sqrt{k1})+k\sqrt{\zeta ^2+4\mu \zeta \sqrt{k1}}}.$$
(34)
The real and imaginary part of the shifted resolvent are shown in Fig. 4.

Now we can turn to the solution of (4) using the expression (34). Since we are looking for the real eigenvalue $`z0`$ we obtain

$$z=\frac{1}{2}\left[(2k)\alpha 4\sqrt{k1}\mu +k\sqrt{\alpha ^2+4\mu ^2}\right].$$
(35)
This solution is valid only above the error threshold: $`\alpha \alpha _c=\mu (k2)/\sqrt{k1}`$. For smaller $`\alpha `$ the equation (4) have no real solution, there is no localised state and the largest eigenvalue of the dynamic matrix is $`z=0`$.

So we have again found that there is for $`k>2`$ the error threshold between the random wandering regime and the adaptive regime at some positive value $`\alpha =\alpha _c`$. For $`k=2`$ the error threshold disappears, $`\alpha _c=0`$. However, this case is of little interest for us now, because the Bethe lattice reduces to a linear chain, hardly being a good approximation for the holey fitness landscape we started with.

Close to the the critical point, we obtain that for $`k>2`$ the largest eigenvalue $`z`$ behaves as a second power of the distance from the error threshold

$$z\frac{\sqrt{k1}}{4\mu }(\alpha \alpha _c)^2,\mathrm{when}\alpha \alpha _c^+.$$
(36)
So we have got the critical exponent $`\beta _b^+=2`$. This reflects the behaviour of the density of states at the band edge, which is $`𝒟(\zeta )|\zeta |^{1/2}`$ for any $`k>2`$. In fact, we obtained the relation between the exponent $`\beta ^+`$ and the behaviour of the density of states at the band edge already in Sec. 4 when we investigated the evolution in ultrametric space.

## 6 Super–Universality of the Critical Exponent $`\gamma ^+`$

In all previous cases the critical exponent $`\beta ^+>0`$ has been found. In this section we show that if we consider power behaviour of the largest eigenvalue $`z`$ near the critical point $`\alpha _c`$, we can find the value of the critical exponent $`\gamma ^+`$.

Assume that the whole sequence space contains $`N`$ distinct points. Then, there are exactly $`N`$ eigenvalues of the matrix $`T`$. It is known that the largest eigenvalue $`0`$ of the unperturbed dynamic matrix $`T^0`$ is non-degenerate. If $`N`$ is very large, we can approximately write for the matrix element of the resolvent $`G(\zeta )\overline{G}(\zeta )`$ corresponding to the diffusion on the lattice with $`N`$ sites that

$$\overline{G}(\zeta )=\frac{N1}{N}G(\zeta )+\frac{1}{N}\frac{1}{\zeta }.$$
(37)
We want to know what is the dependence of the split-off pole $`z`$ on the genome size $`N`$. It is described by the susceptibility defined in (10).

The Koster-Slater condition written as

$$\overline{G}(z)=\frac{1}{\alpha },$$
(38)
which corresponds to (4), can be regarded as an implicit function $`z(N)`$, so we use the formula for the derivative of the implicit function in order to get the susceptibility $`\chi `$ defined in (10). In the thermodynamic limit $`N\mathrm{}`$ the expression simplifies and we find that the susceptibility of the quasispecies to the change of the sequence length is

$$\chi =\frac{G(z)\frac{1}{z}}{\frac{G}{z}}.$$
(39)
As a next step, in order to find $`\gamma ^+`$, we suppose that in the vicinity of the error threshold the highest eigenvalue can be replaced by the term $`zq(\alpha \alpha _c)^{\beta ^+}`$, where $`q`$ is an appropriate coefficient. So, close to the error threshold

$$G(z)\frac{1}{\alpha _c}\frac{1}{\alpha _c^2}\left(\frac{z}{q}\right)^{1/\beta ^+}$$
(40)
hence

$$\chi =\frac{\beta \alpha _c^2}{\alpha \alpha _c},\text{for}\alpha \alpha _c$$
(41)
which corresponds exactly to (11). The value of the critical exponent $`\gamma ^+=1`$ follows, and we see that the exponent does not depend on the model in question. This proves the super-universality of the exponent $`\gamma ^+`$.

## 7 Conclusions

We investigated the transition between adaptive and neutral phase in biological evolution, which occurs at the error–threshold value of the mutation rate. We suggested to take into account the complicated structure of the fitness landscape by effectively choosing non-trivial topology of the genome space. In all cases we observed the formation of a quasispecies around a site with selective advantage.

We compared the evolutionary process in three different topologies. The first one was the hypercube. (We re-derived by a different method the already known results and completed the study by calculating the susceptibility, not studied before.) This corresponds to standard sharply–peaked landscape, with completely flat landscape as a background. Second, inspired by the theory of spin glasses, we modelled the ruggedness of the fitness landscape by taking as the basic space the set of peaks hierarchically organized in an ultrametric space. Third, we took the holey fitness landscape, in which all viable genomes have equal fitness and all lethal genomes are cut off. We approximated it by the Bethe lattice, as usual in various problems dealing with randomly connected networks in high dimensions.

The problem was treated as diffusion in the underlying space, the site with selective advantage acting as a source. The formation of a localized state (quasispecies) around the selected site was studied by observing the behavior of the maximum eigenvalue of the dynamical matrix governing the diffusion. If the maximum eigenvalue is split off from the rest of eigenvalues (the band), the state is localized and we are in the adaptive phase.

The split-off sets on at the critical value of the selective advantage, corresponding to the error threshold. We calculated the value of the critical value in all three topologies. We therefore confirmed that the very existence of the error threshold phenomenon does not depend on what is the structure of the genome space, at least what concerns the types of topology studied here.

As a measure of the critical behavior we chose the approach of the maximum eigenvalue to zero when the selective advantage approaches its critical value. We found the behavior be dominated by a power and calculated the critical exponent. We found that it is always larger or equal to $`1`$. For the hypercube it is $`1`$, for the Bethe lattice it has value $`2`$, while in the ultrametric space it may assume any value larger or equal to $`1`$, depending on details of the geometry.

In order to find the dependence between the largest eigenvalue and the volume of the genome space, we calculated the susceptibility of the quasispecies to the change of the sequence length. The susceptibility diverge in the vicinity of the error threshold. This could be of the partial interest for biologists, since natural viral population seems to have mutation rate very close to the critical point, see . We found that the divergence of the susceptibility at the error threshold is governed by a power-law with super-universal, model-independent exponent $`1`$.

All the topologies which we considered were regular. Our further outlook is to study the influence of the randomness in the topology of the genome space to the evolutionary process.

## Acknowledgement

Authors would like to thank to Jan Mašek for critical reading of the manuscript. Michal Kolář would like to thank to Vladislav Čápek for fruitful discussions and to Anton Markoš for helpful discussions on the topic of biological evolution, too.

## Appendix

We will present the basic facts on ultrametric spaces, $`p`$-adic numbers and $`p`$-adic Fourier transform here .

The ultrametric space is a special kind of a metric space. It is a set $`𝒰`$ of points endowed by a non-negative function expressing the distance between two points in the space. For reasons which will become clear later, in this article we will use the notation $`|ab|_p`$ for the distance between points $`a,b𝒰`$.

As usual in metric spaces, we require that the distance from a point to itself is zero and distance to any other point is positive. The most important property of ultrametric spaces is that instead of the triangle inequality a stronger condition holds,

$$|ab|_p\mathrm{max}\{|ac|_p,|cb|_p\}a,b,c𝒰.$$
(42)
This implies that any triangle is equilateral or isosceles with a small base, that every point inside a ball is its centre, or that the diameter equals to the radius. More detailed presentation can be found in .

There is a class of the ultrametric spaces formed by *p-adic* numbers. These are defined as sums

$$a=\underset{i=0}{\overset{K1}{}}a_ip^i,$$
(43)
where $`p`$ and $`a_i`$, $`0a_ip1`$, are natural integers. Clearly $`a\{0,1,\mathrm{},p^K1\}`$. The distance $`||_p`$ is defined as $`|ab|_p=p^j`$ where $`p^{+j}`$ is the largest divisor of the number $`|ab|`$ in that form ($`j`$ is an integer). Note that $`|a|_p`$ can be considered as an (ultrametric) norm of the number $`a`$.

The integral over all $`p`$-adic numbers is defined as

$$_pT(a)da=\frac{1}{p^K}\underset{a=0}{\overset{p^{K1}}{}}T(a),$$
(44)
where the sum goes over whole set of $`p`$-adic numbers. With so defined integral we have the Fourier transform of the function $`T(a)`$ in the form:

$$\widehat{T}(m)=\frac{1}{p^K}\underset{a=0}{\overset{p^K1}{}}e^{2\pi \mathrm{i}ma}T(a).$$
(45)
$`m`$ is an element of the reciprocal space and it is in the form $`m=_{i=1}^Km_ip^{iK}`$, $`0m_ip1`$. If the function $`T(a)`$ depends only on the ultrametric norm, i. e. $`T(a)=f(|a|_p)`$, and if $`|m|_p=p^k`$, we get (for $`k=1,\mathrm{},K`$):

$$\begin{array}{cc}\hfill \widehat{T}_k\widehat{T}(m)& =\stackrel{~}{T}+\underset{l=k}{\overset{K1}{}}T_lp^{Kl1}(p1)T_{k1}p^{Kk}\hfill \\ \hfill \widehat{T}_0\widehat{T}(0)& =\stackrel{~}{T}+\underset{l=0}{\overset{K1}{}}T_lp^{Kl1}(p1).\hfill \end{array}$$
(46)
In the previous formula we defined, $`\stackrel{~}{T}=T(0)`$, and $`T_l=T(a)`$ when $`|a|_p=p^l`$. On the other side for $`|a|_p=p^l`$ we get the inversion of the Fourier transform in the form (for $`l=0,\mathrm{},K1`$):

$$\begin{array}{cc}\hfill T_lT(a)& =p^K\widehat{T}_0+\underset{m=1}{\overset{l}{}}\widehat{T}_mp^{mK1}(p1)\widehat{T}_{l+1}p^{lK}\hfill \\ \hfill \stackrel{~}{T}T(0)& =p^K\widehat{T}_0+\underset{m=1}{\overset{K}{}}\widehat{T}_mp^{mK1}(p1).\hfill \end{array}$$
(47)


text:
# One-parameter nonrelativistic supersymmetry for microtubules

## References