text:
# Dynamics and nonequilibrium states in the Hamiltonian mean-field model: A closer look

## Abstract

We critically revisit the evidence for the existence of quasistationary states in the globally coupled XY (or Hamiltonian mean-field) model. A slow-relaxation regime at long times is clearly revealed by numerical realizations of the model, but no traces of quasistationarity are found during the earlier stages of the evolution. We point out the nonergodic properties of this system in the short-time range, which makes a standard statistical description unsuitable. New aspects of the evolution during the nonergodic regime, and of the energy distribution function in the final approach to equilibrium, are disclosed.

The globally coupled version of the XY model, given by the Hamiltonian

$$=K+V=\frac{1}{2}\underset{i=1}{\overset{N}{}}p_i^2+\frac{1}{2N}\underset{i,j=1}{\overset{N}{}}[1\mathrm{cos}(\theta _i\theta _j)],$$
(1)
has recently attracted considerable attention RuffoPRE95 ; RuffoPRL98 ; RuffoPHD99 ; RuffoPRL99 ; TamPRL00 ; TsPRE01 ; TsPHA02 ; RuffoCSF02 ; mta . This model, usually referred to as the Hamiltonian mean-field (HMF) model, describes an ensemble of $`N`$ rotators characterized by their angles $`\theta _i`$ and conjugate momenta $`p_i`$. In contrast to the usual XY model, interactions affect all pairs of rotators in the ensemble. The HMF model has been studied in connection with the emergence of collective self-organized states RuffoPRE95 , dynamical behavior near phase transitions RuffoPRL98 , and anomalous diffusion in phase space RuffoPRL99 , among other dynamical and thermodynamical aspects.

The canonical thermodynamical description of the HMF model can be completely carried out RuffoPRE95 ; RuffoPHD99 . In this description, it is useful to introduce the “magnetization” $`𝐌=N^1_i(\mathrm{cos}\theta _i,\mathrm{sin}\theta _i)(M_x,M_y)M\mathrm{exp}(\mathrm{i}\mathrm{\Phi })`$ and the average energy per particle $`U=E/N`$, with $`E`$. The canonical approach predicts that the values of the magnetization and the average energy at equilibrium, $`M_{\mathrm{eq}}`$ and $`U_{\mathrm{eq}}`$, are related by

$$U_{\mathrm{eq}}=\frac{T}{2}+\frac{1}{2}(1M_{\mathrm{eq}}^2),$$
(2)
where $`T`$ is the canonical temperature. The system undergoes a ferromagnetic-like, second-order phase transition at $`T_c=1/2`$, from a state with $`M_{\mathrm{eq}}0`$ at low temperatures to a state of vanishing magnetization at high temperatures. The energy at the transition is $`U_c=3/4`$. These predictions are mostly confirmed by the numerical solution of the HMF equations of motion at fixed total energy, i.e. in a microcanonical scenario. In this situation the equilibrium temperature is defined through the kinetic energy

$$T_{\mathrm{eq}}=\frac{2}{N}K_{\mathrm{eq}},$$
(3)
and $`U=T_{\mathrm{eq}}/2+(1M_{\mathrm{eq}}^2)/2`$ \[cf. Eq. (2)\].

Deviations between the canonical predictions and the microcanonical (numerical) results were however reported in the region of energies just below the critical point, $`0.5<U<U_c`$, where numerical realizations reveal extremely slow relaxation toward the asymptotic state RuffoPRE95 ; RuffoPRL98 . Such deviations, observed when the ensemble of rotators is prepared in specially chosen initial conditions—namely, the so called water-bag initial conditions—have been associated with the existence of long-lived states where, in the thermodynamical limit $`N\mathrm{}`$, the system would spend asymptotically long times. In that limit, canonical equilibrium would never be reached. It has been suggested that the HMF long-lived states—which would replace canonical equilibrium in the thermodynamical limit—would be well described by the equilibrium distributions predicted by Tsallis’s generalization of thermodynamics TsPRE01 ; TsPHA02 ; TsPHA01 . This generalization, in fact, relaxes the assumption of extensivity of the Boltzmann-Gibbs formulation, and is therefore presumably expected to describe the equilibrium statistics of systems with long-range interactions TsPHA95 , as in the Hamiltonian (1).

Water-bag (WB) initial conditions fix $`\theta _i=0`$ for all $`i`$, and the momenta $`p_i`$ are chosen at random from a uniform distribution in such a way that $`_ip_i=0`$ and $`K=_ip_i^2/2=NU`$. For these initial conditions, in fact, the potential energy vanishes. For a typical numerical realization at energies just below the critical point $`U_c`$, the dynamics of the HMF model can be summarized as follows. In the first stage, whose duration is essentially independent of the system size $`N`$ and equals a few time units TsPHA02 , there is a rapid broadening of the distribution in $`\theta `$. This implies an abrupt growth of the potential energy from its initial value $`V(t=0)=0`$, to a value close to that predicted at canonical equilibrium, $`V_{\mathrm{eq}}=N(1M_{\mathrm{eq}}^2)/2`$. The kinetic energy drops accordingly, from $`K(t=0)=NU`$ to a value close to $`K_{\mathrm{eq}}=NT_{\mathrm{eq}}/2`$, and the initial unbalance between $`K`$ and $`V`$ is drastically reduced. As explained below, the details of the evolution in this stage of kinetic relaxation depend strongly on the specific realization of the WB initial condition. Typically, however, the potential energy exceeds $`V_{\mathrm{eq}}`$ and attains a maximum, while $`K`$ reaches a minimum. From then on, the system enters a stage of slow relaxation toward canonical equilibrium.

The slow-relaxation regime of the HMF model was reported in the first systematic study of this system RuffoPRE95 . In Ref. RuffoPHD99,  it was quantitatively characterized by analyzing the relaxation of the nonequilibrium entropy $`S(t)`$ toward its equilibrium value $`S_{\mathrm{eq}}`$ as a function of $`N`$. The difference $`|S(t)S_{\mathrm{eq}}|`$ was shown to approach an asymptotic value of order $`N^{1/2}`$ within a relaxation time $`\tau _rN`$. The same scaling has recently been reported for the relaxation time of the nonequilibrium temperature, $`T(t)=2K(t)/N`$ RuffoCSF02 . The states visited by the system during the slow-relaxation regime have been named quasistationary states, presumably by analogy with the quasiequilibria reported for one-dimensional gravitational systems grav1 ; grav2  and similar long-lived states in many-body systems with long-range interactions Thirring . The name suggests that the time scales associated with the evolution of such states are much longer than those of any other regime during the whole evolution. However, the evolution of $`S(t)`$ toward $`S_{\mathrm{eq}}`$ reveals no dynamical regime other than slow relaxation RuffoPHD99 ; init .

Latora, Rapisarda and Tsallis (LRT), on the other hand, have claimed to demonstrate the entity of quasitationary states (QSS) from the analysis of the evolution of the nonequilibrium temperature $`T(t)`$ TsPRE01 . For clarity in the discussion, we reproduce in the main plot of Fig. 1 their figure 1(b) of Ref. TsPRE01, , where $`T(t)`$ is plotted against time for several values of $`N`$, as obtained from numerical integration of the HMF equations of motion for $`U=0.69`$ LA . This value of the energy corresponds to the largest observed deviation of numerical results with respect to the canonical equilibrium prediction RuffoPRE95 ; RuffoPHD99 . In the plot, the initial stage of kinetic relaxation is disregarded. For short times, we observe a wide plateau where $`T(t)`$ is apparently constant, but substantially differs from its equilibrium value. At long times, $`T(t)`$ has reached $`T_{\mathrm{eq}}`$ and remains at equilibrium. These two plateaus are connected by a crossover stage where $`T(t)`$ exhibits an inflexion point and its variation is seemingly faster. Consequently, LRT conclude that the slow-relaxation regime consists in turn of a first stage where the system is “trapped” in a QSS, and a subsequent crossover toward canonical equilibrium. Since the position $`\tau _C`$ of the crossover shifts to larger times as $`N`$ grows, $`\tau _CN`$, LRT deduce that in the thermodynamical limit the QSS will be observed at all times, and that canonical equilibrium will never be attained. As we argue in the following, however, the presence of the plateaus and, thus, the existence of a QSS in these specific realizations of the HMF model, is an artifact of the peculiar way chosen to plot the data in Fig. 1.

In the first place, in order to exclude the initial stage of kinetic relaxation, LRT choose to plot $`T(t)`$ as a function of $`tt_0`$, with $`t_0=100`$, instead of simply cutting off the time axis at $`t_0`$. Such procedure, which would be innocuous in a linear time scale, has far-reaching consequences in the logarithmic scale of Fig. 1. In fact, the point $`t=t_0`$ becomes shifted to $`\mathrm{}`$ in the logarithmic variable, which considerably enhances the impression of having a plateau at small $`t`$. In the inset of Fig. 1 we plot $`T(t)`$ against its original variable—also in logarithmic scale—showing that the short-time plateau shrinks by one decade. The existence of a crossover state joining the two plateaus is itself an illusion created by the logarithmic scale, even in the original variable $`t`$. Indeed, it can be easily shown that the graph of a function with monotonic first derivative can display an (otherwise inexistent) inflection point by the simple expedient of using a logarithmic scale in the horizontal axis log .

In order to avoid these undesirable effects and to reveal the true nature of the evolution of $`T(t)`$ in the slow-relaxation regime, we opt for using linear time scales. Figure 2 shows, plotted with diamonds, the same data as in Fig. 1 for $`N=500`$ up to $`t=2500`$, while the inset shows the whole temporal range. We see that there are no traces of a short-time plateau, except perhaps for a small interval around $`t=500`$. It might however be argued that the vertical scale in this plot has been exaggeratedly amplified. To decide over this point, we have used a Padé-like approximation to fit the data over the whole time range, shown as a curve in the inset. The quality of the fitting is quite acceptable. Its extension to the short-time range, also shown as a curve in the main plot, proves that in this range the data exhibit the same trend as for larger times. The evolution in the whole range is therefore uniform, and we find no arguments to attribute the condition of quasistationarity to the short-time interval. As for the small interval around $`t=500`$, we also include in Fig. 2 a data set from Ref. TsPHA02, , calculated with higher numerical precision (open circles) LA2 . These data reveal a minimum in $`T(t)`$ at $`t300`$ and confirm the absence of any plateau. Figure 3 shows the data from Fig. 1 for $`N=10^3`$ and $`10^4`$ in linear time scale. Again, we find no evidence of the existence of plateaus or QSS.

Let us mention that, more indirectly, the observation of anomalous diffusion in the slow-relaxation stage of the HMF model RuffoPRL99  confirms that no special regimes exist within that stage. In fact, superdiffusion in the unfolded $`\theta `$-space, $`\theta (\mathrm{},+\mathrm{})`$, which is characterized by an anomalous dispersion law $`\mathrm{\Delta }\theta ^2t^\alpha `$, is recorded in the whole slow-relaxation stage with a constant exponent $`\alpha 1.4`$. The transition to ordinary diffusion ($`\alpha =1`$) occurs at the same time as the system reaches canonical equilibrium. The short-time plateau of Fig. 1 and the crossover to canonical equilibrium are not detected by the transport properties of this system.

The data shown in Figs. 2 and 3 exhibit an inconsistency regarding the behavior of the nonequilibrium temperature in the first stage of the slow-relaxation regime. In fact, while in the results for $`N=500`$ (circles) and for $`N=10^4`$ we observe a well-defined minimum in $`T(t)`$, low-precision data for $`N=500`$ (diamonds) and results for $`N=10^3`$ do not display such minimum. To clarify this aspect of the evolution, we have performed extensive numerical calculations for the same value of the energy, $`U=0.69`$, maintaining the relative error in the total energy conservation around $`\mathrm{\Delta }E/E10^3`$, which makes it possible to average over many realizations of the water-bag initial conditions ($`5\times 10^4`$ for $`N=500`$ to $`10^3`$ for $`N=3000`$). Our results for the nonequilibrium temperature $`T(t)`$ are shown in Fig. 4, as a function of the rescaled time $`t/N`$. Each dot has been obtained as an average in time, over $`0.1N`$ to $`0.2N`$ time units. The arrows indicate, for $`N=1000`$ and $`3000`$, the typical dispersion in the values of $`T(t)`$ after averaging over the whole set of realizations. Fluctuations are indeed large, as expected near a phase transition and due to the highly chaotic HMF dynamics RuffoPRL98 ; RuffoPHD99 , which makes it difficult to obtain reliable numerical results. In any case, our numerical calculations consistently confirm the existence of a minimum in the nonequilibrium temperature, and the subsequent steady growth of $`T(t)`$. The inset shows the position of the minimum, $`t_{\mathrm{min}}`$, as a function of $`N`$, including also the value for $`N=10^4`$ taken from the data of Fig. 3(b). The straight line corresponds to a power-law fitting, $`t_{\mathrm{min}}N^\nu `$, which yields $`\nu =1.7\pm 0.1`$. Though the plot spans less than two decades in $`N`$, the fact that $`t_{\mathrm{min}}`$ grows faster than linearly with the system size is clear.

Meanwhile, no plateau is observed, at least, of the kind shown in Fig. 1. A plateau may of course be compulsively ascribed to the zone around the minimum of $`T(t)`$ at $`t_{\mathrm{min}}`$, but it would be far from coinciding in nature with the QSS reported by LRT. In fact, although the minimum becomes broader for larger $`N`$, it also shifts to longer times, which hardly insures that a stationary (nonequilibrium) temperature would be observed in the limit $`N\mathrm{}`$.

The well-defined minimum of the nonequilibrium temperature $`T(t)`$ at $`t_{\mathrm{min}}`$ establishes a quite natural boundary between short-time and long-time dynamical ranges. It is interesting to analyze in more detail the short-time evolution, with emphasis in the angle and momentum distribution of the ensemble in $`\mu `$-space, where at each time the system is represented by a set of $`N`$ points in the $`(\theta ,p)`$ plane. One-particle distributions of both angles and momenta have received particular attention in previous work RuffoPRE95 ; RuffoPHD99 ; TsPRE01 ; TsPHA02 , so that such representation is useful for comparison. Figure 5 shows four snapshots, at different times, of the $`\mu `$-space distribution in a single (typical) realization of the HMF model with $`N=10^4`$ rotators, at $`U=0.69`$. For clarity, the origin of the angle variable, $`\theta =0`$, has been chosen to coincide in each plot with the phase $`\mathrm{\Phi }`$ of the total magnetization $`𝐌`$. In this way, the center of the distribution coincides approximately with the center of the plot. Take into account, however, that $`\mathrm{\Phi }`$ performs, for short times, diffusive-like motion RuffoPRE95 . The WB initial condition would here correspond to a uniform distribution along a vertical segment at $`\theta =0`$, with $`|p|<\sqrt{6U}2.03`$. As advanced above, we observe a rapid broadening of the distribution in the angle coordinate $`\theta `$, accompanied by a moderate collapse in $`p`$. At $`t=100`$, much of the ensemble has become strongly mixed in $`\mu `$-space, due to the highly chaotic nature of the dynamics for this value of $`U`$ RuffoPRE95 ; RuffoPHD99 . There is however a substantial fraction of the ensemble, represented by the zones where points remain ordered, where the strong correlations imposed by the WB initial condition persist. This fraction is mainly situated near the center of the distribution, in a region of relatively small momentum, but ordered arrays of points are also seen in other zones of $`\mu `$-space. For these times, thus, the system can be thought of as a mixture of two phases, namely, a strongly mixed, “gaseous” phase and a highly correlated, “condensed” phase.

The distribution of the condensed phase in $`\mu `$-space results to strongly depend on the specific realization of the WB initial condition. This sensibility to the initial conditions, illustrated in Fig. 6, is again a consequence of the chaotic dynamics that, as time elapses, amplifies the originally small variations between different realizations. An important byproduct of this property, combined with the persistence of the initial correlations, is that the condensed phase is not ergodic, in the sense that averages over realizations of the WB initial condition yield a poor description of the statistical properties over time of any single realization. The lack of ergodicity of the condensed phase makes in turn the whole system nonergodic in this short-time range. An illustration of this overall lack of ergodicity is provided by the short-time evolution of the nonequilibrium temperature. In any single realization, after its initial drop, $`T(t)`$ displays chaotic oscillations of appreciable amplitude. Due to the sensibility to initial conditions, however, the phases of such oscillations differ between realizations, and the oscillations disappear upon averaging. This observation should warn of the averaging procedures used in previous work TsPHA02  to study statistical properties such as the nonequilibrium temperature and the momentum distribution function in the initial regime of kinetic relaxation. The warning applies also to the average data of $`T(t)`$ shown in our Fig. 4 at the shortest times. Before the nonergodic condensate disappears, in fact, a (standard) statistical approach to the HMF model results to be of limited applicability.

As the evolution proceeds, particles from the condensate “evaporate” into the gaseous phase RuffoPRE95 ; RuffoPRL98 . From Fig. 5 we note that, at $`t=10^4`$, only faint traces of the condensed phase remain in the ensemble. We have studied the dynamics of this evaporation process as follows. First of all, we have numbered the rotators in each realization of the HMF model in such a way that, at the initial time, $`p_i(0)<p_j(0)`$ for $`i<j`$. At subsequent times, we have defined the condensate as formed by the rotators $`i`$ such that the distance to $`i+1`$ in $`\mu `$-space, $`\delta _i=\sqrt{(\theta _{i+1}\theta _i)^2+(p_{i+1}p_i)^2}`$, is lower than a given threshold $`\delta `$. We fix the threshold as the average distance between contiguous pairs at $`t=0`$, $`\delta =2\sqrt{6U}/N`$.

The number $`n`$ of rotators in the condensed phase displays considerable oscillations in time. According to the discussion in the previous paragraph, we choose to measure $`n`$ in a single (typical) realization for $`N=10^4`$ and to average over time intervals. Figure 7(a) shows the evolution of $`n`$ averaged each $`100`$ time units. The condensate size decreases steadily, displaying most of the time a power-law decay with a nontrivial exponent, $`nt^{0.85}`$. This rather smooth decay is abruptly interrupted by a sudden collapse of the condensate, at $`t1.6\times 10^4`$. In Fig. 7(b), we show the evolution of the average energy per particle in the condensate. We see that it varies nonmonotonically with time, though its value is never far from $`u=1`$. Other realizations for the same value of $`N`$ confirm that these results are generic. Moreover, they are robust with respect to the definition of the threshold $`\delta `$. For different values of $`N`$, realizations in the short-time range show that the collapse of the condensed phase occurs at a time preceding, but approximately equal to, $`t_{\mathrm{min}}`$. In other words, when the minimum of the nonequilibrium temperature is attained, the condensate has just evaporated. From then on, the system enters the final stage of its relaxation to canonical equilibrium.

In the long-time range, when the ensemble of rotators has presumably attained an ergodic state, its $`\mu `$-space configuration can be described statistically, in terms of averages over realizations, by means of a distribution function $`f(\theta ,p)`$. Previous results on the $`\mu `$-space distribution of the HMF model refer, in all cases, to reduced distribution functions, $`f_\theta (\theta )`$ and $`f_p(p)`$. In Ref. RuffoPHD99, , the momentum distribution function is used to illustrate, for $`U=0.69`$, slow relaxation for long times. For other values of $`U`$, far from the critical point, both $`f_\theta `$ and $`f_p`$ were compared at moderately long times with the prediction of canonical statistics, showing the relatively fast attainment of canonical equilibrium for such energies. In Refs. TsPRE01,  and TsPHA02, , in contrast, deeper significance is attributed to nonequilibrium distribution functions—always in connection with QSS. There, LRT suggest that $`f_p`$ exhibits the fingerprints of nonextensive thermodynamics, and support this claim showing that the momentum distribution can be acceptably fitted, at low values of $`p`$, by the power-law distribution derived from Tsallis’s nonextensive formulation TsPHA95 . Tsallis’s formulation, however, predicts the one-particle probability distribution for the energy, and not for the momentum. Now, since the one-particle energy RuffoPRE95

$$ϵ=\frac{p^2}{2}+1M_x\mathrm{sin}\theta M_y\mathrm{cos}\theta $$
(4)
depends both on $`\theta `$ and $`p`$, statistical descriptions in terms of $`ϵ`$ and $`p`$ are equivalent only if the $`\mu `$-space distribution is independent of $`\theta `$, i.e. if the energy distribution function $`f_ϵ(ϵ)`$ can be unambiguously given in terms of $`f_p(p)`$. This is not our case for $`U<U_c`$, unfortunately, since the angle distribution is not uniform, as illustrated by our Figs. 5 and 6, and by LRT’s figure 3(a) in Ref. TsPRE01, . Ignoring the angular dependence of the one-particle distribution in $`\mu `$-space invalidates the identification of $`f_p`$ as an equilibrium distribution derived from a generalized thermodynamical formulation like Tsallis’s, which—just as ordinary thermodynamics—yields canonical distributions as functions of the energy.

We have evaluated $`f_ϵ`$ by means of numerical realizations of the HMF model for $`N=10^4`$. From the energy distribution function we can in turn calculate the energy probability density $`P(ϵ)`$, since $`f_ϵ(ϵ)=\rho (ϵ)P(ϵ)`$, where $`\rho (ϵ)`$ is the density of states. The HMF density of states can be expressed in terms of elliptic functions, as

$$\rho (ϵ)=\frac{2\sqrt{2}}{\pi \sqrt{ϵ1+M}}K\left(\frac{2M}{ϵ1+M}\right)$$
(5)
for $`ϵ>1+M`$ and

$$\rho (ϵ)=\frac{2\sqrt{2}}{\pi \sqrt{ϵ1+M}}F(\frac{1}{2}\mathrm{arccos}\frac{1ϵ}{M},\frac{2M}{ϵ1+M})$$
(6)
for $`1M<ϵ<1+M`$, where $`K`$ and $`F`$ are the complete and incomplete elliptic integrals of the first kind, respectively Ab . The energy is bounded from below, at $`ϵ_{\mathrm{min}}=1M`$. Figure 8 shows the results for $`P(ϵ)`$ at three times, plotted as a function of $`ϵϵ_{\mathrm{min}}`$. The straight line represents the canonical prediction, $`P_{\mathrm{eq}}(ϵ)\mathrm{exp}(ϵ/T)`$. Observe that $`t=5\times 10^4`$ corresponds to the minimum of the nonequilibrium temperature, at the very beginning of the long-time range. For this time we note a considerable deviation from canonical equilibrium, with an abrupt cutoff at $`ϵϵ_{\mathrm{min}}1`$ preceded by an overpopulated interval. For very low energies, however, the coincidence with the equilibrium distribution is already quite good. For $`t=2\times 10^5`$ the cutoff persists but it is less abrupt and has shifted to higher energies, while in the low-energy range the distribution is closer to canonical equilibrium. Finally, for $`t=10^6`$, canonical equilibrium is well established in almost the whole range of energies shown in the figure. Deviations of the nonequilibrium temperature with respect to $`T_{\mathrm{eq}}`$ are therefore to be ascribed to the depleted tail at high energies. From these results, we conclude that in the long-time regime the equilibrium energy distribution of the HMF model with WB initial conditions is built up from low energies, with a depleted high-energy tail that recedes as time elapses.

We remark that LRT have assigned a distribution function to states which are clearly in the short-time range. In that range the system is still strongly nonergodic and, as discussed above, such statistical description is therefore of very limited significance. Their best fitting, for $`N=10^5`$ at $`t=1200`$, corresponds in fact to a state where probably most of the ensemble is still in the condensed phase. Moreover, they compare results for different system sizes at the same time \[see their figure 2(c) in Ref. TsPRE01, \], disregarding the scaling laws of time scales with the size, which threatens the validity of their evaluation of the fitting quality as a function of $`N`$ \[their figure 2(d)\]. Finally, we mention that the abrupt drop of $`f_p`$ for large momenta—also seen at larger times in our results for $`P(ϵ)`$, Fig. 8—eludes any description based on Tsallis’s thermodynamics. The artificial cutoff imposed by LRT to their fitting of $`f_p`$, aimed at mimicking its abrupt drop, can hardly find any consistent justification within such formalism. We conclude that no evidence supports the association of nonequilibrium HMF states with generalized thermodynamical equilibria.

In summary, we have reexamined the slow-relaxation dynamics of the Hamiltonian mean-field model, when the evolution starts from water-bag initial conditions. First, we have revised the evidence that would support the existence of quasistationary states in the nonequilibrium dynamics of the model. In that regard, we have shown that the plateaus in the temporal evolution of the nonequilibrium temperature $`T(t)`$ TsPRE01 ; TsPHA02  are a mere artifact of the use of a logarithmic scale for the time axis in the plot of $`T(t)`$. Moreover, a detailed analysis of the evolution of the nonequilibrium temperature shows that the dynamics can be decomposed into two well-defined dynamical regimes: a short-time range characterized by the presence of a condensed phase, in which the system preserves strong memory of the initial conditions, and a long-time range, in which the condensed phase has evaporated and the system slowly relaxes toward equilibrium. A natural boundary between these two regimes is given by the minimum reached by the nonequilibrium temperature, which recedes toward larger times as the system size grows. The numerical evaluation of the probability density function of the one-particle energy in the long-time range demonstrates that the system is slowly relaxing to the Boltzmann distribution, and does not support any connection with the generalized equilibrium distributions derived from Tsallis’s nonextensive formalism.

It must be stressed that neither previous work nor our present results elucidate the nature of the HMF evolution in the thermodynamical limit $`N\mathrm{}`$. The existence of different dynamical regimes, with different scaling laws in terms of $`N`$, makes an extrapolation from the available numerical results extremely uncertain. The possibility of performing significant, extensive numerical realizations of this model with current computational facilities seems to be restricted to ensembles of, at most, $`N10^4`$ elements. Numerical realizations for much larger systems may be necessary in order to attain a meaningful approach to the thermodynamical limit.

The origin of slow dynamics in the Hamiltonian mean-field model, also revealed by the study of aging phenomena and strong memory effects in the long-time range mta , remains to be explained as well. A plausible argument RuffoCSF02  may perhaps arise from the study of the role of the state of zero magnetization, which for $`U<U_c`$ is thermodynamically unstable, and its connection with critical slowing down RuffoPRL98 —a connection which has often been invoked but not yet investigated. Revealing the basic mechanisms of slow relaxation in the Hamiltonian mean-field model should definitely clarify its relation to systems with some degree of structural disorder or frustration, such as glasses, that also exhibit slow-dynamics features. In turn, this would contribute to a unified view on nonequilibrium phenomena in systems with many degrees of freedom.

###### Acknowledgements.

We thank S. Ruffo for his motivating remarks on the HMF model, and C. Tsallis for comments on his own work. D. H. Z. acknowledges hospitality at the Abdus Salam International Centre for Theoretical Physics (Trieste, Italy), where this work was partially carried out.


format: text

text:
# SHORT BASELINE ACCELERATOR-BASED NEUTRINO OSCILLATION SEARCHES

## 1 Introduction

In the last years, tremendous progress has been achieved to firmly establish the nature of neutrino oscillations using neutrinos from the sun  as well as neutrinos produced in the earth’s atmosphere . However, with oscillation parameters accessible to accelerator based experiments, the situation remains unsettled. There is one evidence for oscillations in the appearance mode $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ from the LSND experiment which is ruled out in part by other experiments such as KARMEN, NOMAD and NuTeV.

In this paper, we summarize the status of the search for appearance of $`\overline{\nu }_e`$ ($`\nu _e`$) from $`\overline{\nu }_\mu `$ ($`\nu _\mu `$) sources indicating the flavor oscillation $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ ($`\nu _\mu `$$``$$`\nu _e`$). Restricting the discussion to this oscillation channel, we can simplify the 3-dimensional neutrino mass and mixing scheme to a 2-dimensional one which leads to an oscillation probability $`P`$ for the appearance of $`\overline{\nu }_e`$ in a pure $`\overline{\nu }_\mu `$ beam of

$$P(\overline{\nu }_\mu \overline{\nu }_e)=\mathrm{sin}^2(2\mathrm{\Theta })\mathrm{sin}^2\left(\frac{1.27\mathrm{\Delta }m^2L}{E}\right).$$
(1)
As can be seen from equ. 1, the experimental sensitivity to the oscillation parameter $`\mathrm{\Delta }m^2`$ is given by the distance $`L`$ between source and detector and the neutrino energy $`E`$. The ratio $`L/E`$ in the above formula is typically given in units of m/MeV (see section 2 for LSND and KARMEN or km/GeV (see section 3 for NOMAD and NuTeV).

## 2 Medium-energy beam stop sources

Existing medium-energy high intensity neutrino sources (LANSCE111Los Alamos Neutron Science Center, Los Alamos, USA, ISIS222ISIS at Rutherford Appleton Laboratory, Chilton, UK) are based on accelerators with a proton kinetic energy of 800 MeV and a beam stop or neutrino production target of high-Z material. Due to the massive target (e.g. Copper, Tantalum or Uranium), the resulting neutrino fluxes arise from $`\pi ^+`$ and $`\mu ^+`$ decays at rest (DAR) and hence are well understood. $`\pi ^{}`$ and $`\mu ^{}`$ that stop are readily captured in the high-Z material of the beam stop. The production of kaons or heavier mesons is negligible at these proton energies. The $`\overline{\nu }_e`$ flux is calculated to be only $`<810^4`$ as large as the $`\overline{\nu }_\mu `$ flux, so that the observation of a $`\overline{\nu }_e`$ event rate significantly above the calculated background would be evidence for $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ oscillations.

Figure 1(a) shows the $`\overline{\nu }_\mu `$ spectrum from $`\mu ^+`$ DAR (solid line) together with the $`\overline{\nu }_e`$ contribution. The absolute flux corresponds to the neutrino production of the ISIS accelerator with 200 $`\mu `$A proton current and a Ta–D2O target. If the time structure of the proton pulses is significantly smaller than the muon lifetime of 2.2 $`\mu `$s, the neutrinos show the corresponding exponential decrease (as in Fig. 1(b), an important feature to suppress cosmic induced background). Otherwise, as for LANSCE with its 600 $`\mu `$s long proton pulses, neutrinos are almost constantly produced during beam-on times.

### 2.1 The LSND evidence for $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$

The LSND experiment took data over the years 1993–1998. During this period the LANSCE accelerator delivered 28 896 C of protons on the production target. The LSND detector consisted of an approximately cylindrical tank of $``$210 m3 volume. The center of the detector was 30 m from the beam stop neutrino source. The tank was filled with liquid scintillator of low scintillator concentration which allowed the detection of both Čerenkov and scintillation light. For all details see  and references therein.

A $`\overline{\nu }_e`$ signal from $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ oscillations consists of a spatially correlated delayed (e+,n) sequence from p ( $`\overline{\nu }_e`$ , e+ ) n . This requires a prompt ’electron-like’ event with energy $`E_e>20`$ MeV followed by a low energy $`\gamma `$–event from neutron capture. The information about the delayed event is compressed into a likelihood ratio $`R_\gamma `$: If within one millisecond after the initial event another event is recorded at distance $`\mathrm{\Delta }r250`$ cm, the ratio of likelihood $`R_\gamma `$ in energy, time and distance of being a correlated p ( n,$`\gamma `$ ) over an accidental coincidence is calculated, otherwise $`R_\gamma =0`$.

To determine the oscillation parameters $`\mathrm{sin}^2(2\mathrm{\Theta })`$ and $`\mathrm{\Delta }m^2`$ in a likelihood analysis, only loose cuts are applied. The beam-on event sample then contains 5697 candidates shown in Fig. 2 with their spectral information in visible energy $`E_e`$, spatial distribution along detector axis $`z`$, the angle between the incident neutrino direction and the visible electron-like track $`\mathrm{cos}\theta _\nu `$ and the (e+,n) correlation parameter $`R_\gamma `$.

The event-based lhd function reaches its maximum at $`(\mathrm{sin}^2(2\mathrm{\Theta }),\mathrm{\Delta }m^2)_{max}=(0.003,1.2\mathrm{eV}^2)`$ corresponding to an excess of 89.5 events attributed to oscillations $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$. With a difference of $`lnL_{max}lnL_{noosc.}=14.5`$, the oscillation hypothesis is strongly favored over the no–oscillation case. The contributions in Fig. 2 (from bottom to top) show the result of a likelihood analysis in all 4 parameters: 3664.6 cosmic induced events, 1936.8 events induced by known SM processes and on top (light grey) the about 90 excess events 333The sum of these 3 contributions leaves 6.1 events attributed to the appearance of $`\nu _e`$ from $`\nu _\mu `$$``$$`\nu _e`$ with $`\nu _\mu `$ arising from $`\pi ^+`$ decays in flight. See  for details. from $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$. We will discuss the results in terms of the oscillation parameters later in section 2.3.

### 2.2 The KARMEN limit for $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$

The KARMEN experiment used the ISIS neutrino source. In contrast to the LANSCE source, the protons are extracted from the synchrotron as an intense but narrow double pulse of 525 ns total width. The unique time structure of the ISIS proton pulses allowed a clear separation of $`\nu `$–induced events from any beam unrelated background. The KARMEN detector was a rectangular high resolution liquid scintillation calorimeter, located at a mean distance of 17.7 m from the ISIS target and shielded by a multilayer active veto system and 7000 t steel. The liquid scintillator volume (65 m3) was optically separated into 512 independent central modules. Gadolinium was implemented between the module walls for an efficient detection of thermal neutrons Gd ( n,$`\gamma `$ ) with on average 3 $`\gamma `$’s of energy $`E_\gamma =8`$ MeV.

The KARMEN 2 experiment took data from February 1997 to March 2001. During this time, protons equivalent to a total charge of 9425 Coulombs have been accumulated on the ISIS target.

As for LSND, a $`\overline{\nu }_e`$ signal from $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ oscillations consists of a spatially correlated delayed (e+,n) sequence. The requirements for event sequences in KARMEN are described in detail in  and references therein. Applying all cuts to the data, 15 (e+,n) candidate sequences were finally reduced. Figure 3 shows the remaining sequences in the appropriate energy and time windows.

The background components are also given with their distributions. All components except the intrinsic $`\overline{\nu }_e`$ contamination are measured online in different time and energy windows (see Table 1).

The extracted number of sequences is in excellent agreement with the background expectation, consistent with no oscillation signal. To also include the detailed spectral information of each individual event, an event-based maximum likelihood method is applied. Fig. 4(a) shows the results of this analysis in terms of oscillation events as function of $`\mathrm{\Delta }m^2`$. The best fit line is compatible with zero. Applying a unified frequentist approach  leads to the given 90% C.L. upper limit which can be compared to what one would expect as signal strength in KARMEN taken the evidence from LSND.

### 2.3 LSND and KARMEN: a critical review

Fig. 4 shows a comparison of the LSND and KARMEN results obtained after analysing the final data samples with an event-based maximum likelihood method. In plot 4(a), the LSND evidence with its statistical error is translated into an expected range of oscillation events among the KARMEN data. At large values of $`\mathrm{\Delta }m^2`$, the experimental outcomes are in clear contrast whereas at small $`\mathrm{\Delta }m^2`$, the expected event number is below the 90% C.L. exclusion curve. Fig. 4(b) shows the same results, the KARMEN exclusion curve (for large $`\mathrm{\Delta }m^2`$ $`\mathrm{sin}^2(2\mathrm{\Theta })<1.710^3`$ 90% C.L.) as well as the evidence band of LSND ($`lnL_{max}2.3`$, $`lnL_{max}4.6`$), now as functions of the oscillation parameters $`\mathrm{sin}^2(2\mathrm{\Theta })`$ and $`\mathrm{\Delta }m^2`$.

Since both experiments have similar sensitivity to the oscillation parameters but different central statements about oscillations leading to a partial overlap of exclusion curve and favored region, one has to solve the following problems quantitatively: What is the level of compatibility of both experiments? Assuming statistical compatibility, what are the oscillation parameters accepted by both experiments?

A detailed combined statistical analysis  has been performed to answer these questions. We summarize the main results also shown in Fig. 5. At a combined level of 60 % $`\times `$ 60 % =36 %, the LSND and KARMEN results are incompatible. Taking the rest probability of 64 %, the favored parameter regions are shown in plot 5(a). In Fig. 5(b), the parameter region of 90 % $`\times `$ 64 % =58 % probability of occurrence is shown together with the individual results obtained in the same analysis .

To summarize, LSND and KARMEN are incompatible at a level of 36 % confidence. Assuming statistical compatibility, all parameter combinations with $`\mathrm{\Delta }m^2>1`$ eV2 are excluded apart from a little ’island’ at $`7`$ eV2.

## 3 Experiments with neutrino beams

In section 2, we dealt with neutrino energies of tens of MeV. Neutrinos from high energy proton beams typically have energies of tens up to hundreds of GeV. These higher energies increase the cross section for neutrino interactions with the detector material as well as decrease substantially any beam-uncorrelated background. However, the neutrino flux has to be modeled carefully and cannot –as for the beam stop neutrino source– simply calculated analytically. Detectors are typically at a distance of several hundreds of meters from the neutrino source, leading to a so-called short baseline regime: $`L/E0.5/20`$ km/GeV corresponding to a very high sensitivity in $`\mathrm{\Delta }m^2`$ for values larger than several tens of eV2 mainly limited by the intrinsic ’wrong’ flavor component within the beam, typically of the order of $`\overline{\nu }_e/\overline{\nu }_\mu (\nu _e/\nu _\mu )`$ 1–2 %.

### 3.1 The NOMAD $`\nu _\mu `$$``$$`\nu _e`$ search

The NOMAD experiment took data from 1995–1998 using the CERN wide band neutrino beam from the 450 GeV PS. We summarize here only the search for $`\nu _\mu `$$``$$`\nu _e`$ and refer the reader for all details of the experiment as well as the $`\nu _\mu `$$``$$`\nu _\tau `$ search to reference . The MOMAD detector was a high resolution detector to separate $`\nu `$ induced charged current (CC) events by kinematic criteria with a momentum resolution of 3.5 % ($`p<10`$ GeV). It had a fine grained calorimetry leading to $`\mathrm{\Delta }E/E=3.2\%/\sqrt{E}1`$ % and a particle identification via a TRD resulting in a pion rejection of $`10^3`$ with electron efficiency $`>90`$ %.

Fig. 6 shows a simulation of the $`\nu _\mu `$ and $`\nu _e`$ spectra within a transverse fiducial area of $`260\times 260`$ cm2 of the detector .

The NOMAD $`\nu _\mu `$$``$$`\nu _e`$ search was an appearance experiment based on electron identification with $`\nu _e`$ contamination in the $`\nu _\mu `$ beam of $`\nu _e/\nu _\mu 1.5`$ %. The different energy spectra and radial distributions of $`\nu _e`$ and $`\nu _\mu `$ within the detector were used to calculate an expected ratio $`R_{e\mu }=(\nu _e\mathrm{CC})/(\nu _\mu \mathrm{CC})`$ of $`\nu `$–induced events as function of $`E_\nu `$ and $`r`$ (see Fig. 7(left) for a typical $`\nu _\mu `$–induced CC event). This ratio was then compared with the reconstructed events containing 748000 $`\nu _\mu `$ CC and 8000 $`\nu _e`$ CC candidates. The $`R_{e\mu }`$ distribution in data and MC (Fig. 7) shows no excess of $`\nu _e`$ CC, therefore leading to an exclusion of oscillation parameters, i.e. for large $`\mathrm{\Delta }m^2`$ $`\mathrm{sin}^2(2\mathrm{\Theta })<1.210^3`$ at 90% C.L.  (see also Fig. 10 for the preliminary NOMAD $`\nu _\mu `$$``$$`\nu _e`$ exclusion curve). It is interesting to note that the NOMAD $`\nu _\mu `$$``$$`\nu _e`$ result underlines the KARMEN exclusion at values $`\mathrm{\Delta }m^2>10`$ eV2 and that the 90% C.L. exclusion curve narrowly passes the allowed ’island’ at $`\mathrm{\Delta }m^27`$ eV2.

### 3.2 The NuTeV $`\nu _\mu `$$``$$`\nu _e`$ and $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ searches

The NuTeV experiment at Fermilab took data from 1996–1997, using the 800 GeV primary proton beam from the Tevatron. A sign selected quadrupole train lead to a very effective separation of neutrino and anti-neutrino beams, allowing separate searches of $`\nu _\mu `$$``$$`\nu _e`$ and $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ oscillations by appearance of $`\nu _e`$ or $`\overline{\nu }_e`$, respectively. After a 320 m decay pipe followed by 915 m of shielding, the detector consisted mainly of a sandwich system of steel slabs and scintillator planes. CC events induced by $`\nu _\mu `$ or $`\nu _e`$ can be separated via the energy deposit in the first 3 scintillators compared to the total deposit $`\eta 3=(E_1+E_2+E_3)/E_{\mathrm{cal}}`$: ’short’ tracks indicate electron like events whereas ’long’ tracks point to muons. For all details of the oscillation search see  and references therein for the experimental setup.

A detailed beam simulation predicts the visible energy spectra from various $`\nu `$–induced CC events given in Fig. 8(a) and (b). Fig. 8(c) and (d) show the ratios of the detected over predicted ($`\nu _e`$, $`\overline{\nu }_e`$) event numbers versus visible energy minus 1. The lines correspond to the predictions of various oscillation scenarios, the solid is the 90% C.L. upper limit for $`\mathrm{\Delta }m^2=1000`$ eV2. The exclusion curve for $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ as function of the oscillation parameters relevant for LSND is shown in Fig. 10 excluding the very high $`\mathrm{\Delta }m^2`$ region thereby confirming KARMEN and NOMAD.

### 3.3 The MiniBooNE experiment

The MiniBooNE experiment at Fermilab is designed to cover the entire parameter range favored by LSND with a completely different experimental setup  by looking for $`\nu _\mu `$ disappearance and $`\nu _e`$ appearance from $`\nu _\mu `$$``$$`\nu _e`$ in a $`\nu _\mu `$ beam.

MiniBooNE will use the 8 GeV booster at Fermilab’s main injector to produce a neutrino beam with significant contributions at energies above 500 MeV. The Čerenkov detector is a spherical tank of 40 ft diameter filled with 776 tonnes of mineral oil viewed by 1275 8” PMT’s at a distance of 500 m from the neutrino production point. Again, the neutrino flux must be known precisely by MC simulations and special checks. The expected fluxes as well as a schematic view of the detector are shown in Fig. 9.

After 2 years of data taking corresponding to $`10^{21}`$ protons on target, MiniBooNE expects to detect about $`510^5`$ $`\nu _\mu `$–induced CC events and 3500 electron-like events out of which about 1000 events would arise from a $`\nu _\mu `$$``$$`\nu _e`$ oscillation according to the LSND signal . The expected experimental sensitivity is shown in Fig. 10 as exclusion curve in case there will be no oscillation signal. The detector is already operational, the first neutrino beam is scheduled to be delivered mid of 2002.

## 4 Conclusions

We have discussed the LSND evidence for $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$ oscillations together with the negative searches performed by KARMEN, NOMAD and NuTeV in the same flavor channel. The KARMEN result is the most limiting one confirmed by NOMAD and NuTeV at higher values of $`\mathrm{\Delta }m^2`$. Though there are no further hints than the LSND result for neutrino oscillations with a mass difference $`\mathrm{\Delta }m^2`$ in the order or larger than 0.1 eV2, no existing short baseline experiment can rule it out completely. The upcoming MiniBooNE experiment should be sensitive enough to definitively confirm or refute the LSND result interpreted as appearance of $`\overline{\nu }_e`$ due to oscillations $`\overline{\nu }_\mu `$$``$$`\overline{\nu }_e`$.


text:
# Multiple molecular dynamics time-scales in Hybrid Monte Carlo fermion simulations.

## 1 MULTIPLE MD TIME-SCALES

A scheme for integrating the equations of motion for the MD phase of the HMC algorithm by introducing different time-scales for different segments of the action was introduced in Ref. .

The leap-frog integrator is constructed from the two simple time-evolution operators generated by the the kinetic and potential energy terms. Their effect on the system coordinates, $`\{p,q\}`$ are

$`V_T(\mathrm{\Delta }\tau ):\{p,q\}`$ $``$ $`\{p,q+\mathrm{\Delta }\tau p\}`$
$`V_S(\mathrm{\Delta }\tau ):\{p,q\}`$ $``$ $`\{p\mathrm{\Delta }\tau S,q\}`$ (1)
The simplest reversible leap-frog integrator is then

$$V(\mathrm{\Delta }\tau )=V_S(\frac{\mathrm{\Delta }\tau }{2})V_T(\mathrm{\Delta }\tau )V_S(\frac{\mathrm{\Delta }\tau }{2}).$$
(2)
If the action (and thus the Hamiltonian) is split into two parts,

$`=\underset{_1}{\underset{}{T(p)+S_1(q)}}+\underset{_2}{\underset{}{S_2(q)}},`$ (3)
then the two leap-frog integrators for these two pieces can be constructed as

$`V_1(\mathrm{\Delta }\tau )`$ $`=`$ $`V_{S_1}({\displaystyle \frac{\mathrm{\Delta }\tau }{2}})V_T(\mathrm{\Delta }\tau )V_{S_1}({\displaystyle \frac{\mathrm{\Delta }\tau }{2}})`$
$`V_2(\mathrm{\Delta }\tau )`$ $`=`$ $`V_{S_2}(\mathrm{\Delta }\tau )`$ (4)
and a reversible integrator for the full Hamiltonian can be constructed by combining the two schemes:

$$V(\mathrm{\Delta }\tau )=V_2(\frac{\mathrm{\Delta }\tau }{2})\left[V_1(\frac{\mathrm{\Delta }\tau }{m})\right]^mV_2(\frac{\mathrm{\Delta }\tau }{2}),$$
(5)
with $`mZ`$. This compound integrator effectively introduces two evolution time-scales, $`\mathrm{\Delta }\tau `$ and $`\mathrm{\Delta }\tau /m`$.

We suggest that the multiple time-scale scheme is only helpful when two criteria are fulfilled simultaneously:

1. the force term generated by $`S_1`$ is cheap to compute compared to that of $`S_2`$ and
2. the split captures the high-frequency modes of the system in $`S_1`$ and the low-frequency modes in $`S_2`$.

The popular implementation of the scheme in dynamical fermion simulations with pseudofermions is to split the Hamiltonian into the Yang-Mills term and the pseudofermion action:

$`S_1`$ $`=`$ $`S_G,`$
$`S_2`$ $`=`$ $`\varphi ^{}[M^{}M]^1\varphi `$ (6)
Unfortunately for light fermions, the highest frequency fluctuations are in the pseudofermion action, which also has the more computationally expensive force term, thus the criteria are not met. The key to using the method then is to implement a low-computational cost scheme that separates the high- and low-frequency fermion molecular-dynamics.

## 2 POLYNOMIAL FILTERING

A very low-order polynomial approximation offers a cheap means of mimicking most of the short-distance physics of the fermion interations. Recent interest in polynomial approximations to a matrix inverse was inspired by the multi-boson algorithm of Ref. . This led to the development of the Polynomial HMC (PHMC) algorithm . Following this idea, we write an exact representation of the two-flavour probability measure

$$detM^{}M=𝒟\varphi 𝒟\varphi ^{}𝒟\chi 𝒟\chi ^{}e^{\left\{S_\varphi S_\chi \right\}}$$
(7)
with

$$S_\varphi =\left|[M𝒫(M)]^1\varphi \right|^2\text{ and }S_\chi =\left|𝒫(M)\chi \right|^2$$
(8)
The fields $`\varphi `$ are modified pseudofermions and we term the new fields $`\chi `$, “guide” bosons. The algorithm exactly recovers the probability measure for the two-flavour theory for any choice of polynomial. This scheme is similar to the split-pseudofermion method of Ref.  with one distinction; using polynomials allows us to modify the fermion modes easily and cheaply. The success of the algorithm hinges on the empirical observation that the fermion modes that induce the high-frequency fluctuations are localised and can be handled by very low-order polynomials. This means the force term for the guide bosons, generated by $`S_\chi `$ can be computed cheaply and a multiple time-scale integrator can be built that simultaneously fulfills the two criteria discussed earlier. The time-scale split is then

$$S=\underset{UV,S\text{ cheap}}{\underset{}{S_G+S_\chi }}+\underset{IR,S\text{ costly.}}{\underset{}{S_\varphi }}$$
(9)
Note that in general, $`p`$ time-scales can be introduced straightforwardly by first adding $`p1`$ sets of guide fields, then constructing a heirarchy of leapfrog integrators, with $`m_i`$ leap-frog steps at the $`i^{\mathrm{th}}`$ level. This introduces $`p`$ time-scales, $`\mathrm{\Delta }\tau >\frac{\mathrm{\Delta }\tau }{m_{p1}}\mathrm{}\frac{\mathrm{\Delta }\tau }{m_{p1}\mathrm{}m_1}`$.

## 3 TESTING THE ALGORITHM

The performance of the method is under investigation in simulations of the two-flavour massive Schwinger model ($`2d`$ QED). HMC runs presented here are performed on $`64\times 64`$ lattices, with $`\beta =4.0,\kappa =0.2618`$. Non-hermitian Chebyshev polynomial approximations are used . Measurements of the dependence of the acceptance probability on the polynomial degree and tuning the MD parameters, $`\mathrm{\Delta }\tau `$ and $`m`$ are being made.

Fig. 1 shows the acceptance rate for the standard HMC algorithm (with even-odd preconditioned pseudofermions) along with the polynomial-filtered method. Two different polynomials were used; $`N=4`$ and $`N=16`$. The solid lines are fits to the expected behaviour of HMC as $`\mathrm{\Delta }\tau 0`$, namely

$$P_{\mathrm{acc}}=\text{erfc}\left\{(\mathrm{\Delta }\tau /\tau _0)^2\right\},$$
(10)
where the value of $`\tau _0`$ is determined by the best fit. $`\tau _0`$ is then a characteristic time-scale for the modes encapsulated in the pseudofermionic part of the action, $`S_\varphi `$ in Eqn. 8. Since evaluation of the force term $`S_\varphi `$ dominates in this parameter range, $`\tau _0`$ is a good indicator of algorithm performance.

## 4 TUNING THE ALGORITHM

The algorithm has a number of free parameters, allowing a good deal of scope to optimise the performance. Fig. 2 shows the characteristic time-scale for pseudofermion integration, $`\tau _0`$ as a function of Chebyshev polynomial degree. $`\tau _0`$ rises very rapidly initially as the degree of the filter algorithm is increased, and $`\tau _0`$ for $`N=16`$ is about four times larger than for the standard HMC algorithm. Inverting the pseudofermion matrix requires roughly the same computational cost, and thus the algorithm is four times more efficient than HMC (assuming autocorrelations are the same).

For low values of $`N`$ (the polynomial degree) the computational bottleneck is solver performance while for larger values, the evaluation of the force term arising from interactions between the gauge fields and the guide bosons begins to dominate. Effective and simple strategies for tuning the algorithm are still under investigation. We are also investigating alternative choices of polynomial filters beyond Chebyshev approximation.

MP is grateful to Enterprise-Ireland for support under grant SC/01/306.


text:
# Theories of Everything and Hawking’s Wave Function of the UniverseTo appear in The Future of Theoretical Physics and Cosmology: Stephen Hawking 60th Birthday Symposium, Cambridge University Press with two additional figures.

## 1 Introduction

It is an honor, of course, to participate in this celebration of Stephen’s 60th birthday and to address such a distinguished audience. But for me it is a special pleasure. I first met Stephen over thirty years ago when we were both beginning in science at the then Institute of Theoretical Astronomy here in Cambridge. I have been following his inspiration ever since and have had the privilege of working with him on several directions that he has pioneered. I want to talk to you about one of these today.

## 2 Different Things Fall with the Same
Acceleration in a Gravitational Field

If a cat, a cannonball, and an economics textbook are all dropped from the same height, they fall to the ground with exactly the same acceleration (9.8m/s2) under the influence of gravity. This equality of the gravitational accelerations of different things is one of the most accurately tested laws of physics. The accelerations are known to be equal to an accuracy of a few parts in a thousand billion. That law, however, tells you very little about cats, cannonballs, or economics. That, in a nutshell, is the theme of my talk today (Figure 1).

Were a cat, a cannonball, and an economics text dropped from a real leaning tower, the accelerations would not be equal to an accuracy of a thousand billion. Air resistance, magnetic fields, and a host of other effects mask the purely gravitational one. Rather, the equality of accelerations is revealed to that accuracy only in delicate laboratory experiments shielded from such extraneous effects, or in the study of the accelerations of astronomical bodies for which they are not important. The accuracy record holder at the moment is the lunar laser ranging demonstration that the Earth and the Moon fall with the same acceleration toward the Sun .

The manned and unmanned missions to the Moon in the late ’60s and early ’70s left behind arrays of corner reflectors. A corner reflector is an arrangement of three mirrors making up the sides of a cube that meet in one corner. It has the property that any light ray incident on the reflector bounces from mirror to mirror and is reflected back in exactly the direction from whence it came. Since 1969 a systematic program to monitor the Moon’s orbit using these reflectors has been carried out by the McDonald Observatory in Texas and the Observatoire de Côte d’Azur in Grasse, France . High powered lasers send light pulses to the Moon that bounce off the reflectors. By measuring precisely the time it takes for the light pulse to return (about 1 sec), and dividing that by the velocity of light, the distance to the Moon can be determined very accurately. About a billion, billion photons are sent out ten times a second and one returning photon is detected every few seconds. The returning photons have been detected for over 30 years, and as a consequence of this effort, we know the position of the Moon relative to the Earth to an accuracy of about a few centimeters. Study of its orbit shows it is falling towards the Sun with the same acceleration as the Earth to an accuracy of a few parts in a thousand billion.

This law of equality of accelerations of different things is a cornerstone of Einstein’s general relativity — Stephen’s subject. However, it is not general relativity specifically that I want to talk to you about today, but about such laws of physics more generally.

## 3 The Fundamental Laws of Physics

The equality of gravitational accelerations of different things is an example of a regularity of nature. Everything falls in exactly the same way. The regularity is universal. No exceptions! There are other regularities that may apply to specific systems in similar circumstances. For example, if we drop three cats, they might all go ‘meeeeow’ in roughly the same way on the way down. Or perhaps there are only statistical regularities for cats — 8 out of 10 cats will go ‘meeeeow’ in the average drop. Other systems will do something different according to their own particular regularities. Cannonballs for example do nothing on the way down. But the equality of accelerations of different things is a special kind of law because it applies universally to anything falling in a gravitational field. It is an example of a fundamental law of physics.

Identifying and explaining the regularities of nature is the goal of science. Physics, like other sciences, is concerned with the regularities exhibited by particular systems. Stars, atoms, fluid flows, high temperature super-conductors, black holes, and the elementary particles are just some of the many examples. Studies of these specific systems define the various subfields of physics — astrophysics, atomic physics, fluid mechanics, and so forth. But beyond the regularities exhibited by specific systems, physics has a special charge. This is to find the laws that govern the regularities that are exhibited by all physical systems — without exception, without qualification, and without approximation. The equality of gravitational accelerations of different things is an example. These are usually called the fundamental laws of physics. Taken together they are called informally a “theory of everything”. Stephen has been a leader in the quest for these universal laws . Today, I would like to ask the question: “How much do we know about the world if we have a theory of everything?”

Ideas for the nature of the fundamental laws have changed as experiment and observation have revealed new realms of phenomena and reached new levels of precision. But since they have been studied, it has been thought that the fundamental laws consist of two parts

* Dynamical laws that govern regularities in time. Newton’s laws of motion governing the orderly progression of the planets or the trajectory of a tennis ball are examples, as is the law that different things fall with the same acceleration in a gravitational field and the Einstein equation governing the evolution of the universe.
* Initial conditions that govern how things started out and therefore most often specify regularities in space. The statistical regularities of the large scale distribution of the galaxies in the universe is a possible example (Figure 2).

Laplace is a name associated with this view of a theory of everything. He wrote, famously,

> An intelligence knowing, at any given instant of time, all forces acting in nature, as well as the momentary positions of all things of which the universe consists, would be able to comprehend the motions of the largest bodies of the world and those of the smallest atoms in one single formula. …To it, nothing would be uncertain, both future and past would be present before its eyes.

That was the theory of everything circa 1820 — Newtonian determinism.

Both parts of a theory of everything are needed to make any predictions. Newton’s dynamical laws by themselves don’t predict the trajectory of a tennis ball you might throw. To predict where it goes, you must also specify the position from which you throw it, the direction, and how fast. In short, you must specify the ball’s initial condition. One of Stephen’s most famous achievements is just such an initial condition , but not for tennis balls. Stephen’s no-boundary initial condition is for the whole universe.

The search for a theory of the dynamical laws has been seriously under way since the time of Newton. Classical mechanics, Newtonian gravity, Maxwell’s electrodynamics, special and general relativity, quantum mechanics, quantum field theory, and superstring theory are but some of the milestones in this search. But the search for a theory of the initial condition of the universe has been seriously under way for only the twenty years since Stephen’s pioneering work on the subject. Why this difference? The examples used above to discuss regularities governed by the two parts of a theory of everything hint at the answer. The trajectory of a tennis ball was used to illustrate the regularities of dynamical laws, and the large scale distribution of galaxies was used to illustrate the regularities implied by the law of the initial condition. There is a difference in the kind and scale of regularities that the two laws predict.

Dynamical laws predict regularities in time. It is a fortunate empirical fact that the fundamental dynamical laws are local — both in space and time. The trajectory of a tennis ball depends only on conditions that are nearby both in space and time, and not, for example, either on what is going on in distant parts of the universe or a long time ago. This is fortunate because that means that dynamical laws can be discovered and studied in laboratories on Earth and extrapolated, assuming locality, to the rest of the universe. For example, because it is local, the law that different things fall with the same acceleration in a gravitational field can be discovered by experiments in laboratories here, and indeed all over the galaxy (Figure 4). Without that simplicity of the dynamical laws in the here and now it is possible that we would never have discovered them.

By contrast, the regularities governed by the law of the initial condition of the universe occur on large, cosmological scales. The universe isn’t simple on small spatial scales. Look at the disorder or complexity in the room we’re in right now for example. But the universe is simple on large, cosmological scales — more or less the same in one direction as in any other, more or less the same in one place as in any other.

The mottled map of the cosmic background radiation in Figure 3 illustrates this. In the beginning the universe was opaque, hot, and glowing. As it expanded it got cooler. A few hundred thousand years after the big bang the matter cooled to a temperature at which it was transparent. Light from that time has been propagating to us ever since, cooled to a temperature now of only a few degrees above absolute zero. This is the cosmic background radiation and a map of the temperature of the radiation in different directions, such as that in Figure 3, is the earliest picture we have of the universe.

The mottled map in Figure 3 may not look like evidence that the early universe was smooth, but the temperature difference between the light and dark regions is only 30 millionths of a degree. At resolution of your eye the universe would look completely uniform at this time. The reason serious research on the initial condition at the big bang began only recently is that it is only recently that we have had cosmological observations like this that are detailed enough to reveal the regularities it governs.

## 4 Quantum Mechanics

There is another way in which our vision of the fundamental laws and the nature of a theory of everything has changed since the times of Newton and Laplace. That is quantum mechanics. We don’t yet know the final form the fundamental laws will take. But the inference is inescapable from the physics of the last seventy-five years that they will conform to that subtle framework of prediction we call quantum mechanics.

In quantum mechanics, any system — the universe included — is described by a wave function $`\mathrm{\Psi }`$. There is a dynamical law called the Schrödinger equation which governs how the wave function changes in time

$$i\mathrm{}\frac{d\mathrm{\Psi }(t)}{dt}=H\mathrm{\Psi }(t)(\mathrm{dynamical}\mathrm{law}).$$

Here the operator $`H`$, called the Hamiltonian, summarizes the dynamical theory. There are different forms of $`H`$ for Maxwell’s electrodynamics, for a theory of the strong interactions, etc. Like Newton’s laws of motion, the Schrödinger equation doesn’t make any predictions by itself, it requires an initial condition. This is

$$\mathrm{\Psi }(0),(\mathrm{initial}\mathrm{condition}).$$

When we consider the universe as a quantum mechanical system, this initial condition is Hawking’s wave function of the universe

Probabilities are the key difference between classical and quantum mechanics. Let’s first think about probabilities in classical physics. If I say that there is a 60% chance of hitting an audience member if I toss a ball in this room, I am not expressing a lack of confidence that its trajectory will be governed by the deterministic laws of Newtonian mechanics. Rather the 60% reflects my ignorance of the exact initial speed I’ll impart to the ball, of the influence of air on its motion, and perhaps my ability to do an accurate calculation. If I practice to control the initial condition when it’s thrown, the subsequent evolution of the tennis ball becomes more certain. Probabilities in classical physics result from ignorance.

But in quantum mechanics probabilities are fundamental and uncertainty inevitable. No amount of careful determination of the present state of the tennis ball will achieve certainty for its trajectory. In quantum mechanics there is some probability that a ball will take any trajectory as it leaves my hand. However, in classical situations one trajectory — the one obeying Newton’s laws — is much more probable than all the others. The determinism of classical physics is an approximation, but an approximation on which we can rely in many practical circumstances.

## 5 A Theory of Everything is Not a Theory of Everything

My colleague Murray Gell-Mann used to ask me “If you know the wave function of the universe, why aren’t you rich?” A quantum mechanical theory of the Hamiltonian and the initial state does predict probabilities for every event that might happen in the universe. In that is one sense in which it is a “theory of everything”. However, only a few things are predicted with near certainty. The vast majority of alternatives are predicted with approximately 50%–50% probabilities giving no useful information. That’s the sense in which a “theory of everything” is not a theory of everything.

We hope that the probability is high for the form of the effective dynamical theories that suffice to predict the outcomes of experiments done in every laboratory constructed so far. We also hope that a theory of everything predicts some large scale features of the universe with near certainty such as its approximate smoothness on large scales, the statistics of the distribution of galaxies illustrated in Figure 2, the vast disparity between the age of the universe and the time scales characterizing the fundamental dynamical theory, and so forth.

But it is too much to hope that an interesting probability like that for the FTSE to go up tomorrow is in this category. Since a rise tomorrow is an event in the history of the universe, a quantum mechanical probability could, in principle, be calculated for it from a theory of everything. (Although I suspect that it is well beyond our present powers of computation.) But it’s likely that, after all that work, the predicted probability would be 50% for an upward tick. That is why you can’t get rich knowing the wave function of the universe, and that is why a theory of everything isn’t a theory of everything.

But there is a deeper reason why a theory of everything doesn’t explain everything — it’s too short. Everyone with a PC knows that picture files are typically a lot longer than text files. The files for the undergraduate general relativity text that I am writing, for example, consists of roughly 1Mb of text and 100Mb of pictures. A very rough estimate is that it would take a billion, billion, billion, billion, billion, billion, billion, billion, billion ($`10^{81}`$) compact disks (CDs) to describe the visible universe at a reasonably coarse-grained level of classical description. And that is just at one moment in time! All the matter in the planets, stars, and galaxies in the visible universe would not be enough to make this number of CDs.

However, the description of the universe can be compressed because it exhibits regularities on large scales. Compression is an idea familiar from computation. When a large file exhibits regularities its length can be compressed. For instance, a string of 1000 zeros ‘00000000000000000000…0000’ in a file is a kind of regularity. It can be replaced by a shorter string ‘1000 zeros’. The regularities summarized by the laws of nature similarly permit compression. Rather than saying, ‘cannonballs fall at 9.8m/s2, cats fall at 9.8m/s2, economics textbooks fall at 9.8m/s2 etc.’, the laws of nature permit us to say ‘everything falls at 9.8m/s2 near the surface of the Earth’. That is shorter and therefore more useful. It is possible that everything we see, every detail of every leaf, every event in human history, every thought, is a long string which is compressible to a law implementable by a very short computer program. But there is no evidence that the universe exhibits such regularity, and even Laplace did not propose such a thing.

The initial condition of a deterministic classical theory that did not allow ignorance would have to be as detailed as a present description of the universe. It seems likely therefore that, were world governed by classical laws, it might take so many CDs to write out the law of its initial condition that we could never make the CDs to do it.

Let’s contrast this situation in classical physics with Stephen’s no-boundary quantum wave function of the universe. It’s given by the following simple formula

$$\overline{)\mathrm{\Psi }=\delta g\delta \varphi e^{I[g,\varphi ]}}\left(\genfrac{}{}{0pt}{}{\mathrm{The}\mathrm{no}\mathrm{boundary}\mathrm{wave}}{\mathrm{function}\mathrm{of}\mathrm{the}\mathrm{universe}}\right)$$

That’s short! — maybe 45 key strokes to write out. Furthermore, it’s as complete a description of the initial condition as it is possible to have in quantum mechanics. It implies uncertainty but contains no ignorance. Presumably, the basic equations of the fundamental dynamical theory are similarly short.

You might think this a little misleading because I haven’t said what $`\mathrm{\Psi }`$ means, what $`g`$ means, and what $``$ means, etc. But, even including the lengths of the ten or so texts that physics students read to understand what it all means, a maximum of 10 or 20 Mb would be needed to state the law — easily fitting on one CD. This means the law is the discoverable and implementable, applicable to every happening in the universe, but predicting the near certainty of only a few of its many regularities. A complete discoverable theory of everything is only possible in quantum mechanics where some things are predictable but not everything.

Before I left Santa Barbara to come to Cambridge, my colleague Steve Giddings asked me what I was going to talk about in this lecture. I said I would speak on the question “What do we know if we have a theory of everything?” and answer it “Not all that much.” He replied something like “I hope you’re going to say something more hopeful than that!” But it is a hopeful message. It is only because so little of the complexity of the present universe is predicted by the fundamental universal laws that we can discover them.

## 6 Reduction

Where do all the other regularities in the universe come from, those particular to specific systems — those of the behavior of cats as they fall, those studied by the environmental sciences like biology, geology, economics, and psychology? They are the results of chance events that occur naturally over the history of a quantum mechanical universe. As my colleague, Murray Gell-Mann puts it , they are frozen accidents. “Chance events of which particular outcomes have a multiplicity of long-term consequences, all related by their common ancestry”.

The regularities of cats probably do depend a little bit on the fundamental physical laws, for example, an initial condition that is smooth across the universe, leads to three spatial dimensions, etc. But the origin of most of their regularities can be traced to the chance events of four billion years of biological evolution. Cats behave in similar ways because they have a common ancestry and develop in similar environments. The mechanisms which produce those chance events that led to cats are very much dependent on fundamental biochemistry and ultimately atomic physics. But the particular outcomes of those chances have little to do with the theory of everything.

Do psychology, economics, biology reduce to physics? The answer is YES, because everything considered in those subjects must obey the universal, fundamental laws of physics. Every one of the subjects of study in these sciences — humans, market tables, historical documents, bacteria, cats, etc. — fall with the same acceleration in a gravitational field. The answer is NO, because the regularities of interest in these subjects are not predicted by the universal laws with near certainty even in principle. They are frozen quantum accidents that produce emergent regularities. The answer depends upon what you mean by reduce.

## 7 The Main Points Again

I have always liked the part of the BBC news called “The Main Points Again”. Here are my main points again:

* The fundamental laws of physics constituting a ‘theory of everything’ are those which specify the regularities exhibited by every physical system, without exception, without qualification, and without approximation.
* A theory of everything is not (and cannot be) a theory of everything in a quantum mechanical universe.
* If it’s short enough to be discoverable, it’s too short to predict everything.
* The regularities of human history, personal psychology, economics, biology, geology, etc. are consistent with the fundamental laws of physics, but do not follow from them.

But remember also, especially on this occasion, that all the beautiful regularities that we observe in the universe, certain or not, predictable or not, could be the result of quantum chances following from the fundamental dynamical theory and Hawking’s no-boundary wave function of the universe.

## Acknowledgment

The author’s research was supported in part by the National Science Foundation under grant PHY00-70895.


text:
# On the continuum limit of the Discrete Regge model in 4d

## 1 REGGE QUANTUM GRAVITY

One promising method to quantize the theory of gravitation employs the Euclidean path integral

$$Z=D[q]e^{I(q)},$$
(1)
where the partition function describes a fluctuating space-time manifold. In the Regge approach the (quadratic) link lengths $`q_l`$ represent the dynamical degrees of freedom , deforming continuously a simplicial lattice with fixed incidence matrix, whereas in the somehow complementary approach of dynamical triangulation the incidence matrix is fluctuating with constant link lengths. Any of these two approaches is plagued with various problems; for the Regge approach see .

In “conventional” Regge theory the Regge-Einstein action including a cosmological term ,

$$I_r=\beta \underset{t}{}A_t\delta _t+\lambda \underset{s}{}V_s,$$
(2)
is used. The first sum runs over all products of triangle area $`A_t`$ times corresponding deficit angle $`\delta _t`$ weighted by the bare gravitational coupling $`\beta `$. The second sum extends over the volumes $`V_s`$ of the 4-simplices of the lattice and allows together with the cosmological constant $`\lambda `$ to set an overall scale in the action.

The Discrete Regge model was invented in an attempt to reformulate (1) as the partition function of a spin system . It is defined by restricting the squared link lengths to take on only two values

$$q_l=b_l(1+ϵ\sigma _l),\sigma _lZ_2,$$
(3)
where $`b_l=1,2,3`$, and $`4`$ for edges, face diagonals, body diagonals and the hyperbody diagonal of a hypercube, respectively, is chosen to allow for fluctuations around flat space. The Euclidean triangle inequalities are fulfilled automatically as long as $`ϵ`$ is smaller than a maximum value $`ϵ_{\mathrm{max}}`$. The measure $`D[q]`$ in the quantum gravity path integral is taken to unity for all possible link configurations.

Numerical simulations of the $`Z_2`$ system become extremely efficient by implementing look-up tables and a heat-bath algorithm. In this work typically $`\mathrm{200\hspace{0.17em}000}\mathrm{500\hspace{0.17em}000}`$ iterations have been generated. Calculations have been performed with the parameter $`ϵ=0.0875`$ and the cosmological constant $`\lambda =0`$ because (3) already fixes the average lattice volume.

## 2 RESULTS

The Discrete Regge model – like full Regge theory – exhibits two phase transitions . One is located at a negative value and the other one at a positive value of the bare gravitational coupling. Although earlier work concentrated on the latter transition, there is no reason for favoring a positive value of $`\beta `$ over a negative one. Since a Wick rotation from the Lorentzian to the Euclidean sector of quantum gravity is not feasible in general, the sign in the exponential of the path integral is not fixed a priori.

The transition for positive gravitational coupling has been observed before the one at negative $`\beta `$ . There, histograms, e.g. of $`A_t\delta _t`$ , clearly show a two-peak structure. The two phases coexist and tunneling from one phase to the other and back takes place. The system also exhibits a hysteresis; the transition occurs at a larger value of $`\beta `$ if the simulation is started from a configuration in the well-defined phase than it does if the calculation is started from a “frozen” configuration. Given all these pieces of evidence, we conjecture that the transition at $`\beta >0`$ is of first order.

To determine the order of the transition at negative $`\beta `$ we employed in Ref.  histogram techniques and used the Binder–Challa–Landau (BCL) cumulant criterion . The BCL cumulant is defined as $`B_L:=1\frac{E^4}{3E^2^2}`$, with $`E`$ being the action of the system under consideration. It was evaluated for $`A_t\delta _t`$ on different lattice sizes with $`L=3`$ to $`10`$ vertices per direction, simulating the system at several values of the bare coupling $`\beta `$ with high statistics. For the BCL cumulant a trend towards 2/3 was observed and all histograms showed a clear one-peak structure, cf. .

In the present Monte Carlo simulations we recorded for every run the time series of the energy density $`e=E/N_0`$ and the magnetization density $`m=_l\sigma _l/N_0`$, with the lattice size $`N_0=L^4`$. To obtain results for the various observables $`𝒪`$ at values of the bare gravitational coupling $`\beta `$ in an interval around the simulation point $`\beta _0`$, we applied the reweighting method .

With the help of the time series we can compute the specific heat,

$$C(\beta )=\beta ^2N_0(e^2e^2),$$
(4)
and the (finite lattice) susceptibility,

$$\chi (\beta )=N_0(m^2|m|^2).$$
(5)
Figures 1 and 2 show the finite-size scaling (FSS) of the maxima of the specific heat $`C_{\mathrm{max}}`$ and the susceptibility $`\chi _{\mathrm{max}}`$, respectively. While the behavior of $`C_{\mathrm{max}}`$ could still be explained as critical scaling with a negative exponent $`\alpha `$, the flattening of $`\chi _{\mathrm{max}}`$ is rather unusual for a second-order transition and should be taken as an indication for a cross-over regime.

Another feature of the system can be seen in Figs. 3 and 4 depicting the critical gravitational coupling, as determined from the specific-heat maxima $`C_{\mathrm{max}}`$ and the susceptibility maxima $`\chi _{\mathrm{max}}`$. In the case of a second-order transition the extrapolations of all pseudo-transition points lead to one infinite-volume critical value. This seems to be violated in the four-dimensional Discrete Regge model and thus favors the interpretation as a cross-over phenomenon over a true, thermodynamically defined phase transition.

## 3 CONCLUSION

In the present analysis of the scaling of the maxima of the specific heat and susceptibility we found evidence for a cross-over regime in the four-dimensional Discrete Regge model of quantum gravity in the negative coupling region. If this can be substantiated by further investigations and also for the Regge theory with continuous link lengths, the existence of a continuum limit at negative bare gravitational coupling is questionable. This is of major concern for a continuum theory of quantum gravity with matter fields .


text:
# ITP-Budapest 588, DESY 02-141 A scalable PC-based parallel computer for lattice QCD

## 1 Introduction

The most powerful computers for lattice gauge theory are industrial supercomputers or special purpose parallel computers (see e.g. ). It is more and more accepted that e.g. off-the-shelf PC systems can be used to build parallel computers for lattice QCD simulations .

Single PC hardware has excellent price/performance ratios. For recent review papers see .

We present our experiences and benchmark results on a scalable system, which uses nearest-neighbor communication through Gigabit Ethernet (GigE) cards. The communication is fast enough (consuming 40% of the total time in typical applications) and cheap enough (30% of the total price is spent on the communication). The system can sustain $``$100 Gflops on today’s medium/large lattices. Its price/performance ratio is below $1/Mflops for 32-bit applications (twice as much for 64-bit applications).

Already in 1999 a report was presented  on the PC-based parallel computer project at the Eötvös University, Budapest, Hungary. A machine was constructed with 32 PCs arranged in a three-dimensional 2$`\times `$4$`\times `$4 mesh. Each node had two special, purpose designed communication cards providing communication to the six neighbors. We used the “multimedia extension” instruction set of the processors to increase the performance. The communication bandwidth (16 Mbit/s) at that time was only enough for bosonic simulations. Here we report on a system which needs no hardware development and has two orders of magnitude larger communication bandwidth through GigE cards.

## 2 Hardware

Each node consists of an Intel KD850GB motherboard, Intel-P4-1.7GHz processor, 512 MB RDRAM, 100 Mbit Ethernet card, 20.4 GB IDE HDD and four SMC9452 GigE cards for the two dimensional inter-node communication. The 100 Mbit Ethernet cards with switches are used for job management only.

We have 137 PC nodes. A single machine is used for controlling. Two smaller clusters with 4 nodes are used for development. The remaining 128 machines can be used as one cluster with 128 nodes (or two clusters with 64 nodes or four clusters with 32 nodes) for mass production.

The best number of connected nodes for a given lattice can be determined by optimizing the surface to volume ratio of the local sub-lattice. Changing the node topology needs reconnecting cables, which can be done easily in a few minutes.

In April, 2002 the price of one node including the 100 Mbit Ethernet switches is $380 (see ). The four GigE cards with cables cost additional $160 for each node. The power consumption of the nodes (140W each) requires a cooling system which costs around $13 per node. Thus, the total node (PC+communication+cooling) price is about $553.

The key element of a cost effective design is an appropriate balance between communication and the performance of the nodes. We spent more than twice as much on the bare PC ($393 including cooling) than for the Gigabit communication ($160). These numbers are in strong contrast with Myrinet based PC systems, for which the high price of the Myrinet card exceeds the price of such a PC by a factor of two.

## 3 Software

The main operating system of the cluster is SuSE Linux 7.1, being installed on each node. Job management is done by the “main” computer through the 100 Mbit Ethernet network. We developed a simple job-management utility to distribute jobs on the nodes and collect the results.

To take advantage of the Gigabit communication from applications (e.g. C, C++ or Fortran), a simple C library was written using the standard Linux network interface. Currently, we are using a standard socket based communication with Transmission Control Protocol (TCP) widely used in network applications. The typical bandwidth that we can reach in QCD applications is around 400 Mbit/s in contrast to the theoretical 1000 Mbit/s. Writing a lower level driver for the gigabit cards which would increase the bandwidth is in progress. The functions of the library can be used to open and close the communcation channels and transfer data between the neighbors.

## 4 Performance

Figure 1 gives a summary of our benchmark runs. Single precision is used for local variables (gauge links and vectors, see ) and dot products are accumulated in double precision. Using double precision gauge links and vectors reduces the performance by approximately 50%.

For benchmarking we started from the MILC code . To increase the performance of the code we modified it by three different techniques.

First of all, we used the “multimedia extension” (sse) instructions of the Intel-P4 processors As we pointed out it in 1999  this capability can accelerate the processor by a large factor. We rewrote almost the whole conjugate gradient part of the program to assembly (including also loops over the lattice and not only elementary matrix operations). This way we obtained a speedup factor of $``$2. (see similar results of e.g. Ref. ).

Secondly, an important speedup was obtained by changing the data structure as suggested by S. Gottlieb . In the original “site major” concept all the physical variables of a given site are stored in one structure and the lattice is an array of these structures. Instead of this concept one should use “field major” variables. The set of a given type of variable of the different sites are collected and stored sequentially. This increases the number of cache hits. Similarly to Ref.  this change leads to a speedup factor of $``$2.

The third ingredient of our improvement was a better cache management by extensive use of the “prefetch” instruction.

One of the most obvious features of Figure 1 is the sharp drop of the performance when one turns on the communication. The most economic solution is to turn on the communication only if the memory is insufficient for the single-node mode. Smaller local lattices can be also studied by using the communication (for instance for thermalization or parameter tuning); however, in these cases the communication overhead increases somewhat. The performance as a function of the local lattice volume is shown on Fig. 2 (we also made a direct comparison between our architecture and a Myrinet system at DESY, Hamburg). Two orders of magnitude smaller local lattice means factor of $``$2 drop in the performance. Thus, for local lattices as small as $`8^3\times 4`$ the performance is still acceptable. This indicates clearly the scalability of this architecture. Note that the performance of a Myrinet based system is only 5-10% better than that of the GigE based system.

## 5 Conclusion

We reportde on the status of our PC-based parallel computer project for lattice QCD (for more details see ref. ). Nearest-neighbor communication is implemented in a two-dimensional mesh using Gigabit Ethernet cards. This architecture presents a good compromise between computation and communication. The saturated sustained performances on large lattices are around 0.5(0.8) GFlops/node for staggered (Wilson) fermions, which gives a price/performance ratio better than $1.0(0.7)/Mflop.

Acknowledgment. We thank F. Csikor and Z. Horváth for their continuous help. This work was supported by Hungarian Science Foundation Grants under Contract Nos. OTKA-T37615/T34980/T29803/M37071/OMFB1548/OMMU-708.