text:
# On the Evolution of Time Horizons in Parallel and Grid Simulations

## I Introduction

Parallel and grid computations for simulation of spatially-decomposable models with short-range interactions are arguably the most important topic in traditional applied computer science today. This is because of their application to simulations of many models in science, engineering, social science, manufacturing, and economics. We will first concentrate our discussion on parallel discrete event simulations (PDES), and will show the relationship to all short-ranged grid or parallel simulations in the last sections.

PDES are the execution of a single discrete event simulation program on a parallel computer or on a cluster of computers Fuji90 . It is a challenging area of parallel computing and has numerous applications in physics, computer science, economics and engineering. The number of applications are constantly growing in areas where extensive dynamical processes need to be simulated, especially those requiring a huge memory or wall-clock execution time. For such parallel simulations the system to be simulated is broken spatially into disjoint subsystems, and each processing element (PE) performs simulations on a particular subsystem. The simulated system jumps discontinuously from one state to another, these jumps are called an event. Thus the changes of state occur at discrete points in simulated time, although the time is considered to be continuous. The main challenge is to efficiently process with discrete events without changing the order in which the events are processed, i.e. preserving the causality in the simulation. One technique to help preserve causality is to introduce the idea of a local virtual time Jeff85  on a node or PE, which leads to a surface of local simulated times (LST).

For example, consider a kinetic Monte Carlo simulation for a 2d $`L`$$`\times `$$`L`$ ferromagnetic Ising model on a square lattice. The discrete events are spin flips. If the program is to be simulated using $`N_{\mathrm{PE}}`$ processors, each processor may be allocated an equal spatially-disjoint sublattice of spins. The average interval between two flips of the same spin varies for Metropolis-like algorithms from about 3.3 Monte Carlo steps per spin (MCS) at the critical point $`T_\mathrm{c}2.609`$ to 142.9 MCS at the low temperature of $`T=0.5`$ Over00 . Implementations of such PDES for kinetic Ising models have been performed using both conservative Korn99  and optimistic Over00  methods of preserving causality in the simulations. In the conservative implementation, a PE waits until causality is not violated before it proceeds with its calculations. In the optimistic implementation, if causality may be violated the calculation proceeds using some guess, and if this guess is incorrect the PE must roll back to an earlier state before any causality violations were present.

In this paper we investigate the dynamics of a number of PDES schemes. These parallel schemes are applicable to a wide range of stochastic cellular automata with local dynamics, where the discrete events are Poisson arrivals. We are interested in the evolution of the time horizon formed with the local simulation time (LST) of the nodes. In contrast with the previous work Korn00 ; Kola02 ; Kola03 ; Kola04 , where each PE manages only one node and communication between PEs is implemented according to the conservative scheme Fuji90 , we generalize the scheme so each PE processes a number of nodes which communicate conservatively, whereas nodes from different PEs communicate according to either the conservative or optimistic scenario. The simulated time horizon is analogous to a growing surface. The local time increments of the node correspond to the deposition of some amounts of “material” at the given element of the surface and the efficiency (which is the fraction of non-idling PEs) of the conservative scheme exactly corresponds to the density of local minima in the surface model. It was shown in Korn00 , that the density of the local minima does not vanish when the number of PEs goes to infinity. This remarkable result insures that the simulated time horizon propagates with a nonzero velocity and that the compute phase of the algorithm is asymptotically scalable.

The width of the time horizon in conservative PDES, after an early-time regime and before saturation, diverges with an exponent consistent with the Kardar-Parisi-Zhang one Korn00 . This scaling property is valid for the averages over the ensemble of the surfaces, or for the ensemble of the time horizons in the language of PDES. It is informative however to analyze the dynamics of a single realization of the time horizon as well. A single realization is strongly affected by the surface fluctuations, and the evaluation of a particular surface sometimes loses the full predictability. It is not clear how these fluctuations are connected with the turbulence of the Burgers solitons, although there are some similarity between LST-horizon evolution and the evolution of the surface slope described by noisy Burgers equations Gurbatov ; Chekhlov ; Fog-noisy .

In the original paper Korn00  each PE simulates only one node in a conservative manner. The efficiency of this algorithm is about $`1/4`$, on average, i.e. one PE out of four is working at any given time. Each PE sends messages to its neighbors about once in four time slices. One way to increase the utilization is to have each node contain a large portion of the lattice, then the utilization can be increased Korn99 ; Kola02 ; Kola03 ; Kola04 .

We analyze here a more efficient implementation of the Korniss et al. idea. Namely, we realize each node as a thread. Threads are distributed among processor elements so each PE is responsible for some number of threads. The threads within one PE communicate using the conservative PDES manner. Communications among threads on the same PE do not require the inter-processor communication latency time (one needs only system calls within a PE), while inter-processor communication requires calls to the I/O routines. Processes communicate within a single PE using system calls of the operating system. Inter-processor communication requires some I/O operations involved in addition to system calls. Threads were invented just to accelerate communications of the partially independent parts of the program. Their communication is supported by the kernel of modern operation systems Threads .

We have two choices for the inter-processor communication: either conservative and optimistic Fuji90 .

With the conservative inter-processor communications, the evolution of the time horizon will be exactly the same as in Korn00 , albeit the utilization would be larger than $`1/4`$ and close to unity for large enough $`\mathrm{}`$. Here $`\mathrm{}`$ is the number of nodes on a PE, so the system size simulated is $`L=\mathrm{}N_{\mathrm{PE}}`$. We assume that the check of the local minima condition (to avoid causality violations) between threads (nodes within a PE) is negligibly small in comparison with the time needed to process the event.

With optimistic inter-processor communications, the evolution of the time horizon becomes more complex. The optimistic scenario Fuji90  assumes that messages were not sent during some given time window. Causality is then checked. In the case where some node proceeds with broken causality, anti-messages have to be sent between processors, the corresponding events are rejected. This leads to a roll-back to an earlier time and state.

The possible scenarios can be understood taking into account the mapping of our problem onto surface growth dynamics. The messages sent by one PE to another are nothing but the boundary conditions in our algorithm. In fact, there are three (and not the two!) possible boundary conditions for the most left and the most right nodes for the chain of nodes within a PE: continuous, free and fixed. Continuous boundary conditions imply that the boundary nodes would follow the causality, taking messages from the neighbor node on the neighbor PE. So, this corresponds to the totally conservative case.

A more interesting solution is with fixed boundaries. In this case, a slope develops at boundaries of each PE. The angle of the slope depends on the mean of the event time intervals and a nearly flat top grows according to the conservative rule. This solution can be interpreted as a fixed soliton of the corresponding Burgers equation.

Free boundary conditions lead to an optimistic implementation of the algorithm. The evolution of the time horizon with the free boundaries have mixed features compared with the algorithms with continuous or fixed boundaries. This boundary condition will be analyzed elsewhere.

The algorithm we discuss here is a new PDES implementation scheme. The main purpose of the paper is a detailed analysis of this PDES algorithm scheme with fixed boundary conditions. The paper is organized as follows. In section II we review briefly the main ideas of PDES and the approach of Korniss et al. Korn00 . In section III we introduce our algorithm and discuss its behavior for simple realizations in one dimension, two dimensions, and for the solution of partial differential equations. We discuss in IV the results and provide a more general overview.

## II Conservative PDES

The PDES method is a tool capable of parallelizing any discrete event dynamic algorithm, even those which appear to be intrinsically sequential ones. One example is the development of a kinetic Ising model algorithm by Lubachevsky Luba88 , and its successful implementation Korn99  which preserves the original dynamics of the model.

### II.1 PDES and Kardar-Parisi-Zhang equation

Korniss et al. Korn00  developed an approach for the analysis of such algorithms. They mainly considered the case of one-dimensional systems with only nearest-neighbor interactions and periodic boundary conditions. (See ref. KornL00  for two and three-dimensional cases.) Each site of the Ising model is associated with one PE. Consequently the original model has $`N_{\mathrm{PE}}=L`$ PEs simulating $`L`$ node (or site). The number of nodes per PE are $`\mathrm{}=1`$. Update attempts at each node are independent Poisson processes with the same rate. (In an actual simulation the rate for the kinetic Monte Carlo for the Ising model depends only on the energy change for a single spin flip.) At each PE, the random time interval $`\eta _i`$ between two successive attempts is exponentially distributed.

Let us associate with PE number $`i`$ the value of a local simulated time $`\tau _i(t)`$, where we denote by $`t`$ the discrete time of the parallel steps simultaneously performed by each PE. We start with zero LST on all PEs, i.e. $`\tau _i(0)=0`$, $`i=1,2,\mathrm{},L`$. For $`t1`$ the LST evolves iteratively as

$`\tau _i(t+1)=\tau _i(t)+\eta _i(t)`$ $`\mathrm{if}`$ $`\tau _i(t)\mathrm{min}\{\tau _{i1}(t),\tau _{i+1}(t)\};`$
$`\tau _i(t+1)=\tau _i(t)`$ $`\mathrm{else}`$ , (1)
where $`\eta _i`$ are random exponential variables. Each time the PE number $`i`$ advances in time it sends messages to the right $`(i+1)`$ and left $`(i1)`$ PEs with the time stamp of it’s LST $`\tau _i(t)`$. This insures that the updating process (1) does not violate causality. This category of PDES is called the conservative PDES scenario.

The Korniss et al. algorithm is free of deadlock, since at least the PE with the absolute minimum LST can proceed. The efficiency of the algorithm is the fraction of nonidling PEs and exactly corresponds to the density of local minima of the simulated time horizon.

The iterative process (1) can be rewritten as

$`\tau _i(t+1)=`$ $`\tau _i(t)+\mathrm{\Theta }\left(\tau _{i1}(t)\tau _i(t)\right)`$ (2)
$`\times \mathrm{\Theta }\left(\tau _{i+1}(t)\tau _i(t)\right)\eta _i(t)`$
using the Heaviside step function $`\mathrm{\Theta }`$.

Introducing the local slopes $`\varphi _i=\tau _i\tau _{i1}`$, the density of local minima can be written as

$$u(t)=\frac{1}{L}\underset{i=1}{\overset{L}{}}\mathrm{\Theta }\left(\varphi _i(t)\right)\mathrm{\Theta }\left(\varphi _{i+1}(t)\right)$$
(3)
and its average

$$u(t)=\mathrm{\Theta }\left(\varphi _i(t)\right)\mathrm{\Theta }\left(\varphi _{i+1}(t)\right)$$
(4)
is the mean velocity of the time horizon, equal to $`0.246410(7)`$. Hence, the efficiency of the algorithm (in this worst-case scenareo) is about 25 per cent.

It was argued by Korniss et al Korn00  that the coarse-grained slope of time horizon $`\widehat{\varphi }(x,\widehat{t})`$ in the continuum limit evaluates according to the Burgers equation Kard86

$$\frac{\widehat{\varphi }}{\widehat{t}}=\frac{^2\widehat{\varphi }}{x^2}\lambda \frac{\widehat{\varphi }^2}{x}$$
(5)
and the coarse-grained time-horizon $`\widehat{\tau }`$, $`\left(\widehat{\varphi }=\widehat{t}/x\right)`$ obeys the Kardar-Parisi-Zhang (KPZ) equation

$$\frac{\widehat{\tau }}{\widehat{t}}=\frac{^2\widehat{\tau }}{x^2}\lambda \left(\frac{\widehat{\tau }}{x}\right)^2,$$
(6)
which should be extended with the noise to capture the fluctuations.

### II.2 Time horizon evaluation in conservative PDES

The Monte Carlo simulations of the process (1) showed Korn00  that the average width of time horizon

$$w^2(t)=\frac{1}{L}\underset{i=1}{\overset{L}{}}\left(\tau _i(t)\overline{\tau }(t)\right)^2,$$
(7)
where $`\overline{\tau }(t)=(1/L)\tau _i(t)`$, grows for appropriately chosen times (before the LST saturates) as $`w^2(t)t^{2\beta }`$ with the exponent $`\beta `$ close to the KPZ exponent, $`\beta `$$`=`$$`1/3`$. Running a single PDES on a parallel computer the LST horizon develops as a particular realization of the stochastic growth process, not as the average process.

It is known, that the dynamics of solitons in the noisy Burgers equation develops so-called Burgers turbulence Gurbatov ; Chekhlov ; Fog-noisy . Despite the similarity of the width evolution we have been unsuccessful in finding evidence for Burgers turbulence. The tail of the width distribution (Figure 1a of Korniss et al.) may be due to long-lived fluctuations, rather than due to moving solitons. Figure 1 shows the evolution of the time-horizon width for some realizations of this conservative PDES. A large fluctuation occurs in run number 4 at a time $`t68000`$. The momentary picture of the time horizon profile shown on Fig. 2 looks like a soliton. Nevertheless, a more detailed analysis is needed to claim that it is indeed a soliton. This evidence would be difficult to obtain, since the noise in the PDES conservative iterative process (2) strongly masks the expected soliton-like behavior.

The large fluctuation in the time horizon profile ($`1600`$) is comparable to the system size, $`N_{\mathrm{PE}}=L=10^4`$, for the case shown in Fig. 2. In the next section, we will see that the iterative process (2) has steady-state solutions under some boundary conditions, and the time horizon profile shown in Fig. 2 looks similar to it.

## III Freeze-and-shift algorithm

### III.1 Realizations of PDES.

The conservative PDES algorithm of Korniss et al. is the ideal scheme, in which each PE manages only one process and the LST is associated with the PE. We extend this to a different more general scheme than is present in other publications KornL00 ; Kola02 ; Kola03 ; Kola04 .

Consider the nodes to be vertices of a random graph. The links between the vertices are associated with the possible communications, at least in one direction. Imagine that we can group vertices in clusters, maximizing the connectivity within the cluster and minimizing the number of links between clusters. Let us call those nodes without external links bulk nodes and the rest of the nodes surface nodes. Then, we can map this random graph onto the parallel computer architecture associating one cluster of nodes with the one PE.

Practically, the nodes can be realized as threads Threads . Threads, sometimes called lightweight processes, share the same local memory with the other threads associated with the same PE. Hence, in a practical sense they do not require any extra communication between PEs to communicate with other threads on the same PE. Thus, the communication of the bulk nodes can be organized in the most optimal way, using the conservative PDES implementation.

We have to note that in our approach it is natural to associate the LST with nodes and not with the PEs, in contrast with previous work Korn00 ; KornL00 ; Kola02 ; Kola03 ; Kola04 . All nodes carry its own LST, even those belonging to the same PE.

The communication of surface nodes can be realized in a number of ways.

First, surface nodes can communicate in a conservative manner. The LST horizon will evolve as described by the Korniss et al. scenario discussed in the previous section. Nevertheless, the difference is that the average utilization is not given by $`u`$ (see expr. (4)). Let $`\mathrm{}`$ be the number of nodes per PE, so $`L=N_{\mathrm{PE}}\mathrm{}`$. The probability of having a chosen node not being at a local minimum of the LST is $`1u`$. Assuming complete randomness among the LST of the nodes, i.e. using a mean-field type of argument Kola04  for non-equilibrium properties of the time-evolving surface, gives that the probability when choosing $`\mathrm{}`$ nodes that all of them are not at a local minimum is $`\left(1u\right)^{\mathrm{}}`$. Since a PE can perform an update as long as any of its nodes are at a local minimum, the mean-field argument gives the utilization to be

$$u=1\left(1u\right)^{\mathrm{}}.$$
(8)
Our preliminary simulations show that the average utilization for $`\mathrm{}=1`$ depends on the type of noise (all of which have mean noise per step of unity). For uniformly distributed noise it is $`u_u0.267(4)`$. For Gaussian noise it is $`u_\mathrm{G}0.258(5)`$. These are to be compared with Poisson noise with $`u=0.246410(7)`$. Note that both Gaussian and uniformly distributed noise have an average utilization larger than $`1/4`$.

In the second realization of PDES, the surface nodes postpone messages within a discrete time window interval $`t_\mathrm{w}`$. The system would evaluate in the conservative PDES manner the bulk nodes which lie at a distance from the surface larger than $`d\left(\mathrm{}\sqrt{\mathrm{}^24t}\right)/2`$ as later given by Expr. (11). After the time $`t_1\mathrm{}^2/4`$ (see Expr. (9)) the freezing will reach the inner bulk nodes and the evaluation will stop. The system then could not propagate further without an exchange of the messages between PEs to preserve causality. The whole system will at that time be “frozen”. The second idea to implement this algorithm is the fast exchange of the messages. We propose to do that by redistributing nodes between PEs. A simple realization will be discussed below. We call this algorithm the “freeze-and-shift” algorithm (FaS). Note that the FaS algorithm effectively separates the inter-processor communication from the computations progressing on a PE. Thus computer architectures that are capable of simultaneously performing calculations and inter-processor communication can be used effectively, even for problems with fine-grained parallelism. Another potential application of the FaS algorithm is that it should allow simulations with fine-grained parallelism to be performed on calculations on the grid.

The third possible realization is the optimistic scenario of PDES, in which one assumes that all messages come in an order in which causality is not violated. After a discrete time $`t_\mathrm{w}`$, one has to check causality and send anti-messages to kill those which violate causality. The process of generating anti-messages can lead to “avalanches” of the time-horizon. These seem to finish with a time horizon profile similar to the one generated by the FaS algorithm introduced above. Clearly, the optimistic scenario is more time consuming: first, the anti-message generation is not a short process; and second, some substantial amount of work on the nodes processed will be rejected. As described on a PE with $`\mathrm{}`$ nodes, the optimistic scenario is some generalization of a combination of the first two scenarios. The understanding of the time horizon avalanche process in the optimistic scenareo can be interesting by itself, but will not be explored in this paper.

### III.2 Time horizon evolution for the FaS algorithm

In this subsection, for reasons of clarity, we discuss the simplest case of the FaS algorithm where the nodes effectively form a one-dimensional graph. Despite this simplicity, the one-dimensional case of the FaS algorithm can be applied to the parallel simulation of one-dimensional partial differential equations (see section IV). The one-dimensional FaS algorithm can also be used in simulations of systems that have been organized into disjoint spatial subsystem slices, with the slices forming a one-dimensional graph.

Let us assume there are $`\mathrm{}`$ nodes (or threads) on each of $`N_{\mathrm{PE}}`$ processor elements. For the case of Ising model simulations, each node carries just one spin of the system of $`L=\mathrm{}N_{\mathrm{PE}}`$ spins. Nodes within a PE communicate in the conservative manner. The three possible ways of the inter-PE communication correspond to different boundary conditions: the conservative scenario is associated with continuous boundary conditions; the FaS scenario corresponds to fixed boundary conditions; and the optimistic scenario - to free boundary conditions (and associated LST roll-backs).

As we mentioned already, our implementation of conservative PDES is analogous to that of the Korniss et al $`\mathrm{}=1`$ implementation Korn00  and further extensions to larger $`\mathrm{}`$ Kola02 ; Kola03 ; Kola04 . The only difference, and a major difference for algorithm design, is that in our case the LST is associated with the nodes and not with the PEs.

In the case of the FaS algorithm, the evolution of the LST horizon is quite different.

In Fig. 3 the two first steps of the possible evolution of the time horizon are sketched both for the case of the conservative scenario with $`\mathrm{}=1`$ and for the freeze-and-shift (FaS) scenario with $`\mathrm{}=4`$. Initially, at the time $`t=0`$ all $`\tau _i=0`$. In the first step of PDES all the nodes proceed — the condition (1) is fulfilled for all nodes. The difference between the two scenarious starts just after the first event (time step). In the top of Fig. 3, all the local minima will advance in the conservative scenario. For the frozen boundaries of the FaS algorithm, the left and right boundaries on each PE are fixed, so these minima cannot proceed. As the FaS algorithm progresses, the LST starts to develop a slope with an average angle $`\psi `$. The average value of this angle depends only on the mean of the noise for a single time slice, with it’s tangent being equal to that mean. This is because proceeding from a frozen node the question is: given that the next node can exceed the value of $`\tau `$ on the frozen node in the next time step what is the average distance at which its value of $`\tau `$ freezes ahead of the frozen node’s $`\tau `$? Since we are using a mean noise of unity, the average angle will be $`\psi =\pi /4`$ for all types of noise. This is seen in Fig. 4. Note that this profile would also develop if the algorithm always added unity to each updated node, as would be done in simulations of partial differential equations (see later on).

The process of growth will stop at the moment the time horizon reaches the hill knap. This average time can be calculated by assuming that it is given by the same time as the case where each node update advances $`\tau `$ by one unit. This gives that

$$t_1=2\underset{i=1}{\overset{\mathrm{}/2}{}}i\frac{\mathrm{}^2}{4}.$$
(9)
At that time the LST of each individual PE will look like a saw-tooth (Fig. 4). At time $`t_1`$ the entire time horizon will look like a saw with $`N_{\mathrm{PE}}`$ teeth. It is very important to realize that up to that time all PEs were with high probability busy with an efficiency of one <sup>1</sup><sup>1</sup>1We have ignored the extra book-keeping required for a single PE to decide which of its $`\mathrm{}`$ nodes can still be updated. This will often be the case in actual implementations. However, if this is not the case one expects only a constant difference between serial execution and parallel simulations, with the constant independent of $`N_{\mathrm{PE}}`$.

After a time $`t_1`$ the profile of the time horizon stops developing. Thus the PE can on average perform $`t_1`$ node updates before the LST profile freezes. Figure 4 shows a realization of this frozen steady-state profile. The derivative of the profile $`\varphi _i=\tau _{i+1}\tau _i`$ will represent a kink, the soliton-antisoliton pair. Contrary to the Burgers equation Fog-noisy , these kinks are not moving but are in the steady-state. We have an interesting result: by fixing the boundaries we select the soliton-antisoliton solution of equation (2).

The mean value of the LST-horizon for $`tt_1`$ can also be calculated. Consider the non-random case, starting with a flat distribution for $`\tau =0`$. The time for the LST-horizon to reach a plateau with all middle nodes at the same value of $`\tau `$ is

$$t=\underset{i=1}{\overset{\tau }{}}\left(\mathrm{}\tau \right)\tau \mathrm{}\tau ^2.$$
(10)
This equation can be solved for $`\tau `$ using the quadratic equation and choosing the physical root since $`\tau \mathrm{}/2`$, giving

$$\tau =\left(\mathrm{}\sqrt{\mathrm{}^24t}\right)/2.$$
(11)
For this configuration the average value of the LST is

$$\tau =\frac{\tau (\mathrm{}\tau )}{\mathrm{}}.$$
(12)
Then substituting for $`\tau `$ from Eq. (11) gives an expression for the time dependence for $`\tau `$ for times $`0tt_1`$. In a similar fashion an expression for the time evolution of $`w^2`$ can be derived.

At time $`t_1`$ the LST-horizon is frozen, and the question is how to proceed further with the simulation. The solution we propose is to redistribute the nodes between the PEs. For example, let us shift them by $`\mathrm{}/2`$ to the right (or left), so the tops of the hills (saw-teeth) will be at the boundaries of the processors and fixed for the next time window processing.

Futher evolution of the time horizon is illustrated in Fig. 5. The lowest curve is the LST-horizon depicted on Fig. 4 and cyclically shifted by $`\mathrm{}/2`$. The LST-horizon will grow on average until the time $`t_23t_1=3\mathrm{}^2/4`$ at which time the next steady-state frozen solution will be reached (the highest curve in Fig. 5).

The process can be repeated by alternating freezing of the boundaries for a time window interval not longer than $`2t_1`$ and shifting nodes between PEs by $`\mathrm{}/2`$ for each time window. Figure 6 illustrates this process for $`\mathrm{}=1000`$ nodes and $`N_{\mathrm{PE}}=5`$ for ten shift-cycles with the time window of $`2t_1`$. Consequently there are about $`20t_1`$ time steps (discrete events).

The square width $`w^2(t)`$ of the time horizon (LST) seems to evolve periodically with time between a minimum and a maximum value. The evidence for this is seen in Fig. 7 in which the LST evolution of $`w^2`$ for $`\mathrm{}=10^3`$ and $`N_{\mathrm{PE}}=5`$ is shown. For the time between $`t=2\times 10^4`$ and up to the $`t=2\times 10^5`$ the LST-horizon width grows with the effective exponent $`z=3/2`$, it then stops at the shift moment, and oscillations start. This exponent is not associated with any stochastic process, but reflects the deterministic process of the saw-tooth’s formation according with Expressions (9-12). Probably, the maximum values of the LST width would follow the KPZ exponent, and this is a fruitful subject for future research. Figure 7 seems to show that in the FaS algorithm one can limit (effectively “saturate”) the size of the surface width. Proven ways of saturating the surface width in conservative PDES simulations include imposing a fixed constraint on the width Kola02  and imposing small-world connections between PEs Korn03 . Although the freeze-and-shift method damps out the surface width, much larger studies would be required to see whether or not the governing universality class of the interface is still the KPZ universality class and the FaS algorithm just leads to a coarse-grained length in the KPZ equation.

The time window $`t_\mathrm{w}`$ for the freeze-and-shift algorithm can be chosen as any value in the interval $`1<t_\mathrm{w}t_1`$. We can choose them shorter then $`t_1`$, for example equal to $`t_1/2`$ as shown in Fig. 8 and Fig. 9 for the even and odd shifts, respectively. The difference between the maximum and minimum possible values of the LST will be smaller than in the case of larger $`t_\mathrm{w}`$.

### III.3 FaS on a two-dimensional lattice

The next interesting realization of the freeze-and-shift algorithm is the application to a two-dimensional model. Suppose we have $`N_{\mathrm{PE}}`$ each of which will simulate a $`\mathrm{}\times \mathrm{}`$ block of nodes (spins in kinetic Monte Carlo for the Ising model). The total system size is $`L^2=N_{\mathrm{PE}}\mathrm{}^2`$ where $`N_{\mathrm{PE}}`$ should be a perfect square. We associate $`\mathrm{}^2`$ sites (spins) with $`\mathrm{}^2`$ threads running on each of the $`N_{\mathrm{PE}}`$ PEs. The evolution of the LST horizon is described by the two-dimensional iterative process (compare with expression (2))

$`\tau _{i,j}(t+1)=`$ $`\tau _{i,j}(t)+\mathrm{\Theta }\left(\tau _{i1,j}(t)\tau _{i,j}(t)\right)`$ (13)
$`\times \mathrm{\Theta }\left(\tau _{i+1,j}(t)\tau _{i,j}(t)\right)`$
$`\times \mathrm{\Theta }\left(\tau _{i,j1}(t)\tau _{i,j}(t)\right)`$
$`\times \mathrm{\Theta }\left(\tau _{i,j+1}(t)\tau _{i,j}(t)\right)\eta _i(t)`$
with the product of four theta-functions at the right hand side and $`(i,j)`$ are coordinates of the lattice edges. For $`\mathrm{}=1`$ one expects that the average speed of time horizon for the conservative PDES scenario will be approximately $`u_21/8`$, two times slower than in one-dimensional case. Our computer simulation gives $`u_2=0.1205(2)`$ for a Poisson distribution of PDES arrival times with unit mean value. This value again depends on the distribution law of the random numbers used.

Figure 10 shows the steady-state solution of process (13) which looks like a pyramid with the slope $`1/2`$. This solution is achieved at a time given by the volume of the pyramid, $`t_{2\mathrm{sst}}\mathrm{}^3/12`$.

We have to note that up to the time $`t_{2\mathrm{sst}}`$ the efficiency of the FaS algorithm is practically unity. Shifting the nodes between the PEs can now be accomplished in many different ways. Shifting the nodes by $`\mathrm{}/2`$ in both directions of the lattice is equivalent to the sending of postponed messages in the language of PDES. Then, the LST frozen profile will take a time of double $`t_{2\mathrm{sst}}`$ before reaching the second frozen steady-state position. Repeating this process of freezing and shifting, we can evolve our nodes in parallel as far as we wish. As in the one-dimensional FaS algorithm, the time between shifting can be any value given by $`0<t_\mathrm{w}t_{2\mathrm{sst}}`$. Note that this shifting can usually be implemented very effectively on modern computers which have fast methods of copying whole blocks of memory between processors, and may have the ability to simultaneously perform calculations and message passing.

### III.4 FaS, partial differential equations, grid computing

The discussed algorithms can be effectively applied also to numerical solutions for solving partial differential equations. In this case one has to solve iteratively finite-differences equations defined on the lattice. For simplicity we will discuss partial differential equations of second order in the one-dimensional case. The common form can be written as

$$\psi _i(m)=F(\psi _i(m),\psi _{i1}(m),\psi _{i+1}(m),\psi _i(m1),\psi _{i1}(m1),\psi _{i+1}(m1);\mathrm{\Delta },\delta ),$$
(14)
where the index $`i`$ is associated with the space variable, the space increment is $`\delta `$, and $`m`$ is associated with the time variable for which the increment is $`\mathrm{\Delta }`$.

We have to divide the entire one-dimensional lattice into $`N_{\mathrm{PE}}`$ pieces, each with $`\mathrm{}`$ lattice nodes, and then associate each piece with a single PE. The LST time increment is equal to $`\mathrm{\Delta }`$, and is not random in this case. We freeze the nodes at the boundaries for each PE. The propagation of the algorithm on each PE will be the normal one, except the nodes that do not know the values of their own or neighboring nodes at both times $`m`$ and $`m+\mathrm{\Delta }`$ will be frozen. This condition of the frozen boundaries between neighboring PEs will form a LST horizon on each PE. For each time step $`\mathrm{\Delta }`$ the space coordinates of the LST will be $`\delta `$. If the calculation of the right-hand side of Eq. 14 needs a large time compared with memory shifting between PEs, this algorithm could be very effectively implemented on a parallel computer architecture.

The scheme can be generalized to large dimensions of the lattice, i.e. to the many-dimensional partial differential equations, in the manner demonstrated in the previous subsection for FaS PDES.

Grid computing should enable geographically distributed heterogeneous computations to be performed if appropriate algorithms are available. For previous implementations with fine-grained parallelism, such algorithms are not available. For example, for conservative PDES simulations implemented with one virtual time per PE, the calculation on a particular PE halts at irregular times (whenever the algorithm hits the surface node of the PE). However, conservative PDES implemented with one virtual time per PE allows for an approximately regular and calculable number of time steps $`t`$ that a PE can perform before it needs to halt to preserve causality. Furthermore, the FaS algorithm allows a decoupling of calculation and the communication. Both of these facts can be important for grid computing. One example is since grid communication paths are used by many people the time for communication between grid computers is not constant, the communication between grid computers can take place without the calculations on a computer having to wait for another computer. Furthermore, the communication may be timed so it is performed when others are historically utilizing the communications network the least (say during nights or weekends). Thus the FaS algorithm should allow grid computations to be performed for PDES simulations and for numerical solutions of partial differential equations.

## IV Discussion

We have proposed a new algorithm for parallel discrete event simulations (PDES) which allows for an effective realization on the parallel computers, clusters of computers, and in grid computing. We call this algorithm the FaS (Freeze-and-Shift) algorithm. It effectively decouples the computation phase and communication phase for these parallel computations. This allows, for example, the programmer to utilize fast block-memory-transfer commands. Furthermore, it should allow the efficient simultaneous execution of both computation and communication hardware when they are performed by separate hardware.

There are several essential points of our FaS algorithm. First, we associate a LST with the nodes rather then individual processing elements (PEs). Second, we group the nodes using some rule (see later) and associate each group with a PE (usually the processor running the one process). The nodes are realized as threads, which share the same memory within one process. This allows for a fast realization of the conservative PDES rules within a PE. PEs do not communicate with other PEs for some time interval window $`t_\mathrm{w}`$ (the frozen part of the algorithm), and after that time the message exchange is realized as part-of-the-memory shifting between PEs (the shift part of the algorithm). The last part can be very efficiently realized on some parallel architectures. The width of the LST-horizon, which characterizes the difference in the LST between different shifts is under control and depends on the number of nodes $`\mathrm{}`$ per PE. The FaS algorithm guaranties the efficiency of the simulations and idleness of PEs should be nearly equal to zero for balanced simulations.

The general scheme can be sketched as in Fig. 11, where we have partitioned the nodes into groups (associated with PEs during the freeze part of the FaS algorithm) and we freeze development of those nodes within the interfaces. The efficiency of the algorithm depends on the smallness of the ratio of the interface area to the bulk area, and on the new (alternating) interface area we have to create by the shifting of nodes between PEs (the alternating partitioning of the node manifold).

Analysis of several realizations demonstrate the efficiency of the general idea of the FaS algorithm. Actual implementations on parallel computers, heterogeneous compute clusters, and computer grids for problems of interest will determine the ultimate usefulness of the FaS algorithm.

###### Acknowledgements.

Supported in part by National Science Foundation grants DMR-0113049 and DMR-0120310, and by Russian Foundation for Basic Research. LNS thanks Department of Physics and Astronomy and the ERC Center for Computational Sciences at Mississippi State University for their hospitality.


format: text

text:
# References

Description of pp-interactions
with very high multiplicity at 70 GeV/c

E. Kokoulina

The Gomel State Technical University, Belarus,

## Abstract

The collective behavior of secondary particles in pp - interactions at 70 GeV/c is studied. A Two Stage Gluon Model is offered to describe processes with very high multiplicity. An active role of gluons is shown in multiparticle dynamics. The analysis of multiplicity distributions has revealed a possibility of a thermodynamic interpretation of these interactions. A mechanism of the soft photon production as a signature of the quark-gluon system is proposed.

PACS: 12.38Qk; 12.40Ee; 13.85Rm.

These investigations have been carried out in the framework of the project ”Thermalization”(JINR). This project is aimed at studying the collective behaviour of secondary particles in proton interactions at 70 GeV/c .

According to the present understanding of hadronic structure, based on quantum chromodynamics (QCD) , protons consist of quarks and gluons. We have developed a Two Stage Gluon model  to describe high energy multiparticle production (MP) in proton interactions. At the first stage of the model QCD and thermodynamical approaches are used. At the second stage (hadronization) a phenomenological description is applied.

After the inelastic collision of two protons a part of the energy of these moving particles is converted into the thermal one. Constituents of the proton, quarks and gluons, have enough energy to be described by pQCD, because the strong coupling reaches a small value. Quarks and gluons become liberated. Our model investigations have shown that the quark division of primary protons accompanied by hadron production in pp interactions at 70 GeV is absent. MP is realized by gluons. This study has confirmed the idea of P.Carruthers about a passive role of quarks .

We have used one of the most generally accepted methods to study multiparticle dynamics: the multiplicity distribution (MD) analysis. Two schemes of interactions are offered. They are distinguished only by the quark-gluon (QG) stage. If we want to study the gluon division inside the QG system (QGS), the scheme branch model is used. If we are not interested in what is going inside QGS the second scheme thermodynamic model should be applied. In the both schemes some of gluons (not all of them) leave QGS and convert to real hadrons. Using the thermodynamic interpretation we can say that active gluons are evaporated from hot QGS. After the evaporation they pass the stage of hadronization.

Processes of pp interactions at 70 GeV/c were investigated experimentally . MD of charged particles were limited to 20 secondaries. Among them there are $`n`$ charged mesons ( $`\pi ^+`$ or $`\pi ^{}`$) and two leading protons: $`p+pn\pi +2N`$. In Thermalization project it is planned to observe events with high multiplicities. The production of soft photons could give the information about the early stage of QG interactions .

The choice of the MP scheme is based on comparison with experimental data . Physicists from JINR (Dubna) used generator PYTHIA and obtained MD of charged hadrons . It was shown that PYTHIA generator does not agree with the experimental data at high multiplicities and has a deviation two orders of the magnitude at $`n_{ch}=18`$ equal to (Fig. 1). We have refused from the quark model  (Fig. 1) because of some specific features which we cannot accept and that is why we have built a scheme of hadron interactions for MD description by means of the modified Two Stage Model (TSM) . We consider that at the early stage of pp interactions initial quarks and gluons take part in the formation of QGS. They develop branch processes . Also as in TSM we have used the hypothesis of soft decoloration at the second stage

$$P_n=\underset{m=0}{}P_m^PP_n^H(m),$$
(1)
where $`P_n`$ \- resulting MD of hadrons, $`P_m^P`$ \- MD of partons (quarks and gluons), $`P_n^H(m)`$ \- MD of hadrons (the second stage) from $`m`$ partons.

At the beginning of this study we used the model where some of quarks from protons took part in the formation of hadron jets. But it turned out that parameters of that model had values which were differed a lot from the parameters obtained in $`e^+e^{}`$\- annihilation, and especially from the parameters of hadronization. It was one of the main reasons to give up the scheme with active quarks.

We dwell upon the model where quarks of protons do not take part in the jet production, but remain in the initial protons. All the remained hadrons are formed by gluons. We call these gluons the active ones. It is important to know how many active gluons are in QGS in the moment just after the impact. The simplest MD to describe this gluons is the Poisson distribution:

$$P_k=e^{\overline{k}}\overline{k}^k/k!,$$
(2)
where $`k`$ and $`\overline{k}`$ are the number and mean multiplicity of the active gluons, accordingly.

We begin our MD analysis with the branch scheme of gluons. To describe the MD of k gluons, we have used the Farry distribution

$$P_m^B=\frac{1}{\overline{m}^k}\left(1\frac{1}{\overline{m}}\right)^{mk}\frac{(m1)(m2)\mathrm{}(mk+1)}{(k1)!},k>1,$$
(3)
where m and $`\overline{m}`$ are the number of secondary gluons and their mean multiplicites.

At the second stage some of the active gluons may leave QGS and transform to real hadrons. We call such gluons the evaporated ones. Let us introduce parameter $`\alpha `$ as the ratio of the evaporated gluons, leaving QGS, to all the active gluons, which may be transformed to hadrons. We use binomial distributions for MD of the hadrons from the evaporated gluons at the stage of hadronization

$$P_n^H=C_{\alpha mN}^{n2}\left(\frac{\overline{n}^h}{N}\right)^{n2}\left(1\frac{\overline{n}^h}{N}\right)^{\alpha mN(n2)},$$
(4)
where $`\overline{n}^h`$ and $`N`$ are parameters of hadronization. They mean the average value and maximal possible multiplicities of the hadrons from one active gluon at the second stage. The hadron MD in the process of pp scattering in the Two Stage Gluon Model with branch(TSMB), is:

$$P_n=\underset{k=0}{}\frac{e^{\overline{k}}\overline{k}^k}{k!}\underset{m=k}{}\frac{1}{\overline{m}^k}\frac{(m1)(m2)\mathrm{}(mk+1)}{(k1)!}$$

$$\left(1\frac{1}{\overline{m}}\right)^{mk}C_{\alpha mN}^{n2}\left(\frac{\overline{n}^h}{N}\right)^{n2}\left(1\frac{\overline{n}^h}{N}\right)^{\alpha mN(n2)}.$$
(5)
The parameters are determined from the comparison with experimental data . They are $`N=40`$(fix), $`\overline{m}=2.61\pm .08`$, $`\alpha =.47\pm .01`$, $`\overline{k}=2.53\pm .05`$, $`\overline{n}^h=2.50\pm .29`$. We can conclude that branch processes are absent, since parameters $`\overline{m}`$ and $`\overline{k}`$ are equal to the errors. The fraction of the evaporated gluons is equal to .47. A maximal possible number of hadrons from the gluon looks very much like the number of partons in the glob of cold QG plasma of L.Van Hove . At the fixed parameter of hadronization $`\overline{n}^h`$ equal to $`1.63`$ (see below the thermodynamic model) the fraction of the evaporated gluons is about .73 (Fig. 1).

In the thermodynamic model without branches the active gluons which appeared in the moment of the impact, may leave QGS and fragment to hadron jets. We assume that the active gluons the evaporated from QGS have the Poisson MD as in (2). Using the binomial distribution for hadrons (4) and the idea of convolution of two stages (1), we obtain MD of hadrons in pp-collisions in the framework of the Two Stage Thermodynamic Model (TSTM)

$$P_n=\underset{m=0}{\overset{M}{}}\frac{e^{\overline{m}}\overline{m}^m}{m!}C_{mN}^{n2}\left(\frac{\overline{n}^h}{N}\right)^{n2}\left(1\frac{\overline{n}^h}{N}\right)^{mN(n2)}(n>2)$$
(6)
($`P_2=e^{\overline{m}}`$). The comparison (6) with experimental data  (see Fig. 2), gives the following parameter values: $`N=4.24\pm .13`$, $`\overline{m}=2.48\pm .20`$, $`\overline{n}^h=1.63\pm .12`$. In sum (6) we constrain the maximal possible number of the evaporated gluons equal to $`M=6`$. At the description of experimental data of $`e^+e^{}`$ annihilation the hadronization parameter $`N`$ was found equal to $`45`$ . We can see that our parameter $`N`$ obtained in TSTM coincides with this value. Both TSMB and TSTM describe the data equally well.

The maximal possible number of the charged particles from TSTM is $`26`$. It is interesting to get MD for neutral mesons. For this purpose we have taken the $`\pi ^0`$ mean multiplicity in pp-interactions at 69 GeV/c and used it for normalization. It is equal to $`2.57\pm .13`$ . We have determined the parameter of hadronization $`\overline{n}_0^h`$ for neutral mesons. It is equal to $`1.036\pm .041`$. We are based on TSM  to estimate the probability of the different hadron production from one active gluon.

MD for neutral mesons is shown in Fig. 3. From this distribution we see that the maximal possible number of $`\pi ^0`$’s from TSTM is equal to 16. MD for total multiplicity is shown on Fig. 4. We see that the maximal possible number of total number of particles in this case is equal to42.

The dependence of the mean multiplicities of neutral mesons versus the number of the charged particles may be obtained by means of MD for total multiplicity $`P_{n_{tot}}`$:

$$\overline{n}_0(n_{ch})=\underset{n=n_1}{\overset{n_2}{}}P_{n_{tot}}(nn_{ch})/\underset{n=n_1}{\overset{n_2}{}}P_{n_{tot}},$$
(7)
we take the Bayess theorem into account, $`n_1`$ and $`n_2`$ are determined only by conservation laws. Fig. 5 shows a difference with the data at small multiplicities. The noticeable improvement will be reached if we decrease the top limit at low multiplicities ($`n_{ch}10`$) to $`n_2=2n_{ch}`$ (Fig. 6). We don’t know what is happening in the region of VHMP. This behavior of $`n_1`$ and $`n_2`$ in (7) indicates that Centauro events with a large number of charged particles and practically no accompanying neutrals, may be realized in the region of VHMP. AntiCentauro events with a large number of neutral particles and with very small number charged ones must be absent.

In Two Stage Model with gluon branch it was established that several of active gluons are staying inside of hot QGS and don’t give hadron jets. New formed hadrons are catching up small energetic gluons which were free before this time. These hadrons are excited because they have additional energy at the expense of absorbed gluons. This energy may be thrown down by means of the photon radiation.

In project ”Thermalization” it is planned to investigate of soft photons (SP) . It was shown that measured cross sections of such photons are several times larger than expected from QED. For the explanation of the SP excess the phenomenological glob model  and the modified soft annihilation model Lichard and Thomson  are used.

We want to estimate the number of them. We consider that at certain moment QGS or exited hadrons may be in almost equilibrium state on. That’s why we try to use for the description of the massless bosons (gluons and photons) the black body emission spectrum

$$\frac{d\rho (\nu )}{d\nu }=\frac{8\pi }{c^3}\frac{\nu ^2}{e^{h\nu /T}1},$$
(8)
where $`\nu `$ is the energy of photon. These spectra help us to calculate the number of SP . The gluon density at the deconfinement temperature $`T_c160200`$ MeV can be estimated by comparison with relic one: $`\rho _{gl}(160)=0.13(fm)^3`$ and $`\rho _{gl}(200)=0.25(fm)^3`$. The number of gluons $`N_{gl}`$ in the hot QGS of size $`L^3`$, where $`L=2`$fm, will be: $`N_{gl}(160)1`$, $`N_{gl}(200)2`$.

Using the spectral spatial density of relic photons (8) it is possible to get the number of SP $`N_\gamma `$ in the region of size of our system (new formed hadrons). This size must be bigger than one in the gluon case. If the size of our system about 2 fm and average energy of photons 15-20 MeV/c the number of such SP will be the order of $`10^3`$.

I thank members of Organizing Committee for wonderful possibility to take part in ISMD XXXIII. I also are deeply indebted to V.Nikitin for considerable support of my study. I thank E.Kuraev for the help in the understanding of SP nature.


text:
# Discretisation effects in the topological susceptibility in lattice QCD

## I Introduction

QCD is distinguished from the constituent quark models of hadronic physics by the presence of light sea quarks in the vacuum and one of the principal goals of lattice field theory is to provide an accurate theoretical description of their effects.

The topological charge, $`Q`$, and its associated susceptibility, $`\chi `$, are especially sensitive to the presence and properties of the sea quarks. Chiral perturbation theory describes their variation with the quark mass, $`m_q`$. By comparing lattice Monte Carlo measurements with this expectation, we may understand how near we are to achieving the above aim.

Initial work on unimproved sea quark formulations (Wilson Bali:2001gk  and staggered Bitar:1991wr ; Kuramashi:1992vq ; Kuramashi:1993mv ; Alles:1999kf ; Alles:2000cg ; Hasenfratz:2001wd ) showed only weak evidence for the expected chiral ($`m_q0`$) suppression of $`\chi `$ relative to the quenched ($`m_q\mathrm{}`$) theory at an equivalent lattice spacing, $`a`$ (see Hart:2000wr ; Durr:2001ty ; Hart:2001pj  for reviews). A clearer chiral reduction was seen for improved fermion formulations, probably as a result of better chiral symmetry at finite lattice spacing Hasenfratz:2001wd . $`𝒪(a)`$–improved Sheikholeslami–Wohlert (SW, or “clover”) Wilson fermions showed a strong sea quark effect in a study at quasi–fixed lattice spacing, $`a0.1\mathrm{fm}`$ Hart:1999hy ; Hart:2000wr ; Hart:2000gh ; Allton:2001sk ; Hart:2001pj ; Hart:2001tv ; mean field improved SW quarks also showed an effect AliKhan:2001ym . Improved staggered fermion studies have also shown clear suppression of $`\chi `$ over a range of lattice spacings $`0.09a0.12\mathrm{fm}`$ Hasenfratz:2001wd ; Bernard:2002sa ; Bernard:2003gq . Studies using lattice sea quark formulations with chiral properties more comparable to those of the continuum have so far been limited to two dimensional models Durr:2003xs .

In the light of the improved staggered results, it is interesting to examine the effects of finite lattice spacing on the topological susceptibility from $`𝒪(a)`$–improved SW quarks. This is useful not only from the theoretical viewpoint, but also when comparing the relative costs of producing independent field configurations in Monte Carlo simulations (the topological modes are expected to be amongst the slowest to decorrelate).

In this paper we extend the analysis of two flavour $`𝒪(a)`$–improved SW results at $`a0.1\mathrm{fm}`$ (included in this work with due acknowledgment) to cover a comparable range of lattice spacings to the improved staggered fermion study, using lattice QCD ensembles produced by the UKQCD and QCDSF collaborations. Moving to small quark masses (the lightest yet produced for Wilson–like fermions), we find strong evidence for the chiral suppression of $`\chi `$. Finite volume effects are shown to be under control, and no strong autocorrelation is seen for the topological charge. We perform a joint continuum chiral and continuum extrapolation, and find good agreement with theoretical expectations.

The structure of this paper is as follows: in Section II we review the theoretical expectations for the chiral variation of the topological susceptibility in the continuum and on the lattice. We describe our measurements in Section III, and discuss their variation with $`m_q`$, $`a`$ and the lattice volume. We conclude with a discussion in Section IV.

## II The topological susceptibility

Four–dimensional, Euclidean gauge fields can be classified by an integer–valued topological charge,

$$Q=d^4x\frac{1}{2}\epsilon _{\mu \nu \sigma \tau }F_{\mu \nu }(x)F_{\sigma \tau }(x).$$
(1)
In the absence of a $`\theta `$–term, $`Q=0`$ and the topological susceptibility, $`\chi =Q^2/V`$, has a finite, non–zero limit as the volume $`V\mathrm{}`$.

As the mass of the sea quarks is reduced the topological susceptibility is suppressed, the leading order behaviour being Crewther:1977ce ; DiVecchia:1980ve ; Leutwyler:1992yt

$$\chi (m_\pi ^2)=\frac{(f_\pi m_\pi )^2}{4}$$
(2)
for two degenerate flavours, with corrections at $`𝒪(m_\pi ^4)`$. This relation should hold when $`V`$ is large enough that chiral symmetry is not restored, i.e.

$$x_{\mathrm{LS}}m_q\mathrm{\Sigma }V(f_\pi m_\pi )^2V1,$$
(3)
where $`\mathrm{\Sigma }`$ is the chiral condensate Leutwyler:1992yt  111 By the “pion” we mean the lightest pseudoscalar meson with valence quarks of the same mass as the two degenerate sea quarks flavours The PCAC relation predicts $`m_qm_\pi ^2`$ and we shall use the latter as a measure of the quark mass..

As $`m_q\mathrm{}`$, $`\chi `$ approaches the (constant) quenched value, $`\chi ^{(\mathrm{qu})}`$, from below. Higher order corrections to Eqn. (2) thus introduce a negative curvature at some intermediate quark mass. Various interpolating ansätze between the chiral and quenched regimes have been proposed Leutwyler:1992yt ; Hart:2000wr ; Durr:2001ty . In this paper we shall focus on one particular form, the large-$`N_c`$ form Leutwyler:1992yt , as motivated for SU(3) by Dürr Durr:2001ty :

$$r_0^4\chi =\frac{c_3c_0(r_0m_\pi )^2}{c_3+c_0(r_0m_\pi )^2}\mathrm{with}\begin{array}{c}c_0=(r_0f_\pi )^2/4,\hfill \\ c_3=r_0^4\chi ^{(\mathrm{qu})}.\hfill \end{array}$$
(4)
Assuming the physical values of $`r_0=0.5\mathrm{fm}`$ Sommer:1994ce  and $`f_\pi =93\mathrm{MeV}`$, then $`c_0=0.0139`$. This formula has been seen to work well in QCD Hart:2000wr ; Durr:2001ty ; Hart:2001pj  at fixed lattice spacing.

### II.1 Lattice QCD

Eqn. (2) also holds for lattice QCD at finite lattice spacing for fermions with sufficiently good chiral properties Chandrasekharan:1998wg ; Giusti:2001xh . In general, however, the lattice topological susceptibility is related to the continuum by additive and multiplicative renormalisation factors, $`Z`$ and $`M`$:

$$r_0^4\chi ^{(\mathrm{lat})}=Z^2r_0^4\chi ^{(\mathrm{cont})}+M.$$
(5)
These competing factors are in general functions of the lattice spacing and the quark mass (partly implicitly through the coupling, $`g^2`$). The upshot of this is that discretisation effects act to reduce $`\chi `$ below its continuum limit at large $`m_q`$, but to increase it in the chiral limit. A recent discussion can be found in Bernard:2003gq . Formally, before examining the sea quark mass dependence, we should first perform a continuum extrapolation of $`\chi `$ at fixed $`m_q`$ to remove the effect of these renormalisation factors. This is, however, not possible with the ensembles we have available (especially if we do not admit interpolation in $`m_q`$ of measurements at finite $`a`$). To attempt a full analysis of our data we must assume that the continuum and chiral limits commute sufficiently well that we can either perform the $`m_q0`$ limit first at fixed $`a`$, or that we can carry out a joint chiral and continuum extrapolation using a single formula.

The first approach was used in Hart:2000wr ; Hart:2001pj , fitting Eqns. (2) and (4) to data at $`a0.1\mathrm{fm}`$. Constraints of simulation mean that this study cannot contribute many ensembles at comparably light sea quark mass but different $`a`$. We cannot, therefore, extend the above method by carrying out chiral fits at other lattice spacings so we can look at discretisation effects in the fit parameters. Instead we follow the second approach. Allowing $`c_{0,3}`$ to have a leading order lattice spacing correction (quadratic for the lattice action considered here)

$$c_0(r_0)=c_{00}+\frac{c_{01}}{r_0^2},c_3(r_0)=c_{30}+\frac{c_{31}}{r_0^2},$$
(6)
in Eqn. (4), we obtain:

$$r_0^4\chi =\frac{c_{30}c_{00}(r_0m_\pi )^2}{c_{30}+c_{00}(r_0m_\pi )^2}+\frac{1}{r_0^2}\left(k_1+\frac{(r_0m_\pi )^2}{(c_{30}+c_{00}(r_0m_\pi )^2)^2}\left(c_{01}c_{30}^2+c_{31}c_{00}^2(r_0m_\pi )^2\right)\right).$$
(7)
We have added a general constant, $`k_10`$, to allow $`\chi `$ to have a non-zero chiral limit at finite lattice spacing, as might occur in the absence of exact zero modes in the lattice Dirac operator Hart:2000wr . A check of our results is consistency of the fitted parameters with the results for $`r_05`$ from Hart:2000wr ; Hart:2001pj .

## III Measurements

We have calculated the topological charge and susceptibility using ensembles generated by the UKQCD Allton:2001sk ; Irving:2002fx ; Ukqcd:2003  and QCDSF Booth:2001qp ; Khan:2003cu  collaborations. Details of these ensembles are given in Table 1 (taken in the main from Khan:2003cu ). The lattice spacings lie in the range $`0.08a0.12\mathrm{fm}`$. The SU(3) gauge fields are governed by the Wilson plaquette action, with two flavours of SW fermions. The improvement parameter, $`c_{\text{sw}}`$, has been chosen so that the leading order discretisation errors vary quadratically with the lattice spacing Jansen:1998mx . The exact Hybrid Monte Carlo (HMC) simulation algorithm is used, which avoids finite step-size errors.

We measure $`Q`$ using the method of Hart:2000wr : ten cooling sweeps are applied using the Wilson gauge action. The cooling action has relatively little effect on the topological susceptibility Bernard:2003gq . Ten cools strikes a good balance between adequate suppression of these ultraviolet dislocations and excessive destruction of the long range topological structure Hart:2000wr ; Bernard:2003gq . The choice of operator does not matter in the continuum limit, at least for the quenched theory Teper:1999wp , and a reflection–symmetrised “twisted plaquette” lattice topological charge operator is used. In Fig. 1 we show the variation of the topological charge as a function of Monte Carlo simulation time for three ensembles over a range of sea quark mass. The histograms of $`Q`$ (using unit width bins, although in the rest of the analysis the topological charge is not rounded to the nearest integer value) show good agreement with the expected Gaussian form. Although $`Q=0`$ within statistical errors, we opt (for consistency with earlier studies) to subtract terms in $`Q`$ from $`\chi `$. Both this and the choice of whether to round $`Q`$ to integer values affect $`\chi `$ by much less than one standard deviation.

This algorithm for measuring $`\chi `$ has already been studied in the quenched theory (the Wilson gauge action), and continuum extrapolation yields $`c_{30}=0.065(3)`$ and $`c_{31}=0.28(4)`$ Hart:2001pj .

The decorrelation of $`Q`$ is good (see Fig. 1, and Hart:2000wr ; Hart:2001pj ; Ukqcd:2003  for details of some other ensembles). More quantitatively, we measure the integrated autocorrelation, $`\tau _{\mathrm{int}}(Q)`$, using an implementation of the sliding window method Madras:1988ei , with results shown in Table 2 (scaled to units of HMC trajectories) 222Autocorrelations are notoriously hard to estimate in realistic ensembles, and for consistency of analysis here there are slight differences from earlier studies of the intermediate quark mass ensembles Hart:2000wr ; Hart:2001pj . With this in mind, the results presented here have been compared to those from a number of other algorithms for determining $`\tau _{\mathrm{int}}`$, and show a good degree of robustness.. With the exceptions of ensembles $`q_{8,12,15}`$ and $`u_{10}`$, there is weak evidence for a slow increase in $`\tau _{\mathrm{int}}`$ from $`𝒪(10)`$ HMC trajectories for the least chiral ensembles to $`𝒪(30)`$ trajectories for the most chiral. We plot these data for the $`16^332`$ ensembles in Fig. 2 (omitting the least chiral point). The trend is no stronger for any particular gauge coupling, $`\beta `$, but the statistical errors and available range of $`(r_0m_\pi )^2`$ make definitive statements difficult. Of the outliers, we note that ensembles $`q_{8,15}`$ have only a limited number of trajectories and, as we discuss later, $`q_{12}`$ is the smallest physical volume studied.

Barring these possible exceptions, we have confidence that the topological susceptibility may be estimated free from autocorrelation effects and the statistical errors quoted are from a jack–knife analysis using 10 bins. No increase in the statistical error estimates for $`\chi `$ was found with larger bins, as expected when comparing $`\tau _{\mathrm{int}}`$ to the bin size. Results for the topological susceptibility are given in Table 2 and plotted against quark mass in Fig. 3.

### III.1 Discretisation effects

There is a significant suppression of the topological susceptibility in the chiral limit. In Fig. 3 we show the leading order chiral expectation of Eqn. (2), using the physical value for $`f_\pi `$. There is approximate agreement for $`(r_0m_\pi )^23`$, although we note that discretisation effects tend to increase $`\chi `$ and exaggerate the range of agreement. We do not see clear evidence for $`\chi `$ deviating from a rapid decrease with $`m_q`$ at finite lattice spacing as in Bernard:2003gq , although critical slowing prevents exploration of sea quarks comparably light with those simulated in that study.

To quantify the effects of lattice artifacts, it is useful to consider $`\chi `$ over the full range of sea quark masses. Assuming that Eqn. (7) holds, we can use it to perform a chiral extrapolation of $`\chi `$ whilst allowing for (leading order) discretisation effects. We exclude ensemble $`u_{10}`$, whose statistical error is probably unreliable, and ensemble $`q_{12}`$ which has a small volume (see on). Ensembles $`q_{8,15}`$ also have large autocorrelation, and we remark that excluding them from the fitted data had no statistically significant effect. We fit only the dynamical results, with the quenched parameters fixed. In principle we could avoid this and fit the quenched data as well, assigning all such points an arbitrarily large value of $`(r_0m_\pi )^2`$. Such fits do not converge well, however. This was probably as a result of the data being arranged as two distinct clusters in $`(r_0m_\pi )^2`$.

We show the fits to Eqn. (7) in Table 3, with the $`k_1=0`$ results (Fit 1) plotted in Fig. 3 for a representative sample of lattice spacings, and in Fig. 4.

It is interesting to note that even for the heaviest sea quarks simulated, the discretisation effects act to increase $`\chi `$. To see a crossover to the quenched behaviour (where a decrease is seen), the fits suggest we need $`(r_0m_\pi )^225`$. It is thus clear that for all the ensembles the vacuum has been qualitatively altered from that of the quenched theory by the presence of the sea quarks. As expected, topological degrees of freedom appear to be more sensitive to such changes than most hadronic observables.

For $`r_0=5.0\pm 0.2`$ (the range studied in Hart:2000wr ; Hart:2001pj ), these fits suggest $`c_0=0.024(3)(1)(1)`$, where the errors arise from those on $`c_{00}`$, $`c_{01}`$ and the above range of $`r_0`$ respectively. This figure agrees very closely with $`c_0=0.22(7)`$ obtained in the previous work for Eqn. (4), and gives confidence in the consistency of our fits.

In addition to the errors arising from the statistical variation in the data, there is also a systematic error arising from the uncertainty in the quenched behaviour. We attempt to account for this by repeating the fits with $`c_{30}`$ and $`c_{31}`$ varying by one standard deviation in each direction: $`c_{00}`$ varies by $`\pm 0.003`$ under this, and it is clear that this is a sub–leading effect that can be ignored.

The fitted $`c_{00}=0.0135(26)`$ corresponds to $`f_\pi =91.6(8.8)\mathrm{MeV}`$. This agrees surprisingly well both with the experimental (and 2+1–flavour) QCD measurement ($`93\mathrm{MeV}`$) and measurements by other methods on $`N_f=2`$ SW actions ($`97.0(0.3)\mathrm{MeV}`$ in these conventions, rescaled from Aoki:2002uc ). The leading order lattice spacing correction to this increases $`f_\pi `$, as expected. The size of this correction is consistent with those obtained for other methods (see Fig. 35 of Aoki:2002uc , for instance, although the different gauge actions and operators will, of course, lead to different lattice artifacts).

It is useful here to consider whether there is any evidence for $`\chi `$ having a non–zero chiral limit at finite $`a`$. In Fit 2 of Table 3 we show the result lifting the restriction $`k_1=0`$. The fits are extremely unstable; the first column is the result of starting with all free parameters being zero, whilst the second takes the results of Fit 1 as the initial values. This instability probably arises from the most chiral data points being clustered around $`r_0=5`$. We conclude that whilst we cannot rule out a non–zero chiral limit, the data appears to be described more consistently when such a term is not present. For the remainder of the paper we therefore concentrate on the results from Fit 1.

### III.2 Finite volume effects

In Bernard:2003gq  the lattice volumes were sufficiently large that the topological charge on subvolumes (of size $`L^4`$) varied independently; indeed a reduction in statistical error on $`\chi `$ was achieved by exploiting this self–averaging. We did not find the same in the this study, and this prompts us to consider more carefully whether the topological susceptibility exhibits finite volume contamination.

In general we expect finite volume effects to lead to a reduction in $`\chi `$: a small lattice excludes large instantons. In the quenched theory this restriction is negligible once $`L/r_02.5`$ Lucini:2001ej . With dynamical sea quarks present, small volumes can, in addition, see a restoration of chiral symmetry with the leading order chiral variation now being $`\chi (m_q)^{N_f}(m_\pi )^{2N_f}.`$ Avoiding this requires $`x_{\mathrm{LS}}1`$, although the precise limit is not known. What is also not clear is how the behaviour interpolates between this and Eqn. (3) as we vary $`x_{\mathrm{LS}}`$; we do not address this in this study.

Finally, we need the lattice to be large in units of the pion correlation length, $`Lm_\pi 1`$. Spectroscopy suggests $`Lm_\pi 5.7`$ is necessary to render the low-lying hadron states free of finite volume effects Allton:1998gi . This limit is, however, associated more with pion exchange than phase structure and is probably unduly stringent for this study.

In Table 1 we show all three measures of FVEs. All the above limits seem well satisfied and we do not expect, or see, any sign of chiral symmetry restoration or significant finite volume effects. Certainly we do not see the smaller volumes lying consistently below the fitted curve.

The ensembles $`q_{12,13,14}`$ vary only in volume, and no strong suppression of $`\chi `$ is seen as $`V`$ is reduced. Nonetheless, we choose to exclude $`q_{12}`$, the smallest lattice, from the scaling analysis as the volume is particularly small.

We can attempt to quantify the effects of finite volume by repeating the fit with $`k_1=0`$, progressively excluding the ensembles with the smallest $`x_{\mathrm{LS}}`$ ($`u_5`$, $`u_4`$…) until the fits become unreliable. The fit parameters remained consistent within statistical errors. We conclude that finite volume effects are not significant at this level of statistical accuracy if $`x_{\mathrm{LS}}10`$.

## IV Summary

We have carried out a systematic study of the topological susceptibility, $`\chi `$, in lattice QCD with two light quark flavours over a range of lattice spacings and sea quark masses using ensembles generated by the UKQCD and QCDSF collaborations. We find compelling evidence for suppression of $`\chi `$ with decreasing quark mass. Finite volume contamination of the results was not seen, and $`x_{\mathrm{LS}}10`$ is necessary to avoid such effects on the topological degrees of freedom.

We observe a slight (and, in truth, statistically insignificant) increase in the integrated autocorrelation of the topological charge, $`Q`$, as the sea quark mass is reduced, although this effect is not marked despite a critical slowing of the HMC algorithm Irving:2002fx ; Ukqcd:2003 . The autocorrelation times are, in general, small enough not to bias the jack–knife estimates of the statistical errors on our data.

There have been suggestions that at $`\beta 5.2`$ the $`N_f=2`$ $`𝒪(a)`$–improved SW action exhibits evidence for a nearby phase transition Hart:2001fp ; Hart:2001tv ; Sommer:2003ne . This transition affects most strongly observables with the quantum numbers of the vacuum, and we do not expect (or find) a strong effect on the topological susceptibility. It may, however, induce additional mass dependence in $`r_0`$ Sommer:2003ne , although we see no signal in $`r_0^4\chi `$ (whose fourth power would presumably magnify any such effect, distorting the data away from Eqn. (7)). It may also lead to large autocorrelation Sommer:2003ne . Whilst a large integrated autocorrelation is seen for one ensemble at large sea quark mass, this is at $`\beta =5.29`$ which does not immediately fit the hypothesis of Sommer:2003ne . No unexpected trends were seen for $`\beta =5.2`$.

The chiral variation of the data was compared to a theoretically motivated ansatz for the behaviour of $`\chi `$ over the full range of sea quark masses. This was extended to allow for variation of $`\chi `$ with the lattice spacing, and a joint chiral/continuum extrapolation was carried out for $`0.08a0.12\mathrm{fm}`$ (over a factor of 2 in $`a^2`$). The behaviour was as expected: the continuum limit of $`f_\pi =91.6(8.8)\mathrm{MeV}`$ was remarkably consistent with physical expectations, with discretisation effects tending to increase its value (as seen in other determinations AliKhan:2001tx ).

The discretisation effects tended to increase the topological susceptibility on all the ensembles. This agrees with expectations for the chiral limit and is in contrast to the quenched theory, where they suppress $`\chi `$. This suggests that, relative to the quenched theory, there are fundamental differences in the vacuum due to sea quarks, even for the largest $`m_q`$. The topological susceptibility is more sensitive to these differences than, for example, most quantities in the light hadron spectrum.

In summary, the topological susceptibility is a sensitive probe of the vacuum and chiral properties of the lattice action. Although the $`𝒪(a)`$–improved SW action formally breaks the chiral symmetry at finite lattice spacing, the behaviour of $`\chi `$ suggests that the effects are not significant at the lattice spacings simulated. The chiral behaviour of this action is comparably good to that of other actions currently used in large–scale simulations, and makes it a suitable laboratory for studying topologically sensitive states, including the $`\eta ^{}`$ meson.

###### Acknowledgements.

The author is grateful to G. Schierholz and M. Teper for useful discussions; to D. Hepburn, A.C. Irving and D. Pleiter for providing estimates of quantities for ensembles $`u_5`$ and $`q_{15}`$ in Table 1; and to the Royal Society for financial support.