text:
# 𝐽/𝜓 Radiative Decays

## 1 Introduction

One of the distinctive features of QCD as a non-Abelian gauge theory is the interaction of quarks and gluons, which predicts the existence of other types of hadrons with explicit gluonic degrees of freedom – glueballs and hybrids. The indirect evidence for gluon-gluon interactions has been obtained at high energies. However, glueballs, the bound states of gluons, and hybrids predicted by QCD, have not been confirmed yet. Therefore, the observation of glueballs and hybrids is, to some extent, a direct test of QCD, and the study of the glueball and hybrid spectroscopy, together with the meson and baryon spectroscopy will be a good laboratory for the study of the strong interactions in the strongly coupled non-perturbative regime.

Many QCD-based models and calculations, for example, bag models , flux-tube models , QCD sum rules  and lattice QCD  are developed to make predictions to the properties of glueballs and hybrids. Of them, lattice QCD is considered as the most relevant since it originated from the first principle of QCD.

According to the calculations of different lattice QCD groups, the lightest glueball is found to be a $`0^{++}`$ state with the mass in the region of 1.5-1.7 GeV, while the next lightest glueball is a $`2^{++}`$ with the mass around 2.3GeV.

$`J/\psi `$ decay is a good lab. for the study of hadron spectroscopy, as well as glueball and hybrid search. Fig.1 shows the Feymann diagrams of some $`J/\psi `$ decays. A naive estimation of the production rate of various particles based simply on counting powers of the electromagnetic and strong coupling constants yields,

$$\mathrm{\Gamma }(J/\psi \gamma G)O(\alpha \alpha _s^2),\mathrm{\Gamma }(J/\psi \gamma H)O(\alpha \alpha _s^3)$$

$$\mathrm{\Gamma }(J/\psi \gamma M)O(\alpha \alpha _s^4),\mathrm{\Gamma }(J/\psi \gamma F)O(\alpha \alpha _s^4)$$

where, G, H, M and F represent glueball, hybrid state, meson and four-quark state, respectively. Apparently, glueballs have enhanced production rates in $`J/\psi `$ radiative decays.

## 2 Previous results on $`J/\psi `$ radiative decays

Many experiments at the $`e^+e^{}`$ storage ring facilities, such as Crystal Ball, MARKIII, DM2 and BES have been dedicating to the study of the hadron spectroscopy and the search for non-$`q\overline{q}`$ states through $`J/\psi `$ radiative decays.

DASP, MARKII, MARKIII, DM2, Crystal Ball and BES studied the radiative production of $`q\overline{q}`$ mesons, such as tensors, scalars, pseudoscalars and axialvectors, and measured the resonance parameters, decay branching ratios and polarization parameters. In search for glueballs and new resonances, these experiments provided much information on $`f_0(1500)`$, $`f_0(1710)`$, $`\iota /\eta (1440)`$ and $`\xi (2230)`$.

## 3 Preliminary results from BESII $`J/\psi `$ data

BES, shown in Fig. 2, is a large general purpose solenoidal detector at the Beijing Electron Positron Collider (BEPC). The beam energy range is from 1.0 to 2.8 GeV and the luminosity at $`J/\psi `$ peak is around $`5\times 10^{30}cm^2s^1`$. The details of BESI are described in ref. .

The upgrades of BESI to BESII include the replacement of the central drift chamber with a vertex chamber composed of 12 tracking layers, the installation of a new barrel time-of-flight counter (BTOF) with a time resolution of 180ps and the installation of a new main drift chamber (MDC), which has 10 tracking layers and provides a $`dE/dx`$ resolution of $`\sigma _{dE/dx}=8.4\%`$ for particle identification and $`\sigma _p/p=1.7\%\sqrt{1+p^2}`$ ($`p`$ in GeV) momentum resolution for charged tracks. The barrel shower counter (BSC), which covers $`80\%`$ of $`4\pi `$ solid angle, has an energy resolution of $`\sigma _E/E=22\%/\sqrt{E}`$ ($`E`$ in GeV) and a spatial resolution of 7.9 mrad in $`\varphi `$ and 2.3 cm in z, is located outside the TOF. Outermost is a $`\mu `$ identification system, which consists of three double layers of proportional tubes interspersed in the iron flux return of the magnet.

With the upgraded BESII detector, till the summer of 2001, about $`5.8\times 10^7J/\psi `$ events have been accumulated. This is the largest $`J/\psi `$ data sample in the world. Fig. 3 shows the inclusive $`K_s^0`$ and $`\varphi `$ signals, which indicates a good data quality and a large data sample.

### 3.1 Partial Wave Analysis of $`J/\psi \gamma K\overline{K}`$ (preliminary)

The spin-parity of the structure around 1.7 GeV has a long history of uncertainty. Partial wave analyses (PWA) are applied to $`J/\psi \gamma K^+K^{}`$ and $`\gamma K_s^0K_s^0`$ channels, based on BESII $`5.8\times 10^7J/\psi `$ events. Fig. 4 shows the invariant mass spectra of $`K^+K^{}`$ and $`K_s^0K_s^0`$, where $`f_2^{}(1525)`$ and the structure around 1.7 GeV are clearly seen in both cases.

The amplitudes are fitted to the relativistic covariant tensor expressions, and the maximum likelihood method is applied in the fit. Global and slice fits are performed to $`KK`$ mass region of 1-2 GeV. For slice fit, 40MeV bin width is adopted. Fig. 5 shows $`0^{++}`$ and $`2^{++}`$ intensity distributions as a function of the invariant mass of $`K^+K^{}`$ and $`K_s^0K_s^0`$, where, dots with error bars are the efficiency corrected data points, solid curves stand for the coherent superposition of the individual Breit-Wigner resonance fits and dashed histograms represent the global fit results. Both $`K^+K^{}`$ and $`K_s^0K_s^0`$ D wave intensity shows a clear $`f_2^{}(1525)`$ with the mass and width being $`M=1518\pm 6`$ MeV and $`\mathrm{\Gamma }=84_{24}^{+28}`$ MeV. There is an evidence for $`f_2(1270)`$ and a weak $`2^{++}`$ intensity in the mass region around 1.7 GeV. While $`0^{++}`$ is found to be dominant in 1.7 GeV region, which is well determined by a Breit-Wigner resonance. The mass and width of the $`0^{++}`$ component are:

$$M=1703_{10}^{+8}MeV,\mathrm{\Gamma }=163_{22}^{+27}MeV$$

### 3.2 Partial Wave Analysis of $`J/\psi \gamma \pi ^+\pi ^{}`$ (preliminary)

Fig. 6 is the invariant mass spectrum of $`\pi ^+\pi ^{}`$ in $`J/\psi \gamma \pi ^+\pi ^{}`$.

Except for the well known $`f_2(1270)`$, a shoulder in the high mass region of $`f_2(1270)`$ and a bump near 1.7 GeV are seen. PWA is performed to $`\pi ^+\pi ^{}`$ mass in 1-2 GeV region. $`0^{++}`$ and $`2^{++}`$ intensity distributions are shown in Fig. 7, where, the same as in $`K\overline{K}`$ case, dots with error bars are the efficiency corrected data points, solid curves stand for the coherent superposition of the individual Breit-Wigner resonance fit and dashed histograms represent the global fit results.

A strong $`2^{++}`$ is well determined at around 1.27 GeV, and two $`0^{++}`$’s are observed at around 1.4 and 1.7 GeV. If the $`0^{++}`$ at around 1.7 GeV is considered as the same one as in $`K\overline{K}`$ and the mass and width are fixed to be the values obtained from $`K\overline{K}`$ analyses, i.e.

$$M=1703_{10}^{+8}MeV,\mathrm{\Gamma }=163_{22}^{+27}MeV,$$

the fit results give the masses and widths of $`2^{++}`$ and another $`0^{++}`$ as:

$$M=1266\pm 6MeV,\mathrm{\Gamma }=170_{30}^{+37}MeV$$

$$M=1383_{18}^{+20}MeV,\mathrm{\Gamma }=274_{128}^{+59}MeV$$

The preliminary results indicate a relative $`\pi \pi `$ to $`KK`$ branching ratio of around 30%.

### 3.3 Measurement of $`\eta _c`$ mass

A precise knowledge of the mass difference between $`J/\psi (1^{})`$ and $`\eta _c(0^+)`$ charmonium states is useful to the determination of the strength of spin-spin interaction term in non-relativistic potential model. The mass of $`J/\psi `$ has been measured with high accuracy, while the $`\eta _c`$ mass measurements listed on PDG2000 are very different and the fit to the measurements has a confidence level of only 0.001. Based on BESII $`5.8\times 10^7J/\psi `$ events, we measure the mass of $`\eta _c`$ from 6 channels: $`J/\psi \gamma \eta _c`$, $`\eta _cK^+K^{}\pi ^+\pi ^{}`$, $`\pi ^+\pi ^{}\pi ^+\pi ^{}`$, $`K^\pm K_S^0\pi ^{}`$, $`\varphi \varphi `$, $`K^+K^{}\pi ^0`$ and $`p\overline{p}`$. Fig. 8 shows $`\eta _c`$ signals in above 6 channels.

The fit values of the number of events and $`\eta _c`$ mass in the individual channels are listed in Table 1 (the errors are statistical only). In the fit, the width of $`\eta _c`$ is fixed at 16.5 MeV, which is the weighted average of PDG2000, BESI(2000) and recent CLEO two-photon collision experiment . Combine the weighted average with the results from the 5 channels, listed in Table 1, the mass of $`\eta _c`$ is given as $`M_{\eta _c}=2977.6\pm 0.8`$(stat.) MeV.

### 3.4 Near threshold structure in $`J/\psi \gamma p\overline{p}`$

There is an accumulation of evidence for anomalous behavior in the $`p\overline{p}`$ system near $`2m_p`$ mass threshold. A narrow dip-like structure at a center of mass energy of $`2m_pc^2`$ was observed in $`e^+e^{}hadrons`$. A narrow dip structure which is just above $`2m_p`$ was also observed in diffractive photoproduction of $`3\pi ^+3\pi ^{}`$ final states. Based on BESII $`J/\psi `$ data, we analyze $`J/\psi \gamma p\overline{p}`$. The invariant mass spectrum of $`p\overline{p}`$ is shown in Fig. 9. Except for $`\eta _c`$ signal, there is a clear enhancement near $`2m_p`$.

An S-wave Breit-Wigner function, as shown below, is used to fit the near threshold enhancement.

$$BW\frac{M_0\mathrm{\Gamma }_0(q/q_0)}{(M^2M_0^2)^2+(M_0\mathrm{\Gamma }_0(q/q_0))^2}$$

Where, $`q`$ is the momentum of daughter particle, $`q_0`$ is the momentum of daughter particle at peak. Fig. 10 shows the fit curve. The preliminary fit gives (statistical errors only):

$$M=1894.1\pm 1.4MeV,\mathrm{\Gamma }=48.0\pm 5.9MeV$$

## 4 Summary

$`J/\psi `$ radiative decay is a good laboratory in searches for glueball, the study of light hadron spectroscopy, as well as the glueball spectroscopy. In this talk, we mainly reported the PWA results from $`J/\psi \gamma K\overline{K}`$, $`\gamma \pi ^+\pi ^{}`$, the measurement of $`\eta _c`$ and the search for interesting states based on BESII $`5.8\times 10^7`$ $`J/\psi `$ data. More analyses are going on and more results will be produced very soon with this data set. CLEO is going to reduce the beam energy to collect a large amount of $`J/\psi `$ data in a few years. An experimental plan of upgrading BESII/BEPC to BESIII/BEPCII, with a luminosity of $`10^{33}`$ and a good detector is underway now. More results on the study of the light hadron spectroscopy and the search of the new form of hadronic states are expected in the near future.

## 5 Acknowledgements

We acknowledge the staff of the BEPC accelerator and IHEP computing center for their efforts. The work was supported in part by the National Natural Science Foundation of China under Contracts No. 19991480, No. 19825116 and No. 19605007, and by the Department of Energy of US under Contracts No. DE-FG03-93ER40788 (Colorado State University), No. DE-AC03-76SF00515 (SLAC), No. DE-FG03-94ER40833 (University of Hawaii) and No. DE-FG03-95ER40925 (University of Texas at Dallas).


format: text

text:
# Particle ratios at RHIC: Effective hadron masses and chemical freeze-out

## I Introduction

Thermodynamical equilibrium calculations of particle production in high energy particle- and nuclear collisions have been carried out for a long time ferm50 ; land53 ; west76 ; stoe78 ; cser86 ; hahn86 ; brau96 ; Let00 ; raf01 ; bec00 ; brau01 . Recently hadron abundances and particle ratios have been measured in heavy-ion collisions from SIS, AGS, SPS to RHIC energies. These data have revived the interest in the extraction of temperatures and chemical potentials from thermal equilibrium ”chemical” model analyses. The experimentally determined hadron ratios can be fitted well with straightforward noninteracting gas model calculations brau96 ; Let00 ; brau98 ; cley99 ; raf01 ; bec00 ; brau01 , if a sudden breakup of a thermalized source is assumed and once the subsequent feeding of the various channels by the strongly decaying resonances is taken into account. From the $`\chi ^2`$ freeze-out fits one has constructed a quite narrow band of freeze-out values in the $`T\mu _B`$ plane (see e.g.brau98 ; cley99 ). The extracted freeze-out parameters are fairly close to the phase transition curve for SPS and RHIC energies. However, when we are indeed so close to the phase transition or to a crossover as suggested by the data for $`T`$ and $`\mu _B`$, we can not afford to neglect the very in-medium effects we are after - and which, after all, do produce the phase transition. Thus, since noninteracting gas models neglect any kind of possible in-medium modifications they can not yield information about the phase transition.

Therefore, we will employ below a relativistic selfconsistent chiral model of hadrons and hadron matter developed in paper2 ; paper3 ; springer . This model can be used as a thermodynamically consistent effective theory or as a toy model, which embodies the restoration of chiral symmetry at high temperatures or densities. Therefore the model predicts temperature and density dependent hadronic masses and effective chemical potentials, which have already been proposed and considered in stoe78 ; thei83 ; scha91 ; springer ; BroRho ; mich01 . Thus, using the chiral $`SU(3)`$ model we can investigate, whether the freeze-out in fact takes place close to the phase transition boundary (if it exists) and if the extracted $`T,\mu _B`$ parameters are strongly model dependent. Depending on the chosen parameters and degrees of freedom different scenarios for the chiral phase change are predicted by the model: Strong or weak first order phase transition or a crossover. The transitions take place around $`T_c=155\text{ MeV}`$ zsch01 ; springer , which is in qualitative agreement with lattice predictions Karsch:2000kv  for the critical temperature for the onset of a deconfined phase which coincides with that of a chirally restored phase kar98 .

## II Model description

The chiral $`SU(3)`$ model is presented in detail in paper3 ; springer . We will briefly introduce the model here: We consider a relativistic field theoretical model of baryons and mesons built on chiral symmetry and broken scale invariance. The general form of the Lagrangean looks as follows:

$$=_{\mathrm{kin}}+\underset{W=X,Y,V,𝒜,u}{}_{\mathrm{BW}}+_{\mathrm{VP}}+_{\mathrm{vec}}+_0+_{\mathrm{SB}}.$$

$`_{\mathrm{kin}}`$ is the kinetic energy term, $`_{\mathrm{BW}}`$ includes the interaction terms of the different baryons with the various spin-0 and spin-1 mesons (see paper3  for details). The baryon masses are generated by both, the nonstrange $`\sigma `$ ($`<q\overline{q}>`$) and the strange $`\zeta `$ ($`<s\overline{s}>`$) scalar condensate. $`_{\mathrm{VP}}`$ contains the interaction terms of vector mesons with pseudoscalar mesons. $`_{\mathrm{vec}}`$ generates the masses of the spin-1 mesons through interactions with spin-0 fields, and $`_0`$ gives the meson-meson interaction terms which induce the spontaneous breaking of chiral symmetry. It also includes a scale-invariance breaking logarithmic potential. Finally, $`_{\mathrm{SB}}`$ introduces an explicit symmetry breaking of the U(1)A, the SU(3)V, and the chiral symmetry. All these terms have been discussed in detail in springer ; paper3 .
The hadronic matter properties at finite density and temperature are studied in the mean-field approximation serot97 . Then the Lagrangean (II) becomes

$`_{BX}+_{BV}`$ $`=`$ $`{\displaystyle \underset{i}{}}\overline{\psi _i}[g_{i\omega }\gamma _0\omega ^0+g_{i\varphi }\gamma _0\varphi ^0+m_i^{}]\psi _i`$ (1)
$`_{vec}`$ $`=`$ $`{\displaystyle \frac{1}{2}}m_\omega ^2{\displaystyle \frac{\chi ^2}{\chi _0^2}}\omega ^2+{\displaystyle \frac{1}{2}}m_\varphi ^2{\displaystyle \frac{\chi ^2}{\chi _0^2}}\varphi ^2+g_4^4(\omega ^4+2\varphi ^4)`$
$`𝒱_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}k_0\chi ^2(\sigma ^2+\zeta ^2)k_1(\sigma ^2+\zeta ^2)^2k_2({\displaystyle \frac{\sigma ^4}{2}}+\zeta ^4)k_3\chi \sigma ^2\zeta `$ (2)
$`+`$ $`k_4\chi ^4+{\displaystyle \frac{1}{4}}\chi ^4\mathrm{ln}{\displaystyle \frac{\chi ^4}{\chi _0^4}}{\displaystyle \frac{\delta }{3}}\chi ^4\mathrm{ln}{\displaystyle \frac{\sigma ^2\zeta }{\sigma _0^2\zeta _0}}`$ (3)
$`𝒱_{SB}`$ $`=`$ $`\left({\displaystyle \frac{\chi }{\chi _0}}\right)^2\left[m_\pi ^2f_\pi \sigma +(\sqrt{2}m_K^2f_K{\displaystyle \frac{1}{\sqrt{2}}}m_\pi ^2f_\pi )\zeta \right],`$
where $`m_i`$ is the effective mass of the hadron species $`i`$. $`\sigma `$ and $`\zeta `$ correspond to the scalar condensates, $`\omega `$ and $`\varphi `$ represent the non-strange and the strange vector field respectively, and $`\chi `$ is the scalar-isoscalar dilaton field, which mimics the effects of the gluon condensate sche80 . Only the scalar ($`_{BX}`$) and the vector meson terms ($`_{BV}`$) contribute to the baryon-meson interaction, since for all other mesons the expectation value vanishes in the mean-field approximation. The grand canonical potential $`\mathrm{\Omega }`$ per volume $`V`$ as a function of chemical potential $`\mu `$ and temperature $`T`$ can be written as:

$$\frac{\mathrm{\Omega }}{V}=_{vec}_0_{SB}𝒱_{vac}T\underset{i}{}\frac{\gamma _i}{(2\pi )^3}d^3k\left[\mathrm{ln}\left(1\pm e^{\frac{1}{T}[E_i^{}(k)\mu _i^{}]}\right)\right],$$
(4)
with the baryons (top sign) and mesons (bottom sign). The vacuum energy $`𝒱_{vac}`$ (the potential at $`\rho _B=0,T=0`$) has been subtracted in order to get a vanishing vacuum energy. $`\gamma _i`$ denote the hadronic spin-isospin degeneracy factors. The single particle energies are $`E_i^{}(k)=\sqrt{k_i^2+m_{i}^{}{}_{}{}^{2}}`$ and the effective chemical potentials read $`\mu _i^{}=\mu _ig_{i\omega }\omega g_{\varphi i}\varphi `$.
The mesonic fields are determined by extremizing $`\frac{\mathrm{\Omega }}{V}(\mu ,T=0)`$. The density of particle $`i`$ can be calculated by differentiating $`\mathrm{\Omega }`$ with respect to the corresponding chemical potential $`\mu _i`$. This yields:

$`\rho _i=\gamma _i{\displaystyle \frac{d^3k}{(2\pi )^3}\left[\frac{1}{\mathrm{exp}[(E_i^{}\mu _i^{})/T]\pm 1}\right]}.`$ (5)
All other thermodynamic quantities can also be obtained from the grand canonical potential. In the present calculation the lowest lying baryonic octet and decuplet and the lowest lying mesonic nonets are coupled to the relativistic mean fields. Depending on the coupling of the baryon resonances (the decuplet) to the field equations, the model shows a first order phase transition or a crossover (for details see zsch01 ). We will use three different parameter sets: Parameter set CI treats the members of the baryon decuplet as free particles, which yields a crossover behaviour. Parameter sets CII and CIII include also the (anti)-baryon decuplet as sources for the meson field equations. They differ by an additional explicit symmetry breaking for the baryon resonances along the hypercharge direction, as described in paper3  for the baryon octet. This is included in CII and not used in CIII. This leads to a weak first order phase transition at $`\mu =0`$ for CII and two first order phase transitions for CIII, which can be viewed as one strong first order phase transition. Heavier resonances up to $`m=2`$ GeV are always included as free particles. The resulting baryon masses for CI and CIII are shown in fig. 1. We observe a continous decrease of the baryon masses for CI starting at $`T150\text{ MeV}`$. In contrast, CIII shows two jumps around $`T=155\text{ MeV}`$.

The critical energy densities, the entropy densities and the transition temperatures for $`\mu _q=\mu _s=0`$ ($`\mu _q=\mu _B/3,\mu _s=\mu _B/3\mu _S`$) are specified in table 1.

## III Particle ratios in the chiral $`SU(3)\times SU(3)`$ model

Since the chiral SU(3) model predicts density and temperature dependent hadronic masses and effective potentials, in contrast to noninteracting models, the resulting particle ratios and therefore the deduced freeze-out temperatures and baryon chemical potentials are expected to change zsch00 . Hence in the following, we identify combinations of temperatures and chemical potentials that fit the observed particle ratios in the chiral model. In all calculations the value of the strange chemical potential $`\mu _S`$ is chosen such that the net strangeness $`f_s=0`$. We are looking for minima of $`\chi ^2`$ with

$$\chi ^2=\underset{i}{}\frac{\left(r_i^{exp}r_i^{model}\right)^2}{\sigma _i^2}.$$
(6)
Here $`r_i^{exp}`$ is the experimental ratio, $`r_i^{model}`$ is the ratio calculated in the model and $`\sigma _i`$ represents the error in the experimental data points. We use the same ratios as in brau01 : $`\overline{p}/p`$, $`\overline{\mathrm{\Lambda }}/\mathrm{\Lambda }`$, $`\overline{\mathrm{\Xi }}/\mathrm{\Xi }`$, $`\pi ^{}/\pi ^+`$, $`K^{}/K^+`$, $`K^{}/\pi ^{}`$, $`K_0^{}/h^{}`$, $`\overline{K_0^{}}/h^{}`$.

Even though the only parameters in a thermal and chemical equilibrium approach on first sight are the temperature and the baryon chemical potential, there exist further unknowns: On the one hand, some decays of high mass resonances are not well known and on the other hand the effect of weak decays in the experiments strongly depends on the detector geometry and on the reconstruction efficiency of the experiments. The feeding correction from the strong and electromagnetic decays of the hadronic resonances used here employs the procedure used in the UrQMD model bass98 ; urqmdhome . Weak decays are not considered here. We rather focus on the principal question whether an interacting chiral SU(3) approach with $`m^{}m_{vac}`$ can at all describe the particle yields at RHIC. Fine tuning of the $`\chi ^2`$ by adjustment of the weak decay scheme is not our intention. Even though it has been shown mich\_dr  that $`\chi ^2`$ values may be improved by including weak decays.

To compare the quality of the fits obtained in the chiral model with those obtained from the noninteracting gas approach, we set all masses and chemical potentials contained in the chiral model to their vacuum values and again use the same UrQMD feeding procedure as for the interacting model. This yields the ideal gas denoted $`ig_{FFM}`$. We find that the resulting ideal gas ratios are not identical but comparable to those obtained in the literature zsch00 ; mich\_dr ; brau01 ; flor01 . The differences should only result from a different treatment of weak interactions and from the uncertainty in the decay scheme of high mass resonances.

## IV Results for Au + Au Collisions at RHIC

First, we find that a reasonable fit of the measured particle ratios at RHIC is possible in all three phase transition scenarios of the chiral model and the ideal gas case with comparable quality.

Second, the resulting freeze-out values depend on the model employed, i.e. crossover, weak first order, strong first order or free thermal gas.

Third, a reasonable description of the data is impossible above $`T_c`$ in the models showing a first order phase transition. This shows that no direct freeze-out from the restored phase is observed.

Figure 2 shows the value of $`\chi ^2`$ in the $`T\mu _B`$ plane for the crossover case and for the strong first order phase transition. We see that the best fit $`T\mu _B`$ values differ in both models. Furthermore, in the crossover case $`\chi ^2`$ is well behaved as a function of $`T`$ and $`\mu _B`$. In contrast, the model with a strong first order phase transition shows a very steep increase of $`\chi ^2`$ at the phase transition boundary: the quality of the fit decreases drastically due to the jump of the effective masses at the phase transition boundary. Above $`T_c`$ the $`\chi ^2`$ values are inacceptable, $`\chi ^2>500`$.

The resulting best-fit particle ratios, $`\chi ^2`$-values and thermodynamic quantities are shown in table 2 and figure 3.

The $`\chi ^2`$ values for the chiral model are: $`\chi _{CI}^2=5.50`$, $`\chi _{CII}^2=5.73`$ and $`\chi _{CIII}^2=5.40`$. Thus, all three parameter sets describe the data equally well. Furthermore, the agreement is as good as in the noninteracting gas calculation ($`\chi _{ig}^2=5.72`$ brau01 , $`\chi _{FFM}^2=5.66`$). The best fit $`T\mu _B`$ parameters vary quite considerably between the different models. The noninteracting gas calculation yields $`T=187.6\text{ MeV}`$ and $`\mu _B=44.1\text{ MeV}`$. These freeze-out values can be compared to those obtained in other ideal gas calculations: $`T=174\text{ MeV},\mu _B=46\text{ MeV}`$ in brau01 , $`T=165\text{ MeV},\mu _B=41\text{ MeV}`$ in flor01  and $`T=190\text{ MeV},\mu _B=45\text{ MeV}`$ in xu01 . The crossover case in the interacting chiral model (CI) yields $`T=170.8\text{ MeV},\mu _B=48.3\text{ MeV}`$. Very strong deviations are found for the models with a first order phase transition (CII,CIII): The freeze-out temperatures are $`T=155`$ MeV (CII) and $`T=153.3`$ MeV (CIII), more than $`30`$ MeV lower than for $`ig_{FFM}`$. The fitted baryon chemical potentials $`\mu _B`$ increase by about $`710\text{ MeV}`$. These $`T\mu _B`$ pairs are very close to the phase boundary (CII) or even right on it (CIII) and are about $`10`$ MeV higher than the values obtained at SPS-energies zsch02 . Mainly due to the different freeze-out temperatures the values of the corresponding thermodynamic quantities vary between the different approaches. However, the energy per particle $`E/A`$ is approximately $`1`$ GeV in all cases. This ’unified freeze-out condition’ has already been proposed in cley98 .

The fact that the freeze-out appears right at the phase boundary or at crossover implies that there are large in-medium corrections, in particular for the effective masses, a phenomenon observed already in thei83 . The effective masses shown in figure 4 are shifted up to $`15\%`$ from their vacuum values. However, all the interacting models show similar values for the effective mass of a given hadron. The strongest in-medium modifications are observed for the nonstrange baryons ($`\mathrm{\Delta }m_i^{}/m_i15\%`$). Mesons and strange baryons show smaller changes of the effective masses, e.g. about $`10\%`$ for $`\mathrm{\Lambda },\pi ,K^{}`$, about $`5\%`$ for the Kaons and nearly no change for the $`\mathrm{\Xi }`$s.

These results, together with the steep $`\chi ^2`$ contours from Fig. 2, suggest that the relative particle abundances “freeze” shortly after the spontaneous breaking of chiral symmetry. The success of our fit suggests extremely rapid chemical equilibration (through abundance-changing reactions) in the state with broken symmetry. Fig. 2 shows that the chemical composition of the hadronic system has to change substantially within a small temperature interval, just before freeze-out, even for the crossover transition (i.e. parameter set CI); for reference, we have indicated the dynamical path in the $`T\mu _B`$ plane corresponding to the expansion of a perfect fluid (i.e. with constant entropy per net baryon sub86 ). While $`2n`$ reactions are perhaps too slow to explain such rapid chemical equilibration bass\_ch ; bass00 , $`mn`$ processes with several particles in the initial state may be important as well rapp00 ; grei01 ; grei02 ; grei02b .

Alternatively, the appearence of chemical equilibrium right after the phase transition (or the crossover) to the state of broken chiral symmetry might just be the outcome of the dynamical symmetry breaking process itself scav02 , with statistical occupation of the various hadronic channels according to phase space beca97 ; beca98 ; stock99a ; stock99b . If so, number-changing reactions in the broken phase need not proceed at a high rate. To test this picture experimentally, it might be useful to consider central collisions of small ions like protons or deuterons, at similar energy and particle densities in the central region as for central Au+Au. For systems of transverse extent comparable to the correlation lengths of the chiral condensates, the dynamical symmetry breaking process should be different from that in large systems (for example, the mean field approximation should not apply). The correlation lengths $`\xi _{\sigma ,\zeta }`$ are given by

$$\xi _\sigma ^2=\frac{^2\left(\mathrm{\Omega }/V\right)}{\sigma ^2},$$
(7)
and accordingly for $`\xi _\zeta `$. We evaluate the curvature of the thermodynamical potential at the global minimum and for $`T`$, $`\mu _B`$, $`\mu _S`$ at the freeze-out point. For parameter sets CI, CII, CIII we obtain $`\xi _\sigma =0.37`$ fm, 0.41 fm, 0.40 fm, respectively. For the correlation length of the strange condensate we obtain $`\xi _\zeta =0.20`$ fm in all three cases. The correlation lengths are not very much smaller than, say, the radius of a proton. Thus, even if the freeze-out point for high-energy $`pp`$ collisions happens to be close to that for Au+Au collisions at RHIC energies, the transition from the symmetry restored to the broken phase might be different. Finally, we also note that the correlation lengths obtained from our effective potential are not larger than the thermal correlation length $`1/T`$ at freeze-out, and so corrections beyond the mean-field approximation employed here should be analyzed in the future.

## V conclusion

Particle ratios as calculated in a chiral $`SU(3)\sigma \omega `$ model are compared with RHIC data for Au+Au at $`\sqrt{130}\text{ AGeV}`$ and with noninteracting gas calculations. Since different versions of the chiral model show qualitatively different phase transition scenarios, we investigate whether the particle production, i.e. the chemistry of the system, is sensitive to the phase transition behaviour. Since we have shown that the current data are described by all three different phase transition scenarios and the ideal gas model, we can so far not favour or rule out any one scenario. In all interacting models the effective masses at freeze-out are shifted up to $`15\%`$ from their vacuum values. The fitted chemical freeze-out temperatures and chemical potentials depend on the order of the phase transition. The crossover case yields $`15`$ MeV shifted $`T`$ values as compared to the noninteracting gas model while the models with a first order phase transition yield more than $`30`$ MeV lower temperatures. Furthermore, the fitted freeze-out points are located practically right on the phase transition boundary, in the first order phase transition scenarios, but $`T`$ is always $`T_c`$. *This suggests that at RHIC the system emerges after the chiral chiral phase transition. This of course is only true if a first order phase transition does actually occur in QCD at small chemical potentials and high $`\mathrm{T}`$.* ”Freezing” of the relative abundances of various hadrons in the symmetric phase (at $`T>T_c`$) is excluded.

###### Acknowledgements.

The authors are grateful to P. Braun-Munzinger, A. Dumitru, L. Gerland, M. Gorenstein, I. Mishustin, and K. Paech for fruitful discussions. This work is supported by Deutsche Forschungsgemeinschaft (DFG), Gesellschaft für Schwerionenforschung (GSI), Bundesministerium für Bildung und Forschung (BMBF), the Graduiertenkolleg Theoretische und Experimentelle Schwerionenphysik and by the U.S. Department of Energy, Nuclear Physics Division (Contract No. W-31-109-Eng-38).


text:
# Topology and correlations in structured scale-free networks

## I Introduction

Recently a major scientific effort has been devoted to the characterization and modeling of a wide range of social and natural systems that can be described as networks Albert and Barabási (2002); Dorogovtsev and Mendes (2002). Systems such as the Internet Faloutsos et al. (1999); Caldarelli et al. (2000); Pastor-Satorras et al. (2001); Vázquez et al. (2002) or the World-Wide-Web Albert et al. (1999), social communities Strogatz (2001), food-webs Montoya and Solé (2002), and biological interacting networks Wagner (2001); Jeong et al. (2001); Solé et al. (2002); Vázquez et al. (2001) can be represented as a graph Chartrand and Lesniak (1986), in which nodes represent the population individuals and links the physical interactions among them. Strikingly, many of these networks have complex topological properties and dynamical features that cannot be accounted for by classical graph modeling Erdös and Rényi (1960). In particular, small-world properties Watts and Strogatz (1998) and scale-free degree distributions Barabási and Albert (1999) (where the degree or connectivity of a node is defined as the number of other nodes to which it is attached) seem to emerge frequently as dominant features governing the topology of real-world networks. These global properties imply a large connectivity heterogeneity and a short average distance between nodes, which have considerable impact on the behavior of physical processes taking place on top of the network. For instance, scale-free (SF) networks have been shown to be resilient to random damage (absence of a percolation transition) Albert et al. (2000); Callaway et al. (2000); Cohen et al. (2001) and prone to epidemic spreading (null epidemic threshold) Pastor-Satorras and Vespignani (2001a, b); May and Lloyd (2001); Moreno et al. (2002).

The detailed scrutiny of the topological properties of networks has pointed out that small-world and scale-free properties come often along with non-trivial degree correlations and clustering properties. Recently, an interesting class of networks has been introduced by Klemm and Eguíluz by proposing a growing model in which nodes are progressively deactivated with a probability inversely proportional to their connectivity Klemm and Eguíluz (2002a). Analytical arguments and numerical simulations have lead to the claim that, under general conditions, the deactivation model, allowing a core of $`m`$ active nodes, generates a network with average degree $`k=2m`$ and degree probability distribution $`P(k)=2m^2k^3`$. Interestingly, the scale-free properties are associated to a high clustering coefficient. For this reason the deactivation model has been used to study how clustering can alter the picture obtained for the resilience to damage and epidemic spreading in SF networks Eguíluz and Klemm (2002); Crucitti et al. (2002).

In this paper we revisit the analysis of the deactivation model. We find an analytical solution in the case of minimal values of active nodes $`m`$ (low average connectivity). In addition, large scale numerical simulations exhibits a noticeable variability of the degree distribution with $`m`$. In particular, the degree exponent strongly depends on $`m`$ for the general case considered in Ref. Klemm and Eguíluz (2002a). The model topology is also susceptible to several details of the construction algorithm. By means of large-scale numerical simulations we study the deactivation model topology in the whole range of $`m`$ and for different algorithm parameters. We calculate analytically the clustering coefficient and connectivity correlation functions. Also in this case a variability with respect to the model parameters is found. Extensive numerical simulations confirm the analytical picture presented here.

In the generated networks, we also report the lack of small-world properties. In the whole parameters range we find a network diameter increasing linearly with the number of nodes forming the network Klemm and Eguíluz (2002b). The networks’ topology is therefore similar to an array of connected star-shaped graphs. The networks are thus similar to a one dimensional lattice in what concerns their physical properties. In particular, diffusion and spreading processes might be heavily affected by the increasing average distance among nodes that make the system alike to a one-dimensional chain. In this perspective, we discuss the properties of epidemic spreading and resilience to damage in networks generated with the deactivation model.

## II Deactivation model

The deactivation model introduced by Klemm and Eguíluz Klemm and Eguíluz (2002a) is defined as follows: Consider a network with directed links. Each node can be in two states, either active or inactive. The model starts from a completely connected graph of $`m`$ active nodes and proceeds by adding new nodes one by one. Each time a node is added:

1. It is connected to all active nodes in the network;
2. One of the active nodes is selected and set inactive with probability

$$p_d(k_i^{\mathrm{in}})=\frac{\left[_{j𝒜}(a+k_j^{\mathrm{in}})^1\right]^1}{a+k_i^{\mathrm{in}}};$$
(1)
3. The new node is set active.

The sum in Eq. (1) runs over the set of active nodes $`𝒜`$, $`a`$ is a model parameter, and $`k_i^{\mathrm{in}}`$ denotes the in-degree of the $`i`$-th node.

As we shall show below, this model is quite sensitive to the order in which steps 2 and 3 are performed and, therefore, it is better to discriminate the following cases.

* Model A: step 2 is performed before step 3.
* Model B: step 2 is performed after step 3.

For $`m\mathrm{}`$ both models can be solved analytically in the continuous $`k^{\mathrm{in}}`$ approximation, after introducing the probability density that an active node has in-degree $`k^{\mathrm{in}}`$ Klemm and Eguíluz (2002a). Moreover in this limit the order of steps 2 and 3 is irrelevant, obtaining the same in-degree distribution $`P(k^{\mathrm{in}})(a+k^{\mathrm{in}})^\gamma `$ with

$$\gamma =2+\frac{a}{m}.$$
(2)
The model is usually simulated by using $`a=m`$. In this way the deactivation probability is inversely proportional to the total connectivity of the nodes $`(m+k^{\mathrm{in}})^1`$ and the connectivity distribution results to be $`P(k)=2m^2k^3`$. Interestingly, due to the deactivation mechanism the networks show a high clustering coefficient that approaches a constant value in the infinite size limit Klemm and Eguíluz (2002a).

At lower values of $`m`$, it has been claimed that finite size effects set in and the connectivity distribution shows deviations from the predicted behavior. We shall see in the next section that for $`a=m10`$ the model presents a very different analytical solution which yields a connectivity distribution very far from the $`m\mathrm{}`$ limit. In addition, the deactivation model topology is very sensible to changes in the details of the growing algorithms.

## III Degree distribution

### III.1 Model A

Let us first focus on model A with $`m=2`$, that is the smallest value of $`m`$ for which the model is non-trivial. In this case, after adding a new node we have only two nodes at the deactivation step. One of them will be set inactive and replaced by the new added node, that has in-degree $`0`$. In the worst case, the other node will have in-degree $`0+1`$, the $`0`$ coming from its initial in-degree and the $`1`$ from the connection to the newly added node, and in general it will have in-degree larger than or equal to $`2`$. Later on, at the next deactivation step, the in-degrees of both nodes will have increased by one resulting in one active node with in-degree $`1`$ and another with in-degree $`K2`$, where $`K`$ is the in-degree of the active node with largest in-degree, that coincides with the oldest node. Then, following Eq. (1), one of them will be deactivated with probability

$$p_d(K)=\frac{1+a}{1+2a+K},p_d(1)=1p_d(K).$$
(3)
Each time the oldest node is not deactivated its in-degree increases by one and, therefore, the probability that the oldest node has in-degree $`K`$ is just the probability that it is not deactivated in $`K2`$ steps, with running in-degree $`2,3,\mathrm{}K1`$. Thus, the probability $`\stackrel{~}{P}(K)`$ of creating a deactivated node of in-degree $`K`$ is equal to the probability that the largest node is not deactivated in $`K2`$ steps and is deactivated in the last step, i.e.:

$`\stackrel{~}{P}(K)`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=2}{\overset{K1}{}}}\left[1p_d(\mathrm{})\right]p_d(K)`$ (4)
$`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }(3+2a)}{\mathrm{\Gamma }(1+a)}}{\displaystyle \frac{\mathrm{\Gamma }(a+K)}{\mathrm{\Gamma }(2+2a+K)}},`$
where $`\mathrm{\Gamma }(x)`$ is the standard Gamma function Abramowitz and Stegun (1972). On the other hand, every time that the oldest node is not deactivated, the other, with in-degree $`1`$, is deactivated. Hence, in the $`K1`$ deactivation steps leading to the generation of a node with in-degree $`K`$, $`K2`$ nodes with in-degree $`1`$ are created. The average number of nodes with in-degree $`1`$ created in the process is then

$$\stackrel{~}{P}_1=\underset{K=2}{\overset{\mathrm{}}{}}(K2)\stackrel{~}{P}(K)=\frac{2+a}{a}.$$
(5)
Therefore, the in-degree distribution will be given by

$$P(k^{\mathrm{in}})=C^1\times \{\begin{array}{cc}\stackrel{~}{P}_1& k^{\mathrm{in}}=1\hfill \\ & \\ \stackrel{~}{P}(k^{\mathrm{in}})& k^{\mathrm{in}}>1\hfill \end{array},$$
(6)
where $`C`$ is a constant, obtained from the normalization condition $`_{k^{\mathrm{in}}}P(k^{\mathrm{in}})=1`$, that has the value

$$C=\stackrel{~}{P}_1+\underset{K=2}{\overset{\mathrm{}}{}}\stackrel{~}{P}(K)=\frac{2a+2}{a}.$$
(7)
From this equation, we obtain the analytic expression for the in-degree distribution

$$P(k^{\mathrm{in}})=\{\begin{array}{cc}\frac{2+a}{2+2a}& k^{\mathrm{in}}=1\hfill \\ & \\ \frac{\mathrm{\Gamma }(2+2a)}{\mathrm{\Gamma }(a)}\frac{\mathrm{\Gamma }(a+k^{\mathrm{in}})}{\mathrm{\Gamma }(2+2a+k^{\mathrm{in}})}& k^{\mathrm{in}}>1\hfill \end{array}.$$
(8)
For large $`k^{\mathrm{in}}`$ we can expand the previous expression using Stirling’s approximation, to obtain that the in-degree distribution follows the asymptotic behavior

$$P(k^{\mathrm{in}})k_{}^{\mathrm{in}}{}_{}{}^{\gamma },\gamma =2+a.$$
(9)
Moreover, since the out-degree of all nodes is $`m`$ then the degree $`k`$ of a node (in-degree plus out-degree) is $`m+k^{\mathrm{in}}`$ and will follow the same distribution shifted by $`m`$. For the particular case $`a=m=2`$, the degree distribution takes the form

$$P(k)=\{\begin{array}{cc}\frac{2}{3}& k=3\hfill \\ & \\ \frac{120}{k(k+1)(k+2)(k+3)}& k>3\hfill \end{array}$$
(10)
In Fig. 1 we plot the degree distribution obtained from numerical simulations of model $`A`$ for $`a=m`$. For $`m=2`$ the numerical points are in very good agreement with the exact distribution given in Eq. (10), with a power law decay with exponent $`\gamma =2+a=4`$. In the limiting case of large $`m`$ the continuous approach predicts the exponent 3 Klemm and Eguíluz (2002a) (see Eq. (2)), giving us a lower bound. Hence,

$$\text{model A with }a=m3<\gamma 4$$
(11)
and, therefore, the degree distribution has always a bounded second moment. For larger $`m`$ the distribution follows a power law decay but with an exponent $`\gamma `$ that depends on $`m`$. In order to show that the degree distribution approaches for each $`m`$ an asymptotic power law behavior with $`\gamma >3`$ we performed large scale simulations of networks with $`N=10^7`$ nodes. In Fig. 1 we report the behavior of the exponent $`\gamma `$ as a function of $`m`$. For all values of $`m<10`$ the degree exponents strongly deviates from the $`m\mathrm{}`$ limit.

### III.2 Model B

Using similar arguments we can compute the degree distribution of model B for $`m=1`$. In this case we also have two nodes at the deactivation process, the one just added and the one surviving from the previous deactivation step. The former has in-degree $`0`$ while the latter (the oldest) has in-degree $`K1`$, and one of them is deactivated with probability

$$p_d(K)=\frac{a}{2a+K},p_d(0)=1p_d(K).$$
(12)
The probability that when the oldest node is deactivated it has degree $`K`$ is given by

$`\stackrel{~}{P}(K)`$ $`=`$ $`{\displaystyle \underset{\mathrm{}=1}{\overset{K1}{}}}\left[1p_d(\mathrm{})\right]p_d(K)`$ (13)
$`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }(1+2a)}{\mathrm{\Gamma }(a)}}{\displaystyle \frac{\mathrm{\Gamma }(a+K)}{\mathrm{\Gamma }(1+2a+K)}}.`$
In the process of creating a node of in-degree $`K`$, $`K1`$ nodes of in-degree $`0`$ have been created. The average number of nodes with in-degree $`0`$ created is

$$\stackrel{~}{P}_0=\underset{K=1}{\overset{\mathrm{}}{}}(K1)\stackrel{~}{P}(K)=\frac{a+1}{a1}.$$
(14)
Thus, the analytic expression for the normalized in-degree distribution is given by

$$P(k^{\mathrm{in}})=C^1\times \{\begin{array}{cc}\stackrel{~}{P}_0& k^{\mathrm{in}}=0\hfill \\ & \\ \stackrel{~}{P}(k^{\mathrm{in}})& k^{\mathrm{in}}>0\hfill \end{array},$$
(15)
with the normalization constant

$$C=\stackrel{~}{P}_0+\underset{k=1}{\overset{\mathrm{}}{}}\stackrel{~}{P}(K)=\frac{2a}{a1}.$$
(16)
From here follows the expression for the degree distribution (where $`k=m+k^{\mathrm{in}}`$)

$$P(k)=\{\begin{array}{cc}\frac{1+a}{2a}& k=1\hfill \\ & \\ \frac{\mathrm{\Gamma }(2a)}{\mathrm{\Gamma }(a1)}\frac{\mathrm{\Gamma }(a+k1)}{\mathrm{\Gamma }(2a+k)}& k>1\hfill \end{array}.$$
(17)
For large $`k`$ the degree distribution follows the asymptotic behavior

$$P(k)k^\gamma ,\gamma =1+a.$$
(18)
In Fig. 2 we show the degree distribution obtained from numerical simulations of model B with $`a=m`$. For $`a=m=1`$ we recover the predicted exponent $`\gamma =2`$. Also in this case, we provide large scale numerical simulations ($`N=10^7`$) of networks with larger values of $`m`$. The obtained distributions still follow a power law decay but with an exponent $`\gamma `$ that is a continuously increasing function of $`m`$. It is worth remarking that for $`m<10`$ the degree exponent is stable and strongly differs from the value $`\gamma =3`$.

It is worth noticing that for $`a=m=1`$ the analytic solution Eq. (17) is singular, as can be readily seen from the $`\mathrm{\Gamma }(a1)`$ factor in the denominator. In fact, the solution in this case is $`P(k)=\delta _{k,1}`$, that is, in the thermodynamic limit (infinitely large network), the weight of the nodes with degree $`1`$ is overwhelming with respect to the nodes with different connectivity. This singularity is rooted in the fact that the distribution, with exponent $`\gamma =2`$, lacks a finite first moment in the thermodynamic limit, while we know that, by definition, model $`B`$ has average connectivity $`k=2`$. This necessarily implies that there must be an implicit dependence on the network size $`N`$ in the degree distribution for $`a=m=1`$, dependence that cannot be assessed by our analytic solution since we are already working in the infinite network limit. We can nevertheless estimate the functional form of the degree distribution for a finite network composed by $`N`$ nodes, which has a maximum connectivity $`k_c`$, such that there are no nodes with degree larger that $`k_c`$. Assuming that the distribution for $`k>1`$ follows the same functional form as Eq. (17), we have that for $`a=1`$

$$P_N(k)=\{\begin{array}{cc}C_1& k=1\hfill \\ & \\ \frac{C_2}{k(k+1)}& 1<kk_c\hfill \end{array},$$
(19)
where $`C_1`$ and $`C_2`$ are constants to be determined by the normalization conditions $`_{k=1}^{\mathrm{}}P_N(k)=1`$ and $`_{k=1}^{k_c}kP_N(k)=2`$ (the upper limit in the first normalization condition can be taken to be infinite, since the corrections stemming from $`k_c`$ are of lower order). From this two conditions we obtain, in the continuous $`k`$ approximation that replaces sums by integrals,

$$C_1=1\frac{2\mathrm{ln}(3/2)}{\mathrm{ln}\left({\displaystyle \frac{1+k_c}{2}}\right)},C_2=\frac{2}{\mathrm{ln}\left({\displaystyle \frac{1+k_c}{2}}\right)}.$$
(20)
For finite SF networks with degree distribution $`P(k)k^\gamma `$, the maximum degree $`k_c`$ scales with the number of nodes as $`k_cN^{1/(\gamma 1)}`$ Dorogovtsev and Mendes (2002). In the present case we have $`k_cN`$, and thus, for large $`N`$,

$$1C_1\frac{1}{\mathrm{ln}N},C_2\frac{1}{\mathrm{ln}N}.$$
(21)
Therefore, in the limit $`N\mathrm{}`$, we recover a singular degree distribution with $`C_11`$ and $`C_20`$. We can check numerically this result by noticing that, from Eqs. (19) and (21), the degree distribution at fixed $`k`$ should scale as

$$1P_N(1)\frac{1}{\mathrm{ln}N},P_N(k)\frac{1}{\mathrm{ln}N}k>1.$$
(22)
We have verified this scaling form in Fig. 3. Therefore, in model B with $`a=m=1`$ we obtain a degree distribution that decays as $`k^2`$, but with a normalization constant for $`k>1`$ that decays with the network size as $`1/\mathrm{ln}N`$. Finally, it is worth mentioning that the second moment of the distribution is diverging as $`k^2N/\mathrm{ln}N`$. Despite this singular behavior for $`a=m=1`$, however, Eq. (17) remains exact for any value of $`a1`$.

From the results of Fig. 2, together with the upper bound $`\gamma =3`$ obtained from the large $`m`$ approximation Klemm and Eguíluz (2002a), we have that

$$\text{model B with }a=m2\gamma <3,$$
(23)
and, therefore, the degree distribution has a divergent second moment.

The analysis made above has shown that the deactivation model is quite sensitive to the order in which steps 2 and 3 are performed, yielding degree distributions with a finite or divergent second moment, depending on the order. In addition, the exponent $`\gamma `$ is rather sensible to the value of $`a=m`$, showing a wide range of variation. This fact has not been noticed in previous works where this model has been considered Klemm and Eguíluz (2002a); Eguíluz and Klemm (2002); Crucitti et al. (2002), prompting that some of the conclusions obtained in those works should be reconsidered in this perspective.

## IV Clustering coefficient

We can go beyond the degree distribution and compute the clustering coefficient $`c(k)`$, as a function of the node degree $`k`$ Vázquez et al. (2002); Klemm and Eguíluz (2002b). For this quantity we can perform an analytic calculation for any value of $`a`$ and $`m`$ and for both models A and B. In order to compute the clustering coefficient, we will consider the network as undirected and denote by $`k_i=k_i^{\mathrm{in}}+m`$ the total degree of the node $`i`$.

The clustering coefficient of the node $`i`$ is defined by Watts and Strogatz (1998)

$$c_i=\frac{2e_i}{k_i(k_i1)},$$
(24)
where $`e_i`$ is the number of edges between the neighbors of node $`i`$ and it is divided by its maximum possible value $`k_i(k_i1)/2`$. In the deactivation model new edges are created between the active nodes and the added node. Hence, $`e_i`$ remains constant for inactive nodes and increases only for the active ones. Moreover, all the active nodes are connected. Hence, each time we add a node the degree $`k_i`$ of each active node $`i`$ increases by one and $`e_i`$ increases by $`m1`$, where $`m1`$ are just the new links between the new neighbor of $`i`$ (the added node) and the remaining active nodes. Therefore, the dynamics of $`e_i`$ is given by

$$\frac{e_i}{t}=(m1),$$
(25)
while the connectivity obeys the relation $`k_i(t)=m+t`$. Here $`t=0`$ corresponds to the time at which the node $`i`$ was created. Besides, when the node is added it has degree $`m`$, thus $`e_i(0)=m(m1)/2`$ and, therefore, $`c_i(0)=1`$. Integrating Eq. (25) with this initial condition and substituting the result in Eq. (24), taking into account that $`t=k_im`$, we obtain

$`c(k)`$ $`=`$ $`{\displaystyle \frac{m(m1)}{k(k1)}}+{\displaystyle \frac{2(m1)(km)}{k(k1)}}`$ (26)
$`=`$ $`{\displaystyle \frac{2(m1)}{k}}{\displaystyle \frac{(m1)(m2)}{k(k1)}},`$
where the last expression in Eq. (26) is obtained after some algebraic manipulations. Eq. (26) recover the result previously obtained in Ref.Klemm and Eguíluz (2002b). For $`m=1`$ the network is a tree, and therefore we obviously recover $`c(k)=0`$. For $`m=2`$ we obtain the exact behavior $`c(k)=2/k`$. For $`m>2`$, the asymptotic behavior for large $`k`$ is $`c(k)1/k`$ Klemm and Eguíluz (2002a). Interestingly, we recover in this model the same behavior of $`c(k)`$ found in other systems in Ref. Ravasz and Barabási (2002).

In Fig. 4 we plot the clustering coefficient as a function of the node degree obtained for model A and B and different values of $`m`$ from numerical simulations. As it can be seen the numerical dependency coincides with the analytical expression in Eq. (26).

## V Degree correlation function

Degree correlations can be characterized by analyzing the nearest neighbor average degree introduced in Refs. Pastor-Satorras et al. (2001); Vázquez et al. (2002), defined as

$$k_{nn,i}=\frac{D_i}{k_i},$$
(27)
where $`D_i`$ is the sum of the degrees of the neighbors of node $`i`$. In uncorrelated networks, the quantity $`k_{nn,i}`$ does not show any dependence on the degree of the node $`i`$. This is not the case when degree correlations are present. In this case $`k_{nn,i}`$ is a function of the degree of the node whose nearest neighbors are analyzed. In particular, we can face two possible kinds of correlation. In the first situation, nodes with high connectivity will connect more preferably to highly connected nodes; a property referred to as “assortative mixing”. On the opposite side, it is possible to have “dissortative mixing”; i.e. highly connected nodes are preferably connected to nodes with low connectivity Newman (2002).

In the deactivation model, when the node is added it has degree $`m`$ and $`D_i=mk_𝒜`$, where $`k_𝒜`$ is the average degree among active neighbors. Then, if the node $`i`$ is active, it is, by construction, neighbor of the $`m1`$ remaining active nodes. Thus, every time a new node is added, $`k_i`$ increases by one and $`D_i`$ increases by $`(m1)+m`$, the $`m1`$ because the degree of the remaining $`m1`$ neighbors have also increase by one and the $`m`$ because the new neighbor has degree $`m`$. Hence

$$\frac{D_i}{t}=(2m1),$$
(28)
for each active node $`i`$. Integrating this equation, taking into account the initial condition $`D_i(0)=mk_𝒜`$ and the relation $`t=k_im`$, we obtain that

$$D_i^{}=(2m1)(k_im)+mk_𝒜$$
(29)
when the node $`i`$ is deactivated. Now, when an active node becomes inactive its degree remains fixed but the degree of its active neighbor nodes will still increase until they get deactivated. Therefore, in the infinite time limit, we have

$$D_i=D_i^{}+\mathrm{\Delta }D_i,$$
(30)
where $`\mathrm{\Delta }D_i`$ is the increase of $`D_i`$, since node $`i`$ was set inactive until all its neighbors are set inactive.

Hence, from Eqs. (27), (29), and (30), it follows that

$$k_{nn,i}=2m1+\frac{mk_𝒜+\mathrm{\Delta }D_im(2m1)}{k_i}.$$
(31)
It remains now the task to assess the possible dependence of $`\mathrm{\Delta }D_i`$ on the connectivity $`k_i`$ (it is clear that the long time average of $`k_𝒜`$ must be independent of the connectivity of any deactivated node). For the minimum $`m`$ ($`m=2`$ for model A and $`m=1`$ for model B) the degree of an active node set inactive is not correlated with the degree of the remaining active nodes, since those remaining nodes have always degrees $`2`$ and $`3`$ in model A with $`m=2`$, and degree $`1`$ in model B with $`m=1`$, independently of the degree of the last deactivated node. Therefore in this case $`\mathrm{\Delta }D_i`$ cannot depend on $`k_i`$. This lack of correlations is also clear for $`m1`$ where the sum $`_{j𝒜}(a+k_j)^1`$ in Eq. (1) is a constant Klemm and Eguíluz (2002a) and, therefore, the degree of the active nodes in not correlated with the degree of the inactive nodes. For intermediate values of $`m`$, however, the degree of the active nodes may be correlated in such a way that $`\mathrm{\Delta }D_i`$ depends on $`k_i`$.

In Fig. 5 we plot the dependency of the average nearest neighbors degree $`\overline{k}_{nn}(k)`$ as a function of the degree $`k`$ for models A and B and different values of $`m`$. In the case of model A, $`\overline{k}_{nn}(k)(2m1)1/k`$ even for $`m2`$, in agreement with Eq. (31). In the case of model B, $`\overline{k}_{nn}(k)(2m1)1/k`$ for the $`m=1`$ and $`m=10`$ but decays faster for intermediate values of $`m`$. Thus, in this case the correlations between the active node degrees introduce stronger deviations for intermediate values of $`m`$. In all cases, however, we find that correlations in the deactivation model are of “disassortative” nature; i.e. highly connected nodes are preferably connected with poorly connected nodes. It is also worth stressing that the results for model B with $`m=1`$ must be taken with a grain of salt, given the singular nature of the model exposed in Sec. III.2.

In the deactivation model, either A or B, for a fixed network size $`N`$ and assuming that $`\mathrm{\Delta }D_i`$ does not grow faster than $`k_i`$, we have that in the limit $`k_i\mathrm{}`$, $`k_{nn,i}2m1`$. That is, the average nearest neighbor degree of the hubs (nodes with largest $`k_i`$) equals $`k1`$, as previously pointed out in Ref. Eguíluz and Klemm (2002). However, this fact does not necessarily imply that $`\mathrm{\Delta }D_i`$ is independent of $`N`$. One way to check this point is to compute the average of $`k_{nn,i}`$ over all nodes, $`\overline{k}_{nn}_N=_kP(k)\overline{k}_{nn}(k)`$. Let us assume that $`\overline{k}_{nn}k1+\alpha /k`$, where $`\alpha `$ is depending on $`\mathrm{\Delta }D_i`$. If $`\mathrm{\Delta }D_i`$ is approaching a constant value, we should obtain $`\overline{k}_{nn}_N`$ const., independently of $`N`$. In Fig. 6 we show how $`\overline{k}_{nn}_N`$ behaves with increasing $`N`$ for $`a=m`$. For model A, where $`3<\gamma 4`$, it approaches a stationary value for $`N1`$. Moreover, the asymptotic limit of $`\overline{k}_{nn}_N`$ increases with increasing $`m`$. In fact, with increasing $`m`$ the exponent $`\gamma `$ decreases approaching the limit $`\gamma =3`$ for $`m1`$, where $`k^{nn}`$ diverges logarithmically with $`N`$. On the contrary, for model B, where $`2\gamma <3`$, $`\overline{k}_{nn}_N`$ is growing with $`N`$ following a power law. This implies that $`\mathrm{\Delta }D_i`$ is a diverging function of $`N`$ and that in the thermodynamic limit (in which we perform first the limit $`N\mathrm{}`$) the average nearest neighbor connectivity curve is progressively shifting to larger and larger values. This finally points out that the average nearest neighbor connectivity of hubs is not a well-defined quantity since the $`k_i\mathrm{}`$ limit must be performed only after the $`N\mathrm{}`$ limit. The divergence of $`\overline{k}_{nn}_N`$ with $`N`$ is related to a general property of SF networks with diverging connectivity fluctuations and it is dictated by the detailed balance of connectivity Boguñá and Pastor-Satorras (2002); Boguñá et al. (2002).

## VI Diameter and shortest path length

Another fundamental topological feature of complex networks is identified by the scaling of the average path length among nodes and the network’s diameter. The minimum path between two nodes is given by the minimum number of intermediate nodes that must be traversed to go from node to node. The average minimum path length $`d`$ is thus defined as the minimum path distance averaged over all the possible pairs of nodes in the network. Similarly, the network diameter is defined as the largest among the shortest paths between any two nodes in the network.

While regular networks (for instance hypercubic lattices) have a diameter scaling with the size $`N`$ as the inverse of the Euclidean dimension, many complex networks show striking small-world properties; i.e. in average one can go from one node to any other node in the system by passing through a very small number of intermediate nodes Watts and Strogatz (1998). In this case the graph diameter grows logarithmically, or even slower, with the system’s number of nodes $`N`$.

In Ref. Klemm and Eguíluz (2002b), it has been noticed that for large $`m`$ values $`d`$ is scaling linearly with the network size $`N`$. In the deactivation model (A and B) we measured both the diameter and the average minimum path distance $`d`$ as a function of $`N`$ for values of $`a=m`$ ranging from $`1`$ to $`4`$. In all cases we find that after a small size transient both metrics approach a linear scaling with $`N`$. In Fig. 7 we report the results obtained in the case of the deactivation model with rule B. This evidence implies that the topology of the generated networks is approaching those of a one dimensional lattice. In other words, the deactivation model does not exhibit small-world properties.

In order to provide a visual representation of the deactivation model topology, we report in Fig. 8 the illustration of a network generated with model B and $`a=m=3`$. The linear structure of the network is evident. In particular, we find a long chain of sequentially connected star-like structures. The star-like structures correspond to groups of nodes connected to a node which has been active for longer times and has had the possibility to develop a high number of connections. Once these hubs are deactivated, they do not receive any further connection. The network grows by adding bridge nodes which are rapidly deactivated, until a new star-like structure is developed by a node that is active long enough. The growth mechanism, however, does not allow the formation of shortcuts between the deactivated region of the network and the new active nodes, hindering the development of small-world properties. The linear chain is therefore reflecting the time evolution of the structure: recently added nodes are separated from the original core of active nodes by a sequence of deactivated nodes that increases proportionally to the network size. By inspecting networks with larger $`m`$ we find very similar structures, with an increasing size of the star-like structures forming the linear chain. As we shall discuss in the last section, the absence of small-world properties might have a relevant effect in many physical properties of the network.

## VII Discussion and conclusions

In the present work we have provided a detailed analysis of the deactivation model introduced in Ref. Klemm and Eguíluz (2002a). The model shows a rich behavior, being very sensible to the value of the parameters used in the model and slight variations of the growing algorithm. The most striking result is that the degree distribution is depending on the value of the number of simultaneously active nodes $`m`$ also in the case in which $`a=m`$;  i.e. when the deactivation probability is related to the nodes’ total degree. The degree exponent is asymptotically approaching the value $`\gamma =3`$ only for $`m\mathrm{}`$, and the SF properties of networks suffer large variations in the range $`1m10`$. Along with the high clustering observed in previous works, we find that the model exhibits interesting degree correlation properties. In particular, we find marked disassortative mixing properties; i.e. highly connected nodes link preferably to poorly connected nodes. The analytical expression for the degree correlation is obtained and recovered by numerical simulations. Strikingly, the SF and correlation properties are not associated with small-world properties. The numerical analysis shows that for all values of $`m`$ the network diameter is increasing linearly with the number of nodes. The network thus approaches a linear structure, lacking long-range shortcuts.

One of the most interesting issues related to SF networks is the effects of their complex topological features on the dynamics of spreading phenomena Pastor-Satorras and Vespignani (2001a, b); Lloyd and May (2001); Moreno et al. (2002) and the onset of percolation transitions Albert et al. (2000); Callaway et al. (2000); Cohen et al. (2001). In the case of random SF networks, where degree correlations are absent, it has been found that the epidemic threshold is proportional to $`k/k^2`$ Pastor-Satorras and Vespignani (2001a, b). Uncorrelated SF networks allow the onset of large epidemics whatever the spreading rate of the infection. This is a noticeable result that has a large impact in immunization as well as control and design policies in real networks Pastor-Satorras and Vespignani (2001c); Dezsö and Barabási (2002). On the other hand, most real networks show non-trivial degree correlations and clustering properties as it is the case in the present deactivation model. Similarly, the random removal of nodes does not destroy the connectivity of SF networks with $`\gamma 3`$. In other words the percolation transition is absent, and the networks are extremely robust to random damages Albert et al. (2000); Callaway et al. (2000); Cohen et al. (2001). A natural question is to know whether or not the clustering properties of SF networks plus their correlations alter the general results obtained for uncorrelated networks. For this reason, several recent works have addressed the effect of such correlations in the epidemic spreading occurring on these networks Eguíluz and Klemm (2002); Crucitti et al. (2002); Boguñá and Pastor-Satorras (2002); Warren et al. (2002). In particular, in Ref. Eguíluz and Klemm (2002) it has been claimed the existence of an epidemic threshold in the case of the deactivation model for rule B.

The presence of a finite threshold in the deactivation model has been traced back to the high clustering coefficient and the finite limit of the average nearest neighbor connectivity of the largest hubs Eguíluz and Klemm (2002). On the other hand, we have shown here that the average nearest neighbor connectivity in the system is diverging with the system size. What appears as more fundamental for the properties of spreading in the deactivation model is its linear structure, with a diameter that increases with $`N`$. In a coarse grained picture, the epidemic spreading is dominated by the diffusion of the disease on a linear chain. In order to check this point, we have simulated a standard random walk in the B model with $`m=3`$. In Fig. 9 we plot the mean-square displacement of the random walker, $`R^2(t)^{1/2}=[r(t)r(0)]^2^{1/2}`$, where the brackets denote an average over $`250`$ realizations of the random walk on $`250`$ different networks. For a purely diffusive system, as would be the case of a one dimensional lattice, we would expect a scaling $`R^2(t)^{1/2}t^{1/2}`$. For the deactivation model we observe a slightly subdiffusive behavior, with a mean-square displacement scaling as $`R^2(t)^{1/2}t^{0.46}`$. We thus conclude that dynamics on the deactivation model is almost purely diffusive, as expected from its non small-world character. The analysis of spreading and percolation properties in this network cannot therefore be performed at the mean-field level Pastor-Satorras and Vespignani (2001a, b), but must include diffusion and most probably fluctuations, leading to a much more complex formalism based in a field theory Marro and Dickman (1999). For the sake of comparison, we have also plotted in Fig. 9 the mean-square displacement of a random walker on a Barabási-Albert network Barabási and Albert (1999) and on a Internet snapshot map from 1999, collected by the National Laboratory for Applied Network Research (NLANR) nla . As we observe, in these last two networks $`R^2(t)^{1/2}`$ saturates very quickly to a constant value, proportional to the network’s diameter, indicating the presence of a strong small-world component. The essential difference of the diffusive properties between the Internet and the deactivation model does not allow to extend the conclusions obtained from the model to the spreading in the real system.

The same applies to percolation properties, that naturally exhibit a finite threshold in this case. The fact that spreading and percolation properties on the deactivation model are similar to those of regular lattices because of the absence of small-world features is corroborated by the analysis of Ref. Crucitti et al. (2002), that shows how the introduction of a small amount of shortcuts restores the usual absence of a percolation threshold. In this perspective, it would be extremely interesting to have a detailed study of the epidemic spreading properties in the case of the deactivation model with random rewiring Klemm and Eguíluz (2002b), in order to assess the effect of clustering and degree correlations in spreading processes in SF networks with small-world properties.

## VIII Acknowledgments

This work has been partially supported by the European commission FET Open project COSIN IST-2001-33555. R.P.-S. acknowledges financial support from the Ministerio de Ciencia y Tecnología (Spain).


text:
# Logarithmic roughening in a growth process with edge evaporation

## I Introduction

A large variety of models for surface growth is known to display universal scaling laws Meakin ; KrugSpohn ; Dietrich ; Barabasi ; Krug . In most cases one observes simple power-law scaling, i.e., the width $`W`$ of an initially flat interface grows with time as

$$Wt^\gamma ,$$
(1)
where $`\gamma `$ is the so-called growth exponent. In a finite system the width eventually saturates at a finite value $`W_{\mathrm{sat}}`$, which grows with the linear system size $`L`$ as

$$W_{\mathrm{sat}}L^\alpha ,$$
(2)
where $`\alpha `$ is the roughness exponent. Both asymptotic power laws can be combined in a single scaling form

$$W(L,t)=L^\alpha f(t/L^z),$$
(3)
where $`z=\alpha /\gamma `$ is the dynamic exponent and $`f`$ is a scaling function with an appropriate asymptotic behavior. This type of power-law scaling, also known as Family-Vicsek scaling FamilyVicsek , is generic for a large number of self-similar growth processes. However, in some cases the interface was found to roughen logarithmically with time. This may happen, for example, when a system undergoes a roughening transition. The purpose of this paper is to study a simple model with such a logarithmic scaling behavior in more detail. To this end we consider a class of solid-on-solid growth models introduced a few years ago by Alon et al. Alon ; Alon2  which have the special feature that atoms can only evaporate at the edges of terraces. The models exhibit a robust phase transition from a non-moving smooth phase to a moving phase, where the interface roughens continuously. Remarkably, this roughening transition takes place even in one spatial dimension. At the critical point the interface width was found to increase logarithmically with time. However, so far it was impossible to unfold the complete scaling picture since previous numerical simulations were not accurate enough Alon2 . Here we close this gap by performing high-precision simulations based on a parallelized code. As a main result we find that Eq. (3) has to be replaced by a logarithmic scaling law of the form

$$W^2(L,t)=\lambda \mathrm{ln}L+f(t/L^z),$$
(4)
Moreover, we present a simple approximation which explains why the interface roughens logarithmically.

## II Growth model with evaporation at the edges of plateaus

The growth models investigated in Ref. Alon  are defined as solid-on-solid deposition-evaporation processes with the special property that atoms can only evaporate at the edges of terraces. The models are defined on a $`d`$-dimensional lattice with $`L`$ sites and periodic boundary conditions, where each site $`i`$ is associated with an integer variable $`h_i`$. In Ref. Alon  two variants of models have been considered. Here we will focus on the physically motivated restricted model in one dimension, where the heights at neighboring sites obey the inequality

$$|h_ih_{i+1}|1.$$
(5)
The model evolves by random-sequential updates, i.e., a site $`i`$ is randomly selected and one of the following moves is carried out (see Fig. 1). With probability $`q`$ an atom is deposited at site $`i`$, increasing $`h_i`$ by $`1`$. Otherwise one of the nearest neighbors is randomly selected. If the selected neighbor is at a lower height, indicating the edge of a terrace, an atom is removed from site $`i`$. In both cases a move is abandoned if the resulting configuration would violate the condition (5). Each attempted update corresponds to an average time increment $`\mathrm{\Delta }t=1/L`$.

In the limit $`L\mathrm{}`$ the model has the following phenomenological properties. If the growth rate $`q`$ is small, the interface is flat and pinned to a spontaneously selected bottom layer. Small islands will occasionally grow on top of this bottom layer but will quickly be eliminated by desorption at the islands edges so that the system eventually reaches a stationary state characterized by a finite width. As $`q`$ is increased, more and more islands on top of the bottom layer are produced, until above $`q_c`$, the critical value of $`q`$, the islands merge and the interface starts moving. Approaching $`q_c`$ from below, the stationary width is found to diverge logarithmically.

In Refs. Alon ; Alon2  it was shown that the roughening transition in this class of models is driven by a directed percolation (DP) transition Hinrichsen00a  at the spontaneously selected bottom layer. This means that the density of exposed sites at the bottom layer $`n_0`$ can be interpreted as the density of active sites in a DP process. Therefore, at criticality this density is expected to decay as

$$n_0t^\delta ,$$
(6)
where $`\delta =\beta /\nu _{}0.1595`$ is one of the dynamic critical exponents of DP in 1+1 dimensions. Similarly the critical behavior at the first few layers above the bottom layer can be described in terms of unidirectionally coupled DP processes Coupled . Roughly speaking each layer can be associated with a DP process which is coupled to the layer below. Effective couplings in opposite direction, which certainly exist in the model defined above, turn out to be irrelevant. Field-theoretic renormalization group studies revealed that the scaling exponents $`\nu _{}`$ and $`\nu _{}`$ (and hence their ratio $`z=\nu _{}/\nu _{}`$) are the same at all levels, while the exponents associated with the order parameters decrease with increasing height above the bottom layer. Thus the concept of unidirectionally coupled DP successfully explains the critical behavior at the first few layers above the bottom layer. However, it cannot predict the critical properties of the interface as a whole, especially the scaling of the interface width.

## III Numerical simulations

### III.1 Scaling properties at the critical point

Performing high-precision simulations on a parallel computer we first used Eq. (6) to estimate the critical threshold. As shown in Fig. 2 we obtain the result

$$q_c=1.8861(1)$$
(7)
which is much more accurate than the previous estimate reported in Ref. Alon . Moreover, our simulations confirm that the dynamics at the bottom layer is indeed driven by a DP process.

At the critical point $`q=q_c`$ the average height

$$H=\frac{1}{L}\underset{i=1}{\overset{L}{}}h_i$$
(8)
and the squared interface width

$$W^2=\frac{1}{L}\underset{i=1}{\overset{L}{}}(h_iH)^2$$
(9)
are found to increase logarithmically with time as

$`W^2(t)`$ $``$ $`\tau \mathrm{ln}t,`$ (10)
$`H(t)`$ $``$ $`\sigma \mathrm{ln}t,`$ (11)
where $`\tau =0.102(3)`$ and $`\sigma =0.133(3)`$ (see Fig. 3). Thus, apart from different prefactors, the squared width and the mean height show the same type of logarithmic behavior. We note that in a previous work the width was erroneously conjectured to increase as $`W(t)(\mathrm{ln}t)^{0.43}`$ Alon2 .

### III.2 Finite-size simulations

In a finite system the critical behavior depends on the type of the boundary condition. In what follows we consider (a) periodic and (b) Dirichlet boundary conditions, where the sites at the boundary are fixed at zero height. As shown in Fig. 4, the average height $`H(L,t)`$ and the squared width $`W^2(L,t)`$ first increase logarithmically until finite-size effects become relevant and the system crosses over to a different regime, where the width saturates. Measuring the saturation levels we find that

$$W_{\mathrm{sat}}^2(L)A_{p,f}+\lambda _{p,f}\mathrm{ln}L,$$
(12)
where the indices $`p,f`$ stand for periodic and fixed boundary conditions, respectively. Our best estimates are:

$$\begin{array}{cc}\hfill A_p& =0.16(1),\lambda _p=0.161(2)\hfill \\ \hfill A_f& =0.21(1),\lambda _f=0.163(2)\hfill \end{array}$$
(13)
Apparently the amplitude $`\lambda `$ is universal in the sense that it does not depend on the choice of the boundary condition.

Unlike the width, the interface height does only saturate if fixed boundary conditions are used. Here we find a similar formula

$$H_{\mathrm{sat}}(L)B_f+\mu _f\mathrm{ln}L,$$
(14)
where

$$B_f=0.24(1),\mu _f=0.212(10).$$
(15)
In order to collapse the curves shown in Fig. 4, recall that the roughening transition in this model is driven by a 1+1-dimensional DP process at the bottom layer, which is characterized by a dynamic exponent

$$z=\nu _{}/\nu _{}1.5807.$$
(16)
Moreover it has been shown that the critical dynamics at the first few layers is the same as in a sequence of unidirectionally coupled DP processes, which are all characterized by the same dynamic exponent $`z`$. Therefore, it is reasonable to assume that the roughening transition as a whole is characterized by the dynamic exponent of DP so that finite-size scaling functions can only depend on the scaling-invariant ratio $`t/L^z`$. Thus, the expected scaling forms read

$$\begin{array}{c}\hfill W^2(L,t)\lambda \mathrm{ln}L+f(t/L^z),\\ \hfill H(L,t)\mu \mathrm{ln}L+g(t/L^z),\end{array}$$
(17)
where $`f(\xi )`$ and $`g(\xi )`$ are scaling functions with the asymptotic behavior

$$f(\xi )=\frac{\lambda }{z}\mathrm{ln}\xi ,g(\xi )=\frac{\mu }{z}\mathrm{ln}\xi $$
(18)
for $`\xi 0`$ and $`f(\xi )=A`$ and $`g(\xi )=B`$ for $`\xi \mathrm{}`$. In order to verify these scaling forms we plot $`W^2(L,t)\lambda \mathrm{ln}L`$ and $`H(L,t)\mu \mathrm{ln}L`$ versus $`t/L^z`$. As shown in Fig. 5, we obtain excellent data collapses.

### III.3 Off-critical simulations

A similar situation is encountered in off-critical simulations. In the subcritical regime $`q<q_c`$ the width and the average height saturate at

$$\begin{array}{cc}\hfill W_{\mathrm{sat}}^2(ϵ)& =C\eta \mathrm{ln}|ϵ|,\hfill \\ \hfill H_{\mathrm{sat}}(ϵ)& =D\varphi \mathrm{ln}|ϵ|,\hfill \end{array}$$
(19)
where $`ϵ=qq_c`$ denotes the distance from criticality. Following the same arguments as in the previous subsection, the conjectured scaling forms read

$$\begin{array}{cc}\hfill W^2(ϵ,t)& =\eta \mathrm{ln}|ϵ|+F(t|ϵ|^\nu _{}),\hfill \\ \hfill H(ϵ,t)& =\varphi \mathrm{ln}|ϵ|+G(t|ϵ|^\nu _{}),\hfill \end{array}$$
(20)
where $`\nu _{}1.7338`$ is the temporal scaling exponent of DP in 1+1 dimensions. These scaling forms are expected to hold not only below but also above criticality. In fact, plotting $`W^2(ϵ,t)+\eta \mathrm{ln}|ϵ|`$ and $`H(ϵ,t)+\varphi \mathrm{ln}|ϵ|`$versus $`t|ϵ|^\nu _{}`$ (see Fig. 6), the best data collapses are obtained for

$$\eta =0.173(10),\varphi =0.230(10),$$
(21)
while the constants are estimated by

$$C=0.40(1),D=0.54(1).$$
(22)
### III.4 Relations among the amplitudes

Obviously the scaling forms (17) and (19) suggest that the exponentiated quantities $`\mathrm{exp}(W^2)`$ and $`\mathrm{exp}(H)`$ display ordinary power-law scaling. This would imply that the results for finite-size and off-critical scaling can be combined in a single scaling form, e.g.,

$$\begin{array}{cc}\hfill \mathrm{exp}\left[W^2(ϵ,L,t)\right]& =L^\lambda \stackrel{~}{F}(ϵL^{1/\nu _{}},t/L^z),\hfill \\ \hfill \mathrm{exp}\left[H(ϵ,L,t)\right]& =L^\mu \stackrel{~}{G}(ϵL^{1/\nu _{}},t/L^z).\hfill \end{array}$$
(23)
This in turn would imply that the six amplitudes $`\lambda ,\mu ,\tau ,\sigma ,\eta ,\varphi `$ introduced above, which now play the role of critical exponents, are no longer independent. Only two of them, say $`\eta `$ and $`\varphi `$, are independent while the others can be expressed as

$$\begin{array}{cc}\hfill \lambda & =\eta /\nu _{},\hfill \\ \hfill \mu & =\varphi /\nu _{},\hfill \\ \hfill \tau & =\eta /\nu _{},\hfill \\ \hfill \sigma & =\varphi /\nu _{},\hfill \end{array}$$
(24)
where $`\nu _{}1.0968`$ and $`\nu _{}1.7338`$. As shown in Table 1 these relations are compatible with the numerical results.

## IV Simple approximation

Why does the interface roughen logarithmically at the transition? To answer this question we present a simple approximation for the average interface height $`H(L)`$ in a finite 1+1-dimensional system at criticality. We consider fixed boundary conditions, i.e., the interface is pinned to the bottom layer at the two boundary sites. In between the interface may return to the bottom layer several times, dividing the chain into intervals of different lengths $`\mathrm{}_1,\mathrm{}_2,\mathrm{}`$, as sketched in Fig. 7. The approximation relies on the following assumptions:

1. The probability distribution for the interval sizes $`\mathrm{}`$ is the same as the size distribution of inactive islands in a DP process, which is known to decrease as $`P(\mathrm{})\mathrm{}^{\beta /\nu _{}2}\mathrm{}^{0.25}`$.
2. As we are using fixed boundary conditions, we will assume that this probability distribution is cut off at the system size, i.e., $`P(\mathrm{})=0`$ for $`\mathrm{}L`$.
3. The sites, where the interface touches the bottom layer, can be regarded as fixed boundary conditions for the dynamics taking place at the following layer.

While the first two assumptions can be justified by numerical checks, the third assumption is crucial since the sites at the bottom layer evolve in time and can only be regarded as fixed in an approximate sense. However, using this approximation we can derive a recursion relation since each island (for example the island encircled by the dotted line in Fig. 7) can be regarded as an independent growth process taking place in a finite system of size $`\mathrm{}`$ with fixed boundary conditions at the following layer. Obviously, the average height of such an island will be $`1+H(\mathrm{})`$. In order to calculate the average height of the entire interface $`H(L)`$, one has to add up all these contributions weighted by the island size $`\mathrm{}`$. Replacing sums by integrals, this leads to the integral equation

$$H(L)=\frac{_0^L𝑑\mathrm{}\mathrm{}P(\mathrm{})\left[1+H(\mathrm{})\right]}{_0^L𝑑\mathrm{}\mathrm{}P(\mathrm{})}.$$
(25)
Shifting the denominator to the l.h.s. and differentiating with respect to $`L`$ we obtain

$$\frac{dH(L)}{dL}_0^L𝑑\mathrm{}\mathrm{}P(\mathrm{})=LP(L),$$
(26)
or equivalently

$$\frac{dH(L)}{dL}=\frac{d}{dL}\mathrm{ln}_0^L𝑑\mathrm{}\mathrm{}P(\mathrm{}),$$
(27)
leading to the solution

$$\begin{array}{cc}\hfill H(L)& =const+\mathrm{ln}_0^L𝑑\mathrm{}\mathrm{}P(\mathrm{})\hfill \\ & =const+\frac{\beta }{\nu _{}}\mathrm{ln}L.\hfill \end{array}$$
(28)
Thus the approximation correctly reproduces the logarithmic increase of the average height in a finite critical system. However, the predicted amplitude $`\mu =\beta /\nu _{}0.25`$ differs from the measured value $`\mu 0.21`$.

A similar argument can be used to explain the logarithmic increase of the squared width. To this end we first calculate the second moment of the height

$$M_2(L)=\frac{1}{L}\underset{i=1}{\overset{L}{}}h_i^2.$$
(29)
Using the same assumptions as before, each island of size $`\mathrm{}`$ contributes to $`M_2(L)`$ with a term

$$M_2(\mathrm{})|_{hh+1}=1+2H(\mathrm{})+M_2(\mathrm{}).$$
(30)
Replacing sums by integrals we obtain the integral equation

$$M_2(L)=\frac{_0^L𝑑\mathrm{}\mathrm{}P(\mathrm{})\left[1+2H(\mathrm{})+M_2(\mathrm{})\right]}{_0^L𝑑\mathrm{}\mathrm{}P(\mathrm{})},$$
(31)
leading to the differential equation

$$\frac{dM_2(L)}{dL}_0^L𝑑\mathrm{}\mathrm{}P(\mathrm{})=LP(L)+2LP(L)H(L).$$
(32)
In order to remove the mixed term on the r.h.s. we differentiate the squared width $`W^2(L)=M_2(L)H^2(L)`$ with respect to $`L`$ and combine the result with Eq. (26):

$$\frac{d[W^2(L)M_2(L)]}{dL}_0^L𝑑\mathrm{}\mathrm{}P(\mathrm{})=2LP(L)H(L).$$
(33)
Inserting Eq. (32) the mixed term drops out and we obtain the differential equation

$$\frac{dW^2(L)}{dL}_0^L𝑑\mathrm{}\mathrm{}P(\mathrm{})=LP(L),$$
(34)
which has exactly the same form as Eq. (26). Therefore, $`W^2(L)`$ and $`H(L)`$ will only differ by a constant. Again the approximation correctly reproduces the logarithmic increase but cannot predict the precise value of the amplitude.

## V Conclusions

In the present paper we have analyzed the scaling behavior at the roughening transition of a restricted solid-on-solid growth process with evaporation at the edges of terraces. Restricting to the case of a one-dimensional substrate we have carried out extensive Monte Carlo simulations, exceeding the temporal range of previous numerical studies by three decades. The numerical results confirm that the transition is indeed driven by a directed percolation process at the bottom layer. Furthermore, they clearly show that squared interface width $`W^2(t)`$ and the average height above the bottom layer $`H(t)`$ increase logarithmically with time. Analyzing numerical data from finite-size and off-critical simulations we have postulated appropriate scaling forms which generate excellent data collapses. These scaling forms can be interpreted in such a way that the exponentiated quantities $`\mathrm{exp}(W^2)`$ and $`\mathrm{exp}(H)`$ obey ordinary power-law scaling, which allows one to derive various relations among the amplitudes. Moreover, we have presented a simple approximation which explains why $`W^2`$ and $`H`$ grow logarithmically.

In contrast to the expectations expressed in Alon2 , we cannot find evidence for several competing length scales. Moreover, there is no evidence for logarithmic scaling based on local scale invariance Tang , rather the system shows ordinary scaling after an appropriate redefinition (exponentiation) of the order parameters.

Our results suggest that any roughening transition, which is driven by a phase transition from an active into an absorbing state at the bottom layer, should exhibit this type of logarithmic roughening. An important example is a recently introduced model for dimer adsorption and desorption Dimer , where the transition is driven by a parity-conserving transition at the bottom layer.

Acknowledgements: The simulations were partly performed on the ALiCE parallel computer at the IAI in Wuppertal. I would like to thank B. Orth and G. Arnold for technical support.


text:
# Static hybrid quarkonium potential with improved staggered quarks

## 1 INTRODUCTION

Simulations with dynamical quarks have found that light quarks modify the heavy quark-antiquark potential in a number of ways . At large distances they decrease the string tension in units of the Sommer $`r_0`$ and $`r_1`$ parameters (defined by $`r^2F(r)=1.65`$ and 1.00, respectively) and lead eventually to string-breaking. At shorter distances they modify the running of the coupling constant, deepening the Coulomb well and increasing the ratio $`r_0/r_1`$. In this work, we extend these studies to some of the potentials with excited flux tubes. Of particular interest to quarkonium spectroscopy are the $`\mathrm{\Pi }_u`$ excitations leading to exotic $`Q\overline{Q}g`$ hybrids .

We report results of a study in which our sources and sinks are optimized to create and annihilate a flux-tube state. In the presence of dynamical quarks, string breaking is expected. It is known that in the conventional $`\mathrm{\Sigma }_g^+`$ channel, transitions to the open two-meson channel are exceedingly weak, qualitatively consistent with the small widths of quarkonium states above the heavy-light meson thresholds . Since at present we do not include the open two-meson channel we do not expect to observe string breaking here.

## 2 MEASUREMENTS

We have measured the heavy quark potential on an ensemble of $`28^3\times 96`$ ($`a0.09`$ fm) gauge configurations generated in the presence of $`2+1`$ flavors of Asqtad dynamical quarks of varying masses and a one-loop Symanzik gauge action. The strange quark mass is set approximately to its physical value. Here we compare results from our 358-configuration quenched ensemble with our 495-configuration dynamical quark ensemble for which $`(m_\pi r_0)^21.3`$.

The configurations are first smoothed using a single hypercubic (HYP) blocking pass , a technique that improves significantly the signal-to-noise ratio . The blocking procedure involves replacing all gauge links (timelike as well as spacelike) with an SU(3)-projected average over paths confined to adjacent hypercubes. Thus distortions in the result are local and expected to be confined to distances smaller than about $`2a`$ . After HYP blocking the spacelike links are further smoothed via five cycles of APE smearing with SU(3) projection.

On the thus smoothed lattices we measure the expectation value of the standard $`R\times T`$ Wilson loop on axis and along three different off-axis directions. These measurements yield the conventional ground state $`\mathrm{\Sigma }_g^+`$ and excited state $`\mathrm{\Sigma }_g^+`$ potentials. For the $`\mathrm{\Pi }_u`$ excited state, we measured the expectation value of a bent loop formed by replacing the source and sink flux tubes of length $`R`$ by a superposition of large “staples” of sides $`(2a,R,2a)`$. For example one such loop replaces each on-axis spacelike flux path ($`R\widehat{x}`$) by paths of the form ($`2a\widehat{y}`$, $`R\widehat{x}`$, $`2a\widehat{y}`$) minus its reflection in the $`xz`$ plane .

For the standard Wilson loop we extracted the usual $`\mathrm{\Sigma }_g^+`$ potential $`V_{\mathrm{\Sigma }g+}`$ and its excited state $`V_{\mathrm{\Sigma }g+}^{}`$ by doing a blocked, correlated, double-exponential fit to the Wilson loop data:

$`W(R,T)`$ $`=`$ $`C_{\mathrm{\Sigma }g+}(R)e^{V_{\mathrm{\Sigma }g+}(R)T}`$ (1)
$`+`$ $`C_{\mathrm{\Sigma }g+}^{}(R)e^{V_{\mathrm{\Sigma }g+}^{}(R)T}.`$ (2)
For the $`\mathrm{\Pi }_u`$ potential we did only a single-exponential fit.

In all cases we use the same fit ranges for both quenched and dynamical lattices to reduce possible systematic errors.

## 3 RESULTS

In Fig. 1 we compare the ground state potential on the quenched ensemble and the 2+1 flavor ensemble. Both the distance scale and the potential are plotted in units of $`r_0`$, and a constant has been subtracted from the potential so that it is zero at $`r_0`$. Since $`r_0`$ was determined from this potential, the fits are tangent at this point. Away from $`r_0`$, the potentials have different shapes, namely, the Coulomb attraction is slightly stronger in the light quark ensemble and the string tension is slightly weaker in units of $`r_0`$, confirming earlier findings . A softening of the Coulomb well is also evident. This is an expected consequence of HYP smoothing.

Similarly in Fig. 2 we show light quark effects in the $`\mathrm{\Sigma }_g^{}`$ excitation potential. The potentials are plotted relative to the zero determined in the fit to the ground state potentials. From Fig. 2 we see that the excited state potential $`\mathrm{\Sigma }_g^{}`$ is slightly steeper than that of quenched QCD.

In Fig. 3 we plot the $`\mathrm{\Pi }_u^{}`$ potential. This hybrid potential is weakly repulsive at short range, as would be expected from the Coulomb interaction in a color octet quark-antiquark system. This effect is softened by HYP smoothing. At the level of our statistical errors there are no apparent differences at long range, but better statistics would certainly be of interest.

## 4 CONCLUSIONS

Our measurements at $`a=0.09`$ fm confirm the shape changes in the ground-state potential, seen previously at $`a=0.13`$ fm. In units of $`r_0`$ we find, further, that adding $`2+1`$ flavors of dynamical quarks makes the $`\mathrm{\Sigma }_g^+`$ excited state potential slightly steeper and the $`\mathrm{\Pi }_u^{}`$ slightly more repulsive at short range. We find no clear evidence for a flattening of the potentials that would signal string breaking.

Computations were performed at LANL, NERSC, NCSA, ORNL, PSC, SDSC, FNAL, and the CHPC (Utah). This work is supported by the U.S. NSF and DOE.