text:
# Are there really phase transitions in 1-d heat conduction models?

## Abstract

Recently, it has been claimed (O. V. Gendelman and A. V. Savin, Phys. Rev. Lett. 84, 2381 (2000); A.V.Savin and O.V.Gendelman, arXiv: cond-mat/0204631 (2002)) that two nonlinear classical 1-d lattice models show transitions, at finite temperatures, where the heat conduction changes from being finite to being infinite. These are the well known Frenkel-Kontorova (FK) model and a model for coupled rotators. For the FK model we give strong theoretical arguments why such a phase transition is not to be expected. For both models we show numerically that the effects observed by Gendelman et al. are not true phase transitions but are rather the expected cross-overs associated to the conductivity divergence as $`T0`$ and (for the FK model) $`T\mathrm{}`$.

Heat conduction in classical one-dimensional lattices has recently been investigated by many authors (see lepri  for a recent review). A large class of 1-d systems can be described by the general Hamiltonian

$$H=\underset{i=1}{\overset{N}{}}\left(\frac{p_i^2}{2m}+U(q_{i+1}q_i)+V(q_i)\right)$$
(1)
where $`N`$ is the system size, $`p_i`$ is the momentum of the $`i`$th particle, $`q_i`$ is its displacement from the equilibrium position, $`U`$ is the internal potential, and $`V`$ is an external potential. In the following we shall always use $`m=1`$, without loss of generality. Real isolated systems must have $`V=0`$, since any external potential would itself involve a “scaffolding” which is not rigid and would thus also contribute to heat conduction. Thus, for real systems the support of the external potential must be included in the description of the system, and $`V`$ has to be replaced by a contribution to the internal potential $`U`$. Nevertheless, we shall in the following keep the ansatz (1), understanding that we are calculating only part of the complete heat conduction when $`V0`$. In the following we shall also assume that the system has no frozen disorder.

It is well known that heat conduction is infinite, if all potentials are harmonic harmonic  or if the system is integrable. In this case phonons resp. solitons are not scattered. Thus they propagate ballistically, given a constant heat flux $`J`$ (independent of $`N`$) when a finite constant temperature difference $`\mathrm{\Delta }T`$ is applied to the two ends of a chain of length $`N`$. Thus formally, the conductivity $`\kappa =JN/\mathrm{\Delta }T`$ is proportional to $`N`$. For non-integrable models with an acoustic phonon branch, i.e. without an external potential, one expects ballistics transport in the infinite wave length limit (since phonon scattering in general decreases with energy), which leads to a power behavior $`\kappa N^\alpha `$ with $`0<\alpha <1`$ foot1 . This is the case, e.g., for the Fermi-Pasta-Ulam (FPU) model FPU\_1d  and for the diatomic Toda lattice Toda\_1d .

An exception to this seems to happen for coupled rotators R\_livi ; R\_savin ,

$$U(q_{i+1}q_i)=\mathrm{cos}(q_{i+1}q_i),V(q_i)=0,$$
(2)
which seem to have finite $`\kappa `$ R\_livi . This is explained R2\_savin  by the fact that single highly excited rotators essentially decouple from neighboring rotators, acting thus as barriers for the propagation of any phonons, even soft ones.

Soft acoustic phonons are essentially Goldstone modes due to Galilei invariance. Finite heat conduction (in $`d=1`$) is therefore expected when $`V0`$, since then translation invariance is broken and an acoustic branch does not exist in the phonon spectrum. This is e.g. the case for the Frenkel-Kontorova (FK) model

$$U(q_{i+1}q_i)=\frac{1}{2}(q_{i+1}q_i)^2,V(q_i)=\epsilon \mathrm{cos}q_i.$$
(3)
(notice that this is the commensurate FK model, where the harmonic part of the potential leads to the same particle distance in the ground state as the cosine potential; we shall only discuss this case in the following). In this model, all phonons have a finite mean free path, bounded from above by a finite constant which is independent of the wave number $`k`$, but which diverges for $`T0`$ and for $`T\mathrm{}`$. The latter follows from the fact that $`V`$ given by Eq.(3) effectively becomes negligible when $`T\mathrm{}`$, and effectively becomes a sum of harmonic potentials when $`T0`$ . Thus there is no ballistic transport, and no obvious mechanism which could lead to an infinite conductivity for any finite $`T`$, while one expects $`\kappa `$ to diverge when $`T0`$ or $`T\mathrm{}`$ FK\_bishop .

It was thus very surprising when Savin and Gendelman FK\_savin  claimed to have clear evidence for phase transitions in the FK model, at which the conductivity changed from finite to infinite. They claimed that $`\kappa `$ is finite only in an interval $`T_{c1}<T<T_{c2}`$, with $`T_{c1}`$ and $`T_{c2}`$ dependent on $`\epsilon `$, while $`\kappa =\mathrm{}`$ outside this interval.

Indeed, the same authors had also claimed that there is a phase transition in the rotator model R\_savin ; R2\_savin . There, the density of highly excited rotators should of course go to zero for $`T0`$. Thus one expects that soft phonons exist in this limit, and $`\kappa \mathrm{}`$ for $`T0`$. Instead of this, it was claimed in R\_savin ; R2\_savin  that $`\kappa =\mathrm{}`$ in an entire interval $`0TT_c`$.

In the present paper we want to test these claims by performing simulations on larger lattices and with higher precision than in R\_savin ; R2\_savin .

In order to mimic the simulations of FK\_savin ; R\_savin ; R2\_savin  as close as possible, we also used Langevin thermostats (we do not agree with these authors that Nosé-Hoover thermostats would be unsuitable, but we just don’t want any discussion about this point). More precisely, we simulated a chain of $`N_0+N+N_0`$ oscillators. The central $`N`$ oscillators follow their Hamiltonian equations of motion, while the outer $`2N_0`$ ones satisfy $`\ddot{q}_n=\frac{\varphi }{q_n}\gamma \dot{q}_n+\xi _n`$ and with $`\varphi =U+V`$ being the total potential, $`\xi _n`$ being white Gaussian noises, $`\xi _n(t)\xi _k(t^{})=2\gamma T_n\delta _{nk}\delta (tt^{})`$ and with $`T_n=T_{\mathrm{high}}`$ for $`N_0n0`$ and $`T_n=T_{\mathrm{low}}`$ for $`Nn<N+N_0`$. We used $`N_0=40`$ and $`\gamma =0.1`$, as in FK\_savin ; R\_savin ; R2\_savin . The temperature difference $`\mathrm{\Delta }TT_{\mathrm{high}}T_{\mathrm{low}}`$ was chosen between 10% and 20%.

For the integration we used a simple leap frog LeapFrog . On the one hand this is symplectic and thus more suited for the central region than, say, a Runge-Kutta integrator. On the other hand it should be more robust than higher order symplectic integrators in the boundary regions which are not Hamiltonian. Step size was 0.05, and total integration times were typically $`10^710^8`$ units (i.e. $`10^810^9`$ steps), with some runs going up to $`5\times 10^8`$ units. We checked that this was sufficient to reach a steady state and that the time-averaged heat flux $`J`$ was independent of the site.

We verified also that the temperature profiles were roughly linear in the central region $`0nN`$, but we could not verify the absence of temperature jumps at its boundaries claimed in FK\_savin ; R\_savin ; R2\_savin . More precisely, such jumps were absent only for very large lattices and small conductivities, i.e. if the heat flux was small. Otherwise, for small lattices and/or large conductivities, there were very large jumps, mostly located in the boundary regions $`40<n<0`$ and $`N<n<N+40`$ (see Fig. 1). Thus while the profile in the central region was essentially linear (in contrast to simulations with Nose-Hoover thermostats coupled to single particles, as e.g. in FPU\_yg ; FPU\_1d , it would be very wrong to estimate the conductivity simply by dividing the flux by the nominal imposed $`\mathrm{\Delta }T`$. It seems that this was done in several cases in FK\_savin ; R\_savin ; R2\_savin , which explains some – but not all – of the differences between their results and those of the present paper. In other cases the authors of FK\_savin ; R\_savin ; R2\_savin  must have taken into account boundary jumps, otherwise their results would disagree much more with ours than they actually do. Unless otherwise said, we will estimate $`\kappa `$ by dividing the flux by the temperature drop over the inner half of the central region.

Notice that a similar behavior was found also for the rotator model and for the discrete $`\varphi ^4`$ model. The latter is given by the Hamiltonian Fai4 zhao

$$H=\frac{1}{2}\underset{i}{}\left(p_i^2+(q_{i+1}q_i)^2+\alpha q_i^2+\frac{1}{2}\beta q_i^4\right).$$
(4)
Conductivities for the $`\varphi ^4`$ model with $`T=2`$ and $`\beta =1\alpha `$ are shown in Fig. 2 for various values of $`\beta [0,1]`$. The $`\varphi ^4`$ model is a prototype model with finite conductivities. Indeed we see that all measured values of $`\kappa `$ are not only finite but are independent of $`N`$. This would not have been the case if we had not taken the temperature jumps into account and would have used the nominal value of $`\mathrm{\Delta }T`$ when estimating $`\kappa `$ from $`J/\mathrm{\Delta }T`$. The values of $`\kappa `$ shown in Fig. 2 are in very good agreement with those of Fai4 zhao .

Frenkel-Kontorova model: In all simulations, $`N`$ ranged from 32 to 2048. We made simulations for $`\epsilon =1.0,3.0`$ and 10.0. According to FK\_savin , $`T_{c1}`$ is only weakly dependent on $`\epsilon `$: $`T_{c1}=2.6,2.3`$, and $`2.0`$ for $`\epsilon =1.0,3.0,10.0`$. On the other hand, $`T_{c2}`$ should strongly increase with $`\epsilon `$, $`T_{c2}=3.3,15`$, and 150 for the above three values.

Conductivities for $`\epsilon =1.0`$ are shown in Fig. 3. Obviously, they are finite for $`T=2`$ and $`T=8`$, showing that the temperature range with finite conductivities is underestimated in FK\_savin . For $`T=0.5`$ and $`T=32`$ we see a slow increase of $`\kappa `$ with $`N`$, over a wide range of the latter, but it seems to stop for the very largest lattices ($`N>1000`$). Finally, for $`T=0.125`$ and $`T=128`$ there is a slow increase for all $`N>300`$. The latter could be taken as an indication that $`\kappa `$ diverges for these temperatures, but we think that this would be wrong. On the one hand, the increase with $`N`$ is very slow, much slower than in Fig. 4 of FK\_savin . We would get a similarly fast increase as in FK\_savin  if we would use the nominal temperature difference, i.e. if we would disregard the jumps seen in Fig. 1. On the other hand, the data for $`0.5T32`$ show us that the saturation of $`\kappa `$ happens at larger and larger lattice sizes as we go away from the central energy region, just as we have expected. Thus we must expect in any case an increase of $`\kappa `$ for all reachable lattice sizes, and observing it does not give any relevant information.

Analogous results for $`\epsilon =3.0`$ are given in Fig. 4. There we show data for the temperature range $`[0.75,75]`$. This is again much larger than the range where convergent conductivities were found in FK\_savin . This time all curves become horizontal for large $`N`$, i.e. the conductivity is finite in the entire range. It of course depends strongly on $`T`$, see Fig. 5. It diverges both for $`T0`$ and for $`T\mathrm{}`$, since the problem effectively becomes harmonic in both limits.

Finally, our last simulations for the FK model, for $`\epsilon =10`$, are summarized in Fig. 6. There we only show results for low temperatures, $`0.5T2.75`$. Except for the last temperature, they are all in the regime where the authors of FK\_savin  have found divergent $`\kappa `$. In contrast, all our curves are either horizontal for all $`N`$ or become horizontal for large $`N`$, suggesting that $`\kappa `$ is finite for all finite $`T`$.

Rotator model: For the rotator model, we simulated larger systems, with $`N`$ ranging from 32 to 8192. Conductivities are plotted in Fig. 7 against $`N`$ for $`T=0.6,0.45,0.3,`$ and 0.2 (from bottom to top). According to R\_savin , the phase transition from a high-$`T`$ phase with finite conduction to a low-$`T`$ phase with infinite conduction occurs at some $`T_c`$ between $`0.2`$ and $`0.3`$. According to that, the lowest three curves in Fig. 7 should become flat for $`N\mathrm{}`$, while the uppermost should continue to grow. This is not what is found, although our values for $`T=0.2`$ and 0.3 agree numerically quite well with those of R\_savin . But while the curve for $`T=0.2`$ increases with the same average slope as in Fig. 2 of R\_savin , it is definitely S-shaped and stops to rise for the largest values of $`N`$.

Conclusion: In this paper we studied the size dependence of the effective finite size conductivity of nonlinear 1D lattices, as a function of temperature. We used straightforward but high statistics simulations to show that there are no indications of the phase transitions suggested in R\_savin ; R2\_savin ; FK\_savin  on the basis of similar simulations. For the FK model, this is in agreement with expectations, since phonons should have a finite free path in this model for all finite temperatures. For the rotator model it is less obvious. It suggests that the blocking of the propagation of soft phonons by localized excitations R\_savin  is effective at all finite temperatures. It becomes of course less and less important as $`T0`$, since the density of such excitation decreases exponentially with $`1/T`$. But it is present at all finite $`T`$, and it becomes dominant for $`N\mathrm{}`$.

P. G. wants to thank Roberto Livi and Antonio Politi for very helpful discussions.


format: text

text:
# References

Structure in the early afterglow lightcurve of the $`\gamma `$-ray burst of 29 March 2003

Makoto Uemura, Taichi Kato, Ryoko Ishioka, Hitoshi Yamaoka, Berto Monard, Daisaku Nogami§, Hiroyuki Maehara, Atsushi Sugie, and Susumu Takahashi
Department of Astronomy, Kyoto University, Kyoto 606-8502, Japan
Faculty of Science, Kyushu University, Fukuoka 810-8560
Bronberg Observatory, PO Box 11426, Tiegerpoort 0056, South Africa
§ Hida Observatory, Kyoto University, Kamitakara, Gifu 506-1314
VSOLJ, Namiki 1-13-4, Kawaguchi, Saitama 332-0034
Dynic Astronomical Observatory, Taga 283-1, Taga-chou, Inugami, Shiga, 522-0341 ,Japan

## Abstract

Gamma-ray bursts (GRBs) are energetic explosions that for 0.01–100 s are the brightest $`\gamma `$-ray sources in the sky.  Observations of the early evolution of afterglows we expected to provide clues about the nature of the bursts, but their rapid fading has hampered such studies; some recent rapid localizations   of bursts have improved the situation. Here we report on an early detection of the very bright afterglow of the burst of 29 March 2003 (GRB030329). Our data show that, even early in the aferglow phase, the light curve shows unexpectedly complicated structures superimposed on the fading background.

GRB030329 was detected by the HETE-2 satellite on 29 March 2003 at 11:37:14.67 UT. A prompt identification of its burst position, which was notified 0.051 d after the burst, enabled observers to make observations during an early phase of the afterglow of the GRB. A new very bright optical source of magnitude 12 was discovered within the error circle of the GRB position, and then continuously monitored by a number of observers.   Optical spectroscopic observations determined its redshift to be $`z=0.168`$, and using this, its isotropic released energy was calculated to be $`E_{\mathrm{iso}}=8\times 10^{51}\mathrm{erg}\mathrm{s}^1`$ (30–400 keV).  In the 30-years of studying GRBs, this GRB afterglow is the closest except for the suspected supernova associated with GRB980425.

We initiated time-series CCD photometric observations on 29 March 2003 at 12:53:41 UT, about 0.053 days after the burst, at Kyoto University, with 25- and 30-cm telescopes. After the observation at Kyoto, observations at Saitama, the Dynic Astronomical Observatory, and the Bronberg Observatory started with 20-cm, 60-cm, and 30-cm telescopes, respectively. Follow-up observations were performed on 30 March 30, 3, 5, and 6 April at Kyoto, on 30 March at Hida, and on 30 and 31 March at Bronberg. A dark-current image was subtracted from obtained CCD images, and then flat-fielding was performed. We calculated the magnitudes of the afterglow with neighbour comparison stars (GSC1434.239, GSC1434.129, and GCS1434.192), whose constancy was better than 0.04 mag during our observation. Our unfiltered CCD observation yields a magnitude system near that of the $`R_\mathrm{c}`$-system, since the sensitivity peak of the camera is near to the peak of the $`R_\mathrm{c}`$-system, and the spectra of early afterglows have a smooth continuum without strong emission lines or absorption edges in the optical range. The difference between our unfiltered CCD and the $`R_\mathrm{c}`$-system is less than 0.02 mag when the spectrum is described with a simple power law ($`f(\nu )\nu ^p`$) with index, $`2.0<p<0.3`$, as expected in GRB optical afterglows.

Our obtained light curve is shown in Fig. 1. In Fig. 1a, the abscissa and the ordinate denote the time from the burst in days and the $`R_\mathrm{c}`$-magnitude, respectively. To date, the light curve of afterglows, in particular during the early phase of $`0.1`$ days, has been sampled only sparsely, which has been insufficient to observationally verify theoretical models. As can be seen in this figure, we succeeded in obtaining a completely continuous light curve from 0.053 to 0.500 days, or for 11 hours in GRB030329, which reveals the presence of unambiguous detailed structures in the light curve of early afterglows.

The light curve of the GRB030329 afterglow exhibits clear, and repeating deviations from a canonical power-law model. The dotted line in Fig. 1a is the best fitted power-law model using our light curve and the GRB Coordinates Network data between 0.5 and 10 days. The average decay index, $`\alpha `$ ($`ft^\alpha `$) is calculated to be $`1.115\pm 0.006`$. This is a standard value for early afterglows of GRBs. As can be seen in the figure, the fading of the afterglow cannot be described only by this simple model. The residual from the single power-law model is shown in Fig. 1b. During the first day, the light curve can be described with four segments, that is, two sets of ascending and descending branches of bumps. Their transitions were not smooth, but occurred in short timespans. A relatively rapid fading presumably follows these bumps. At the end of the first day, the afterglow ceased fading, and then, experienced a rebrightening. While they were sampled only partially, the light curve shows that the bumps after 1 d had a rising time-scale of 0.1 day. It is interesting to note that an early bump around 0.4 day in Fig. 1 had a timescale of the same order as the later bumps. This similarity might imply that these large (amplitude $`0.4`$$`0.5`$ mag) bumps were of the same nature.

The afterglow of GRBs is now widely believed to be synchrotron emission from a forward shock region in an expanding jet colliding with ambient medium, whereas the GRB itself is from an internal shock region between shells.   In a number of afterglows, it is well known that light curves show a sudden increase of the fading rate due to the jet geometry, typically a few days after the burst.  In the case of GRB030329, this jet break is expected to appear later than 60 days from its isotropic energy and redshift.  This indicates that the jet break is expected not to appear in Fig. 1, and hence, it is difficult to explain the variation in Fig. 1 with the jet break. As well as the jet break, another well-known phenomenon, the supernova bump is evidently difficult to explain the multiple modulations. In GRB030329, the feature of a supernova spectrum appeared 7 days after the burst. The optical flux before 7 days was hence dominated by the afterglow emission itself.

In some GRB afterglows, similar deviations from a general fading component have been reported.   The recent object, GRB021004 is the first object exhibiting multiple fluctuations 1 day after the burst. In GRB030329, we found that the large amplitude bumps appeared even in the first day after the burst. The mechanism of these bumps has not been established yet, but the continuous light curve of the bumps of GRB030329 should help us to understand them.

Acknowledgements

This work is partly supported by a grant-in aid from the Japanese Ministry of Education, Culture, Sports, Science and Technology. Part of this work is supported by a Research Fellowship of the Japan Society for the Promotion of Science for Young Scientists (MU and RI).

Figure 1. Light curves of the optical afterglow of GRB030329. The abscissa denotes the time from the burst (29 March 2003, 11:37:14.67 UT in days, shown in logarithmic scale. The ordinate denotes the $`R_\mathrm{c}`$-magnitude in a and the residual magnitude from a power-law decay in b. The dotted line in a is the best-fitted power-law model, which was used to calculate the residual magnitude in b. Our observations are indicated with filled circles with error bars, indicating standard errors. The exposure times of each frame were 30, 40, 60, 30-60and 45 s for Kyoto, Saitama, the Dynic Astronomical Observatory, the Hida Observatory and the Bronberg observations, respectively. The points in the figure are binned points of the observations with $`\mathrm{\Delta }\mathrm{log}t_{\mathrm{day}}=0.002`$. Using a smoothly broken power-law model ($`f(t)=(f_1(t)^n+f_2(t)^n)^{1/n}`$ with $`f_i(t)=k_it^{\alpha _i}`$), we determined decay indices ($`\alpha _i`$) to be $`0.74\pm 0.02(0.053<t<0.085)`$, $`0.95\pm 0.01(0.085<t0.163)`$, $`0.65\pm 0.04(0.163<t<0.227)`$, and $`1.16\pm 0.01(0.227<t<0.492)`$. The break times are $`0.085\pm 0.028`$, $`0.163\pm 0.060`$, and $`0.227\pm 0.043`$. Note that the 0.16-day break time was calculated with a broken power-law model without the smoothness parameter. The smoothness parameters, $`n`$, are $`94\pm 51`$ for the earlier break and $`69\pm 89`$ for the later break, indicating rapid state transitions.


text:
# Scalar field perturbations in Fresh inflation

## Abstract

The model of fresh inflation with increasing cosmological parameter provides sufficient e-folds to solve the flatness/horizon problem and the density fluctuations agree with experimental values. In this model the temperature increases during fresh inflation and reach its maximum value when inflation ends where the $`N`$-number of e-folds is $`N(t_e)60`$. The most important characteristic of this model is that provides a natural transition between the end of inflation and the epoch when the universe is radiation dominated.

Pacs numbers: 98.80.Cq

Inflation is one of the most reliable concepts in modern cosmology. The first model of inflation was proposed by A. Starobinsky. A much simpler inflationary model with a clear motivation was developed by A. Guth, in order to solve some of the shortcomings of the big bang theory, and in particular, to explain the extraordinary homogeneity of the observable universe. However, the universe after inflation in this scenario becomes very inhomogeneous. These problems were sorted out by A. Linde in 1983 with the introduction of chaotic inflation. In this scenario inflation can occur in theories with potentials such as $`V(\varphi )\varphi ^n`$. It may begin in the absence of thermal equilibrium in the early universe, and may start even at the Planckian density, in which case the problem of initial conditions for inflation can be easily solved. Although a justification from first principles for dissipative effects has not been firmly achieved in the framework of inflation, such effects should not be ruled out on the basis of readiness alone. Much work can be done on phenomenological grounds as, for instance, by applying nonequilibrium thermodynamic techniques to the problem or even studying particular models with dissipation. An interesting example of the latter case is the warm inflationary picture. As in new inflation, a phase transition driving the universe to an inflationary period dominated by the scalar field potential is assumed. A standard phenomenological frictionlike term $`\mathrm{\Gamma }\dot{\varphi }`$ is inserted into the scalar field equation of motion to represent a continuous energy transfered from $`\varphi `$ to the radiation field. This persistent thermal contact during inflation is so finely adjusted that the scalar field evolves all the time in a damped regime generating an isothermal expansion. As a consequence, the subsequent reheating mechanism is not needed and thermal fluctuations produce the primordial spectrum of density perturbations. More recently was demonstrated that isentropic and warm pictures are just extremme cases of an infinite two-parametric family of possible inflationary scenarios. Very recently, a dissipative mechanism which emerges in generic interacting quantum field systems and leads to robust warm inflation was developed.

Two years ago, a new scenario called fresh inflation was introduced. It can be viewed as a unification of both chaotic and warm inflation scenarios with constant and increasing $`F`$-cosmolgical parameter. As in chaotic inflation the universe begins from an unstable primordial matter field perturbation with energy density nearly $`\mathrm{M}_p^4`$ ($`\mathrm{M}_p=1.2\times GeV`$ is the Planckian mass) and chaotic initial conditions. Furthermore, initially the universe there is not thermalized so that the radiation energy density when inflation starts is zero \[$`\rho _r(t=t_0)=0`$\]. It must to be point out that quantum field theory first principle calculations of radiation generated from a zero temperature state has been done in, but this formalism is realized in an non-expanding Minkowski spacetime.

As initial time we understand the Planckian time $`G^{1/2}`$, where $`G`$ is the gravitational constant. Later, the universe will describe a second-order phase transition. In other words, the inflaton rolls down towards the minimum energetic configuration. Particle production and thermalization occur together during the rapid expansion of the universe, so that the radiation energy density grows during fresh inflation ($`\dot{\rho }_r>0`$). The interaction between the inflation field and the particles produced during inflation provides slow-rolling of the inflaton field. So, in the fresh inflationary model (also in warm inflation), the slow-roll conditions are physically well justified. The decay width of the $`\varphi `$-field grows with time, so when the inflaton approaches to the minimum of the potential there is no oscillation around the minimum energetic configuration. Hence, the reheating period does not happen in fresh inflation. This model attempts to build a bridge between the standard and warm inflationary models, beginning from chaotic initial conditions which provides naturality. In fresh inflation with constant cosmological parameter the universe can be seen as inflating in a $`4D`$ Friedmann-Robertson-Walker (FRW) metric embedding in a $`5D`$ metric in apparent vacuum when the fifth coordinate of the metric is non-compact and remains constant. Furthermore, can be demonstrated the possibility that violation of baryon number conservation and ultralight particle creation can occur during the period out-of-equilibrium, i.e., before the thermal equilibrium to be restored at the end of fresh inflation. Gauge-invariant metric fluctuations in a fresh inflationary model with increasing cosmological parameter was considered in. The main aim of this paper is the study of the evolution of density perturbations in such that model.

We describe fresh inflation with a Lagrangian for a $`\varphi `$-scalar field minimally coupled to gravity, which also interacts with another $`\psi `$-scalar field by means of $`_{int}=\mathrm{g}^2\varphi ^2\psi ^2`$, is

$$=\sqrt{g}\left[\frac{R}{16\pi G}\frac{1}{2}\left(\varphi \right)^2V(\varphi )+_{int}\right],$$
(1)
where $`R=6(a\ddot{a}+\dot{a}^2)/a^2`$ is the scalar curvature, $`a`$ is the scale factor of the universe and $`g`$ is the determinant of the metric tensor $`g^{\mu \nu }`$ with $`\mu ,\nu =0,1,2,3`$. In this paper I consider a FRW metric for a spatially flat, isotropic and homogeneous universe described by the line element $`ds^2=dt^2+a^2dr^2`$. If $`\delta =\dot{\rho }_r+4H\rho _r`$ describes the interaction between the inflaton and the bath for a $`\gamma =4/3`$-fluid which expands with a Hubble parameter $`H=\dot{a}/a`$ and radiation energy density $`\rho _r`$, hence the equations of motion for $`\varphi `$ and radiation energy density are

$$\ddot{\varphi }+3H\dot{\varphi }+V^{}(\varphi )+\frac{\delta }{\dot{\varphi }}=0,\dot{\rho }_r+4H\rho _r\delta =0.$$
(2)
Here, $`\delta =\mathrm{\Gamma }(\theta )\dot{\varphi }^2`$ describes a Yukawa interaction and the $`\varphi `$-decay width is $`\mathrm{\Gamma }(\theta )=[g_{\mathrm{eff}}^4/(192\pi )]\theta `$. Furthermore, $`\theta \rho _r^{1/4}`$ is the temperature of the bath. The cosmological parameter $`F=(p_t+\rho _t)/\rho _t`$ describes the evolution of the universe during inflation

$$F=\frac{2\dot{H}}{3H^2}=\frac{\dot{\varphi }^2+\frac{4}{3}\rho _r}{\rho _r+\frac{\dot{\varphi }^2}{2}+V},$$
(3)
where the total pressure and energy density are given respectively by

$$p_t=\frac{\dot{\varphi }^2}{2}+\frac{\rho _r}{3}V(\varphi ),\rho _t=\rho _r+\frac{\dot{\varphi }^2}{2}+V(\varphi ).$$
(4)
In the previous works only was considered the case where the cosmological parameter $`F`$ is a constant. However, as we can see in eq. (3), during inflation the potential energy density decreases, so that the radiation energy density becomes more important in $`F`$. This means that $`F`$ must be increasing during fresh inflation, but of course, always remaining below $`4/3`$, which corresponds to a radiation dominated universe. We can write $`\rho _r`$ and $`V(\varphi )`$ as a function of $`\varphi `$

$$\rho _r=\left(\frac{3F}{43F}\right)V\frac{27}{8}\left(\frac{H^2}{H^{}}\right)^2\frac{F^2(2F)}{(43F)},V(\varphi )=\frac{3}{8\pi G}\left[\left(\frac{43F}{4}\right)H^2+\frac{3\pi G}{2}F^2\left(\frac{H^2}{H^{}}\right)^2\right],$$
(5)
where $`F`$ is a function of $`\varphi `$ and the time evolution of $`\varphi `$ is described by the equation

$$\dot{\varphi }=\frac{3H^2}{2H^{}}F(\varphi ).$$
(6)
We consider in the second equation of (5) the potential $`V(\varphi )=[^2(0)/2]\varphi ^2+[\lambda ^2/4]\varphi ^4`$, where $`\mathrm{G}=M_p^2`$ is the gravitational constant and $`M_p=\mathrm{1.2\; 10}^{19}\mathrm{GeV}`$ is the Planckian mass. The inflaton field is really an effective field described by $`\varphi =(\varphi _i\varphi _i)^{1/2}`$. Furthermore, $`^2(0)`$ is given by $`_0^2`$ plus renormalization counterterms in the initial potential $`\frac{1}{2}_0^2(\varphi _i\varphi _i)+\frac{\lambda ^2}{4}(\varphi _i\varphi _i)^2`$, the effective potential is $`V_{\mathrm{eff}}(\varphi ,\theta )=[^2(\theta )/2]\varphi ^2+[\lambda ^2/4]\varphi ^4`$. Here, $`\theta `$ is the temperature and $`^2(\theta )=^2(0)+\frac{(n+2)}{12}\lambda ^2\theta ^2`$, such that $`V_{\mathrm{eff}}(\varphi ,\theta )=V(\varphi )+\rho _r(\theta ,\varphi )`$. The temperature increases with the expansion of the universe because the inflaton transfers radiation energy density to the bath with a rate larger than the expansion of the universe. So, the number of created particles $`n`$ \[for $`\rho _r=(\pi ^2/30)g_{\mathrm{eff}}\theta ^4=\left(^2(\theta )^2(0)\right)\varphi ^2`$\], is given by

$$(n+2)=\frac{2\pi ^2}{5\lambda ^2}g_{\mathrm{eff}}\frac{\theta ^2}{\varphi ^2},$$
(7)
where $`g_{\mathrm{eff}}`$ denotes the effective degrees of freedom of the particles and it is assumed that $`\varphi `$ has no self-interaction. Note that eq. (7) only depends on the choice of $`V_{\mathrm{eff}}`$ and $`^2(\theta )`$, and not on $`H`$ and $`F`$. (For a more detailed explanation of eq. (7) the reader can see the papers .) The density fluctuacions are given by

$$\frac{\delta \rho _t}{\rho _t}\mathrm{\Delta }(\varphi )\delta \varphi ^2^{1/2},$$
(8)
where $`\delta \varphi ^2`$ are the squared fluctuations of the inflaton field and $`\mathrm{\Delta }(\varphi )=\frac{V_{eff}^{}}{\dot{\varphi }^2+4/3\rho _r}`$. The notation $`<\mathrm{}>`$ denotes the averaging with respect to a gaussian distribution. Since the temperature of the thermal bath is $`\theta `$, we can consider the Langevin equation for the field fluctuations $`\delta \varphi `$, with a white and gaussian noise $`\kappa (t)`$, due to the stochastic interaction of the field $`\varphi `$ with the thermalyzed environment

$$\dot{\delta }\varphi \frac{H^2}{3H+\mathrm{\Gamma }}\delta \varphi +\kappa (t).$$
(9)
Here, we have made use of the fluctuation dissipation theorem, and the fact that the slow-roll condition holds. This fact implies that $`H^2V^{\prime \prime }`$. The diffusion coefficient for the fluctuations with a wavelenght given by the Hubble radius $`H^1`$, is (see the first paper in ) $`D=\frac{3}{2\pi }\frac{\theta H^3}{3H+\mathrm{\Gamma }}`$, such that the squared mean fluctuations of the field are

$$\delta \varphi ^2\sqrt{\mathrm{\Gamma }H}\theta .$$
(10)
This relation is only valid when the system is thermalyzed, i.e., in the last stages of fresh inflation for which the relation $`\mathrm{\Gamma }>H`$ holds.

As an example we consider the particular case where the Hubble and cosmological parameters are $`H(\varphi )=A\varphi ^2`$ and $`F(\varphi )=B\varphi ^2`$, being $`A`$ and $`B`$ $`\mathrm{G}^{1/2}`$ and $`\mathrm{G}^1`$ dimensional constants which can be obtained by replacing $`H`$ and $`F`$ in eq. (5)

$$A^2=\frac{\lambda ^2}{4}\frac{8\pi G}{3},B=\frac{3\lambda +\sqrt{9\lambda ^2+48\pi G^2(0)}}{3\lambda \pi G}.$$
(11)
From eq. (6) we obtain the time evolution for the inflaton field $`\varphi (t)=\varphi ^{(s)}e^{\frac{3AB}{2}t}`$, where $`\varphi ^{(s)}`$ is its initial value. It must be sufficiently large to assures at least, the $`5060`$ e-folds needed during inflation before transcurred ($`10^{10}10^{12}`$) Planckian times. Since $`H=\dot{a}/a`$, the time evolution of the scale factor will be $`a=a_0e^{\frac{\varphi ^2(t)}{3B}}`$, which increases with time due to the decreasing of $`\varphi (t)`$. Furthermore, the temperature, written as a function of the field is obtained from the second equation in (2)

$$\theta (\varphi )=\left\{4\pi \varphi (8\varphi 3B)\left[16^2(0)+\varphi ^2\left(8\lambda ^218A^2B\right)+9A^2B^2\right]\right\}\left[BAg_{\mathrm{eff}}^4\left(4\varphi ^23B\right)\right]^1,$$
(12)
which is zero when fresh inflation starts. The $`\varphi `$-value at this time is $`\varphi ^{(s)}=3B/8`$. The function $`\mathrm{\Delta }(\varphi )`$ in eq. (8) is

$$\mathrm{\Delta }(\varphi )=\frac{3\left(\pi GB2\right)}{4\varphi }.$$
(13)
Hence, from eqs. (13) and (10) one obtains the density fluctuations $`\frac{<\delta \rho _t>}{\rho _t}`$ in eq. (8)

$$\frac{\delta \rho _t}{\rho _t}\frac{3(\pi BG2)}{4\varphi }\left[\frac{\pi \varphi ^3(3B8\varphi )}{g_{\mathrm{eff}}^4B(4\varphi ^23B)}\right]^{1/2}\left[2\varphi ^2(9A^2B4\lambda ^2)16^2(0)9A^2B^2\right]^{1/2}.$$
(14)
Figures (1), (2) and (3) show respectively the evolution of the temperature $`\theta (t)`$, density fluctuations $`<\delta \rho _t>/\rho _t`$ (as a function of the number of $`e`$-folds) and the number of created particles in eq. (7) (as a function of time). These graphics were made using the values $`(0)=1.5\times 10^9\mathrm{G}^{1/2}`$, $`\lambda =4\times 10^{15}`$ and $`g_{\mathrm{eff}}=25`$. With these parameter values we obtain respectively $`A=0.578\times 10^{14}\mathrm{G}^{1/2}`$ and $`B=488726.728\mathrm{G}^1`$ in eqs. (11). Notice that fresh inflation ends at $`t_e2\times 10^8\mathrm{G}^{1/2}`$, when $`\theta (t_e)6\times 10^8\mathrm{G}^{1/2}`$ \[see Fig. (1)\] and $`<\delta \rho _t>/\rho _t|_{N(t_e)}10^5`$ \[see Fig. (2)\] reach their maximum values. The $`\varphi `$-value when fresh inflation ends is given by $`\varphi ^{(e)}\varphi (t_e)8\times 10^5\mathrm{G}^{1/2}`$. Note that $`N(t_e)60`$ at this moment so that the model here studied solve the flatness/horizon problems. Furthermore, $`\varphi `$ takes transplanckian values during fresh inflation with increasing cosmological parameter. It could be a problem for this model because there is no clear the validity of conventional quantum field theory in the transplanckian regime.

Finally, Fig. (3) shows the evolution of the number of created particles, which grows to reach its asymptotic value $`n9\times 10^6`$, at the end of the inflationary period.

To summarize, the conditions on $`\rho _r`$ lead to two different types of thermodynamic regimes of inflation. At the begining $`\rho _r0`$ and the expansion is isentropic during inflation, but the temperature increases because $`\mathrm{\Gamma }`$ is sufficietly strong to make $`\dot{\rho }_r>0`$. This is the case of fresh inflation, in which this first regime is followed by a non-isentropic period in which $`\rho _v>\rho _r>0`$, so that the temperature of the universe may still be sizable. In this regime, conversion of vacuum energy into radiation energy occurs during the inflationary period and the inflationary regime smoothly terminates into a radiation dominated regime without an intermediate reheating period. The toy model here studied shows that fresh inflation can be a feasible alternative model to standard inflation. The most important characteristic of this model is that shows a natural transition between the end of inflation and the epoch when the universe is radiation dominated. To conclude, a more exaustive calculation in the inflation/radiation-dominated interphase should include cold dark matter (CDM) in this model of fresh inflation with increasing cosmological parameter.

Fig. 1: Evolution of the temperature $`\theta (t)`$ during inflation. The maximum, with parameters $`(0)=1.5\times 10^9\mathrm{G}^{1/2}`$, $`\lambda =4\times 10^{15}`$ and $`g_{\mathrm{eff}}=25`$, occurs at the end of fresh inflation (i.e., around $`t_e2\times 10^8\mathrm{G}^{1/2}`$).

Fig. 2: Evolution of density fluctuations $`\frac{<\delta \rho _t>}{\rho _t}`$ as a function of the number of $`e`$-folds, where $`N(t_e)60`$. The maximum at the end of inflation is of the order of $`10^5`$.

Fig. 3: Temporal evolution of $`n`$ which shows that the number of created particles grows with time to reach the value of $`n9\times 10^6`$ at the end of inflation.


text:
# Thermodynamic properties of the two-dimensional 𝑆=1/2 Heisenberg antiferromagnet coupled to bond phonons

## I Introduction

During the last years, there has been considerable interest in low-dimensional spin systems with spin phonon coupling. While one-dimensional (1D) models including this mechanism have been studied extensively, little is known about two-dimensional (2D) systems with spin phonon coupling.

In 1D systems, the mechanism of spin phonon coupling is closely connected with the phenomenon of the spin Peierls transition. Theoretical understanding of this phase transition goes back to the work of Pytte,Pytte (1974) who showed that a three-dimensional system consisting of uniform antiferromagnetic $`S=1/2`$ Heisenberg chains is unstable towards dimerization of the chains if coupled to three-dimensional lattice vibrations. Additionally, he showed that in the adiabatic limit of small phonon frequencies the dimerized phase can be described by a statically dimerized spin model with temperature-dependent dimerization. Cross and FisherCross and Fisher (1979) improved on the calculation of Ref. Pytte, 1974 by treating the spin part of the Hamiltonian in continuum field theory.Luther and Peschel (1975) Their calculation yielded a convincing description of the phonon softening in the limit of small spin phonon couplings. In the case of CuGeO3, however, which was the first inorganic spin Peierls compound discovered,Hase et al. (1993) this treatment is not sufficient. For this reason various 1D dynamical models have been investigated, taking into account the coupling of magnetic and phononic degrees of freedom explicitly. Most publications on this issue refer to two different models, both showing a quantum phase transition between a dimerized and a Néel ordered phase. While in the so-called difference coupling modelUhrig (1998); Bursill et al. (1999); Weiße et al. (1999); Sun et al. (2000); Raas et al. (2001) the magnetic interaction between nearest neighbours depends on the distance between neighboured sites, the bond coupling modelSandvik et al. (1997); Weiße et al. (1999); Kühne and Löw (1999); Sandvik and Campbell (1999); Trebst et al. (2001); Raas et al. (2002) is considered to be more realistic for describing the spin phonon coupling mechanism in CuGeO3.Werner et al. (1999); Feldkemper and Weber (2000)

In the case of the 2D $`S=1/2`$ Heisenberg antiferromagnet with spin phonon coupling discussions of statically dimerized models are found in the literature.Tang and Hirsch (1988); Koga et al. (1998, 1999); Al-Omari (2000); Sirker et al. (2002) As in one dimension, these models are thought to describe the dimerized phase of dynamical models in the adiabatic limit. In contrast to the 1D case though it is not clear how to place alternating magnetic couplings in both spatial directions on a square lattice. By comparing ground state energies of three statically dimerized models Sirker *et. al.*Sirker et al. (2002) conclude that a stair-like distortion of the lattice is the energetically favoured dimerization pattern, contradicting an older result by Tang and HirschTang and Hirsch (1988) who find a plaquette-like distortion.

The aim of this paper is to investigate thermodynamic properties of the 2D Heisenberg model coupled to bond phonons. Such a model takes into account the elastic energy due to lattice distortions which is not included in statically dimerized models. Furthermore, the whole range of phonon frequencies is accessible in our treatment.

The paper is organized as follows. In Sect. II we introduce the model Hamiltonian and describe our quantum Monte Carlo algorithm. In Sect. III we discuss the temperature dependence of the magnetic susceptibility, mean phonon occupation numbers and the specific heat. An analysis of the spin correlation function is found in Sect. IV. Section V concludes with a summary.

## II Model Hamiltonian and Quantum Monte Carlo method

The model we consider is a generalization of the 1D bond coupling model from Refs. Kühne and Löw, 1999; Raas et al., 2002. The Hamiltonian readsAits (2002)

$`H`$ $`={\displaystyle \frac{J}{2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}(\stackrel{}{\sigma }_{ij}\stackrel{}{\sigma }_{i+1,j}1)(1+g[a_{ij}+a_{ij}^{}])`$
$`+{\displaystyle \frac{J}{2}}{\displaystyle \underset{i,j=1}{\overset{N}{}}}(\stackrel{}{\sigma }_{ij}\stackrel{}{\sigma }_{i,j+1}1)(1+g[b_{ij}+b_{ij}^{}])`$ (1)
$`+\omega {\displaystyle \underset{i,j=1}{\overset{N}{}}}(a_{ij}^{}a_{ij}+b_{ij}^{}b_{ij}).`$
Here $`\stackrel{}{\sigma }_{ij}`$ denote Pauli spin operators at lattice site $`(i,j)`$ on a square lattice, while $`a_{ij}`$ and $`a_{ij}^{}`$ ($`b_{ij}`$ and $`b_{ij}^{}`$) are phonon annihilation and creation operators on the bond between site $`(i,j)`$ and site $`(i+1,j)`$ (between sites $`(i,j)`$ and $`(i,j+1)`$). Note that we assume periodic boundary conditions.

By shifting the phonon operators according to

$$a_{ij}a_{ij}+\frac{gJ}{2\omega },b_{ij}b_{ij}+\frac{gJ}{2\omega }$$
(2)
and neglecting a constant energy contribution model (II) can be mapped onto the phenomenologically more realistic Hamiltonian

$`\stackrel{~}{H}`$ $`={\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}(J^{}+g^{}[a_{ij}+a_{ij}^{}])\stackrel{}{\sigma }_{ij}\stackrel{}{\sigma }_{i+1,j}`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}(J^{}+g^{}[b_{ij}+b_{ij}^{}])\stackrel{}{\sigma }_{ij}\stackrel{}{\sigma }_{i,j+1}`$ (3)
$`+\omega {\displaystyle \underset{ij}{}}(a_{ij}^{}a_{ij}+b_{ij}^{}b_{ij}),`$
with rescaled coupling constants $`J^{}J(1+g^2J/\omega )`$ and $`g^{}gJ`$. Note that Eq. (II) differs from (II) by the absence of the (unphysical) static terms $`Jg/2(a_{ij}^{}+a_{ij})`$,$`Jg/2(b_{ij}^{}+b_{ij})`$. Model (II), however, cannot be analyzed directly in a quantum Monte Carlo study because a sign problem occurs. For this reason we choose Hamiltonian (II) as a starting point for our analysis and shift the numerical results if needed.

To study the properties of model (II), we developed a quantum Monte CarloSuzuki et al. (1977) algorithm similar to the algorithm described in Ref. Kühne and Löw, 1999. First, the partition function of the 2D quantum system (II) is mapped onto a three-dimensional classical system. Technically, this is done by means of a Trotter Suzuki decomposition.Kühne and Löw (1999) We then apply an update procedure which treats spin and phononic degrees of freedom separately. For the spin updates, we make use of a modified loop algorithmEvertz et al. (1993); Evertz (1997) for quantum spin systems. The main advantage of the loop algorithm is that it allows global spin updates, substantially reducing autocorrelation times. Furthermore, so-called improved estimators can be used in the evaluation of the magnetic susceptibility and spin correlations. To modify the phonon occupation numbers we apply local heat bath updates. By building clusters of phonons in imaginary time direction we extended the algorithm from Ref. Kühne and Löw, 1999, diminishing autocorrelation effects even more. Obviously the detailed balance condition is fulfilled for both steps separately and thus for the whole procedure.

For the phonon updates we had to introduce a cutoff, allowing occupation numbers up to 40 phonons per bond. The effect of such a truncation of the Hilbert space is negligible if the measured mean phonon occupation numbers are more than an order of magnitude smaller than the cutoff. In Sect. III.2 we show explicitly that this condition is fulfilled. In order to take into account the high dimension of the phonon subspace of the Hilbert space, we employed the importance sampling technique and made 30 phonon updates per spin update, using only the last configuration for the evaluation of expectation values. Additionally, for each temperature the first 25% of the sweeps were skipped for thermalization.

For Monte Carlo simulations based on a Trotter Suzuki decomposition the estimates of thermodynamic quantities depend on the inverse Trotter number squared.Suzuki (1985) In the following sections we give the explicit value for the Trotter number $`M`$. With the values for $`M`$ chosen, we find the statistical fluctuations of our results larger than the effect of the finite Trotter number.

Before discussing the results in detail we add one further remark concerning statistical errors. In our calculations we neglected autocorrelation effects and – as a rough estimate – calculated root-mean-squared errors only. If no error bars in the plots of this paper are shown, the errors are smaller than the symbol size used.

## III Thermodynamic properties

In this section we discuss the finite temperature properties of model (II). We expect that the knowledge of how a non-vanishing spin phonon coupling influences the thermodynamic properties will be of importance for the interpretation of experiments. This might be of particular interest for substances which display e. g. acoustic anomalies or for which it is known that the exchange integral depends sensitively on the positions of the ions.

Here, we confine ourselves to temperatures $`0.5JT4J`$. In this temperature range, we find that the dependence of measured quantities on the system size is negligible if we consider linear system sizes $`N12`$. All results presented in this section were calculated on a lattice with 12$`\times `$12 sites, providing statements about system properties in the thermodynamic limit. If not stated differently, at each temperature 105 spin updates were executed. For the Trotter number a value of $`M=80`$ was chosen.

### III.1 Magnetic susceptibility

We start with a discussion of the magnetic susceptibility per site $`\chi `$ for vanishing magnetic fields. Figure 1 shows the dependence of the susceptibility on the spin phonon coupling $`g`$ for fixed phonon frequency $`\omega `$, Fig. 2 the dependence on $`\omega `$ for fixed $`g`$. In both figures Monte Carlo results for the susceptibility of the 2D Heisenberg model are included.

The results can be summarized as follows. For fixed $`\omega `$, the overall height of the susceptibility is diminished with increasing spin phonon coupling. As in the 1D case,Kühne and Löw (1999) a large spin phonon coupling tends to reduce the magnetic response of the system. On the other hand, for fixed $`g`$ the susceptibility is growing with increasing phonon frequency. In the antiadiabatic limit phononic degrees of freedom are suppressed, yielding the results of the Heisenberg model. The curves show a broad maximum which is typical for antiferromagnetic spin models. This maximum is shifted to higher temperatures with increasing $`gJ/\omega `$.

We find that both the shift of the position of the maximum and the reduction of magnetic response with increasing $`\frac{gJ}{\omega }`$ is mainly due to the static terms in (II). In units of the rescaled magnetic coupling $`J^{}`$ from the transformed Hamiltonian (II) the shift of the maximum is reduced significantly (see Fig. 3), and the reduction of the magnetic response due to the spin phonon coupling is not very strong. The same behaviour has already been reported for the 1D bond coupling model in Ref. Raas et al., 2002.

For comparison between $`d=1`$ and $`d=2`$ we return to model (II) without phonon shift. Figure 4 shows the susceptibilities for $`g=0.5`$, $`\omega =2J`$. Compared to the 1D case, the overall height of the susceptibility of the 2D model is diminished. This effect is explained by the larger coordination number in the 2D case, reducing the response of the system to an external magnetic field.

Qualitatively the influence of the spin phonon coupling is similar in $`d=1`$ and $`d=2`$. In Fig. 4 the exact result from Refs. Eggert et al., 1994; Klümper, 1998 for the 1D and the Monte Carlo results for the 2D Heisenberg model are shown. Both in 1D and 2D we find a significant shift of the maximum and a strong reduction of the maximum height as compared to the Heisenberg model.

### III.2 Mean phonon occupation numbers

Further insight into the model can be gained by studying the influence of the spin phonon coupling on the mean phonon occupation numbers

$`n`$ $`={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{ij}{}}a_{ij}^{}a_{ij}`$ (4)
$`m`$ $`={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{ij}{}}b_{ij}^{}b_{ij}.`$ (5)
as compared to the free phonon case. These numbers can be viewed as a measure of the strength of lattice vibrations and therefore allow to analyze how the lattice is influenced by the spin degrees of freedom.

As expected, we find no difference in the mean occupation numbers $`n`$ and $`m`$, and therefore restrict the following discussion to the behaviour of $`n`$. Figure 5 shows the Monte Carlo results for $`n`$ for different values of $`g`$ and $`\omega `$ in a plot vs. $`T/\omega `$. The data are compared to the Bose distribution for free Einstein phonons

$$n_{\text{free}}(T)=\frac{1}{e^{\omega /T}1},$$
(6)
which is also shown in Fig. 5. Again we find a striking similarity to the 1D bond coupling model. In $`d=1`$ it has been found that the mean phonon occupation numbers obey the relationKühne and Löw (1999)

$$n(T)=n_0+n_{\text{free}}(T)$$
(7)
with a temperature-independent constant $`n_0`$. As a good approximation, this relation is valid in a temperature range $`0.5JT3J`$ in the 2D case as well. This can be seen most clearly in Fig. 6, which shows the same results as Fig. 5 with relation (6) subtracted from the Monte Carlo data.

In order to derive an expression for $`n_0`$, we averaged the differences $`nn_{\text{free}}`$ over the temperature range $`0.5JT3J`$ and plotted these values vs. $`g^2J^2/\omega ^2`$ (see Fig. 7). By applying linear regression we find that in $`d=2`$ the shift obeys the relation

$$n_0(1.375\pm 0.003)\left(\frac{gJ}{\omega }\right)^2.$$
(8)
Thus the only difference between the relations for the mean phonon occupation numbers in the 1D and 2D case is given by the numerical prefactor in (8), the value being 1.375 in $`d=2`$ and 2 in $`d=1`$.Kühne and Löw (1999)

We close this section with a technical remark concerning our choice for the cutoff for the phonon occupation numbers. As can be seen in Fig. 5, in the whole temperature range the measured mean occupation numbers are more than an order of magnitude smaller than the cutoff 40. In retrospect our choice therefore is justified. With (7) and (8) we have also found an expression that might be important for other numerical methods (e. g. exact diagonalization) which depend crucially on a (low) cutoff in the phonon numbers.

### III.3 Specific heat

Another important thermodynamic quantity is the specific heat per site $`C`$. Although in principle this observable is directly accessible in the experiment, it is dominated by lattice vibrations, making it difficult to extract its magnetic part. Even for a simple model as given by Hamiltonian (II), we find this behaviour confirmed. Figure 8 shows Monte Carlo data for the specific heat in a system with $`g=2`$ and $`\omega =2J`$ and the exact result

$$C_{\text{free}}(T)=2\left(\frac{\omega }{T}\right)^2\frac{e^{\omega /T}}{(e^{\omega /T}1)^2}$$
(9)
for free phonons of the same frequency (the factor two accounts for two phonons per lattice site). There is only a small difference in the overall height of $`C`$ and $`C_{\text{free}}`$. At high temperatures, both curves approach the same constant value, yielding the Dulong-Petit rule.

The spin phonon coupling influences the specific heat significantly though. As can be seen in Fig. 8 as well, the curve for the system with spin phonon coupling differs significantly from the sum of Monte Carlo results for the 2D Heisenberg model and the contribution (9) for free Einstein phonons. Note that both the strong fluctuations and the divergency of the data for $`T0`$ are due to difficulties in evaluating the specific heat within Monte Carlo procedures as discussed in Ref. Fye and Scalettar, 1987.

## IV Spin correlation function and ground state properties

We now turn our attention towards ground state properties of model (II). In principle the Monte Carlo method is only applicable at finite temperatures. By analyzing the behaviour of the spin correlation function at low temperatures, however, it is possible to make statements about system properties at $`T=0`$.

The argument is as follows. Suggest that we choose the coupling constants in (II) such that the system is Heisenberg-like, showing long range Néel order in the ground state. Then for reasons of universality we expect that the spin correlation function

$$G(\stackrel{}{d})=\frac{1}{N^2}\underset{\stackrel{}{r}}{}\stackrel{}{\sigma }_\stackrel{}{r}\stackrel{}{\sigma }_{\stackrel{}{r}+\stackrel{}{d}}$$
(10)
obeys the result known for the 2D Heisenberg modelChakravarty et al. (1989); Ding and Makivić (1990); Barnes (1991)

$$G(\stackrel{}{d})(1)^{d_1+d_2}|\stackrel{}{d}|^\lambda e^{|\stackrel{}{d}|/\xi (T)},$$
(11)
with the algebraic exponent $`\lambda `$ closeDing and Makivić (1990) to the classical Ornstein-Zernike value of $`\frac{1}{2}`$. Here $`\xi (T)`$ is the spin correlation length which can be interpreted as the mean size of domains with antiferromagnetic order. At $`T=0`$ these domains get macroscopic, because for $`T0`$ the correlation length diverges exponentially.Chakravarty et al. (1989); Ding and Makivić (1990); Barnes (1991); Beard et al. (1998) Assuming $`\lambda =\frac{1}{2}`$ in (11), this means that the static structure factor

$$S(\stackrel{}{q})=\underset{\stackrel{}{d}}{}e^{i\stackrel{}{q}\stackrel{}{d}}G(\stackrel{}{d})$$
(12)
diverges for momentum111In the following we set $`a1`$ for the lattice constant. $`\stackrel{}{q}=(\pi ,\pi )`$ with the linear system size like $`N^{\frac{3}{2}}`$. As long as the correlation length in the infinite system stays significantly larger than the system sizes considered, this behaviour should be visible at low temperatures.

We first illustrate this in case of the 2D Heisenberg model ($`g=\omega =0`$ in (II)) at $`T=0.1J`$, taking $`10^5`$ spin updates and choosing a Trotter number of $`M=160`$ for the evaluation of spin correlations. At this temperature, the correlation length in the infinite system is of the order of $`10^9`$ lattice spacings.Chakravarty et al. (1989); Ding and Makivić (1990); Barnes (1991); Beard et al. (1998) As can be seen in the upper plot of Fig. 9, the static structure factor shows a pronounced maximum at $`\stackrel{}{q}=(\pi ,\pi )`$, and the peak height roughly scales with the system size like $`N^{\frac{3}{2}}`$. Note that deviations from the $`N^{\frac{3}{2}}`$-dependence might indicate that in the quantum system $`\lambda `$ differs slightly from the value of $`\frac{1}{2}`$ chosen above.

We now return to model (II) with spin phonon coupling and discuss our results for two different choices of coupling constants. To compare our results for different values of $`g`$ and $`\omega `$ we drop the unphysical terms $`Jg/2(a_{ij}^{}+a_{ij}),Jg/2(b_{ij}^{}+b_{ij})`$ in (II). We therefore measure the temperature in units of the rescaled magnetic coupling $`J^{}`$ of the effective Hamiltonian (II). We calculated spin correlations at $`T=0.1J^{}`$, taking $`5\times 10^5`$ spin updates and $`M=160`$ in our calculations. First, we consider a system with a small value for $`\frac{gJ}{\omega }`$ ($`g=0.1`$, $`\omega =8J`$) where no dimerization is expected. Again the static structure factor shows a pronounced peak for $`\stackrel{}{q}=(\pi ,\pi )`$, and the peak height scales with the system size like $`N^{\frac{3}{2}}`$ (middle of Fig. 9), indicating Heisenberg-like behaviour as anticipated.

For the second system we choose $`g=2`$ and $`\omega =2J`$. In the 1D case this choice corresponds to a system which strongly dimerizes in the ground state.Raas et al. (2002) Even here we find a pronounced peak of $`S`$ for $`(\pi ,\pi )`$ (bottom of Fig. 9), the maximum height scaling like $`N^{\frac{3}{2}}`$. The interpretation is that even in the regime of large values for $`\frac{gJ}{\omega }`$ the system shows antiferromagnetic order in the ground state. This is a striking difference to the 1D bond coupling model. Compared to the case with $`g=0.1`$ and $`\omega =8J`$, however, we find that the spin phonon coupling counteracts the tendency of the system to order antiferromagnetically in the ground state. As can be seen in Fig. 9, for fixed system size the peak heights $`S(\pi ,\pi )`$ in the case of strong spin phonon coupling are slightly diminished as compared to the weak coupling regime.

Our results can be confirmed by a direct analysis of the temperature dependence of the spin correlation length. For both systems and at various temperatures we extracted finite system correlation lengths $`\xi _N`$ by fitting the function

$$f(\stackrel{}{d})=a(1)^{d_1+d_2}\left(\frac{e^{|\stackrel{}{d}|/\xi }}{\sqrt{|\stackrel{}{d}|}}+\frac{e^{(N|\stackrel{}{d}|)/\xi }}{\sqrt{N|\stackrel{}{d}|}}\right)$$
(13)
with two free parameters $`a,\xi `$ to our data for system sizes $`N=10,12,14,20,24`$. For $`g=0.1`$ and $`\omega =8J`$, we selected values $`M=120`$ for temperatures $`0.5JT0.9J`$ and $`M=80`$ for $`TJ`$, taking $`10^5`$ spin updates (for $`N=24`$ we chose $`M=120`$ for all temperatures). For $`g=2`$ and $`\omega =2J`$ the choice was $`M=160`$ for $`1.5JT1.9J`$, $`M=120`$ for $`2JT4J`$ and $`M=80`$ for $`T5J`$, again averaged over $`10^5`$ Monte Carlo sweeps (for $`N=20`$ we took $`1.5\times 10^5`$, for $`N=24`$ and $`1.5JT1.9J`$ we took $`2\times 10^5`$ spin updates). The values for the Trotter number are sufficiently large to avoid finite size effects in Trotter direction.

As has been said above, in case of the Heisenberg model in leading order $`\xi `$ behaves like $`e^{\frac{J}{T}}`$ at low temperatures.Chakravarty et al. (1989); Ding and Makivić (1990); Barnes (1991); Beard et al. (1998) In the upper panel of Fig. 10 the natural logarithm of $`\xi _N`$ is plotted vs. $`\frac{1}{T}`$ for $`g=0.1`$ and $`\omega =8J`$. At high temperatures, no dependence on the system size is visible, and as expected we find the same behaviour as in the Heisenberg-model. At low temperatures finite size effects become important and the curves branch off from the asymptotic linear behaviour. In the case of strong spin phonon coupling the graph shows very similar features (bottom of Fig. 10). We therefore find Heisenberg-like behaviour of $`\xi `$ at finite temperatures also in the regime of strong spin phonon coupling. The main difference between the two cases is that at the same effective temperature the correlation length for $`g=2`$, $`\omega =2J`$ is significantly smaller than for $`g=0.1`$, $`\omega =8J`$. The analysis of $`\xi `$ therefore also implies that a strong spin phonon coupling weakens antiferromagnetic order. Both observations are consistent with the conclusions drawn from the analysis of $`S(\stackrel{}{q})`$. Note that in both plots of Fig. 10 the temperatures are given in units of $`J^{}`$.

By means of scaling arguments our analysis can be extended to make direct statements about ground state properties. Suppose the system shows long range Néel order in the ground state. In terms of the renormalization group this means that there is a critical fixed point at $`T=0`$ which controls the system properties at low temperatures. In this case a finite size scaling ansatzCaracciolo et al. (1995, 1995); Beard et al. (1998)

$$\frac{\xi _{2N}(T)}{\xi _N(T)}=F\left(\frac{\xi _N(T)}{N}\right)$$
(14)
holds, where $`F`$ is a universal scaling function. With the data from Fig. 10 it is possible to test the scaling prediction (14). Plotting $`\frac{\xi _{2N}}{\xi _N}`$ vs. $`\frac{\xi _N}{N}`$ with $`N=10,12`$ we find that both for $`g=0.1`$, $`\omega =8J`$ (top of Fig. 11) and $`g=2`$, $`\omega =2J`$ (bottom of Fig. 11) the data lie on one curve. The shape of the scaling function $`F`$ in (14), however, depends on the choice of the coupling constants. The interpretation is that in the weak and in the strong coupling regime model (II) shows long range Néel order in the ground state, strongly confirming the conclusions drawn from our analysis of the static structure factor at low temperatures.

It is instructive to compare these results to the low temperature behaviour of spin correlations in $`d=1`$. Our argumentation is completely analogous to the 2D case. As has been discussed in Sect. I, the 1D model shows a quantum phase transition between a Néel ordered and a dimerized phase.Weiße et al. (1999); Kühne and Löw (1999); Raas et al. (2002) Though only the approximate shape of the phase separation line in coupling constant space has been determined, it is known that for small values $`\frac{gJ}{\omega }`$ the system is Heisenberg-like, showing quasi long range Néel order in the ground state. This means that the 1D correlation function

$$G^{\text{1D}}(d)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\stackrel{}{\sigma }_i\stackrel{}{\sigma }_{i+d}$$
(15)
decays exponentially at finite temperatures, with a rate given by the correlation length $`\xi ^{\text{1D}}\frac{1}{T}`$.Takada and Kubo (1986); Suzuki et al. (1990); Nomura and Yamada (1991); Suzuki et al. (1992); Klümper (1993) At $`T=0`$, there is a crossover to an algebraic decayLuther and Peschel (1975); Fogedby (1978); Korepin and Izergin (1984); Bogoliubov et al. (1986)

$$G^{\text{1D}}(d)\frac{(1)^d}{d},$$
(16)
and at $`T=0`$ the static structure factor

$$S^{\text{1D}}(q)=\underset{d=1}{\overset{N}{}}e^{iqd}G^{\text{1D}}(d)$$
(17)
diverges for momentum $`q=\pi `$ like the $`N`$th partial sum of the harmonic series with the system size. As in $`d=2`$, signs of this divergence should be visible at low temperatures.

For large values $`\frac{gJ}{\omega }`$, on the other hand, the chain dimerizes. This means that long range dimer order is established in the ground state. Therefore the spin correlation length $`\xi ^{\text{1D}}`$ stays finite at $`T=0`$, and we expect no dependence of $`S^{\text{1D}}`$ on $`N`$ in the low temperature regime.

In Fig. 12 the static structure factors $`S^{\text{1D}}`$ for two systems with the same choice of coupling constants as in $`d=2`$ are plotted. As in $`d=2`$, we measure the temperature in units of the rescaled magnetic coupling $`J^{}`$ of the 1D counterpart of (II). We selected $`T=0.1J^{}`$ and $`M=160`$, executing $`5\times 10^5`$ spin updates. In both systems $`S^{\text{1D}}`$ shows a maximum at $`q=\pi `$. In the Heisenberg-like system the maximum is more pronounced though, and for small system sizes the peak height depends on $`N`$. For larger system sizes, however, such a behaviour is not visible. This leads to the conclusion that in contrast to the 2D case the correlation length in the infinite system is not larger than the system sizes in consideration. In the system with dimerization in the ground state we find the expected behaviour: The values $`S^{\text{1D}}(\pi )`$ do not depend on the system size, indicating that $`\xi ^{\text{1D}}`$ is very small.

By analyzing the temperature dependence of $`\xi ^{\text{1D}}`$ it is possible to distinguish more clearly between the two regimes. Similar to $`d=2`$ we extracted correlation lengths $`\xi ^{\text{1D}}`$ by fitting an exponential decay with two free parameters to our data. For $`g=0.1`$ and $`\omega =8J`$, we executed $`5\times 10^5`$ spin updates and selected $`M=160`$, which is large enough to avoid effects by the finite Trotter number. Furthermore, the system sizes were chosen that large that finite size effects are negligible ($`N=500`$ for $`T=0.025J,0.05J`$, $`N=400`$ for $`T=0.075J,0.125J`$ and $`N=300`$ for $`T=0.1J,0.15J`$). This can be seen in the inset of Fig. 13, where $`\xi _N^{\text{1D}}`$ is plotted vs. $`N`$ at $`T=0.05J0.05J^{}`$. For $`g=2`$ and $`\omega =2J`$ also $`5\times 10^5`$ spin updates were made. The correlation lengths are that small that a chain length of $`N=200`$ is sufficient to make statements about the thermodynamic limit. The effect of the finite Trotter number is more important in this system. For $`M=400`$ at the lowest temperatures, however, the effect is smaller than the error which enters our analysis during the fitting procedure (we chose $`M=400`$ for $`T=0.05J,0.075J,0.1J`$, $`M=360`$ for $`T=0.15J`$ and $`M=160`$ for $`T0.2J`$).

Figure 13 shows the correlation lengths in a plot vs. the inverse temperature in units of the rescaled effective coupling $`J^{}`$. In both systems we find the expected behaviour: In the Heisenberg-like system the correlation length grows linearly with the inverse temperature, while in the dimerized system $`\xi ^{\text{1D}}`$ takes small values and shows no such dependence in the temperature range shown.

We close this section with a final remark concerning the 2D model (II). The results from this section need not mean that the model shows no lattice distortion. In both the statically dimerized stair and plaquette models e. g. a phase with coexisting dimerization and long range antiferromagnetic order is known.Sirker et al. (2002) For model (II) it is therefore conceivable that small lattice distortions appear which cannot be detected by analyzing spin correlations at low temperatures. Even a finite temperature phase transition seems possible, because due to the Mermin-Wagner theoremMermin and Wagner (1966) a breaking of the discrete lattice symmetry at finite temperatures cannot be excluded in a 2D system. However, we find no hints on a finite temperature phase transition in the behaviour of thermodynamic properties in the temperature range discussed in Sect. III. Further investigations of the model seem appropriate to clarify whether dimerization appears and which dimerization pattern is realized.

## V Summary

By combining loop updates for spin and cluster updates for phononic degrees of freedom we have developed a quantum Monte Carlo algorithm to study the properties of the 2D antiferromagnetic Heisenberg model coupled to bond phonons.

As thermodynamic quantities are concerned, we studied the susceptibility, mean phonon occupation numbers and specific heat in the temperature range $`0.5JT4J`$. The properties of the model at finite temperatures are similar to the 1D case.

For temperatures $`0.5JT3J`$, we derived an expression for the mean phonon occupation numbers which is of practical value for further investigations of the model.

We investigated the temperature dependence of the spin correlation length for two choices of coupling constants. Our analysis indicates that the model shows long range Néel order in the ground state both in the regime of weak and strong spin phonon coupling.

###### Acknowledgements.

The authors would like to thank M. Braden, A. Klümper, E. Müller-Hartmann, C. Raas, W. Weber and R. Werner for valuable discussions. We acknowledge funding by the DFG in the Sonderforschungsbereich 608.


text:
# Continuum in the spin excitation spectrum of a Haldane chain (erratum).

\[

\]

In Ref.  we reported inelastic neutron scattering study of the spin excitation spectrum in the S=1 quasi-one-dimensional antiferromagnet CsNiCl3. A cross-over from the coherent, single-particle magnon excitation at wavevectors close to $`q=\pi `$, to a continuum dominated by the two-particle states at small wavevectors, $`q0.5\pi `$, was clearly demonstrated by the measured intensities shown in the Fig. 2 of that paper. For the sake of comparison, the two-magnon lineshapes predicted from the nonlinear $`\sigma `$-model by Affleck and Weston in Ref.  were also shown, by the dotted lines, in the figure.

It recently came to our attention that theoretical curves shown in Fig. 2 of our paper  are slightly incorrect. This is clear if one compares the onset of the dotted line in Fig. 2 with the 2-magnon boundary shown in Fig. 3 of the paper, which is at lower energy. The error in Fig. 2 came about because an incorrect value of the Haldane gap, about twice the expected value of $`0.41J`$, was, by mistake, used for the single-particle dispersion in the formulas of Ref. . As a result, the discrepancy between the prediction of the nonlinear $`\sigma `$-model and the experimental results looks larger than it actually is.

Figure 1 here is the corrected version of the Fig. 2 of Ref. , and is now in agreement with the Fig. 3 of that paper. Clearly, there is also much better agreement between the data and the $`\sigma `$-model prediction now (note the absence of an arbitrary intensity scaling). Nevertheless, the observed continuum still lies noticeably below the lowest energy of the two non-interacting purely 1D magnons, supporting our previous conclusion that there is an attraction between the quasiparticles which form the continuum states . Such attraction is probably a result of the inter-chain coupling which leads to a 3D ordering in CsNiCl3. On the mean field level it can be understood in a manner similar to the spinon attraction in the coupled S=1/2 chains, if the excitations are the “spin-zero defects” of Ref. , as shown in Fig. 2 below.

We thank Ian Affleck for pointing out this error, and F. Essler for encouraging discussions. This work was carried out under Contract DE-AC02-98CH10886, Division of Materials Sciences, US DOE. The work on SPINS was supported by NSF through DMR-9986442.


text:
# The apeNEXT project

## I INTRODUCTION

APE is one of several projects (for a review see 1 ) in the theoretical physics community that have developed massively parallel, high-performance computer architectures. The driving force why physicists develop and build computers by themselves is the success of numerical simulations in understanding the interactions of elementary particles, in particular their strong interactions described by quantum chromodynamics (QCD). In the absence of closed-form analytical solutions for theories, like QCD, one of the most interesting approximation schemes is a reformulation of the theory on a discrete lattice (see 2  for a short introduction, or 3 ). The original theory is recovered as the lattice spacing $`a`$ goes to zero. This approach, pioneered by K. Wilson more than 25 years ago 4 , is the starting point for Lattice Gauge Theory (LGT). This discrete and computer-friendly formulation of quantum field theory has triggered an immense activity.

The phenomena investigated with such simulations range from the permanent confinement of quarks inside hadrons to the cosmological phase transition that occurred in the early phases of the universe or in matter under extreme conditions as produced in heavy-ion collision experiments. Within the framework of LGT, fundamental parameters of QCD, like the masses of quarks or the strength of the running strong coupling constants, have been computed from first principles. Also, theoretical concepts such as spontaneous chiral symmetry breaking and even the mathematical structure of the theory itself can be tested with modern simulation techniques. One of the big challenges is the determination of weak matrix elements of hadronic states to understand the interplay between weak and strong interactions. Problems like the $`\mathrm{\Delta }I=1/2`$ rule or the violation of CP symmetry are still open. The study of the heavy quark semileptonic decays is crucial for the determination of the Cabibbo-Kobayashi -Maskawa angles which are basic parameters of the Standard Model. A further example is the non-leptonic decay $`K\pi \pi `$, relevant to understand CP violation 7 .

Many of the current LGT projects focus on the simulation of QCD with dynamical fermions. Because of the limits in available computing power one is forced in many cases to apply the so-called quenched approximation, where the effects of vacuum fermion loops are neglected. Although the currently available computing resources allow to relieve this approximation, it will be extremely hard to lower the masses of the dynamical quarks towards their physical values. It will be even more difficult to reduce the lattice spacing and to do simulations closer to the continuum limit. A tremendous amount of computer power is required to overcome these limitations. A panel of the European Committee for Future Accelerator (ECFA), which proposed an ambitious research program for the coming years, estimates that European research groups would need well over 10 TFlops of compute power 8 .

In order to make these computing resources available at a reasonable price, various research groups have engaged in the development of supercomputers which are specifically optimized for their applications. In this paper we describe the Array Processor Experiment (APE) project, which was started in the mid eighties by the Istituto Nazionale di Fisica Nucleare (INFN) and is now carried out within the framework of a European collaboration with DESY and the University of Paris Sud.

The structure of this paper is as follows: in the next section we briefly cover the older members of the APE family. We then describe in some detail APEmille, the APE generation currently used in physics production simulations. Subsequently, we discuss the architecture of apeNEXT, the new generation of APE systems. This is the most important part of our paper, followed by a short discussion of the apeNEXT software environment. The paper ends with some concluding remarks.

## II THE FAMILY OF APE MACHINE

The evolution over more than one decade of APE systems is briefly recollected in Table 1.

The first generation of APE computers dates to the mid eighties. In Fig. 1, a picture of the original APE processor, made out of off-the-shelf electronic components is shown as a historical remark. APE100, the second generation of APE supercomputers, had been the leading workhorse of the European lattice community since the middle of the 1990s. Several parts of the APE100 machine are shown in Fig. 2.

Commissioning of APEmille, the third generation of APE systems, started in the year 2000. These machines make a further 2 TFlops of computing power available to the LGT community. A description of the APEmille architecture is given in a later section.

In order to keep up with future and growing requirements, the development of a new generation of a multi-TFlops computer for LGT, apeNEXT, is in progress. The main goal proposal  is the development and commissioning of a supercomputer with a peak performance of more than 5 TFlops and a sustained efficiency of $`O(50\%)`$ for key lattice gauge theory kernels. Aiming for both large scale simulations with dynamical fermions and quenched calculations on very large lattices the architecture should allow for large on-line data storage (of the order of 1 TByte) as well as input/output channels which sustain at least $`O(0.5)`$ MByte per second per GFlops. Finally, the programming environment should allow smooth migration from older APE systems, i.e. support the TAO language, and introduce for the first time a C language compiler.

## III APEmille SYSTEMS

APEmille is a massively parallel computer optimized for simulating QCD. The architecture is single instruction multiple data (SIMD) and all nodes run strictly synchronously at a moderate clock frequency of 66 MHz. The communication network has a three-dimensional topology and offers a bandwidth of 66 MBytes/s/node. The smallest APEmille unit (see Fig. 3) is a processing board with $`2\times 2\times 2`$ nodes. The largest stand-alone systems built until now consist of $`4\times 8\times 8`$ nodes (see Fig. 4).

Three different integrated circuits (ASICs) have been custom-developed for APEmille. Program execution is controlled by a control processor, which also performs the subset of integer arithmetics common to the whole SIMD partition. Computations using local integers and all floating-point operations are done in parallel by all computing nodes. At each clock cycle, the arithmetic processors are able to complete the “normal” operation $`a\times b+c`$, where $`a`$, $`b`$ and $`c`$ are single precision (32 bit) complex operands. This gives a peak performance of 528 MFlops per node. Each node has 32 MBytes of local memory and a very large register file, holding up to 512 data words. Remote communications between the nodes are implemented as direct memory access which is controlled and routed by the communication processors. The processors are controlled by very long instruction words (VLIW). This allows efficient scheduling of the microcode at compile time. Much effort has therefore been put into the development of software tools for generating efficient code (see later for more details). Loading of the executables and all other operating system services are handled via PCs running Linux. One host PC per four boards is directly attached to the APEmille backplane. It uses a PCI bus to communicate with the processing boards. The user interfaces with the system by logging onto an front-end PC, which is also running Linux and allows to spawn the program on APEmille and to monitor its execution. Altogether, a large APEmille installation may be controlled by more than 20 PC’s. These are connected via an Ethernet network.

APEmille machines are installed at several sites all over Europe, as detailed in table 2. They provide a very stable and reliable computing environment, with typical up-times of the order of 85 %.

## IV apeNEXT PROCESSOR AND GLOBAL DESIGN

apeNEXT has been designed with the main goal of having an architecture as close as possible to its previous generation, as perceived at the user level, while improving on performance as much as possible through the use of more advanced technology. As a consequence, there are many similarities with previous generation systems, while a number of new design challenges had to be solved.

For apeNEXT all processor functionalities, including the network devices, are integrated into one single custom chip running at a clock frequency of 200 MHz. Unlike former machines, the nodes will run asynchronously, which means that apeNEXT follows the single program multiple data (SPMD) programming model.

Our performance goal is based on a peak performance of 1.6 GFlops per processor in 64-bit double precision, while the communication bandwidth between neighboring nodes is 200 MByte/s. We envisage large apeNEXT systems with 2000 processing nodes, delivering a peak performance of 3.2 TFlops. The key design parameters are listed in table 3.

Two new key features have been introduced in apeNEXT. First, apeNEXT is a SPMD (as opposed to SIMD) system. Each processing node is a fully independent processor, with a full-fledged flow-control unit and, of course, a number-crunching unit. The node has access to its own memory bank, where both program and data are stored. It executes its own copy of the program at its own pace. Nodes need to be synchronized only when a data-exchange operation is performed. This architecture may be labeled as a distributed-memory parallel computer, in which nodes exchange data through some sort of “message-passing” scheme. Internode communications are started by the program on the sending node that initiates a data communication step. This operation is matched by a corresponding instruction on the destination node, that explicitly receives the data packet. The latency associated to a “message” is extremely short, of the order of 2 to 3 times the latency associated to an access to local memory. For this reason, the actual data rate between nodes is bandwidth-limited (as opposed to latency-limited) even for short packets, so sequences of short accesses can be freely programmed without significant performance losses. This is an important feature in LGT, where the natural size of data packet transferred to remote processors is not larger than about 200 bytes.

A second important architectural enhancement lays on the possibility of routing all read memory accesses (to local or remote nodes) through a receiving queue, which can be later accessed by the processor with zero latency. This feature is mainly used to perform data pre-fetch in critical kernel loops, taking into account that the address patterns are regular and easily predictable. The basic idea here is that all data items needed to perform iteration $`(i+k)`$ ($`k=1,2,\mathrm{}`$) of the loop are pre-fetched during iteration $`i`$ and stored into the queue. When iteration $`(i+k)`$ starts, data will be immediately available to the processor, effectively hiding almost all latency effects. Note that some of the memory accesses will be local, and some will be remote. They are started in sequence, but they may complete in a different order (a remote access may take longer than a local one). However, the queue mechanism automatically ensures that data are delivered to the processor in the same order in which they were requested from (remote or local) memory.

The complete processing element is contained in just one custom-designed integrated circuit, called J&T, connected to a memory bank of 256–1024 MBytes with Double Data Rate (DDR) Dynamic RAM chips. In turn, the processor and its memory is housed on a small piggy-back printed circuit board shown in Fig. 6. This assembly is basically a complete processor, that delivers more than 1.5 GFlops of processing power in double precision with a power consumption of approximately 7 Watt, that is about 10 times less than current generation high-end PCs. The compactness and the low power consumption of this basic building block are key ingredients to build the very compact multi-node system described later on.

The block diagram of the J&T chip is shown in Fig. 7. It is a 64-bit architecture, optimized for floating point performance.

We would like to highlight some selected details of the processor shown in Fig. 7:

* A large register file of 256 registers each containing a pair of 64 bit words. All operands for the arithmetic unit arrive from the register file and all results are written back here.
* An arithmetic box which performs floating point as well as integer operations. The basic operation executed by the floating point unit (FPU) is the “normal” floating-point operation $`a\times b+c`$. At each clock cycle one normal operation can be started to provide a maximal throughput of eight floating-point operations per clock cycle. The operands are either complex values or pairs of real values. All floating-point data is represented in the 64-bit double-precision format of the IEEE standard. The integer unit (INT) operates on pairs of 64-bit integers. The arithmetic box also contains the logical unit (LU). Finally, a “special function” block (LUT) provides initial approximations for iterative or series-expansion evaluation of inverses, square-roots, logarithmic and exponential functions.
* An address-generation unit (AGU) which computes addresses for memory access independently and concurrently with the main arithmetic box. This is an important feature to boost sustained performance.
* A memory controller (MC) supporting a memory bank of 256–1024 MBytes based on standard DDR-SDRAM. The memory is used to store both data and program instructions. A consequence of this organization is that conflicts between data and instruction load-operations are present. Two strategies have been employed to avoid these conflicts. First, the hardware supports compression of the microcode. The compression rate depends on the level of optimization, typically values are in the range of 40–70%. Instruction de-compression is performed on-the-fly by dedicated hardware. Second, an instruction buffer allows pre-fetching of (compressed) instructions. Under complete software control, a section of the instruction buffer can be used to store performance critical kernels for repeated execution.
* A flow-control unit that executes programs specified as a sequence of compiler prepared microcode words, using the VLIW control style.
* A network interface which contains seven LVDS link interfaces. Each link is bi-directional allowing send and receive operations to run concurrently. Once a communication request is queued it is executed independently of the rest of the processor, which is a prerequisite for overlapping network and floating point operations. Each transmitter (TX) is able to send one byte per clock cycle, i.e. the gross bandwidth is 200 MByte per second per link. Due to protocol overhead the effective network bandwidth is $`180`$ MByte per second. The network latency is $`O`$(0.1 $`\mu `$s) and therefore at least one order of magnitude smaller than for today’s commercial high performance network technologies. While six links are used to connect the processing node to its neighbors, the seventh link is used for input/output operations.
* A set of fifos that implements the queue mechanism described earlier in the text:

+ A first set of fifos (the TX ones) hold the data words received from the main memory (or the register file) until the network is able to send them to another node. The main role of this system is to decouple the (fast) memory system from the (slower) network links.
    + A second set of fifos receives data items from the network and stores them until the processor wants to load them into the register file.
    + Finally, a third set of fifos is used to keep control information. These fifos are needed to guarantee the correct order of the data bursts.
* A slow serial interface based on the I2C standard, used for system initialization, debugging and exception handling.

The overall structure of an apeNEXT system is similar to APEmille. We have a three-dimensional array of processors (see Fig. 5). Periodic boundary conditions are applied. Although all nodes are connected to their nearest neighbors only, the hardware allows routing across up to three orthogonal links to all nodes on a cube, i.e. connecting nodes at distance $`(\mathrm{\Delta }_x,\mathrm{\Delta }_y,\mathrm{\Delta }_z)`$ with $`|\mathrm{\Delta }_i|1`$.

Clusters of 16 apeNEXT processors will be assembled onto just one printed-circuit board. The processors are arranged in the configuration of a three-dimensional structure of $`4\times 2\times 2`$ processors (see fig. 8). A set of $`16`$ boards is housed within one system crate. All communication links between these nodes are established through the crate backplane. Larger systems are assembled connecting together several crates using external cables.

The system is completed by a number of host PCs that can be tailored to match user-specific input/output requirements (see fig. 5). The main input/output channel is based on one of the seven data links available on each processor. One such link from each processing board can be connected to a host interface board that follows the PCI specifications and can be plugged into a standard PC (see Fig. 9). The actual number of input/output channels can be tailored to match user-specific input/output requirements. From the apeNEXT point of view, an input/output operation is simply a remote communication with a special remote node. This structure has the main advantage of requiring a minimum of programming and operating system overhead.

## V SOFTWARE AND BENCHMARKS

We will provide both a TAO and a C compiler for apeNEXT. The latter is based on the freely available lcc compiler lcc  and supports most of the ANSI 89 standard with a few language extensions required for a parallel machine. Both compilers generate a high-level assembly. A assembler pre-processor (mpp) is used to translated this into a low-level assembly. For machine specific optimizations at this assembly level, e.g. address arithmetics and register move operations, the software package sofan is under development. Finally, the microcode generator (shaker) optimizes instruction scheduling, which for APE machines is completely done in software. We finally plan to develop a linker which allows to combine several microcode files into one executable. An overview on the compilation chain is shown in Fig. 10.

For most parts of the compiler software stable prototype versions are available and were already used to benchmark the apeNEXT design. For this purpose we considered various typical linear algebra operations, like the product of two complex vectors. To execute this operation iteratively needs reading two complex numbers per iteration, which on apeNEXT takes at least two clock cycles. The arithmetic instructions executed per iteration, however, have a throughput per clock cycle of up to one. Therefore, this operation is on apeNEXT restricted by the memory bandwidth and the theoretical sustained performance therefore limited to 50%. From VHDL simulations that include all machine details the efficiency was found to be 41%. Even higher performance rates can be achieved for operations where the data after being loaded to the register file is re-used several times. This is, e.g., the case when multiplying arrays of SU(3) matrices. For this operation we measured an efficiency of 65%.

In QCD simulations most of the time is spent applying the Dirac operator, e.g. the Wilson-Dirac operator $`M=1\kappa H`$. Therefore, we investigated the multiplication with the so-called hopping term, i.e. $`H\psi `$. This operation involves remote communications and therefore depends on the number of processors involved. The maximum number of processors is limited by the size of the problem, i.e. the lattice volume. In case of dynamical fermion simulations a multiplication with the hopping term has to be done much more often than updating the gauge fields. We therefore kept a local copy of the gauge fields to save network bandwidth. Considering the worst case where the problem is distributed over the maximum number of processors, we found the sustained performance to be 56%. This figure is made possible by extensive use of the pre-fetch features of the processor, which allowed to completely overlap floating point operations and network communications, such that the time when the processor waits for data becomes almost zero. This eventually indicates an excellent scaling behavior of the apeNEXT architecture.

## VI OUTLOOK AND CONCLUSIONS

The hardware design of the next generation of APE custom built computers has been completed. Prototype boards and the communication backplane are already available. A prototype apeNEXT processor is expected out of the foundry in late August 2003. A larger prototype installation is planned to be running by end of 2003. There exists a stable prototype version for all parts of the compiler software. Based on this software we expect that key lattice gauge theory operations will be able to run at a sustained performance of $`O(50\%)`$ or more. We hope that apeNEXT will be a key element of LGT simulations in the next few years.

###### Acknowledgements.

We would like to thank W. Errico for his important work in the early phases of the project, and A. Agarwal, T. Giorgino and M. Lukyanov for their contributions.