text:
# Uncovering CDM halo substructure with tidal streams

## 1 Introduction

A generic prediction of Cold Dark Matter (CDM) cosmology is the existence of a significant amount of substructure in gravitationally collapsed structures such as galaxy clusters, galaxy halos and even dwarf galaxies (Klypin et al., 1999; Moore et al., 1999b, 2001). On the scale of large spiral galaxies like the Milky Way, some 500 dense clumps are expected to orbit in the Halo. These structures have very dense cores, with a very steep radial profile of universal shape (increasing as $`r^1`$ or $`r^{1.5}`$ in the central regions, Navarro, Frenk & White 1997; Moore et al. 1999a), which renders them essentially impervious to Galactic tides.

The current “dark matter crisis” stems from the difference between these predictions and observations of low-mass galaxies. The first problem is that rotation curves of gas-rich low-mass disk galaxies are not as peaked as required by CDM. Van den Bosch et al. (2000) suggested that this disagreement could be due to an error in the analysis of HI observations, in which the beam-smearing effect of extant HI radial velocity curves had not been taken into account. However, this appears not to be the case, as recent (high spatial resolution) $`H_\alpha `$ rotation curves have shown (Côté et al., 2000; de Blok et al., 2001; Marchesini et al., 2001).

A second fundamental problem also exists: the mass of the CDM substructures in a galaxy-sized halo is similar to that inferred for galaxy satellites, yet the number of CDM substructures in a galaxy-sized halo is much larger than the number of observed dwarf galaxies(Klypin et al., 1999; Moore et al., 1999b). It is, of course, possible that these galaxies have a quite different non-cosmological origin (perhaps due to tidal interactions in the early universe). In either case the problem remains: where is the population of 500 CDM substructures orbiting the Milky Way and other large galaxies? Recently, several studies have investigated the effect of reioninzation on the early evolution of small structures (Bullock et al., 2000; Somerville, 2001; Tully et al., 2001). This re-ionizing radiation sets in at a redshift $`6`$$`10`$, and originates from either the first stars or quasars. Indeed, the first direct evidence for reionization (at $`z6`$) has recently been found (Becker et al., 2001). The intergalactic and galactic mediums are ionized by this radiation, and any gas that is not in deep potential wells is lost from the protogalaxies. The limiting mass for maintaining a gas component appears to correspond to a circular velocity of $`30\mathrm{km}\mathrm{s}^1`$, well above the circular velocity of all the dwarf spheroidal galaxies known. (However, it may be possible for low-mass systems to maintain a gas-fraction if the gas managed to cool to dense molecular form before the epoch of reionization). This “squelching” of star-formation in low-mass galaxies neatly solves the CDM satellite over-production problem.

A potential candidate for the missing dwarf systems that are supposed to inhabit the Galactic halo is the large population of high velocity clouds that pepper the sky (see Wakker 2001 for a recent review). With little or no stellar content, these are seen to possess velocities incompatible with Galactic rotation models, although Blitz et al. (1999) have shown that the kinematic characteristics of the velocities of the clouds are compatible with them being members of an accreting population distributed within a megaparsec of the Milky Way. At these distances, the kinematic signatures of the clouds imply a mass-to-light ratio of 10-50 Braun & Burton (2000). Currently, their distances are extremely difficult to determine, with Complex A being the only system to which the distance, of between 4 and 10kpc, has been measured (van Woerden et al., 1999). With the promise of accurate distances from $`H_\alpha `$ (Bland-Hawthorn & Maloney, 1999), a full appraisal of the nature of high velocity clouds and their relation to the formation of the Milky Way can be made; but a preliminary analysis suggests that the clouds are scattered in the Galactic halo (Bland-Hawthorn & Maloney, 2001), and are not of cosmological origin.

Nevertheless, CDM theory predicts that a large population of completely or almost completely dark “galaxies” should be present in the Milky Way. Here we propose an observational test of this scenario. In the early galaxy, many globular clusters are expected to have surrounded the Milky Way. Galactic tidal forces have been destroying these objects (Gnedin & Ostriker, 1997), and of the initial population of a few hundred objects, only $`100`$ remain, and many of these show signs of tidal disruption (Leon, Meylan & Combes, 2000; Odenkirchen et al., 2001).

In a static potential, as a low-mass stellar system, such as a globular cluster, loses progressively more and more mass to tidal forces, it develops long tidal tails that closely follow the orbit of the globular cluster. (In the limit of a zero mass system, the tidal stream exactly follows the progenitor orbit). Here we explore how the presence of the population of dark matter clumps predicted by CDM alters the phase-space structure of a globular cluster tidal stream. For the special case when the global Galactic potential is nearly spherical, this corresponds to a broadening of the stream from a thin great-circle stream into a wide band on the sky. In a companion paper (Ibata et al. 2002; Paper 2), we use these results to motivate the search for tidal streams in the 2 Micron All Sky Survey (2MASS) dataset.

Other ways to probe the dark matter substructure have been proposed. Recently, Metcalf (2001) and Chiba (2001) suggested that substructure in external galaxies may be detected through their gravitational lensing of background quasars, as a lumpy halo will produce image and brightness configurations different to a smooth matter halo. While such an approach may reveal any missing sub-halos, the method is subject to the vagaries of gravitational lens modeling. During the final preparation of this work, Mayer et al. (2001) presented a study investigating the signatures that are imprinted by different dark matter models on tidal streams of dwarf spheroidal galaxies. They find that the structure and kinematics of the outer regions of dwarf galaxies can also allow one to differentiate between dark matter models.

## 2 Numerical simulations

Our aim is to quantify the difference in the structure of the tidal tail of an ancient globular cluster (or even an ancient globular cluster remnant) if the Galactic halo has a smooth mass distribution compared to the case when the Halo also contains significant substructure. We have modelled the smooth components of the Galaxy as fixed potentials, using the Galactic mass model of Dehnen & Binney (1998), which contains a disk, thick disk, interstellar medium, bulge, spheroid and halo components. The parameters of this model are detailed in Ibata et al. (2001). The halo was taken to have a mass normalization that gives a total model circular velocity of $`v_c=210\mathrm{km}\mathrm{s}^1`$ at $`50\mathrm{kpc}`$. The mass flattening of the potential was left as a free parameter $`q_m`$.

We adapted the fast parallel code PKDGRAV to include the forces due to these mass distributions by including a multipole expansion code kindly supplied by W. Dehnen. In all the following simulations we maintained accurate forces by using an opening angle of $`0.75`$ and expanding the cell moments to hexadecapole order. Two body relaxation was suppressed by using a spline softening of 10 pc, such that the inter-particle forces are completely Newtonian at 20 pc. A variable time-step scheme is used based on the local acceleration, $`\mathrm{\Delta }t<\eta \sqrt{|\mathrm{\Phi }|}/a`$, and density, $`\mathrm{\Delta }t<\eta /\sqrt{G\rho }`$, with the accuracy parameter $`\eta =0.03`$ (Quinn et al., 1997). Typically we used 100000 base steps to integrate 12 Gyr. With the variable time-steps, this was equivalent to taking $`5\times 10^7`$ time-steps.

The initial globular cluster was modelled with a King model (King, 1966), populated with $`10^4`$ particles. The total mass of the globular cluster was taken to be $`10^6\mathrm{M}_{}`$, the model had concentration parameter $`c=1.0`$ and central velocity dispersion $`\sigma =4.1\mathrm{km}\mathrm{s}^1`$, yielding a tidal radius of $`r_t=300\mathrm{pc}`$. This model is of approximately the same size and mass as Omega Centauri, but of lower concentration (by a factor of $`1.7`$). Note that we are not aiming at presenting an accurate model of the disruption of a globular cluster — we do not have sufficient resolution for that — instead, here we seek to study the kinematic behaviour of the unbound tidal streams that emerge from disrupting low-mass systems.

The initial position and systemic velocity of the King model were randomly chosen from the halo component of a galaxy model (Boily, Kroupa & Penarrubia-Garrido, 2001) that gives a similar rotation curve to the adopted Dehnen & Binney mass model.

Integrating the globular cluster models for $`10\mathrm{Gyr}`$ in the smooth potential gives rise to very narrow tidal tails. An example of one of these models is shown in Figure 1, which displays the position, heliocentric radial velocity and distance distribution of the stream particles. At its narrowest, the width of this stream is $`100\mathrm{pc}`$, similar in width to the tidal radius of the initial King model, while at its widest, the stream is $`3\mathrm{kpc}`$ wide. The shape of the halo mass distribution for the simulation displayed in Figure 1 is spherical, $`q_m=1`$, which is motivated by recent analyses of the outer halo of the Milky Way (Ibata et al., 2001; Chiba & Beers, 2001).

For each one of these models, we also investigate the effect of the presence of a large number of moving substructures in the halo. However, we first need a realistic way of modeling these substructures. According to Navarro, Frenk & White (1997) (hereafter NFW), and other authors, galaxy halos have a “universal” density profile $`\rho (r)1/(r(1+r)^2)`$, which fits all scales currently probed by the numerical cosmology simulations. The potential corresponding to this density profile is very simple: $`\varphi (r)1\mathrm{ln}(1+r)/r`$, and the forces due to this potential would be straight-forward to add in to our N-body integrator. Unfortunately, however, this density profile is not well-behaved either at small or at large radii (the total mass diverges); to fix these problems, the force near the center of the sub-halo would need to be softened, and the density would have to be made to fall off to zero at some large radius (the tidal radius of the substructure would be a good choice). A simpler solution is to approximate the NFW profiles with a softened point-mass potential. Figure 2 compares the radial acceleration due to an NFW potential (dotted line), a Kepler potential (dashed line), and a softened Kepler potential (solid line) with spline softening of $`3.41`$ times the scale radius $`r_s`$ of the NFW profile. Note that, at all radii, this softened point-mass potential underestimates the forces compared to the NFW potential.

The substructures are added into the Galactic halo as softened point masses, with a distribution of sub-halo masses that follows a relation $`\mathrm{log}_{10}(N_c)=3(14v_c/V_{Global})`$ (our parameterization of Figure 2 of Moore et al. 2001). At the low-mass end of this distribution, there are 435 sub-halos of circular velocity $`v_c>0.03V_{Global}`$ ($`M10^7\mathrm{M}_{}`$). The halo mass fraction in this lumpy component is relatively small, as it accounts for less than $`10`$% of the mass of the smooth halo. The initial positions and velocities of the particles representing the sub-halos are also chosen from the halo component of the Boily, Kroupa & Penarrubia-Garrido (2001) Galaxy model, but with a pericenter cutoff of $`10\mathrm{kpc}`$, since dynamical friction will likely make the orbits of dark satellites with smaller pericenter radii decay quickly.

The effect of these numerous dark matter clumps is marked: in Figure 3 we show the position, velocity and distance distribution of the tidal stream particles. In this simulation, only 1% of the mass of the halo is not in the smooth component, yet the tidal stream has been substantially heated and is substantially wider than in the smooth halo (the stream is so fluffy that it is difficult to measure a width).

A better way to look at the simulations is in the space of the integrals of the motion. In the axisymmetric potentials considered here, the integrals are the total energy per unit mass $`E`$ and the $`z`$ component of angular momentum per unit mass $`L_z`$. It is also useful to inspect the total angular momentum per unit mass, $`L`$, which is an approximate integral of the motion in the relatively spherical potentials we consider here. Figure 4 shows the relation between these three quantities. The top row of Figure 4 is derived from the simulation of the globular cluster model in the smooth spherical halo whose structure on the sky was previously displayed in Figure 1. After $`12\mathrm{Gyr}`$ the particles have remained localised in the space of the integrals of the motion. The bottom panel shows the same globular cluster model on the same orbit, but simulated in the halo model that contains the additional NFW sub-halos (corresponding to the simulation shown in Figure 3). In this lumpier potential, the globular cluster stream is comparatively much more dispersed, especially in the $`L_z`$ parameter. The r.m.s. dispersion in $`L_z`$ is a factor of 5 larger in the lumpy halo simulation.

Inspecting the spatial distribution of particles can be useful for differentiating between a smooth and lumpy halo when the mass distribution is spherical (or close to spherical). However, in a flattened halo, differential precession mimics the effect of heating by CDM clumps, and it becomes impossible to differentiate between the two causes of spatial dispersion of the stream. It is only when viewing the space of the integrals of the motion that the effect of precession decouples from the stochastic heating by the CDM clumps. We have repeated the same experiments as above in a flattened halo, with $`q_m=0.7`$, a flattening value predicted by numerical cosmology (Katz, 1991; Dubinski & Carlberg, 1991; Warren et al., 1992; Katz & Gunn, 1991; Summers, 1993; Dubinski, 1994). The sub-halos are chosen from a halo distribution function which is also flattened to $`q_m=0.7`$. Figure 5 shows the relation between the integrals of the motion $`E`$ and $`L_z`$ and the approximate integral $`L`$. Again, the dispersing effect of the sub-halos is very strong, with the r.m.s. dispersion in $`L_z`$ increasing by a factor of 9 with respect to the smooth halo case.

### 2.1 Suite of simulations

To obtain a statistical understanding of the change in $`L_z`$, we undertook a series of simulations of the globular cluster model with initial position and velocity randomly-drawn from the spherical halo distribution function model mentioned above. Ten simulations were undertaken in the Galaxy model with the smooth spherical halo. A further ten simulations probed the same models with the addition of the NFW sub-halos, for a total of 20 simulations with $`q_m=1.0`$. To probe the situation with a more flattened halo, we undertook a further 20 simulations with $`q_m=0.7`$.

The result of these experiments is shown in Figure 6, where we display the r.m.s. dispersion in $`L_z`$ for simulations in a smooth halo (filled circles) and in in the presence of halo substructures (star symbols). The top panel shows the case in a spherical halo ($`q_m=1.0`$), while the bottom panel is for the flattened halo ($`q_m=0.7`$) simulations.

Both in a spherical and in a flattened halo, the presence of CDM-like substructures increases substantially, on average, the RMS dispersion in $`L_z`$ of the stream stars. On some orbits (such as that of our simulation #9, the effect can be small, however. Our simulations suggest that globular cluster streams with a very dispersed specific angular momentum, say $`\sigma L_z>200\mathrm{kpc}\mathrm{km}\mathrm{s}^1\mathrm{M}_{}^{}{}_{}{}^{1}`$ require that the halo possesses significant substructure. Smaller values, $`\sigma L_z<50\mathrm{kpc}\mathrm{km}\mathrm{s}^1\mathrm{M}_{}^{}{}_{}{}^{1}`$ say, would require a smooth halo.

## 3 Observational considerations

We next consider whether it will be possible to identify a $`10^6\mathrm{M}_{}`$ stream orbiting the Halo at distances up to $`100\mathrm{kpc}`$ ($`mM=20`$) with photometry, radial velocity, parallax and proper motion data from the future astrometric mission GAIA (Perryman et al., 2001). The great advantage of the GAIA dataset is that gives access to the full 6-dimensional phase space information for many stars, and at least 4-dimensional phase-space information (for those stars that are too far to have measured parallaxes and radial velocities). In addition, the photometric information (in 15 bands) can give age and/or metallicity estimates, providing a powerful discriminant between stellar populations.

Developing a full model of the Galaxy as viewed by GAIA is not a trivial task and it is beyond the scope of this paper. Here we simply make an order of magnitude estimate of the contrast of plausible globular cluster streams over the background by estimating the relative phase space density of stream stars over the large-scale Galactic spheroid (or stellar halo) component. We estimate the phase-space size of the stream as follows. The $`10\mathrm{km}\mathrm{s}^1`$ radial velocity uncertainty of GAIA exceeds the intrinsic radial velocity dispersion of the stream, and is about 10% of the halo radial velocity dispersion. The proper motion dispersion at $`100\mathrm{kpc}`$ corresponds to about $`10\mathrm{km}\mathrm{s}^1`$ or 10% of the halo dispersion (two coordinates). The positional width of the streams are $`\stackrel{<}{}3\mathrm{kpc}`$, or about 1% of the width of the halo. With the full 15 photometric bands, the distance uncertainty of photometric parallax measurements is likely to be $`<10`$%. Thus the phase-space volume occupied by the stream is likely to be about one millionth of that of the halo. Since the stellar halo has a mass of only $`10^9\mathrm{M}_{}`$, the contamination will amount to approximately one star in a thousand. With a limiting magnitude of $`V21`$, the GAIA measurements will easily reach down to the horizontal branch at $`100\mathrm{kpc}`$. Placing all of the stream at a distance of $`d=100\mathrm{kpc}`$, we expect $`1600`$ stars for $`V>21`$ per $`10^6\mathrm{M}_{}`$ of progenitor mass (assuming the luminosity function of the Draco by Odenkirchen et al. 2001, corrected to a mass to light ratio of $`M/L=3`$). For radial velocity measurements the limiting magnitude is $`V17`$, and the number of available stream stars from the $`10^6\mathrm{M}_{}`$ population drops to $`180`$. The detection of these streams will be straightforward with GAIA, and the increase in the dispersion of $`L_z`$ due to any halo substructure will be easily resolved.

With results expected towards the end of the second decade of this century, the GAIA mission is not for the impatient. So it is worthwhile to investigate whether any presently available datasets already contain sufficient information to uncover the predicted Galactic halo dark substructures. To this end we undertook an analysis of the 2MASS dataset, which provides positions and infrared CCD photometry for sources brighter than $`K=14.3`$. This gives access to a relatively restricted sample of Galactic halo M-giants, with no information on their kinematics, and only a very poor constraint on their distances. The search for the halo streams in the 2MASS dataset is presented in Paper 2; the analysis reveals a single very strong stream, that due to the ongoing disruption of the Sagittarius dwarf galaxy.

Is it possible to use the Sagittarius stream to constrain the lumpiness of the halo? To investigate this issue we ran two simulations, in the same manner as above, one in a smooth halo, and the other with NFW clumps, using as an initial model for the Sagittarius dwarf the model D2 of Ibata et al. (2001). Figure 7 shows the result of this experiment. This model is much more massive ($`M=5\times 10^8\mathrm{M}_{}`$) than the globular cluster model simulated above, and is initially much more extended (half-mass radius of $`0.9\mathrm{kpc}`$), so it is naturally widely dispersed in $`E`$, $`L_z`$ and $`L`$. Though the stream particles are heated and dispersed by the sub-halos, the effect is subtle, and will be very difficult to disentangle from the normal dynamical evolution of the dwarf galaxy. So low mass streams from globular clusters need to be found to act as probes of the dark matter.

## 4 Conclusions

We have shown that stellar streams of low mass systems such as disrupted globular clusters can be used to constrain the substructure of the dark halo. If current cosmological simulations correctly represent reality, ancient stellar streams from globular clusters should be substantially dispersed over the sky, even if on average, the Galactic halo potential is close to spherical. This effect is readily noticeable, as the streams will have a large dispersion in the z-component of specific angular momentum $`L_z`$. For the particular halo studied in this contribution $`\sigma L_z>200\mathrm{kpc}\mathrm{km}\mathrm{s}^1\mathrm{M}_{}^{}{}_{}{}^{1}`$ implies the presence of substantial dark matter substructure.

We have simplified the analysis in two major ways. Highly softened point-masses were used to model the halo substructures, which leads to an underestimate of the increased dispersion in $`L_z`$ in the presence of the NFW clumps. Our other major assumption, that the Galactic potential is axisymmetric and does not evolve with time, will likely have the opposite effect, namely to underestimate the dispersion in the absence of the NFW clumps. Our study should be refined once the Galactic potential is better constrained (with data from the GAIA and SIM missions).

This work has shown that an all-sky survey with high precision astrometry, such as will be obtained with the GAIA satellite, we will easily be able to detect globular cluster streams, and measure their kinematic properties, especially the dispersion in $`L_z`$.

The heating due to the numerous encounters with the dark substructure must also affect the computed destruction rate of globular clusters and low-mass dwarf galaxies alike. It is possible that the lower limit to the mass of dSph galaxies $`10^7\mathrm{M}_{}`$ may be set by these disruptive encounters.

Discovering and determining the thickness of ancient stellar streams will allow us to probe the invisible dark halo structures predicted by Cold Dark Matter cosmology. Finding a single narrow stream of width comparable to the tidal radius of a globular cluster would stand in strong contradiction with this theory. Data from the forthcoming GAIA mission planned by the European Space Agency will provide the opportunity to undertake this study.


format: text

text:
# The Weakly Pushed Nature of “Pulled” Fronts with a Cutoff

## I Introduction

Pulled fronts are fronts that propagate into an unstable state, for which the propagation dynamics is essentially that they are being pulled along by the growth and spreading of small perturbations about the unstable state, into which the front propagates. Concretely, this means that their asymptotic speed $`v_{as}`$ is equal to the linear spreading speed $`v^{}`$ of perturbations around the unstable state, $`v_{as}=v^{}`$ . Fronts that propagate into an unstable state but for which $`v_{as}>v^{}`$ are often termed “pushed”. The name stems from the intuitive idea  that in this regime, the dynamics in the nonlinear front region or the bulk region behind the front actually drives the front propagation: effectively it pushes the front from behind, and the front moves with a speed that is higher than the natural speed with which small perturbations about the unstable state spread by themselves ahead of the front.

It is clear from the definition that the concept of a pulled front essentially pertains to a continuum formulation of the relevant dynamical variables. The linear spreading speed $`v^{}`$ is defined and calculated in practice by considering perturbations of arbitrarily small amplitude about the unstable state of the dynamical equations; the value of $`v^{}`$ then follows from an asymptotic analysis of the linearized dynamical equations . However, in all cases, in which one cannot ignore the fact that matter is made of discrete particles, one cannot perturb the unstable state by any arbitrary small amount, because this amount must be at least one “quantum” of particle large.

To model this discrete nature of the constituent particles by means of a continuum equation, Brunet and Derrida  studied the nonlinear diffusion equation

$`{\displaystyle \frac{\varphi }{t}}={\displaystyle \frac{^2\varphi }{x^2}}+f(\varphi ),`$ (1)
with a cutoff $`ϵ`$ in the growth term $`f(\varphi )`$,

$`f(\varphi )=\mathrm{\Theta }(\varphi ϵ)[\varphi \varphi ^n],`$ (3)
$`n>1,\text{e.g., }n=\mathrm{\hspace{0.17em}2}\text{or}\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}3}.`$
Without the cutoff, i.e., for $`ϵ=0`$, this equation is the well-known nonlinear diffusion equation, which has been used since long as the simplest model to study front propagation into an unstable state . Brunet and Derrida  found that the asymptotic front speed $`v_{as}`$ goes as

$`v_{as}=v_ϵv^{}{\displaystyle \frac{\pi ^2}{\mathrm{ln}^2ϵ}},`$ (4)
where $`v^{}=2`$ is the asymptotic speed of the corresponding pulled front of Eq. (1) for $`ϵ=0`$. The above formula shows that the front speed $`v_ϵ`$ converges very slowly to the asymptotic speed $`v^{}`$; this illustrates that unlike pushed fronts, pulled fronts are very sensitive to small changes in the dynamics of the phase into which they propagate.

In comparing with stochastic models of particles on a lattice, Brunet and Derrida associated the cutoff $`ϵ`$ with $`1/N`$, where $`N`$ is the average number of particles in a correlation region in the saturation phase behind the front . Although the validity of this identification has been the matter of some debate, it appears that Eq. (4) with $`ϵ=1/N`$ does give the proper asymptotic correction to the front speed even for very large $`N`$. We refer to the literature  for a further discussion of the applicability of these ideas to stochastic models.

It is intuitively clear that as soon as we introduce this cutoff, fronts that are pulled for $`ϵ=0`$, must actually become weakly pushed as soon as $`ϵ>0`$. After all, any perturbation around the value $`\varphi =0`$ does not start to grow until the local $`\varphi `$ value crosses $`ϵ`$, so strictly speaking, the linear spreading velocity $`v^{}(ϵ)`$ of arbitrarily small linear perturbations about the state $`\varphi =0`$ vanishes. As $`v_ϵ>v^{}(ϵ)=0`$, one clearly must have a weakly pushed front. With this idea in mind, it is natural to address the convergence of the front speed to the asymptotic value, since it is well known that the speed of pulled fronts relaxes algebraically slowly to the asymptotic value $`v^{}`$ , while pushed fronts normally have exponential relaxation to their asymptotic speed.

These observations motivate us to investigate here the slowest relaxation modes of the stability spectrum of fronts for the nonlinear diffusion equation (1), with a cutoff $`ϵ`$ in the growth term (3). We calculate these relaxation modes explicitly for small $`ϵ`$, and find that the slowest relaxation times $`\tau _m`$ are given by

$$\tau _m^1\frac{\left[(m+1)^21\right]\pi ^2}{\mathrm{ln}^2ϵ}m=\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2}\mathrm{}$$
(5)
Hence, the relaxation times of the front velocity and profile approach zero as $`ϵ0`$, but only logarithmically slowly. Just like the corrections to the front speed for practical values of $`ϵ`$ are often significant, so is the exponential relaxation, for example, for $`ϵ=10^5`$, the longest relaxation time $`\tau _1`$ is about 4.48. Thus, while in the absence of a cutoff the front speed is approached very slowly, only as $`3/2t`$ where $`t`$ is the time , with a realistic value of $`ϵ`$, the front speed converges relativly quickly to the asymptotic value.

## II Stability analysis of the asymptotic front solution

### A The Stability Operator

The asymptotic shape of the front is a uniformly translating front solution $`\varphi _ϵ(x,t)`$ which is a function of only the comoving coordinate $`\xi =xv_ϵt`$, and which is obtained by solving the ordinary differential equation

$`v_ϵ{\displaystyle \frac{d\varphi _ϵ(\xi )}{d\xi }}={\displaystyle \frac{d^2\varphi _ϵ(\xi )}{d\xi ^2}}+f(\varphi _ϵ(\xi )).`$ (6)
In carrying out the linear stability analysis of this front solution it is convenient to follow the standard route of transforming the linear eigenvalue problem into a Schrödinger eigenvalue problem . We consider a function $`\varphi (x,t)`$, which is infinitesimally different from $`\varphi _ϵ(\xi )\varphi _ϵ(xv_ϵt)`$ in the comoving frame, i.e., $`\varphi (x,t)=\varphi _ϵ(xv_ϵt)+\eta (\xi ,t)`$. Upon linearizing the dynamical equation in the comoving frame, one finds that the function $`\eta (x,t)\eta (\xi ,t)`$ obeys the following equation:

$`{\displaystyle \frac{\eta }{t}}=v_ϵ{\displaystyle \frac{\eta }{\xi }}+{\displaystyle \frac{^2\eta }{\xi ^2}}+{\displaystyle \frac{\delta f(\varphi )}{\delta \varphi }}|_{\varphi =\varphi _ϵ}\eta .`$ (7)
Since this equation is linear in $`\eta `$, the question of stability can be answered by studying the spectrum of the temporal eigenvalues. To this end, we express $`\eta (\xi ,t)`$ as

$`\eta (\xi ,t)=e^{Et}e^{v_ϵ\xi /2}\psi _E(\xi ),`$ (8)
which converts Eq. (7) to the following one-dimensional Schrödinger equation for a particle in a potential (with $`\mathrm{}^2/2m=1`$):

$`\left[{\displaystyle \frac{d^2}{d\xi ^2}}+{\displaystyle \frac{v_ϵ^2}{4}}{\displaystyle \frac{\delta f(\varphi )}{\delta \varphi }}|_{\varphi =\varphi _ϵ}\right]\psi _E(\xi )=E\psi _E(\xi ).`$ (9)
In Eq. (9), the quantity $`V(\xi )=\left[{\displaystyle \frac{v_ϵ^2}{4}}{\displaystyle \frac{\delta f(\varphi )}{\delta \varphi }}|_{\varphi =\varphi _ϵ}\right]`$ plays the role of the potential. If we now denote by $`\xi _0`$ the coordinate of the point where $`\varphi _ϵ(\xi )=ϵ`$, then for the nonlinearity (3) the potential $`V(\xi )`$ is easily seen to have the form

$`V(\xi )=\left[{\displaystyle \frac{v_ϵ^2}{4}}1+n\varphi _ϵ^{n\mathrm{\hspace{0.17em}1}}(\xi )\right]\mathrm{\Theta }(\xi _0\xi ){\displaystyle \frac{1}{v_ϵ}}\delta (\xi \xi _0)`$ (11)
$`+{\displaystyle \frac{v_ϵ^2}{4}}\mathrm{\Theta }(\xi \xi _0).`$
The $`\delta `$-function in Eq. (11) appears from the functional derivative in Eq. (7), since there is a discontinuity of magnitude $`ϵ`$ in $`f(\varphi )`$ at $`\varphi =ϵ`$. This discontinuity contributes an amount equal to

$$\varphi _ϵ\frac{d\mathrm{\Theta }(\varphi _ϵϵ)}{d\varphi _ϵ}=\varphi _ϵ\delta (\varphi _ϵϵ)=\frac{ϵ}{|\varphi _ϵ^{}(\xi _0)|}\delta (\xi \xi _0)$$
(12)
to $`V(\xi )`$. If we combine this with the fact that $`|\varphi _ϵ^{}(\xi _0)|=ϵv_ϵ`$, which follows immediately from the fact that one simply has $`\varphi _ϵ(\xi )=ϵe^{v_ϵ(\xi \xi _0)}`$ for $`\xi \xi _0`$, one obtains the $`\delta `$-function term in the potential given in Eq. (11).

The form of the potential $`V(\xi )`$ is sketched in Fig. 1. Notice that $`\varphi _ϵ(\xi )`$ is a monotonically increasing function from $`ϵ`$ at $`\xi _0`$ towards the left, asymptotically reaching the value $`1`$ as $`\xi \mathrm{}`$. As a result, for $`\xi <\xi _0`$, $`V(\xi )`$ also increases monotonically towards the left, from $`v_ϵ^2/41+nϵ^{n1}\pi ^2/\mathrm{ln}^2ϵ`$ at $`\xi =\xi _0`$, to $`(n\pi ^2/\mathrm{ln}^2ϵ)n`$ as $`\xi \mathrm{}`$. At $`\xi _0`$, there is an attractive $`\delta `$-function potential of strength $`1/v_ϵ1/2`$ and a finite step of height 1. The crucial feature for the stability analysis below is the fact that $`V(\xi )`$ stays remarkably flat at a value $`\pi ^2/\mathrm{ln}^2ϵ`$ over a distance $`(\xi _0\xi _1)|\mathrm{ln}ϵ|`$, and then on the left of $`\xi _1`$, it increases to the value $`n`$, over a distance of order unity. As argued in Sec. II.B, this is a consequence of the nature of the solution $`\varphi _ϵ(\xi )`$.

If there are negative eigenvalues of the above Schrödinger equation, then according to Eq. (8) $`\eta (\xi ,t)`$ grows in time in the comoving frame, i.e., the front solution $`\varphi _ϵ(\xi )`$ is unstable. On the other hand, if there are no negative eigenvalues, then the asymptotic front shape is stable, and the spectrum of the eigenvalues then determines the nature of the relaxation of $`\varphi (x,t)`$ to the solution $`\varphi _ϵ(\xi )`$.

The full spectrum in general depends on the boundary conditions imposed on the eigenfunctions $`\psi _E`$. Here we consider only localized perturbations, for which we need to have $`\eta (\xi ,t)0`$ as $`\xi \pm \mathrm{}`$. Because of the exponential factor in Eq. (8), any eigenfunction $`\psi _E`$ which vanishes as $`\xi \mathrm{}`$ is consistent with vanishing $`\eta `$ towards the right. However, for $`\xi \mathrm{}`$, the eigenfunctions $`\psi _E`$ need to vanish exponentially fast with a sufficiently large exponent, so that when it is combined with the exponentially diverging term $`e^{v_ϵ/2}`$, they are still consistent with the requirement that $`\eta `$ vanishes for $`\xi \mathrm{}`$. For the lowest “energy” eigenvalues, which we will investigate below, we will demonstrate that these requirements are obeyed.

### B Shape of $`\varphi _ϵ(\xi )`$ and the Zero Mode of the Stability Operator

From the form in the potential, it is clear that the lowest “energy” eigenmodes, i.e., the slowest relaxation eigenmodes, are the ones that are confined to the bottom of the potential. This is the region where the nonlinear terms proportional to $`\varphi ^{n1}`$ are negligible, and which is often called the “leading edge” of the front profile. For $`ϵ1`$, the solution of $`\varphi _ϵ(\xi )`$ in this leading edge is given by

$`\varphi _ϵ(\xi ){\displaystyle \frac{|\mathrm{ln}ϵ|}{\pi }}\mathrm{sin}[z_i\xi ]e^{z_r\xi }\text{for}\xi _1\xi \xi _0|\mathrm{ln}ϵ|`$ (14)
$`=ϵe^{v_ϵ(\xi \xi _0)}\text{for}\xi \xi _0.`$
Here, $`z_i\pi /|\mathrm{ln}ϵ|`$ and $`z_r=1+O(ϵ^2)`$. The values of $`\varphi _ϵ(\xi )`$ and $`{\displaystyle \frac{d\varphi _ϵ}{d\xi }}`$ are continuous at $`\xi =\xi _0`$, and $`\varphi _ϵ(\xi _0)=ϵ`$. Although Eqs. (14) and (15) suggest at first sight that the $`\varphi _ϵ(\xi )`$ has a node at $`\xi =0`$, Eq. (14) is only valid in the leading edge, and $`\varphi _ϵ(\xi )`$ crosses over to other behavior around $`\xi _1`$, which makes the front solution $`\varphi _ϵ(\xi )`$ a mononically decreasing function of $`\xi `$. The value of $`\xi _1`$ is set by the criterion that around $`\xi _1`$, the nonlinear terms of $`f(\varphi _ϵ(\xi ))`$ starts to become significant, just like $`\xi _1`$ marks the point where the potential $`V(\xi )`$ crosses over from the asymptotic value on the left to the bottom value. The coordinate $`\xi _1`$, therefore, is more or less fixed; on the other hand, $`\xi _0`$ asymptotically diverges as $`|\mathrm{ln}ϵ|`$ for small $`ϵ`$, making $`(\xi _0\xi _1)`$ also diverge as $`|\mathrm{ln}ϵ|`$. This is an immediate consequence of the overall exponential decay of $`\varphi _ϵ(\xi )`$ in $`\xi `$ at the leading edge.

From the form in the potential, it is clear that the lowest “energy” eigenmodes, i.e., the slowest relaxation eigenmodes, are the ones that are confined to the bottom of the potential. We notice that among these modes, invariably there is a zero mode of the stability operator that is associated with the uniformly translating front solution of a dynamical equation like Eq. (1): since $`\varphi _ϵ(\xi )`$ and $`\varphi _ϵ(\xi +a)`$ are solutions of Eq. (6) for any arbitrary $`a`$, we find by expanding to first order in $`a`$ that $`\psi _0(\xi )=e^{v_ϵ\xi /2}{\displaystyle \frac{d\varphi _ϵ}{d\xi }}`$ is a solution of Eq. (9) with eigenvalue $`E=0`$. From the result (14) for the asymptotic front solution, we then immediately get to dominant order

$$\psi _0\mathrm{sin}z_i\xi ,z_i\pi /|\mathrm{ln}ϵ|,\xi _1\xi \xi _0.$$
(15)
Furthermore, since $`\varphi _ϵ(\xi )`$ is a monotonically decreasing function of $`\xi `$, the solution $`\psi _0(\xi )=e^{v_ϵ\xi /2}{\displaystyle \frac{d\varphi _ϵ}{d\xi }}`$ is nodeless. Since we know from elementary quantum mechanics that the nodeless eigenfunction has the lowest eigenvalue, this implies that all the other eigenvalues of Eq. (9) are positive, i.e., the solution $`\varphi _ϵ(\xi )`$ is stable.

The spectrum of eigenvalues of Eq. (9) for $`E>0`$, therefore, is going to determine the decay property of localized perturbations $`\eta (\xi ,t)`$ in time. We notice that for $`E>v_ϵ^2/41`$, the value of the potential on the far right, the spectrum of eigenvalues will be continuous. However, we are particularly interested in the smallest eigenvalues $`E_m>0`$ for small $`m`$, since these are the eigenmodes that decay the slowest in time. These are the eigenvalues associated with bound states in the potential well.

### C Lowest Eigenmodes and Eigenvalues for $`ϵ1`$

As $`ϵ0`$, the bottom well of the potential becomes very wide: its width diverges as $`|\mathrm{ln}ϵ|`$. As we know from elementary quantum mechanics, the lowest “energy” eigenfunctions then become essentially sine or cosine waves in the potential well with small wave numbers $`k`$ and correspondingly small “energy” eigenvalues.

Based on the fact that the potential $`V(\xi )`$ on the left rises over length scales of order unity, we now make an approximation. In the limit that the bottom well is very wide and the $`k`$ values of the bound state eigenmodes very small, it becomes an increasingly good and an asymptotically correct approximation to view the left wall of the well simply as a steep step, as sketched in Fig. 2 — we thus approximate the potential by

$`V_0(\xi )=n\left[\mathrm{\hspace{0.17em}1}\mathrm{\Theta }(\xi )\right]{\displaystyle \frac{\pi ^2}{\mathrm{ln}^2ϵ}}\mathrm{\Theta }(\xi )\left[\mathrm{\hspace{0.17em}1}\mathrm{\Theta }(\xi \xi _0)\right]`$ (17)
$`{\displaystyle \frac{1}{2}}\delta (\xi \xi _0)+\mathrm{\Theta }(\xi \xi _0).`$
On the right hand side, there is an attractive delta-function potential at the point where the potential shows a step to a value close to 1. It is easy to check that the prefactor of the delta-function of 1/2 is not strong enough to give rise to bound states with $`E<0`$, and as a result, for very small values of $`ϵ`$, the low-lying eigenmodes approach sine waves with nodes at the position of the walls of the potential ,

$$\psi _m\mathrm{sin}\left[k_m(\xi \xi _1)\right].$$
(18)
The condition that these solutions have nodes at the right edge of the well then yields

$$k_m\frac{(m+1)\pi }{\xi _0\xi _1}\frac{(m+1)\pi }{|\mathrm{ln}ϵ|},$$
(19)
implying that the corresponding eigenvalues are given by

$`E_m{\displaystyle \frac{\left[(m+1)^2\mathrm{\hspace{0.17em}1}\right]\pi ^2}{\mathrm{ln}^2ϵ}}m=\mathrm{\hspace{0.17em}0},1,2,\mathrm{}`$ (20)
Here, the first term between square brackets comes from the “kinetic energy” term $`k^2`$, while the second term orgininates from the value of the potential at the bottom.

Note that for $`m=0`$, the eigenmode $`\mathrm{sin}k_0`$ with eigenvalue $`E_0`$ is indeed the same as the zero eigenmode of Eq. (15) with $`k_0=z_i`$, which we calculated from the shape of the front solution $`\varphi _ϵ`$ in the leading edge. Besides verifying the consistency of our approach, this also confirms that there are no corrections to Eq. (20) for $`m=0`$: for $`m=0`$ it will yield an eigenvalue zero to all orders in $`ϵ`$. Therefore, the smallest nonzero eigenvalue, which governs the relaxation of the front velocity and profile to the asymptotic ones is $`E_1`$ with relaxation time $`\tau _1`$ given by

$`\tau _1^1=E_1{\displaystyle \frac{3\pi ^2}{\mathrm{ln}^2ϵ}}.`$ (21)
Equation (20) also confirms that as $`ϵ0`$, the gap between the spectral lines decreases as $`\mathrm{ln}^2ϵ`$, which is consistent with the fact that for a pulled front $`ϵ=0`$ and the spectrum becomes gapless. Also notice that for the eigenvalues in Eq. (20), the corresponding eigenmodes $`\psi _E(\xi )`$ decay as $`\mathrm{exp}[\sqrt{n}|\xi |]`$ for $`\xi \mathrm{}`$ and as $`\mathrm{exp}[v_ϵ\xi /2]`$ for $`\xi \mathrm{}`$, which make $`e^{v_ϵ\xi /2}\psi _E(\xi )`$ to go to zero for $`\xi \pm \mathrm{}`$, satisfying the boundary conditions discussed previously at the end of section II.A.


text:
# On the dominance of non-exotic meson-meson scattering by 𝑠-channel 𝑞𝑞̄ confinement states and the classification of the scalar mesons

## Introduction

Lattice QCD in principle offers the most direct way to link to experiment what we believe to be the fundamental theory of strong interactions . However, in view of the constant evolution of results from sophisticated lattice QCD solutions , it becomes ever more puzzling why mesonic resonances can be described as simple quark-antiquark systems in effective theories . But apparently, such a scenario works!

Seemingly, the perturbative vacuum states of QCD at low energies are not quarks and gluons, but rather confined constituent quarks and residual interactions, a picture that has been successful for several decades by now. What do lattice calculations teach us about constituent quarks? One might think of colour-triplet configurations of quarks, antiquarks, and glue, maybe with admixtures of higher colour multiplets, which mutually feel colour interactions. Can lattice QCD identify such substructures? Moreover, what happens exactly when those substructures suddenly turn colourless and cease to be confined? Does one observe colourless substructures that drift apart on the lattice?

In unitarised meson models one assumes an effective mass for the constituent quark, a confinement force for the remaining colour interactions, and a mechanism for decay . However, not knowing how the separation into constituent quarks, confinement, and decay can be derived from QCD, each model is the result of educated guesses, rather than of rigorous derivations starting from QCD. This frustrating state of the affairs has in the past three decades led to a proliferation of effective models and theories. However, no model exists so far which completely describes the resonances of meson-meson scattering to a satisfactory degree of accuracy. Nevertheless, some educated guesses are less successful than others. For example, $`q\overline{q}`$ models for resonances that do not take meson loops into consideration  can never be in agreement with experiment, since the resulting spectrum consists of zero-width bound states, whereas large widths are measured.

## Non-exotic meson-meson scattering

The assumption that non-exotic meson-meson scattering is dominated by the $`s`$-channel states of the confinement mechanism has been worked out in a series of papers . Here, we will confine our attention to a toy model that we studied in Ref. . There we obtain for the elastic low-energy partial-wave meson-meson scattering matrix,

$$S_{\mathrm{}}(p)=\mathrm{exp}\left(2i\delta _{\mathrm{}}(p)\right),$$
(1)
the relation

$$\text{cotg}\left(\delta _{\mathrm{}}(p)\right)=\frac{n_{\mathrm{}}(pa)}{j_{\mathrm{}}(pa)}\left[2\lambda ^2\mu paj_{\mathrm{}}^2(pa)\underset{n=0}{\overset{\mathrm{}}{}}\frac{\left|_n(a)\right|^2}{EE_n}\right]^1,$$
(2)
where $`p`$ represents the relative momentum in the CM frame of the two mesons, $`\mathrm{}`$ their relative angular momentum, and $`\mu `$ their reduced mass; $`E_n`$ ($`n=0`$, 1, 2, $`\mathrm{}`$) represent the energy eigenvalues of the constituent $`q\overline{q}`$ system to which the meson pair couples, and $`_n`$ the corresponding $`q\overline{q}`$ eigenfunctions; $`n_{\mathrm{}}`$ and $`j_{\mathrm{}}`$ stand for the spherical Bessel and Neuman functions, respectively.

The intensity of the coupling between the meson-meson system and the $`q\overline{q}`$ system is described by the parameter $`\lambda `$, whereas $`a`$ stands for the average distance at which the transitions from one system to the other take place .

For $`\lambda =0`$, we find $`S_{\mathrm{}}(p)=1`$, which describes a system of two non-interacting mesons. For small values of $`\lambda `$, the scattering matrix (1) has poles in the lower half of the complex-energy plane, which can approximately be given by

$$E_{\text{pole}}E_n\left|_n(a)\right|^2\left[\underset{n^{}n}{}\frac{\left|_n(a)\right|^2}{E_nE_n^{}}\frac{i}{2\lambda ^2\mu paj_{\mathrm{}}(pa)h_{\mathrm{}}^{(1)}(pa)}\right]^1,$$
(3)
indicating that to each value of the radial quantum number $`n`$ corresponds one singularity, i.e., one meson-meson scattering resonance. For higher values of $`\lambda `$, one can determine the locations of the poles in the scattering matrix by numerical methods.

The real parts of the singularities roughly correspond to the central resonance positions, $`E_r`$, whereas the moduli of the imaginary parts approximately equal half the resonance widths, $`\mathrm{\Gamma }_r`$. In short,

$$E_{\text{pole}}E_ri\frac{\mathrm{\Gamma }_r}{2}.$$
(4)
Singularities may be located on the real energy axis below threshold, representing stable (with respect to strong decay) mesons (e.g. $`K`$, $`J/\mathrm{\Psi }`$, $`\mathrm{{\rm Y}}`$ ).

For practical purposes, one might truncate the sum in formula (2) and substitute the truncated part by a constant .

For $`\lambda \mathrm{}`$ one finds

$$\text{cotg}\left(\delta _{\mathrm{}}(p)\right)=\frac{n_{\mathrm{}}(pa)}{j_{\mathrm{}}(pa)},$$
(5)
which represents scattering from a hard sphere of radius $`a`$. In this case, the interior of the $`q\overline{q}`$ state becomes unobservable and no resonance spectrum can be deduced from meson-meson scattering.

## Comparison with experiment

In Ref.  we compare the predictions of formula (1) for $`K\pi `$ $`S`$\- and $`P`$-wave scattering in $`I=1/2`$ with the experimental cross sections. For $`P`$ waves, in the region of the $`K^{}`$(892) resonance, we find that the lowest-lying state of the confinement spectrum is at some 945 MeV, whereas the corresponding pole comes out at $`(88727i)`$ MeV.

In a more refined model , which also takes inelasticity into account by considering the coupling to all channels with allowed initial and final states of pseudoscalar and vector mesons, and furthermore employs a more sophisticated mechanism for the coupling of $`q\overline{q}`$ confinement states to meson-meson scattering channels, it is found that the ground state of the confinement spectrum in this case comes out at 1.19 GeV, some 300 MeV above the position of the $`K^{}`$(892) pole. It shows that bare states can be several hundreds of MeVs away from the actual central resonance positions, and, moreover, that such conclusions are model dependent. The toy model of Ref.  yields a shift of only some 60 MeV. In this perspective, the question as to where a bound-state model should find its bare states is hard to be answered.

The latter question becomes even more difficult in the case of $`S`$-wave scattering. There we find, in the toy model of Ref. , that the ground state of the confinement spectrum is at 1.31 GeV and the corresponding pole at $`(1.460.12i)`$ GeV. However, further inspection of the singularity structure of the scattering matrix reveals another pole far below this energy region, namely at $`(714228i)`$ MeV. The latter singularity has no direct relation to any of the bare states. Hence, the $`I=1/2`$, $`J^P=0^+`$ ground state of bound-state models should be at some 1.3 GeV and not below 1 GeV. This result is confirmed by the full model  ( pole at $`(727263i)`$ MeV), in which also the poles belonging to the two isoscalars $`f_0(980)`$ and the rather controversial $`f_0(470208i)`$ ($`\sigma `$ meson), as well as to the isovector $`a_0(980)`$, have no direct relation with the ground states of the corresponding confinement spectrum at about 1.2 GeV.

## A model study of poles

The relation between the $`q\overline{q}`$ bare states and the poles in the corresponding scattering matrix can be found in the coupled-channel model (1) by considering the process of stepwise reducing the coupling constant $`\lambda `$. Such a study has been performed in detail in the toy model of Ref. , with the following result. All poles move towards a corresponding $`q\overline{q}`$ bare state on the real axis, as predicted by formula (3), except for the $`S`$-wave singularity below 1 GeV. The negative imaginary part of this pole grows inversely proportionally to $`\lambda ^2`$, implying that the corresponding “resonance” disappears into the background of $`K\pi `$ scattering. Unfortunately, such processes cannot be tested in experiment, since Nature corresponds to a fixed value for $`\lambda `$.

## Constituent-quark-pair creation

For $`P`$\- and higher-wave meson-meson scattering, we do not find singularities other than those which can be related to a $`q\overline{q}`$ state of the confinement spectrum. The extra poles below 1 GeV, described through pole doubling in Ref. , exclusively appear in $`S`$-wave meson-meson scattering. Hence, we must conclude that the latter “resonances” are a consequence of the mechanism of constituent-quark-pair annihilation and creation, which couples meson-meson initial and final states to the $`q\overline{q}`$ confinement states.

For $`P`$ and higher waves, the centrifugal barrier prevents the formation of such resonances in meson-meson scattering. But in the absence of a centrifugal barrier for $`S`$ waves, resonances are formed that in the cases of the $`f_0(980)`$ and the $`a_0(980)`$ are narrow enough to be clearly observable, but which for the $`f_0(470208i)`$ and $`K_0^{}(727263i)`$ are too broad to be firmly established.

One should note here that, when a pole with a large imaginary part lies close to threshold — close meaning that the distance from threshold to the real part of the singularity is smaller than or of the same order as the imaginary part — then the corresponding cross section has a shape which is very different from a standard Breit-Wigner. In the Argand plot, one finds a resonance motion that rapidly slows down for higher energies. If then, moreover, new thresholds get open and other rapid Breit-Wigner resonances show up, its appearance can hardly be recognised as that of a resonance, within the experimental accuracy.

Nevertheless, whether or not one associates a resonance with the controversial $`J^P=0^+`$ singularities is of no importance. What is crucial in the above observations is the fact that it settles the classification of the $`f_0(980)`$ and $`a_0(980)`$ resonances in a nonet scheme for mesons rather more naturally than other proposals.

## Resonance shapes

The model result that some of the light scalar resonances are broad, while others are narrow, has its origin in the effects of inelasticity. In Ref. , Table 1 of the Appendix, a list of inelasticity channels for the three scalar-meson isomultiplets is presented, as well as the intensities of the relative couplings.

We learn from this table that the isotriplet couples twice as strongly to $`\eta _n\pi `$ as to $`K\overline{K}`$, other thresholds lying at higher or much higher energies, which makes their effect hardly relevant here. However, only a small part of the $`\eta _n\pi `$ channel decays into $`\eta \pi `$, the rest into $`\eta ^{}\pi `$. This implies that, of the lowest channels, the $`K\overline{K}`$ and also the $`\eta ^{}\pi `$ channel are far stronger than the $`\eta \pi `$ channel (see also Ref. ). Elastic $`S`$-wave $`\eta \pi `$ scattering in the absence of inelasticity can be described by the toy model of formula (2). In Ref. , formula (2) has been applied to elastic $`S`$-wave $`K\pi `$ scattering. Now, when we substitute there the $`K`$ mass by the $`\eta `$ mass, then we obtain a toy model for elastic $`S`$-wave $`\eta \pi `$ scattering. With this substitution, we find for the model of formula (2) indeed a pole close to the $`\eta \pi `$ threshold, and with a relatively large imaginary part ($`763199i`$ MeV). Moreover, the related toy-model prediction for the $`\eta \pi `$ elastic $`S`$-wave scattering cross section does not show a clear resonance, exactly as in the case of the $`K_0^{}(727263i)`$ pole in $`K\pi `$ elastic scattering.

Inelasticity, which has been taken into account in Ref. , has two consequences here: first, the pole moves close to the $`K\overline{K}`$ threshold, with a smaller imaginary part ($`96828i`$ MeV), and, second, since the $`\eta \pi `$ threshold is far enough below that pole, its resonance shape turns more Breit-Wigner-like (see e.g. Fig. 2 of Ref. ). Nevertheless, upon reducing the coupling constant, the modulus of the imaginary part of this pole increases in a similar way as does the lower pole, the $`K_0^{}(727263i)`$, in $`K\pi `$ elastic $`S`$-wave scattering. Consequently, the two poles have a similar origin, not directly related to the bare spectrum. It is only through the strong interference of the $`K\overline{K}`$ and the $`\eta ^{}\pi `$ channels that a reasonable Breit-Wigner-like shape appears for the light isotriplet resonance $`a_0(980)`$.

We also observe from the above-referred table of Ref.  that, of the isoscalar complex ($`n\overline{n}`$ coupled to $`s\overline{s}`$), the $`n\overline{n}`$ couples strongly to $`\pi \pi `$, whereas the $`s\overline{s}`$ couples strongly to $`K\overline{K}`$, with, again, other thresholds lying higher or much higher, thus making their influence of little importance here. Furthermore, $`n\overline{n}`$ and $`s\overline{s}`$ are coupled to one another through the $`K\overline{K}`$ channel, which implies that the $`s\overline{s}`$ component of the isoscalar complex also couples to $`\pi \pi `$, but quite weakly. Hence, for the $`s\overline{s}`$ resonance $`f_0(980)`$ we can now repeat the arguments we gave for the isotriplet resonance shape, with $`\eta \pi `$ replaced by $`\pi \pi `$. Moreover, since the $`\pi \pi `$ threshold in the isoscalar case lies much lower than the $`\eta \pi `$ threshold in the isotriplet case, we find a more convincing Breit-Wigner-like shape for the $`f_0(980)`$ in $`\pi \pi `$ scattering  than for the $`a_0(980)`$ in $`\eta \pi `$ scattering.

However, the $`n\overline{n}`$ component of the isoscalar complex, which yields a pole close to the $`\pi \pi `$ threshold, has no further strong-inelasticity channel to allow for a Breit-Wigner-like shape for the corresponding resonance $`f_0(470208i)`$. The same happens to the isodoublet, which, according to the table of coupling constants mentioned before, couples strongly to the $`K\pi `$ channel. No further lower-lying inelasticity channel exists in this case. Consequently, also the $`K_0^{}(727263i)`$ has no Breit-Wigner-like shape.

None of the poles of this nonet of scalar resonances has a direct relation to the bare spectrum. By stepwise reducing the model coupling constant, all nine poles stepwise disappear into the complex plane with increasing negative imaginary part, whereas the corresponding structures in the meson-meson scattering cross sections stepwise disappear into the background.

## The $`𝑲`$-matrix

As one can easily observe from formula (2), we have no poles in the $`K`$-matrix at the energy eigenvalues $`E_n`$ of the confinement spectrum, since the hard-sphere-scattering part in the expression for the cotangent of the phase shift does not vanish at energies $`E_n`$. This contradicts the observation of Sarantsev and collaborators (Ref.  and references therein) that bare states are the singularities of the $`K`$-matrix, so this issue deserves further study.

In the limit of an infinitely strong coupling between the confinement and scattering sectors, formula (2) predicts that no bare spectrum can be observed in meson-meson scattering other than the hard-sphere spectrum. This is reasonable, since in that limit the mesons become impenetrable and thus do not allow the observation of the interior dynamics. However, when the hard-sphere-scattering part of formula (2) is removed, then one just obtains stronger resonances close to the poles of the $`K`$-matrix when the coupling constant $`\lambda `$ is increased.

For small coupling, both models give similar results, except that the real shifts for formula (2) can be much larger when the hard-sphere-scattering part is present. This is probably the reason why the bare states of Ref.  are always close to the central resonance energies.

We may therefore conclude that the behaviour of both models for moderate coupling is very similar, except for the interpretation of the bare states. In Ref. , we studied other consequences of the fact that mesons are not point particles, but finite distributions of constituent quarks.

Nevertheless, the fits of Ref.  to the data are too good for the corresponding model to be totally wrong. It might be that, with a small modification, the latter model would also yield the extra $`J^P=0^+`$ nonet of singularities and no related bare states, without destroying the excellent fits to the data.

## $`𝑱^𝑷\mathbf{=}\mathrm{𝟎}^\mathbf{+}`$ resonances nonets

In Table 1, we classify  the experimentally observed non-exotic scalar mesons into nonets.

The $`f_0(470208i)`$ comes in the Tables of Particle Properties, Ref. , under $`f_0(400\text{}1200)`$, but is not well established as a resonance, while the $`K_0^{}(727263i)`$ is not even mentioned in Ref. , although evidence for the existence of structure in that energy region has been reported . Moreover, a pole in the $`S`$-matrix is not necessarily observable as a clear resonance in meson-meson scattering, as we have argued before. Furthermore, the $`s\overline{s}`$ assignment of the $`f_0(980)`$  hints at the existence of a corresponding, most probably lower-lying, $`n\overline{n}`$ structure in $`\pi \pi `$ scattering .

The $`f_0(1370)`$ and $`f_0(1500)`$ resonances have been studied in many works , with a diversity of explanations as to their nature, out of which the above nonet classification is the most comprehensive.

Of all resonances in Table 1, the $`K_0^{}(1950)`$ does not seem to be well in place: the general level splittings of some 300 – 400 MeV do not agree with the jump of 520 MeV from the ground state to the first radial excitation of the scalar isodoublet. However, in the analysis of Ref.  one finds in Table 2 a set of possible singularities (in the third Riemann sheet) related to the $`K_0^{}(1950)`$ resonance, which all have real parts in the energy region 1.7 – 1.77 GeV. Moreover, in Ref.  a central resonance position of 1.82$`\pm `$0.04 GeV is reported for this resonance .

Furthermore, the $`f_0(1710)`$ is placed at 1.77 GeV in Ref. , which, nevertheless, does not alter the above classification, whereas both the $`f_0(2020)`$ and the $`f_0(2060)`$ need confirmation and might very well represent the same resonance. However, if there really exist two resonances in this energy region, then our classification indicates that one of them must be of a nature other than $`q\overline{q}`$.

In conclusion, one should note that several analyses find too many $`f_0`$ resonances, whereas in our analysis we lack an $`f_0(1840)`$.

Acknowledgements. One of us (EvB) wishes to thank the organisers of this conference for the kind invitation and warm hospitality at Protvino. We wish to thank Joseph Schechter, H. Noya, David Bugg, Carla Göbel, Brian Meadows, Michael Pennington, Teiji Kunihiro, Shin Ishida, Muneyuki Ishida, Bernard Metsch, and Tadayuki Teshima for valuable exchanges of ideas. This work was partly supported by the Fundação para a Ciência e a Tecnologia under contract numbers POCTI/35304/FIS/2000 and CERN/P/FIS/40119/2000.


text:
# CERN-TH/2001-262 October 2001 Lattice QCD on PCs? Plenary talk, XIX International Symposium on Lattice Field Theory, August 19–24, 2001, Berlin, Germany

## 1 INTRODUCTION

Parallel computers built from PC components are being increasingly used in many branches of science. The obvious advantages of such machines are that the hardware is relatively cheap and that the software environment (linux operating system, the gcc compiler suite and an implementation of the MPI communication library) complies with the established standards.

Some doubts have however been raised that PC clusters are good machines for numerical simulations of lattice QCD. According to the benchmark results presented at last year’s lattice conference , PC processors in fact appear to perform rather poorly in these calculations, particularly so when the lattices get reasonably big. Moreover, it remains unclear how well such clusters scale to large numbers of processors, where heat dissipation, component reliability and the network performance become limiting factors.

The main message of this talk is that lattice QCD programs can be accelerated by a large factor if use is made of the vector arithmetic unit and other enhancements of current PC processors. As a consequence (and also for various other reasons) the prospects for doing numerical simulations on PC clusters are now much brighter than they were only a year ago.

## 2 PC PROCESSOR PERFORMANCE

### 2.1 Multimedia extensions

The vector processing capabilites mentioned above have been added to the recent generations of PC processors to speed up multimedia applications. Depending on the brand and type of processor, the associated instruction sets go under the name of MMX, 3DNow!, SSE and SSE2. In most cases these instructions operate on short vectors of data in parallel and require one or two machine cycles to complete.

Evidently the number of arithmetic operations per cycle that can be performed also depends on the rate at which the data can be moved between the memory and the arithmetic units. To reduce the associated latencies, current processors support memory-to-cache prefetch instructions and streaming memory access modes, and they include a second-level cache memory that is clocked at the processor frequency.

### 2.2 Benchmark results

Most of the computer time in numerical simulations of lattice QCD is spent on the solution of the Dirac equation in the presence of a given background gauge field. There are different ways to solve the equation, but in all cases the program that applies the lattice Dirac operator to a prescribed fermion field is the one that dominates the execution time. The rest of the time is used for linear combinations and scalar products of Dirac fields.

In table 1 some benchmark results for these basic programs are reported, for the case of a PC with $`1.4`$ GHz Intel Pentium 4 processor and $`256`$ MB of PC800 RDRAM. The two lines in this table contain the processor times required for the specified task using 32-bit or 64-bit arithmetic. In the second column, for example, the times for the application of the Wilson–Dirac operator $`D_\mathrm{w}+m_0`$ to a given field $`\psi `$ are listed.

All figures quoted in the table are for a $`16^4`$ lattice, but on larger lattices they would be practically the same, because the memory latencies in these programs are effectively masked by the use of memory-to-cache prefetch instructions and a cache-optimized data layout. To accelerate the floating-point arithmetic, extensive use has been made of the SSE2 vector instructions that are supported on the Pentium 4 processor.

Compared with the program for the Dirac operator, the linear algebra routines (columns 3 and 4 in table 1) appear to be rather slow. This can be explained by noting that these programs spend most of the time to move data from the memory to the processing units. In a typical QCD code the linear algebra programs fortunately consume only a fraction of the total time so that their influence on the overall performance of the code is limited. This is illustrated by the figures for a standard conjugate gradient iteration quoted in the last column of table 1.

The bottom line then is that, with the current generation of PC processors and if use is made of their advanced features, it is possible to achieve sustained computational speeds that are about 10 times higher than those quoted last year .

In the remainder of this section, the vector unit of the Pentium 4 processor and its usage are discussed in some detail. Another issue that will be addressed is the cache management. To a large extent the impressive performance of current PC processors is in fact due to improvements of the memory system.

### 2.3 Streaming SIMD extension (SSE)

The vector unit on the recent Intel processors (Pentium III, Pentium 4, Itanium) has 8 registers that are denoted by xmm0,...,xmm7. They are 128 bits wide and can accommodate 2 double-precision or 4 single-precision floating-point numbers. The associated machine instructions operate on these numbers in parallel, i.e. the vector unit is a Single Instruction Multiple Data (SIMD) machine (see fig. 1).

All basic arithmetic operations (parallel addition, subtraction, multiplication and division) are supported and there are further instructions for data moving and shuffling. On the Pentium III only single-precision numbers can be manipulated in this way, but the Pentium 4 supports an extended instruction set, denoted by SSE2, that allows parallel double-precision and integer arithmetic to be performed in the SSE registers. In all cases the IEEE-754 standard is respected, which guarantees proper rounding and exception handling.

### 2.4 Programming example

When an ISO C compliant program is compiled with the gcc compiler, the generated object code will contain general purpose and x87 FPU instructions only. The SSE unit is thus not used. There is, however, a simple way to access the SSE registers from a C program and to operate on them. Let us consider the loop

for (i=0;i\<100;i++)
a\[i\]=b\[i\]+c\[i\];

where a, b, c are assumed to be arrays of single-precision floating-point numbers. Since there are no data dependences in this loop, we may use the addps instruction described above to do four additions at once. The code then becomes

for (i=0;i\<100;i+=4)
\__asm\__ \__volatile\__ (
"movaps %1, %%xmm0 \n\t"
"movaps %2, %%xmm1 \n\t"
"addps %%xmm0, %%xmm1 \n\t"
"movaps %%xmm1, %0"
:
"=m" (a\[i\])
:
"m" (b\[i\]),
"m" (c\[i\]));

The gcc compiler understands these lines, and if they are substituted in any C program, the executable will perform all 100 additions correctly.

The syntax of the inline assembly statement \__asm\__ \__volatile\__ (..) is documented in the manual pages of the gcc compiler. It has three arguments in the present case, a\[i\], b\[i\] and c\[i\], which are labelled from 0 to 2. The first line "movaps ..." moves four consecutive floating-point numbers from the address of argument 1 to the SSE register xmm0. The register then contains the array elements b\[i\],...,b\[i+3\]. Similarly the second line moves c\[i\],...,c\[i+3\] to xmm1. The parallel addition is then performed and the result is finally moved from xmm1 to the array elements a\[i\],...,a\[i+3\].

This example is rather trivial, but more complicated code sequences can be rewritten in essentially the same way, with perhaps several inline assembly statements in a row. Apart from those mentioned above, many more vector instructions are supported by the Pentium 4 processor that make it possible to perform almost any calculation in the SSE registers. The complete list of SSE2 instructions, together with detailed explanations of what precisely they do, can be downloaded from the Intel web pages .

In lattice QCD the SU(3) matrix times vector multiplication is one of the basic operations that is worth being coded for the SSE unit. On the Pentium 4 processor an average throughput of $`1.8`$ single-precision floating-point operations per cycle can been achieved in this case. The performance of the double-precision code is equally impressive and reaches $`0.9`$ operations per cycle.

### 2.5 Cache management

A common limitation of current computer architectures is that the memory system is often unable to deliver the data to the processor at a sufficiently high rate to keep the arithmetic units busy. The cache memory (which is significantly faster than the main memory) serves to enhance the data availability in those cases where the same data are used more than once.

In the program that applies the Wilson–Dirac operator to a given fermion field, for example, each spinor of the input field is processed 9 times when the program runs through the lattice. So if the spinor components that are needed at the current lattice point are already in the cache memory (because they have previously been used), it takes much less time to forward them to the arithmetic units than would otherwise be required.

Cache memories are expensive and usually too small to contain a complete lattice Dirac field. The spinors thus tend to be overwritten before they are reused. It is possible, however, to increase the cache-hit probability by dividing the lattice into small blocks that are visited one after another. This technique (which is referred to as strip-mining) is easy to implement and often results in significant speed-up factors.

In many cases the cache-hit rate can also be enhanced by explicitly moving the data to the cache memory slightly before they are going to be used. This can be achieved by inserting prefetch instructions such as

\__asm\__ \__volatile\__ (
"prefetcht0 %0" :: "m" (\*(address)));

at the appropriate places in the program. The important point to note is that these are executed “out of order”, i.e. while the arithmetic unit is busy on other data. In ideal cases the memory latencies can thus be masked almost completely.

Data are moved from the main memory to the cache memory in blocks that are referred to as cache lines. The cache line sizes and transfer rates for some of the recent PC processors and memory types are listed in table 2. When a prefetch instruction is issued, a cache line containing the byte at the specified address is moved to the cache memory. Dirac spinors occupy one or more cache lines, depending on the line size and the data alignment, and the programmer needs to take this into account to make sure that all components are prefetched.

### 2.6 Program portability

Inline assembly and the SSE2 instructions are system-specific and thus not as portable as an ISO C compliant code. It is conceivable that the C language will be extended to include short vector data types once the vector units become an industry standard, and it would then be left to the compiler to produce assembler code that makes use of any advanced features of the processor.

For the time being, program portability can be preserved by first writing an ISO C program for the specified task. Where this appears worth while, inline assembly code may then be added that is conditionally compiled, i.e. only when the macro SSE2 for example is defined. In this way it can be decided at compilation time whether some parts of the ISO C program are to be replaced by system-specific code or not.

### 2.7 FermiQCD

Programs that make use of the SSE unit along the lines described above are now part of the FermiQCD package . Many macros for SU(3) vector and matrix operations in the SSE registers can be found there, and the programs may also serve to illustrate the use of prefetch instructions. The package is freely accessible and can be downloaded from a webpage at Fermilab .

## 3 NETWORK ISSUES

### 3.1 Commercial PC clusters

In the course of the last two years or so, fully configured rack-mounted clusters of PCs have become a commodity. These machines usually have a switched network where any node can exchange data with any other node at full speed (see fig. 2).

As will become clear below, the bandwidth of the network is a critical parameter for clusters that are intended for numerical simulations of lattice QCD. The currently preferred network for such machines is Myrinet$`^{\mathrm{TM}}`$ , which provides a one-way bandwidth per channel of up to 250 Mbyte/s. Switches for 8 to 128 nodes are available, and these can, in principle, be connected to build clusters of any size.

### 3.2 How much bandwidth is needed?

Instead of developing a general formula for the required bandwidth, it may be more instructive to consider a realistic study case. So let us suppose that we have a machine with 32 nodes, each equipped with 512 MB of memory and a processor that delivers a sustained computational power of 1 Gflop/s in 32-bit arithmetic. As for the network, we assume that it provides a simultaneous bandwidth of 200 Mbyte/s per channel (the associated latencies are unimportant in the present context since the data can be transferred in large packages).

Such a machine is sufficiently big for quenched QCD simulations with Wilson fermions on a $`96\times 48^3`$ lattice. We may, for example, distribute the lattice to the nodes in such a way that each node is associated with a sublattice of size $`12\times 12\times 48^2`$ (see fig. 3). The question is then how much time is required by the communication overhead in the program that applies the Wilson–Dirac operator to a given fermion field.

When the program runs through the lattice, it reads the components of the input fermion field on the points of the local lattice and on all its boundary points. The Dirac spinors residing on the latter are stored in the memory of the logical neighbours of the current node and need to be communicated through the network before the calculation starts.

If 32-bit arithmetic is used, each Dirac spinor occupies 96 bytes of storage, and the total amount of data that must be moved to every node is thus about 11 MB. On our machine the time required for this is

$`t_{\mathrm{com}}`$ $`=2\times {\displaystyle \frac{96\mathrm{byte}}{200\mathrm{Mbyte}/\mathrm{s}}}\times \frac{1}{3}V=0.32\mu \mathrm{s}\times V,`$
where $`V`$ denotes the volume of the local lattice. Once all the data are available, the application of the Wilson–Dirac operator takes about $`1.4\mu `$s per lattice point, and the communication overhead is thus estimated to be 19% of the total time.

The result of this theoretical exercise shows that current network bandwidths and processing speeds are not out of balance. It should be noted, however, that the ratio of the boundary to the bulk volume of the local lattice, and hence the communication time, tends to increase proportionally to the number of nodes that are used for a given lattice. In many cases it may, therefore, be more efficient to simulate several independent lattices in parallel.

### 3.3 Possible improvements

PC processors are rapidly becoming faster, and it is clear that the network performance will have to be improved if the communication overhead is to remain small. Increasing the one-way bandwidth or having simultaneous bi-directional data transfers is certainly an option. With the current specification of the PCI bus (64 bit @ 66 MHz), transfer rates of up to 500 Mbyte/s can perhaps be reached.

An important reduction of the communication overhead might also be achieved if computation and communication can be made to overlap in time. The obvious difficulty here is that the PCI bus and the processor may be unable to access the memory concurrently with sufficient bandwidth, particularly so if the program is highly optimized.

## 4 PC CLUSTERS FOR LATTICE QCD

### 4.1 Recent installations

This year several dedicated machines have been installed or will shortly be delivered (see table 3). All these machines are commercial products that are equipped with MyrinetTM and come with the standard software environment.

The clusters listed in the last two lines of table 3 are being funded through the Scientific Discovery through Advanced Computing (SciDAC) initiative of the US government . This programme is intended to help closing the gap between nominal and actual performance of massively parallel computers and supports “research and development on scientific modeling codes, as well as on the mathematical and computing systems software that underlie these codes”.

### 4.2 Do PC clusters scale to Tflop/s?

Clusters with (say) 128 nodes currently reach sustained computational speeds of about 100–200 Gflop/s in lattice QCD programs. A small farm of such machines thus delivers an integrated computational power in the Tflop/s range. Having a farm of clusters is in fact not such an odd idea, since the total memory of each cluster will in many cases be sufficiently big to accommodate the lattice fields. Statistics can then be accumulated by running the same program (with different sequences of random numbers) on several machines.

Whether much larger numbers of processors can be integrated in a single cluster in a useful way is not obvious, however, and detailed studies will be needed before a solid answer to this question can be given. It seems likely, though, that the network switch will turn out to be a bottleneck when the cluster exceeds a certain critical size.

A development project, recently initiated at the University of Rome II and now supported by INFN, tries to overcome this difficulty by arranging the nodes in a logical $`n\times n`$ matrix. The nodes in each row and in each column are then connected to one another through $`2n`$ independent switches. This network geometry naturally maps to the block divisions of the physical lattice that are usually adopted. The routing of the data packages is thus simplified, and switching latencies are expected to grow only relatively slowly with the number of nodes.

### 4.3 Uses of small machines

When discussing computer performance and the need for ever more powerful machines, an important fact is sometimes forgotten that Giorgio Parisi pointedly described by :

“In general a brute force approach to numerical simulations is not paying and a lot of ingenuity is needed in order to obtain the wanted results.”

Innovation in computational strategies and algorithms is not only essential for progress in lattice QCD, but it also represents an intellectual challenge and is certainly one of the driving forces in this field. PC clusters are well suited to try out new physical and algorithmic ideas, because they are relatively easy to use and are in many cases sufficiently powerful for significant tests to be performed in a reasonable time. Even small machines can be very useful for this kind of work on which the large-scale projects eventually depend.

For illustration, let us consider the Polyakov loop two-point correlation function in the SU(3) gauge theory. When the time extent $`T`$ of the lattice is larger than (say) $`1\mathrm{fm}`$, the theory is in the confinement phase and the correlation function at distance $`r`$ decays exponentially, roughly like

$`P(r)^{}P(0)\mathrm{e}^{\sigma A},A=rT,`$
where $`\sigma 1\mathrm{GeV}/\mathrm{fm}`$ denotes the string tension. Now if the established simulation algorithms are used, the correlation function is obtained with statistical errors that tend to be essentially independent of $`A`$. As a consequence the amount of computer time needed for a specified relative precision is exponentially rising and areas $`A1\mathrm{fm}^2`$ or so are practically unreachable.

It is clear that significant progress in this field can only be made if better simulation algorithms are found that lead to an exponential reduction of the statistical errors. The recently proposed multilevel algorithm  is of this type and works exceedingly well, at least in the case of the Wilson plaquette action.

An impressive demonstration of this is given by the data plotted in fig. 4. At the chosen coupling, the lattice spacing $`a`$ is about $`0.17\mathrm{fm}`$, and the values of $`A`$ covered in the plot thus reach nearly $`3.5\mathrm{fm}^2`$. In spite of the fact that the signal decays over many orders of magnitude in this range, the equivalent of no more than about $`1000`$ hours of processor time on a PC with $`1.4`$ GHz Pentium 4 processor was required to keep the statistical errors below 2% at all values of $`T`$.

## 5 CONCLUSIONS

At this point it is quite clear that PC clusters are going to be widely used for lattice QCD simulations. While smaller machines are ideal for development work and physics projects at an early stage, large clusters or farms of clusters are certainly an option for big installations, where the aim is to maximize the total computational power for a given budget.

The cost per (32-bit sustained) Mflop/s of a commercial PC cluster with Myrinet$`^{\mathrm{TM}}`$ network is currently $3–4. In terms of cost-effectiveness, these machines are hence very competitive. How well PC clusters scale to large numbers of nodes, or whether it will be more useful to have a farm of smaller machines, remains to be seen. Further experience with both the hardware and the software will be needed, before it can be decided which is the best way to go.

Computers become more powerful at a breathtaking rate, but it seems unlikely that the hard problems in lattice QCD will gradually go away as a result of this alone. Improved numerical methods (such as those reviewed by Mike Peardon ) are in fact often crucial for the success of a numerical simulation project.

To explore a new approach can be fairly painful, however, if there is no reasonably fast computer around that can be freely used. PC clusters provide a good solution to this problem, and the fact that even small groups of lattice physicists may be able to afford them opens the way for more people to contribute to the development of lattice QCD.

I am indebted to Roberto Petronzio and Jarda Pech for many informative discussions on current network technologies. Thanks also go to Filippo Palombi for communicating his benchmark results for the Athlon processors and to Robert Edwards, Karl Jansen, Paul Mackenzie, Mike Peardon, Massimo Di Pierro and Peter Wegner for helpful correspondence.


text:
# 1 The LO parton–level cross section of the lightest Higgs boson as a function of the phase Φ for m_H₁=80 GeV (solid line), 90 GeV (dashed line), 100 GeV (dotted line), 110 GeV (dash–dotted line), 115 GeV (thick dashed line), and 120 GeV (thick solid line). We take the parameter set () with κ=1.6 (upper) and κ=2.0 (lower) and two values of tanβ=4 (left) and 10 (right).

## Acknowledgments

The work of SYC was supported by the Korea Science and Engineering Foundation through the KOSEF–DFG exchange program (Grant No. 20015–111–02–2) and the work of JSL was supported by the Japan Society for the Promotion of Science (JSPS).


text:
# Ordered phase and non-equilibrium fluctuation in stock market

## Abstract

We analyze the statistics of daily price change of stock market in the framework of a statistical physics model for the collective fluctuation of stock portfolio. In this model the time series of price changes are coded into the sequences of up and down spins, and the Hamiltonian of the system is expressed by spin-spin interactions as in spin glass models of disordered magnetic systems. Through the analysis of Dow-Jones industrial portfolio consisting of 30 stock issues by this model, we find a non-equilibrium fluctuation mode on the point slightly below the boundary between ordered and disordered phases. The remaining 29 modes are still in disordered phase and well described by Gibbs distribution. The variance of the fluctuation is outlined by the theoretical curve and peculiarly large in the non-equilibrium mode compared with those in the other modes remaining in ordinary phase.

The number of days in which prices of all stock issues in Dow-Jones industrial portfolio moved in the same direction is 80 within 3636 trading days in the period from 9/Jul/86 to 22/Nov/00, while its probability is $`2^{29}`$ when we assume Bernoulli trials. How can we explain the factor $`10^7`$ in the difference between these values? The methods and the concepts, as scaling and criticality, developed in material science have been applied to the study of financial markets Mandelbrot (1982); Mantegna and Stanley (1995, 2000). Recently a framework based on spin glass model to study the collective price changes of stock portfolios was proposed Maskawa ; H.Takayasu (in press). The application to 1-minute price changes in D-J portfolio made clear that the concept of energy works even in financial markets as well as the above physical concepts. Here we study the properties of daily price changes in D-J portfolio based on the same model and attempt an explanation of the factor $`10^7`$ by Gibbs factor of canonical distribution. Through this study we clarify the applicability of the spin glass picture to the price fluctuations in a wide range of time scale and find a novel feature of stock market. That is a non-equilibrium fluctuation mode on the point close to the boundary between ordered and disordered phases. D-J portfolio has been quenched into the unstable region, but does not reach equilibrium. The variance of the fluctuation, which is physically equal to susceptibility and is called (the square of) risk in financial economy, is outlined by the theoretical curve and peculiarly large in the non-equilibrium mode compared with those in the other modes remaining in ordinary phase.

In this paper, we concentrate on the statistics of the sign of price change Maskawa ; H.Takayasu (in press). The time series of price changes are coded into the sequences of up and down spins $`S_i=\pm 1`$ (i=1 to portfolio size N) and the Hamiltonian is expressed by long-range spin-spin interactions as Sherrington-Kirkpatrick model of spin glass Sherrington and Kirkpatrick (1975), which is given by

$$H[S,h]=\underset{<i,j>}{}J_{ij}S_iS_j\underset{i}{}h_iS_i,$$
(1)
We consider a portfolio as a subset of the whole stock market, and the complement of the subset works as heat reservoir. Various observable quantities are obtained as the statistical averages with Gibbs weight assigned to each configuration in equilibrium. The interaction coefficients $`J_{ij}`$ are constant jij  but not fixed yet. The external field $`h_i`$ is set to be zero in the analysis of actual data ext . We use Thouless-Anderson-Palmer (TAP) mean field theory D.J.Thouless et al. (1977); M.Mezard and M.A.Virasoro (1987); K.H.Fischer and J.A.Hertz (1991) whose fundamental equation is given by

$$m_i=\mathrm{tanh}(\underset{j}{}J_{ij}m_j+h_i\underset{j}{}J_{ij}^2(1m_j^2)m_i).$$
(2)
of $`m_i=<S_i>`$ to determine those coefficients. The partial differentiation of TAP equation (2) with respect to $`h_i`$ yields the equation

$$\underset{k}{}A_{ik}[𝐦]\chi _{kj}=\delta _{ij}$$
(3)
where $`\chi _{ij}=m_i/h_j`$ is susceptibility and

$$A_{ij}[𝐦]=J_{ij}2J_{ij}^2m_im_j+\delta _{ij}[\underset{k}{}J_{ik}^2(1m_k^2)+\frac{1}{1m_i^2}].$$
(4)
On the other hand, fluctuation-response theorem relates the susceptibility $`\chi _{ij}`$ and the covariance $`C_{ij}=<(S_im_i)(S_jm_j)>`$ as

$$\chi _{ij}=C_{ij}.$$
(5)
Substituting the equation (5) into (3), we can derive the relation between $`J_{ij}`$ and $`C_{ij}`$. Interpreting $`C_{ij}`$ as the time average of empirical data over the observation time, $`J_{ij}`$ are phenomenologically determined by the equation

$$J_{ij}=\frac{1+\sqrt{18C_{ij}^1m_im_j}}{4m_im_j}$$
(6)
where $`C_{ij}^1`$ is the (i, j)-element of the inverse of the covariance matrix. In the case with $`m_i0`$ the approximation $`J_{ij}C_{ij}^1`$ is applicable.

We investigate three datasets a, b and c of D-J portfolio consisting of N=30 stock issues in the framework of the model. The dataset a is the time series of stock prices sampled at 1-minute intervals in the period from 16-May-2000 to 21-Jun-2000, b is sampled at 10-minutes intervals in the period from 1-Dec-1999 to 11-Nov-2000, and c is the daily price in the period from 9-Jul-86 to 22-Nov-2000. The time evolutions of Intel’s stock price on each time scale are shown in Fig. 1 for an example. The time series of price changes are coded into the sequences of up and down spins. The covariance $`C_{ij}`$ for 435 pairs of i and j are derived from the coded data. Then the interaction coefficients $`J_{ij}`$ are calculated by the equation (6) fer . The energy spectra of the system (portfolio energy) are defined as the eigenvalues of the Hamiltonian $`H[S,0]`$. The probability density of portfolio energy for each dataset is empirically obtained from the relative frequency during the observation times as

$$p(E)=P(E\frac{\mathrm{\Delta }E}{2}H[S,0]E+\frac{\mathrm{\Delta }E}{2})/\mathrm{\Delta }E.$$
(7)
In our formalism that is given by the equation $`p(E)=n(E)e^E/Z`$ with the density of states $`n(E)`$ and the partition function $`Z`$ if the system is in equilibrium. The empirical probability weight $`p(E)n(0)/p(0)n(E)`$ is plotted in Fig. 2 and compared with Gibbs weight as the theoretical prediction. In the panels a and b we see a remarkable fit of equilibrium line to empirical data. From these figures, the price changes sampled at 1 and 10-minutes intervals are well described by the model and are close to equilibrium. On the other hand, in the panel c some deviation of data from theoretical line is observed. Monte Carlo simulation of the probability distribution of the system magnetization $`m=(1/N)_{i=1}^NS_i`$, which describes the degree of the alignment of price changes in portfolio, make more clear this difference in the nature of price changes on the different time scales. The result is shown in Fig. 3. We see the theory explains the empirical data very well in the panels a and b, while the dataset c shows entirely different profile from the theoretical prediction. Does this difference indicate the break down of the theory on the time scale of day? We investigate the theory in more detail by the mode analysis of TAP equation (2), and intend to explain this phenomenon in the framework of the theory.

The stability of the ordinary phase is analyzed by the linear analysis of the equation (2). In our case, that is

$$\underset{j}{}A_{ij}[\mathrm{𝟎}]m_j=h_i.$$
(8)
The diagonalization of symmetric matrix $`A_{ij}`$ of the equation (4) solve the equation (8) as $`m_\lambda =h_\lambda /A_\lambda `$, where $`A_\lambda `$ are eigenvalues of $`A_{ij}`$, $`m_\lambda =<\lambda i>m_i`$, and $`h_\lambda =<\lambda i>h_i`$ with the real orthogonal eigenvectors $`<\lambda i>`$. The critical temperature $`T_c`$ is determined by the equation

$$A_{\lambda _{min}}=0$$
(9)
$`A_{\lambda _{min}}`$ is the minimum eigenvalue of $`A_{ij}`$, which correspond to the maximum eigenvalue of the matrix $`J_{ij}`$ in the linear approximation. However $`\beta =1/T`$ is included into $`J_{ij}`$ in our case, so $`A_{\lambda _{min}}`$ itself fills the role of the measure of the distance from the critical point. In our systems a, b and c, $`A_{\lambda _{min}}=`$ 0.353, 0.195 and -0.018 respectively. The systems c is in the position slightly below the critical point, while the systems a and b remain in ordinary phase. In order to visualize the statement, TAP free energy D.J.Thouless et al. (1977) of the system c is shown in Fig. 4, in which the parameter space $`\{m_i\}`$ is projected on the mode $`m_\lambda `$. In Fig. 4, the upper and lower curves represent TAP free energy as the function of the modes $`m_\lambda `$ belonging to the maximum and the minimum eigenvalues of $`A_{ij}`$ respectively. The system c is in the unstable region, because the lower curve of the free energy is concave at $`m_{\lambda _{min}}=`$-0.067. The behavior of the price fluctuation belonging to this mode is entirely different from the remaining 29 modes. The probability of the magnetization fluctuation $`S_\lambda =<\lambda i>S_i`$ of the unstable mode and that of a ordinary mode belonging to the maximum eigenvalue $`A_{\lambda _{max}}`$ are plotted in Fig. 5 by filled circles and filled boxes respectively. The fluctuation distribution of the modes in ordinary phase are narrow and shows good agreement with the equilibrium distribution, while that of unstable mode is broad and does not reach equilibrium distribution.

The variance of each mode in the system a, b and c is plotted against eigenvalue $`A_\lambda `$ in Fig. 6. Those are compared with the theoretical curve $`1/A_\lambda `$ for infinite equilibrium system. The variance of the fluctuation, which is physically equal to susceptibility and is called (the square of) risk in financial economy, is outlined by the theoretical curve and peculiarly large in the non-equilibrium mode compared with those in the other modes remaining in ordinary phase. Closing to the critical point $`A_\lambda =0`$, however, we find a systematic deviation of data from the theoretical curve, which is due to the non-equilibrium property and the finite size effect of the systems.

The factor $`10^7`$ problem prompting us to introduce the interaction energy between stocks in the analysis of its collective fluctuation was partially explained by Gibbs weight in our model besides the single non-equilibrium mode. The physical picture of financial markets given here will be useful in the risk management and the selection of portfolios. The applicability of spin glass model to economical system seems to suggest the ubiquity of the applications to other fields such as competitive or frustrated ecological systems.


text:
# Efficient quantum computation of high harmonics of the Liouville density distribution

\[

## Abstract

We show explicitly that high harmonics of the classical Liouville density distribution in the chaotic régime can be obtained efficiently on a quantum computer . As was stated in , this provides information unaccessible for classical computer simulations, and replies to the questions raised in .

\]

In our Letter  we showed on the example of Arnold cat map that classical chaotic dynamics of exponentially many orbits can be simulated in polynomial time on a quantum computer. The Liouville density distribution $`P(x_i,y_j)`$ is encoded on a discretized lattice ($`2^{n_q}\times 2^{n_q}`$) using $`3n_q1`$ qubits organized in three registers. After each map iteration, the distribution is coded in the quantum state $`_{i,j}a_{ij}|x_i>|y_j>|0>`$ with $`P(x_i,y_j)=a_{ij}`$ and $`x_i=i/N`$, $`y_j=j/N`$, $`N=2^{n_q}`$ . One measurement of qubits in this basis gives one point in the phase space and therefore the distribution $`P(x_i,y_j)`$ can be obtained approximately in polynomial number of measurements. However the same information can be obtained via classical Monte Carlo simulation with a polynomial number of orbits, as it was discussed by us in  and later repeated in . Based on this observation, the comment  makes a general claim that no new information can be extracted efficiently from quantum computation of such classical maps (paragraph 3 in ) (see also the comment ). Here we show that this statement is incorrect. Indeed, the quantum Fourier transform (QFT) of $`a_{ij}`$ provides nondiagonal observables , namely the Fourier components $`\stackrel{~}{P}(k_x,k_y)=_{i,j}\mathrm{exp}(i2\pi (k_xx_i+k_yy_j))a_{ij}/N`$. They obviously contain important information relevant for the classical dynamics, and require $`O(n_q^2)`$ operations including measurement. On the contrary, all known classical algorithms will require exponential number of operations to obtain correct probabilities at high harmonics $`k_{x,y}N`$. Such harmonics are important since due to chaos a significant part of total probability is transfered to wavevectors with $`k\mathrm{exp}(ht)`$ where $`h`$ is the Kolmogorov-Sinai entropy, and $`t`$ is the number of iterations (see ). We note that the claim of  applies equally to the Shor algorithm, where all information is also encoded only in squared moduli of amplitudes, but where the QFT produces classically unaccessible information.

In Fig.1 we present the probability $`\stackrel{~}{P}(k_x,k_y)`$ in Fourier space for different times $`t`$. We note that a 2-dimensional (2d) QFT can be efficiently implemented by application of usual QFT to each register consecutively. The results show that $`\stackrel{~}{P}`$ is composed of well-pronounced peaks, most of which move with time to high wavevectors $`k`$. They remain stable in presence of noise in the quantum gates (e.g. top right panel in Fig.1 is unchanged if $`1\%`$ noise is added in each gate). The location and amplitude of these peaks can be extracted from a polynomial number of measurements of qubits after the 2d QFT.

For the Arnold cat map the dynamics in $`(k_x,k_y)`$ space is especially simple, given by $`\overline{k_x}=k_xk_y,\overline{k_y}=2k_yk_x\text{(mod}\text{N)}`$. However, generally this dynamics is very complicated. To exemplify this, we simulated the perturbed cat map $`\overline{y}=y+x+x^2,\overline{x}=x+\overline{y}\text{(mod}\text{1)}`$ (Fig. 1). It can be iterated in $`O(n_q^2)`$ operations on a quantum computer using modular multiplication. In this case, main peaks can be seen directly for short times, while for larger times a polynomial number of measurements of the first $`n_f`$ qubits  gives a coarse-grained image of $`|\stackrel{~}{P}(k_x,k_y)|^2`$, including very high harmonics, which are unaccessible to classical computation .

As concerns the issue of errors in the computation, raised in , it should be stressed that the main aim of  was to compare the errors natural for classical and quantum computers. In fact, we showed that the natural/minimal (last bit) errors for classical computer grow exponentially with number of map iterations while the natural errors in quantum gates do not destroy the time reversibility (Fig. 1 of ). As a further example, Fig. 1 of  clearly shows that the errors of relative precision $`10^8`$ (comparable with ordinary precision on Pentium III) completely destroy the reversibility of classical dynamics. At the same time, the quantum errors of relative precision $`10^2`$ in operations of quantum gates preserve the reversibility. Of course, since the quantum algorithm simulates the classical discretized map exactly, the last bit errors made on a quantum computer lead to exponential divergence of nearby trajectories and exponential drop of fidelity. This is clearly illustrated by Fig. 3 in . Thus exponential instability of classical chaos is preserved in quantum simulations. The quantum computer has enormous capacity in memory and precision which grows exponentially with the number of qubits; thus in the quantum case the last bit errors are much less important than for the classical computer with its limited memory space. However, the quantum computer has its own natural errors related to imperfections in gate operations. The results presented in  show that the quantum simulation is stable with respect to its natural errors in quantum gates operations while the classical computation is unstable with respect to its natural errors in the last bit of dynamical variables.

This work was supported in part by the NSA and ARDA under ARO contract number DAAD19-01-1-0553.