text:
# Δ-scaling and Information Entropy in Ultra-Relativistic Nucleus-Nucleus Collisions

thanks:  Supported by the Major State Basic Research Development Program under Contract No G200077400, the National Natural Science Foundation of China under Grant No 10135030, and the National Natural Science Foundation of China for Distinguished Young Scholars under Grant No 19725521thanks: Corresponding author. Email: ygma@sinr.ac.cn

## Abstract

The $`\mathrm{\Delta }`$-scaling method has been applied to ultra-relativistic p+p, C+C and Pb+Pb collision data simulated using a high energy Monte Carlo package, LUCIAE 3.0. The $`\mathrm{\Delta }`$-scaling is found to be valid for some physical variables, such as charged particle multiplicity, strange particle multiplicity and number of binary nucleon-nucleon collisions from these simulated nucleus-nucleus collisions over an extended energy ranging from $`E_{lab}`$ = 20 to 200 A GeV. In addition we derived information entropy from the multiplicity distribution as a function of beam energy for these collisions.

Liquid gas phase transition, critical fluctuation, fragment topological structure

Ultra-relativistic heavy ion collisions provide a unique means to search for a new state of matter, the Quark-Gluon Plasma (QGP), where quarks and gluons over an extended volume are de-confined \[1-3\] and a phase transition between hadrons and the QGP will occur \[4-6\]. There could be a discontinuity in global features, in particular, the entropy from particle multiplicities, in nuclear collisions associated with the onset of the phase transition. In this Letter we report the first application of the $`\mathrm{\Delta }`$-scaling method and information entropy calculation in ultra-relativistic heavy ion collisions.

In intermediate energy heavy ion collisions, $`\mathrm{\Delta }`$-scaling was proposed by Botet and Ploszajczak . They applied the $`\mathrm{\Delta }`$-scaling law to the INDRA data in intermediate energy heavy ion collisions (Xe+Sn, 25-100 MeV/nucleon) by using $`Z_{max}`$ (the maximum of charge in reactions) as order parameter  and found that the distributions of $`Z_{max}`$ obey the $`\mathrm{\Delta }`$ = 1/2 scaling law below 32 MeV/nucleon while they obey the $`\mathrm{\Delta }`$ = 1 scaling law above 32 MeV/nucleon. This indicates that a transition from an order phase to the maximum fluctuation phase (disorder phase) occurs around 32 MeV/nucleon. Recently, Ma et al. analysed the multi-fragmentation data for mass around $``$ 36 light nuclei systems. They found a similar change of $`\mathrm{\Delta }`$-scaling occurs when the excitation energy of the system around 5.6 MeV/nucleon. Combination of analysis with other probes of phase transition (eg., Zipf law  and maximum fluctuations ) reflects a transition of matter phases with the change of $`\mathrm{\Delta }`$-scaling . Since $`\mathrm{\Delta }`$-scaling is a useful tool to identify the phase transition, we try to make the similar analysis in ultra-relativistic heavy-ion collisions.

Botet and Ploszajczak proposed $`\mathrm{\Delta }`$-scaling to identify the transition in intermediate energy heavy ion collisions . The $`\mathrm{\Delta }`$-scaling law is observed when two or more probability distributions $`P[m]`$ of the stochastic observable $`m`$ collapse onto a single scaling curve $`\mathrm{\Phi }`$(z) if a new scaling observable is defined by:

$$z=\frac{(mm^{})}{m^\mathrm{\Delta }}$$
(1)
This curve is:

$$m^\mathrm{\Delta }P[m]=\mathrm{\Phi }(z)=\mathrm{\Phi }[\frac{mm^{}}{m^\mathrm{\Delta }}]$$
(2)
where $`\mathrm{\Delta }`$ is a scaling parameter, $`m^{}`$ is the most probable value of $`m`$, and $`m`$ is the mean of $`m`$. When $`\mathrm{\Delta }`$ = 1, this kind of scaling law is called the first scaling law that is caused by self-similarity of system. This self-similarity means that, if these distributions with different $`m`$ by a new kind of variable $`z`$, they entirely collapse on the same curve. In fact, the famous KNO scaling  is the special case that the $`\mathrm{\Delta }`$=1 scaling law holds with a stochastic observable of multiplicity of particles. The INDRA data was explored using the $`\mathrm{\Delta }`$ = 1/2 and 1 scaling laws and a phase transition was observed . If we assume that $`P[m]`$ is a Gaussian distribution, we have:

$$P[m]=\frac{1}{\sigma \sqrt{2\pi }}exp[\frac{1}{2}(\frac{m\mu }{\sigma })^2]$$
(3)
where $`\mu =m=m^{}`$ , $`\sigma `$ is the width of the Gaussian distribution, they both depend on incident energy. If this Gaussian distribution $`P[m]`$ obeys $`\mathrm{\Delta }`$-scaling law , we should have:

$$\mu ^\mathrm{\Delta }\sigma $$
(4)
On the other hand, the information entropy is observable such that it was proposed to describe the fluctuation and disorder of a system by Shannon . It is defined by:

$$H=P(m)lnP(m)𝑑m$$
(5)
where $`P(m)𝑑m`$ = 1. This indicates that the quantity of fluctuation of system and depend on the distribution $`P[m]`$. If $`P[m]`$ is the Gaussian distribution, we have :

$$H=ln\sigma +(1+ln2\pi )/2ln\sigma +1.419$$
(6)
In this work, we investigate the $`\mathrm{\Delta }`$-scaling law and the information entropy for p+p, C+C and Pb+Pb in high-energy collisions (20-200 AGeV) with help of LUCIAE 3.0 of Sa and Tai . The head-on collisions are simulated in this work. The LUCIAE is of a Monte Carlo model and is an extension of the FRITIOF . Here the nucleus - nucleus collision is described as the sum of nucleon-nucleon collisions. LUCIAE was improved in the following three aspects: (1) Re-scattering of final hadrons, spectator and participant nucleons are considered in the LUCIAE, because the role of re-scattering  can not be neglected in high-energy domain. (2) The LUCIAE implements the Firecracker Model that includes collective multi-gluon emission from the color fields of interacting strings in the early stage of relativistic ion collisions. (3) LUCIAE introduces the suppression of strange quarks and effective string tensor to make some parameters related to product of strange quarks and string tensor in the JETSET depend on incident energy, size of system, centrality etc. Generally the LUCIAE can explain many data very well .

First, we choose the particle multiplicity produced in each event as stochastic observable $`m`$ and simulated 10000 events for p+p at every energy point to obtain normalized multiplicity distributions correspondingly. Then we performed $`\mathrm{\Delta }`$-scaling. As a result, we find that the multiplicity distributions are approximately satisfied with the $`\mathrm{\Delta }`$ = 1 scaling law, as shown in Figure 1.

We calculated $`m=m_iP(m_i)`$ and $`m_{RMS}=\sqrt{(m_i^2(m^{})^2)P(m_i)}`$ in terms of these normalized multiplicity distributions directly. While we fitted these distributions with Gaussian distributions to obtain $`\mu `$ and $`\sigma `$. By comparing $`\mu `$ with $`m`$ and $`\sigma `$ with $`m_{RMS}`$, we found that these distributions are basically Gaussian. Table 1 shows the fitting parameters and information entropy ($`H_{direct}`$ or $`H_{gauss}`$, which was calculated directly or by Gaussian parameter, respectively (see details in Ref.  and in the following).

In order to better investigate the $`\mathrm{\Delta }`$-scaling, we defined a coefficient

$$L=\frac{m^\mathrm{\Delta }}{m_{RMS}},$$
(7)
which characterizes the validity of the $`\mathrm{\Delta }`$-scaling. We investigat its dependence of incident energy with different $`\mathrm{\Delta }`$ values. Figure 2 shows that the $`L`$-dependence on incident energy from $`\mathrm{\Delta }`$ = 0.5 to $`\mathrm{\Delta }`$ = 1.6 . As a result, we find that $`L`$ is nearly a constant in whole investigated range (20-200 AGeV) when $`\mathrm{\Delta }`$=1.35, i.e., the system obeys the $`\mathrm{\Delta }`$-scaling law most suitably with $`\mathrm{\Delta }`$ = 1.35 (see Fig.3a).

We also deal with another two systems of C+C (5000 events) and Pb+Pb (1000 events). Similarly, their multiplicity distributions basically obey Gaussian distributions. In the same way, the best $`\mathrm{\Delta }`$-scaling is obtained with $`\mathrm{\Delta }`$ = 1.00 for C+C and $`\mathrm{\Delta }`$ = 0.80 for Pb + Pb, respectively. Figure 3b and 3c illustrate these results.

Through $`\mathrm{\Delta }`$-scaling for different systems, we found that the distributions of particle multiplicity obey Gaussian distributions approximately, especially for C + C and Pb + Pb systems and they obey the $`\mathrm{\Delta }`$-scaling law in a wide energy range for a given system, which indicates that no phase transition and no change of reaction mechanism exist. This is expected for the simulated data because the underlying particle production dynamics in the LUCIAE is a smooth function of beam energy without a phase transition.

We calculate the respective information entropy in terms of these distributions of particle multiplicity with the method proposed by Ma in Ref. . Figure 4 shows the dependences of the information entropy on incident energy for the p + p, C + C and Pb + Pb systems. The information entropy increases with the incident energy and with the sizes of the system monotonously. Also, the calculated values of $`H_{direct}`$ are consistent with the value $`H_{gauss}`$ obtained from Eq. (6) in the Gaussian distribution. Again, no indication of phase transition exists.

In addition, we use the distribution of strange particles multiplicity  and the distribution of the number of binary collisions  as two stochastic observables to make the $`\mathrm{\Delta }`$-scaling plots and extract the information entropy in the same way. The similar $`\mathrm{\Delta }`$-scaling was obtained and a monotonous information entropy was also observed.

In summary, we have demonstrated, for the first time to our knowledge, the $`\mathrm{\Delta }`$-scaling of charged particle multiplicity, strange particle multiplicity and the number of binary collisions using simulated p + p, C + C and Pb + Pb collisions from $`E_{lab}`$ = 20 to $`E_{lab}`$ = 200 AGeV. The LUCIAE 3.0, which includes the final state interactions, was used for the simulation. The $`\mathrm{\Delta }`$-scaling means that these observables obey a certain kind of universal laws, regardless of beam energy and collision system. We found that the scaling values $`\mathrm{\Delta }`$ for charged particle multiplicity distributions are 1.35, 1.00 and 0.80 for p + p, C + C and Pb + Pb collisions, respectively. Moreover, the information entropy calculated from charged multiplicity distributions increases with the beam energy and with the colliding system size monotonously. Both the $`\mathrm{\Delta }`$-scaling and the entropy values show no dis-continuity as a function of beam energy as expected because the LUCIAE has no change of particle production dynamics, while they are associated with a phase transition in the simulated data. We expect that the $`\mathrm{\Delta }`$-scaling and the entropy variable can be a valuable tool to search for possible dis-continuities in nucleus-nucleus collisions associated with the onset of a QCD phase transition. Further checks for the models which incorporate the QGP phase transition are in progress.

Acknowledgement: YGM is grateful to Prof. J.B. Natowitz for calling his attention to the $`\mathrm{\Delta }`$-scaling in intermediate energy heavy ion collision.


format: text

text:
# The QCDOC supercomputer: hardware, software, and performance

## I Introduction

Lattice QCD directly performs the path integral for the QCD Lagrangian by Monte-Carlo integration on a computer. Space-time is discretised on a 4-torus, and a large number of snapshots of typical vacuum configurations is used to evaluate hadronic correlation functions non-perturbatively. The numerical integration scheme introduces finite-volume, discretisation, and statistical errors that can be removed with sufficient compute power: lattice QCD is systematically improvable.

With improved actions and extrapolation in the lattice spacing, fairly modest lattice sizes, such as $`32^3\times 64`$, are believed to be adequate for controlling the discretisation and finite volume effects on most hadronic observables.

Algorithms for including the effects of quark loops in the selection of typical vacuum configurations are numerically very expensive. The expense is believed to grow with a large power of the inverse quark mass for current best algorithms, and has thus far not been paid in full by lattice simulations, either choosing to ignore such effects (quenching) or simulating with artificially large quark masses. The goal for lattice QCD is to simulate in a region where the dynamical quark mass at least reliably connects to chiral perturbation theory, if not to the physical masses, and requires at least many tens of Teraflop years of computer power.

The notable characteristic of this problem is the need to focus ever more compute power on reducing the quark mass for a fixed problem size, rather than scaling up the problem size as more compute power becomes available. To this end QCDOC has been designed to allow efficient distribution of a single lattice QCD simulation over a very large multi-dimensional grid of a few tens of thousands of compute nodes.

At such extreme scalability the sparse matrix inversions involved require both global summation and nearest neighbour communication at a much shorter timescale than on smaller MPPs. Consequently both global summation time and nearest neighbour latency bite much harder, and both the QCDOC hardware and software have been designed to address these issues very effectively, with an order of magnitude improvement over traditional cluster technology on these key operations.

## II QCDOC hardware

### II.1 Overview

Continuing advances in the microelectronics industry have made it possible to integrate almost all components that make up a computer system on a single chip. This is known as system-on-a-chip technology. Using this technology, the individual processing nodes in a massively parallel computer can be greatly simplified. This is the idea behind the design of the QCDOC supercomputer: the processing elements consist of a single application-specific integrated circuit (ASIC) and an industry-standard DDR memory module. Large machines can then be build by simply adding many such processing elements.

The main ingredients of the QCDOC ASIC are

* 500 MHz, 32-bit PowerPC 440 processor core
* 64-bit, 1 GFlops floating-point unit
* 4 MBytes of on-chip memory (embedded DRAM)
* controllers for embedded and external memory
* nearest-neighbor serial communications unit with aggregate bandwidth of 12 Gbit/s in 12 independent directions
* other components such as Ethernet controller, interrupt controller, etc.

In the following subsections, the QCDOC hardware is described in more detail (see also Ref. lat01 ). The main advantages of the QCDOC design are

* high scalability: 50% sustained performance for typical applications on machines with several 10,000 nodes
* low price-performance ratio of $1 per sustained MFlops
* low power consumption
* high reliability

### II.2 The QCDOC ASIC

IBM is the leading vendor in the system-on-a-chip industry. The QCDOC ASIC was designed in close collaboration with IBM Research (Yorktown Heights) and is shown schematically in Fig. 1.

It contains standard cores from the IBM library as well as custom-designed components for faster memory access and high-speed off-chip communication.

The on-chip components are linked by three busses from IBM’s CoreConnectTM technology PLB :

* Processor local bus (PLB). This is a 128-bit wide, fully synchronous bus running at 1/3 of the CPU frequency. It contains independent read and write data busses. The PLB protocol (implemented by the PLB arbiter) is rather sophisticated, allowing for pipelining, split transactions, burst transfers (fixed and variable length, early termination possible), DMA transfers, programmable arbitration, and other features such as timeout, abort, etc.
* On-chip Peripheral Bus (OPB). This is a 32-bit wide, fully synchronous bus running at 1/6 of the CPU frequency. Its basic purpose is to off-load slower devices from the PLB bus.
* Device Control Register Bus (DCR). This is a very simple bus which is used to read and write various control registers in the on-chip devices. The 440 is the only master on this bus.

The IBM library components in the ASIC are

* The PowerPC-based 440 CPU core 440  with attached 64-bit IEEE floating point unit FPU . The 440 is a Book E compliant 32-bit processor with a 32 kByte prefetching instruction cache and a 32 kByte data cache with flexible cache control options (64-way associative, partitionable, lockable). The processor includes memory management with a 64-entry translation-lookaside-buffer which supports variable page sizes from 1 kByte to 256 MByte. It also features dynamic branch prediction and a 7-stage, dual issue pipeline. The target frequency of the 440 is 500 MHz, i.e. the peak performance is 1 GFlops.
The 440 also features a JTAG interface. JTAG (Joint Test Action Group) is an industry-standard protocol that allows an external device to take complete control of the processor. This functionality will be used for booting and debugging, see below.
* 4 Mbytes of embedded DRAM (or EDRAM) which is accessed with low latency and high bandwidth through a custom-designed controller, the PEC, see below.
* The PLB arbiter provides programmable arbitration for the up to eight allowed masters that can control PLB transfers. We are using six masters: the 440 instruction read, data read, and data write interfaces (the last two are channeled through the PEC), the EDRAM DMA, the SCU DMA, and the MAL DMA used by the Ethernet controller.
* The PLB-OPB bridge is used to transfer data between the two busses. It is the only master on the OPB and a slave on the PLB.
* The Universal Interrupt Controller (UIC) processes the interrupts that are generated on- and off-chip and provides critical and non-critical interrupt signals to the 440.
* The DDR controller is a slave on the PLB, capable of transferring data to/from external DDR (double data rate) SDRAM at a peak bandwidth of 2.6 GBytes/s. It supports an address space of 2 GBytes and provides error detection, error correction, and refresh of the off-chip SDRAM.
* The Ethernet media access controller (EMAC) provides a 100 Mbit/s Ethernet interface (it also supports Gbit-Ethernet, but we are not making use of this capability). The media-independent interface (MII) signals at the ASIC boundary are connected to a physical layer chip on the daughterboard. The EMAC is a slave on the OPB and has sideband signals to the MCMAL on the PLB.
* The DMA-capable Memory Access Layer (MCMAL) loads/unloads the EMAC through the sideband signals. It is a master on the PLB.
* The inter-integrated circuit (I2C) controller is a slave on the OPB. It is used to communicate with off-chip devices supporting the I2C protocol, such as the serial presence detect EPROM on the DDR DIMM or voltage and temperature sensors on the motherboard.
* The general purpose I/O (GPIO) unit is another slave on the OPB whose 32-bit wide data bus is taken out to the ASIC boundary. It is used, e.g., to drive LEDs, to control the global interrupt tree, and to receive interrupts from off-chip devices.
* The high-speed serial links (HSSL) used by the SCU each contain eight independent ports (four each for send/receive) through which bits are clocked into/out of the ASIC at 500 Mbit/s per port. The bits are converted to bytes (or vice versa) in the HSSL. The HSSL input clocks are phase-aligned by another IBM macro, the phase-locked loop (PLL).

In addition to the components provided by IBM, the QCDOC ASIC contains custom-designed components. Most importantly, these components provide essential support for the high-speed communications required in a massively parallel machine.

* Serial Communications Unit (SCU). The task of the SCU is to reliably manage the exchange of data between neighboring nodes with minimum latency and maximum bandwidth. The design takes into account the particular communication requirements of lattice QCD simulations. A schematic picture of the SCU is shown in Fig. 2.

The custom protocol governing the data transfers defines packets that contain a 64-bit data word and an 8-bit header containing control and parity bits. When the receive unit receives a packet, it first interprets the header, buffers the bytes from the HSSL, and assembles the 64-bit word. It then transfers the word to the receive register or passes it on to the send unit. The receive buffer can store three 64-bit words so that the send unit (in a neighboring node) can send three words before an acknowledgment is received. The functionality of the send unit is essentially the inverse of that of the receive unit. Send and receive operations can proceed simultaneously. Each send or receive unit is controlled by a DMA engine which then transfers the data between memory and a send/receive register. Each DMA engine is controlled by block-strided-move instructions stored in SRAM in the SCU itself.

A low-latency passthrough mode is provided for global operations. Because of the latencies associated with the HSSL, the most efficient method to perform global sums is “shift-and-add”, using a store-and-forward capability built into the SCU. The main advantage of this scheme is that the software latency of about 300 ns is paid only once per dimension, rather than for each node in this dimension.

The total end-to-end latency is estimated to be about 350 ns for supervisor transfers and about 550 ns for normal transfers. This is at least an order of magnitude lower than the latency associated with Myrinet. Since a write instruction from the 440 can initiate many independent transfers on any subset of the 24 send or receive channels, the latencies associated with multiple transfers can be overlapped to some degree. The total off-chip bandwidth using all 24 HSSL ports is 12 Gbit/s. In a 4-dimensional physics application only 16 of the 24 HSSL ports will be used, resulting in a total bandwidth of 8 Gbit/s. This provides a good match for the communications requirements of typical applications. Concrete performance figures are given in Sec. III.4 below.
* Prefetching EDRAM Controller (PEC). The PEC is designed to provide the 440 with high-bandwidth access to the EDRAM. It interfaces to the 440 data read and data write busses via a fast version of the PLB that runs at the CPU frequency and that we call processor direct bus (PDB). The PEC also contains a PLB slave interface to allow for read and write operations from/to any master on the PLB as well as a DMA engine to transfer data between EDRAM and the external DDR memory.

The access to EDRAM (which is memory-mapped) proceeds at 8 GBytes/s. ECC is built into the PEC, with 1-bit error correct and 2-bit error detect functionality. The PEC also refreshes the EDRAM. The latency of the PDB itself is 1-2 CPU cycles. This very low latency eliminates the need for pipelining. The maximum PDB bandwidth is 8 GBytes/s for read and write. However, due to internal latencies in the 440, the maximum sustained bandwidth is 3.2 GBytes/s.

The read data prefetch from EDRAM occurs in two 1024-bit lines. Three read ports (PDB, PLB slave, DMA) arbitrate for the common EDRAM. The coherency between PDB, PLB slave, and DMA is maintained within the PEC. Each read port has four 1024-bit registers that are paired in two sets to allow for ping-ponging between different memory locations. There are also two 1024-bit write buffer registers each for the PDB/PLB slave/DMA write interfaces.
* Ethernet-JTAG interface. As mentioned above, the 440 core has a JTAG interface over which one can take complete control of the processor. In particular, this interface can be used to load boot code into the instruction cache and start execution. This completely eliminates the need for boot ROM. The question is how the JTAG instructions should be loaded into the 440. (There are special tools that use the JTAG interface, but it would be impractical to connect one tool per ASIC for booting.) A solution to this question has already been developed at IBM Research, implemented using a field-programmable gate array (FPGA), that converts special Ethernet packets to JTAG commands and vice versa. This logic is now part of the QCDOC ASIC and will be used not only for booting but also to access the CPU for diagnostics/debugging at run time. The unique MAC address of each ASIC is provided to the Ethernet-JTAG component by location pins on the ASIC, and the IP address is then derived from the MAC address.

The QCDOC ASIC is manufactured using CMOS 7SF technology (0.18 $`\mu `$m lithography process).

### II.3 Mechanical design

The mechanical design is indicated in Figs. 35.

Two QCDOC ASICs will be mounted on a daughterboard, together with two industry-standard DDR SDRAM modules (one per ASIC) whose capacity will be determined by the price of memory when the machine is assembled. A maximum of 2 GBytes per ASIC are supported. The daughterboard also contains four physical layer chips (two per ASIC for their Ethernet and Ethernet-JTAG interfaces) and a 4:1 Ethernet hub so that a single 100 Mbit/s Ethernet signal is taken off the daughterboard.

32 daughterboards are mounted on a motherboard. The motherboard also contains eight 4:1 Ethernet hubs so that the total Ethernet bandwidth off a motherboard is 800 Mbit/s. Furthermore, the motherboard contains power transformers as well as temperature and voltage sensors.

8 motherboards are mounted in a crate with a single backplane. The final machine then consists of a certain number of such crates connected by cables.

### II.4 Networks

There are three separate networks: the high-speed physics network, an Ethernet-based auxiliary network, and a global interrupt tree.

The physics network consists of high-speed serial links between nearest neighbors with a bandwidth of 2$`\times `$500 Mbits/s per link. Transfers across these links are managed by the serial communications unit in the QCDOC ASIC as described above. The nodes are arranged in a 6-dimensional torus which allows for an efficient partitioning of the machine in software, described in more detail below. On a motherboard, the node topology is $`2^6`$, with three dimensions open and three dimensions closed on the motherboard (one of which is closed on the daughterboard).

The auxiliary network is used for booting, diagnostics, and I/O over Ethernet, with an Ethernet controller integrated on the ASIC. The Ethernet traffic to and from the ASIC will run at 100 Mbit/s. As mentioned above, hubs on the daughter- and motherboards provide a total bandwidth of 800 Mbit/s off a motherboard to commercial Gbit-Ethernet switches, a parallel disk system, and the host workstation (a standard Unix SMP with multiple Gbit-Ethernet cards).

The global interrupt tree consists of three separate interrupt lines that are visible across all partitions. An interrupt asserted by any ASIC is first transmitted to the top of the tree and then propagated down the tree to all other ASICs in the full machine. It will remain asserted until cleared by the ASIC from which it originated.

## III System Software

The QCDOC system software is minimally required to boot and manage the machine, load and run application code, and service application, communication and I/O requests. The system can be thought of as composed of three major parts.

* The front-end operating system, known as the qdaemon.
* The node operating system, known as the run-kernel.
* The run-time support libraries used by applications.

In this section we discuss the qdaemon and run-kernel and defer discussion of the run-time support libraries till a later Sec. IV.

### III.1 qdaemon

The qdaemon runs like a normal Unix daemon and is responsible for booting, managing and partitioning the ”back-end” grid of QCDOC nodes. The qdaemon is the sole point of access to the back-end for users and communicates with many nodes concurrently via RPC over the multiple gigabit Ethernet links.

The qdaemon is contacted either via a PBS based queuing system, or directly via a client program called the ”qcsh”. The qcsh is a modified shell, with additional built-in commands for sending requests to the qdaemon to perform operations on the QCDOC, such as running code on a partition.

The qdaemon is aggressively multi-threaded and supports many partitions and connected users simultaneously.

### III.2 Run-Kernel

The node operating system is a simple (non-preemptive) kernel. The overall design goal for the run-kernel is to run one compute process, and run it well.

Thus the kernel uses the 440’s MMU for memory protection, but not translation. This allows for protection of the O/S from errant user code, and for zero-copy communication with simple (i.e. non-virtual memory aware) hardware. Further, the entire memory map can be covered by the 64 entry TLB, so that TLB misses (which are a source of significant performance loss in many HPC machines) cannot occur on QCDOC.

The kernel does not implement scheduling so that the application is guaranteed 100% of the CPU. This is very important since a more traditional kernel on such a large and *very* tightly coupled machine (QCDOC will self-synchronise every 22 microseconds in some codes) would be impacted by a few nodes running their scheduler during any given dslash application.

The run kernel does, however, service both hardware and software interrupts. The features of the standard PowerPC architecture allow user code errors to be trapped cleanly and robustly from software interrupts. The kernel also services system call requests from user code, both to access the communication hardware through a very lean software layer, and to implement the standard C run-time environment (Cygnus newlib).

The kernel includes an Ethernet driver allowing for host-node communication, and an NFS client has been implemented allowing for file I/O from the nodes, both to the front end and to a parallel disk system composed from standard network attached storage.

Both run-kernel and application circular print buffers are maintained on the nodes. This allows for post-mortem readout of each node’s output, or optionally any subset of nodes can be configured at run-time to output directly to the console.

### III.3 Partitioning

As discussed, QCDOC is based on a six-dimensional toroidal grid of nodes.

A partition is a rectangular subvolume of the total machine and can be defined by the 6-coordinates of two nodes in the machine grid, namely the bottom-left, upper-right pair in six dimensions.

Each node has its own mapping from application directions to machine directions. By changing this mapping with node coordinate, we can change the topology and dimensionality seen by application code. This is done by successively folding machine axes together into a single application axis.

A simple 2d example of mapping an in principle non-periodic $`4\times 4`$ square into a 1d-torus of length 16 is shown in Fig. 6.

This process may be repeated, and an example of remapping a $`2^6`$ single motherboard is shown in Fig. 7. For any given machine it is possible to configure QCDOC for any application dimensionality from 1 to 6. From the point of view of the communications software, all partitions are six dimensional, but some variable number of these dimensions may be trivial. Calls in trivial application dimensions implement local copies from a node to itself, while calls in non-trivial application dimensions are mapped to the appropriate machine directions.

### III.4 Performance of optimised code

In this section we present the performance of optimised code on cycle accurate simulation of QCDOC. Table 1 shows the performance of key single-node assembler kernels. Most of these assembler kernels were in fact generated via a C++ program which will also output Alpha and Sparc assembler. This allows a very fair comparison with other contemporary RISC chips.

Anecdotally, the phenomenal bandwidth from EDRAM allows the FPU performance of QCDOC on double precision floating point to match that of high end workstations from L2 cache.

The very low latency communication on QCDOC (circa 550 ns) allows for extreme scalability of QCD sparse matrix multiplies. Performance is very acceptable down to local volumes as small as $`2^4`$, corresponding to a very spread out lattice. These benchmarks are shown in Table 2.

Finally, we display estimates of the scalability of a CG Wilson solver on large machines, including the 10 microsecond global summation time, in Fig. 8.

## IV User software and performance

User software can access the hardware capabilities of QCDOC by linking to the run-time library provided. Optimized lattice QCD computation kernels as well as communication routines will be available. The software infrastructure, including run-time library and data management and execution environment, will be in compliance with the SciDAC lattice QCD software effort to ensure greater portability and efficiency for a wide range of codes not necessarily optimized for the QCDOC hardware.

### IV.1 SciDAC software Program

The SciDAC National Computational Infrastructure for Lattice Gauge Theory initiative SciDAC  is a multi-institutional effort to create the software and hardware infrastructure needed for Teraflop scale lattice QCD simulations.

The goal of the software effort is to create a unified programming environment that will enable the US lattice community to achieve high efficiency on diverse multi-Teraflop scale hardware platforms. More specifically, the SciDAC software effort aims at providing a run-time environment for existing lattice QCD codes such as the Columbia Physics System (CPS) CPS , MILC code MILC , or the SZIN software system SZIN  so that they can be run and achieve high performance on various hardware platforms by adopting the programming environment provided by the SciDAC Software effort. A major part of this activity is to provide C and C++ API routines for lattice QCD (QCD-API) on which we will focus in this report. The performance of MILC routines with QMP on the QCDOC simulator is also presented.

The organization of the QCD-API is as follows.

* Level 1:

+ QCD Message Passing API (QMP) : Inter-node communication
    + QCD Linear Algebra API (QLA): Single node linear Algebra
* Level 2: QCD Data Parallel API (QDP): Lattice-wide operations with both communication and computation.
Example: parallel transport
$`\chi ^{}(x)=U_\mu (x)\chi (x+\widehat{\mu })`$ for all $`x`$ and fixed $`\mu `$, where $`\chi (x)`$ is the fermion field at position $`x`$ and $`U_\mu (x)`$ is the SU(3) gauge field connecting $`x`$ and $`x+\widehat{\mu }`$.
* Level 3: Highly optimized routines for lattice QCD

+ Dirac matrix inverter

    - Wilson, Clover
        - Staggered
        - Asqtad (Improved staggered)
    + Hybrid Monte Carlo (HMC) Asqtad force term

So far, many Level 3 QCD-API routines as well as QMP have been optimized and implemented for QCDOC.

### IV.2 Brief description of QMP

QMP aims to provide portable, low-latency, high-bandwidth communication routines suitable for lattice QCD. Here we describe some of the features of QMP. (For a more detailed description as well as for the C and C++ binding of the QMP routines, see Ref. QMP\_doc .)

* Point-to-point communication

+ Nonblocking (computation and communication can be overlapped)
    + Simultaneous, multi-directional transfer capability
    + Chained block/strided transfer capability
    + Separate routines for initialization and commencement of transfers $``$ opened communication channels can be reused to minimize overhead for repeated transfers
* Global operation

+ Global reduction
    Global Sum, Maximum, Minimum operations for integer, single and double precision numbers are available, as well as general binary reduction
    + Broadcast
    + Barrier

Since QCDOC has optimized hardware for lattice QCD, it is crucial that the implementation does not introduce an excessive software overhead. For the nearest neighbor communication, the QMP specifications are close to QCDOC native calls, and the additional software overhead is small. Both C and C++ bindings are implemented over QCDOC system calls. This also applies to QMP global operations. These are implemented using the store-and-forward capability of QCDOC so that the communication time for global operations grows only as $`L\times N_d`$, where $`N_d`$ is the number of dimensions and $`L`$ is the number of sites per dimension. Currently, the implementation of QMP is complete except for non-nearest neighbor communication, which is discussed in the next subsection.

### IV.3  Non-Nearest Neighbor Communication (NNC) on QCDOC

While communication patterns in lattice QCD are often to and from nearest neighbors only, a well optimized NNC QCD-API can be used for an implementation of routines with more complicated communication patterns, such as improved discretizations of the continuum QCD actions (Asqtad action asqtad ) and the MPI implementation on QMP (which is under way MPI\_QMP ). However, additional software routines are needed to manage NNC on top of the QCDOC hardware where all high-bandwidth connections are only to nearest neighbors.

Some of the desired features for NNC are as follows:

* Minimal buffer copying.
* Minimal overhead over QCDOC native communication call when destination is nearest neighbor.
* Capability to co-exist with nearest neighbor communication. Namely, both nearest and non-nearest communication can be open simultaneously and proceed without explicit user intervention or excessive delay on either communication channel.

The implementation strategy of NNC on QCDOC is as follows. After the outgoing data is packetized in the source node, the header is transferred and acknowledged via the interrupting communication channel (supervisor). The user data is then transferred using the non-interrupting physics network. If the receive node is not the destination of the packet, the packet is pushed onto the send/intermediate packets queue which sends the packets to the next node on the path until it reaches the destination. To avoid possible lock-ups and excessive delays, static, Manhattan-style routing is employed. It should be noted that in most lattice QCD routines, the communication pattern is uniform relative to the source, in which case Manhattan-style is close to being optimal. Figure 9 shows the data flow between QCDOC nodes in NNC. The implementation of NNC is under way MPI\_QMP .

### IV.4  Performance of MILC + QMP on QCDOC

MILC code MILC  is a body of C codes for SU(3) gauge theory, developed and maintained by the MILC collaboration. It has been ported to various parallel computers, workstations, and communication protocols, including MPI and QMP. Implementation of the QMP back-end for the MILC code was done by James Osborne (osborn@physics.utah.edu). We ran several lattice QCD kernels in MILC code with QMP on the QCDOC simulator and measured the performance.

As shown in Table 3, 15 $``$ 20% of peak performance is observed on the QCDOC simulator for many computationally intensive MILC C routines in single precision compiled by the XLC compiler XLC , which generated significantly more efficient assembly compared to GCC for PowerPC.

### IV.5 Improved action (Asqtad) force term

The Asqtad action asqtad  is one of the improved discretizations of continuum QCD fermion actions which exhibits smaller lattice spacing errors and flavor mixing. It uses many different paths connecting $`\chi (x)`$ and $`\chi (x+\widehat{\mu })`$ or $`\chi (x)`$ and $`\chi (x+3\widehat{\mu })`$ in contrast to the Wilson or staggered actions where only one gauge link $`U_\mu (x)`$ is used for each pair of neighbors. A diagram of the paths in the Asqtad action is shown in Fig. 10.

Calculation of the force term for the gauge link $`U_\mu (x)`$, $`S_{\mathrm{Asqtad}}/U_\mu (x)`$, mostly consists of parallel transport $`\chi ^{}(x)=U_\mu (x)\chi (x+\widehat{\mu })`$ and SU(3) vector outer products $`P(x)=\chi (x)\chi _{}^{}{}_{}{}^{}(x)`$. Because of the complexity of the gauge paths used in the action, the number of floating point operations for the Asqtad force term is much larger than for the application of the Dirac operator and the force term for simpler actions. The number of floating point operations per site is $``$ 250,000. Considering that the number of Flops for the Asqtad Dirac operator is only 1146 per site, and the typical number of conjugate gradient iterations is about 1000 for the mass used in typical lattice simulations, the force term calculation is a significant portion of the Hybrid Monte Carlo simulation, where a Dirac matrix inversion is alternated with a force term calculation, and thus should be optimized as well as the Dirac operator.

We found the performance of the Asqtad force term was initially rather low, only 3% of peak for $`2^4`$ and 6% for $`4^4`$ local lattice volumes. Upon examining the source code, we found that the communication channels are being created and destroyed for each parallel transport. This amounts to 30$``$60% of the total number of cycles for parallel transport. After modifying the parallel transport to reuse communication channels and combining small routines used in the outer product routine to eliminate function call overhead, the overall performance increased significantly to 12% for a $`2^4`$ and 13% for a $`4^4`$ local lattice.

Further optimization was done by eliminating function calls for computation routines defined for each lattice site within parallel transport and outer product routines. Together with loop unrolling, this made it possible to preload data into the L1 cache and registers to avoid cache misses. This improved the performance of computation routines by a factor of 1.5$``$1.7. The overall performance increased to 14% for $`2^4`$ and 20% for $`4^4`$ local volume, a 300$``$400% increase over the original performance and quite acceptable for a C routine.

## V Status and conclusions

QCDOC is a massively parallel computing architecture optimized for lattice QCD. The node has been designed to balance the floating point performance, memory bandwidth, and communication performance such that for QCD no single subsystem limits the performance. The design improves the network and memory subsystem performance relative to the floating point peak when compared with current commercial MPPs.

Simulation measurements of the nearest neighbour latency and global summation suggests at least an order of magnitude improvement over current commercial machines. This has enabled us to demonstrate very efficient use of the floating point unit that is maintained on remarkably small local volumes, such as $`2^4`$, and estimate scalability on typical lattices to machines as large as several tens of thousands of nodes.

The six dimensional mesh network is dealt with transparently by the operating system, and the machine can be dynamically partitioned to run multiple applications of any dimension from one through six. This should enable applications with similar characteristics to QCD (local communication on a regular multi-dimensional mesh) to make efficient use of the machine.

A run-time environment compliant with the SciDAC software effort will be available. The message passing interface defined by the SciDAC effort, QMP, provides portable, efficient communication routines for lattice QCD. The QMP implementation on QCDOC is complete except for non-nearest neighbor communication, which is in the process of being implemented. Performance numbers for lattice QCD routines from MILC code with QMP were presented. To take full advantage of the hardware capabilities of QCDOC provided via QMP, the user programs should be written in a way that minimizes repetitive overheads, as shown by the performance numbers for the Asqtad force term.

The first batch of QCDOC ASICs has been manufactured and is being assembled for initial testing. While we thank the editors for their patience, we regret having been unable to further delay the submission of these proceedings until after the hardware went online. Performance figures obtained from real hardware will be available in the near future.

###### Acknowledgements.

This work was supported in part by U.S. Department of Energy contracts DE-AC02-98CH10886 (CJ) and DE-FG02-91ER40608 (TW), by the RIKEN-BNL Research Center (TW), and by PPARC (PAB).


text:
# Fluctuations in network dynamics

## Abstract

Most complex networks serve as conduits for various dynamical processes, ranging from mass transfer by chemical reactions in the cell to packet transfer on the Internet. We collected data on the time dependent activity of five natural and technological networks, finding that for each the coupling of the flux fluctuations with the total flux on individual nodes obeys a unique scaling law. We show that the observed scaling can explain the competition between the system’s internal collective dynamics and changes in the external environment, allowing us to predict the relevant scaling exponents.

Recent advances in uncovering the mechanisms shaping the topology of complex networks networks  are overshadowed by our lack of understanding of common organizing principles governing network dynamics. In particular, we are far from understanding how the collective behavior of often millions of nodes contribute to the observable dynamical features of a given system, prompting us to continue the search for dynamical organizing principles that are common to a wide range of complex systems. To make advances in this direction we need to complement the available network maps with data on the time resolved activity of each node and link.

Traditional approaches to complex dynamical systems focus on the long time behavior of at most a few dynamical variables, characterizing either a single node or the system’s average behavior. To simultaneously characterize the dynamics of thousands of nodes we investigate the coupling between the average flux and fluctuations. Our measurements indicate that in complex networks there is a characteristic coupling between the average flux $`f_i`$ and dispersion $`\sigma _i`$ of individual nodes (Fig. $`1`$). To quantify this observation we plot $`\sigma _i`$ for each node $`i`$ in function of the average flux $`f_i`$ of the same node (Figs. 2 & 3). We find that for five systems for which extensive dynamical data is available the dispersion depends on the average flux as

$$\sigma f^\alpha .$$
(1)
Most intriguing, however, is the finding that the dynamical exponent $`\alpha `$ is in the vicinity of two distinct values, $`\alpha =1/2`$ (Fig. 2) and $`\alpha =1`$ (Fig. 3), suggesting that diverse real systems can display two distinct dynamical universality classes.

The $`\alpha 1/2`$ systems (Fig. 2): The Internet, viewed as a network of routers linked by physical connections, serves as a transportation network for information, carried in form of packets vazquez . Daily traffic measurements of $`374`$ geographically distinct routers indicate that the relationship between traffic and dispersion follows (1) for close to seven orders of magnitude with $`\alpha ^I=1/2`$ (Fig. 2a). In a microprocessor, in which the connections between logic gates generate a static network sole , information is carried by electric currents. At each clock cycle a certain subset of connections $`i`$ are active, the relevant dynamical variable $`f_i(t)`$ taking two possible values, $`0`$ or $`1`$. The activity during $`8,862`$ clock cycles on $`462`$ nodes of the Simple12 microprocessor indicates that the average flux and fluctuations follow (1), with $`\alpha ^m=1/2`$ (Fig. 2b).

The $`\alpha 1`$ systems (Fig. 3): The WWW, an extensive information depository, is a network of documents linked by URLs lawrence . As many websites record individual visits, surfers collectively contribute to a dynamical variable $`f_i^w(t)`$ that represents the number of visits site $`i`$ receives during day $`t`$. We studied the daily breakdown of visitation for $`30`$ days for $`3,000`$ sites scattered over three continents, determining for each node $`i`$ the average $`f_i^w`$ and dispersion $`\sigma _i^w`$. As Fig. 3a shows, $`\sigma _i^w`$ and $`f_i^w`$ follow (1) over three orders of magnitude with dynamical exponent $`\alpha ^w=1`$. The highway system is an example of a transportation network, the relevant dynamical variable being the traffic at different locations. We analyzed the daily breakdown of traffic measurements at $`127`$ locations on Colorado and Vermont highways. The results, shown in Fig. 3b, again document scaling spanning over five orders of magnitude with $`\alpha ^h=1`$. Finally, the river network is a natural transportation system banavar , whose dynamics is probed via time resolved measurements on the stream of several US rivers on $`3,495`$ different locations. While these fluctuations are driven by weather patterns, the relationship between the average stream and its fluctuations again follows (1) with $`\alpha ^r=1`$ (Fig. 3c).

To understand the origin of the observed dynamical scaling law (1) we study a simple dynamical model that incorporates some key elements of the studied systems. While the topology of these systems vary widely, from a tree (rivers) to a scale-free network (WWW, Internet), a common feature of the studied systems is the existence of a transportation network that channels the flux toward selected nodes. Therefore, we start with a network of $`N`$ nodes and $`L`$ links, described by an adjacency matrix $`M_{ij}`$, which we choose to describe either a scale-free or a random network networks . As the dynamics of the studied systems varies widely, we study two different dynamical rules. Model 1 considers the random diffusion of $`W`$ walkers on the network, such that each walker that reaches a node $`i`$ departs in the next time step along one of the links the node has. Originally each walker is placed on the network at a randomly chosen location and removed after it performs $`M`$ steps, mimicking in a highly simplified fashion a human browser surfing the Web for information. To probe the collective transport dynamics counters attached to each node record the number of visits by various walkers. To capture the day to day fluctuations on individual nodes we repeat independently $`D`$ times the diffusion of $`W`$ walkers on the same fixed network and denote by $`f_i(t)`$ the number of visits to node $`i`$ on day $`t=1,\mathrm{},D`$. As Fig. 4a indicates, the average flux and fluctuations follow (1) with $`\alpha =1/2`$. In Model 2 we replaced the diffusive dynamics with a directed flow process. In this case each day $`t`$ we pick $`W`$ randomly selected pairs of nodes, designating one node as a sender and the other as a recipient, and send a message between them along the shortest path. Counters placed on every node count the number of messages passing through. This dynamics mimics, in a highly schematic fashion, the low density traffic between two nodes on the Internet. As Fig. 4d shows, we find that Model 2 also predicts $`\alpha =1/2`$, indicating that the $`\alpha =1/2`$ exponent is not a particular property of the random diffusion model, but it is shared by several dynamical rules.

We can understand the origin of the $`\alpha =1/2`$ exponent if we inspect the nature of fluctuations in Model 1. In the $`M=1`$ limit walkers arrive to randomly selected nodes but fail to diffuse further, reducing the dynamics to random deposition, a well known model of surface roughening barabasi-book . Therefore, the average visitation on each node grows linearly with time, $`ft`$, and the dispersion increases as $`\sigma t^{1/2}`$, providing $`\alpha =1/2`$ barabasi-book . While for $`M>1`$ diffusion generates correlations between the nodes, we find that the fluctuations on the individual nodes, $`\sigma _i^{int}`$, continue to be dominated by the internal randomness of the walker arrival and diffusion process, following the $`\alpha =1/2`$ dynamical exponent poisson .

To understand the origin of the second ($`\alpha =1`$) universality class, we note that in real systems the fluctuations on a given node are determined not only by the system’s internal dynamics, but also by changes in the external environment. To incorporate externally induced fluctuations we allow $`W`$ (the number of walkers and messages in Models 1 and 2), to vary from one day to the other. Assuming that the day to day variations of $`W(t)`$ define a dynamic variable chosen from an uniform distribution in the interval $`[W\mathrm{\Delta }W,W+\mathrm{\Delta }W]`$, for $`\mathrm{\Delta }W=0`$ we recover $`\alpha =1/2`$. However, when $`\mathrm{\Delta }W`$ exceeds a certain threshold, in both models the dynamical exponent changes to $`\alpha =1`$ (Fig. 4b and e).

To understand the origin of the $`\alpha =1`$ exponent we notice that on each node the observed day to day fluctuations have two sources. For $`\mathrm{\Delta }W=0`$ we have only internal fluctuations, coming from the fact that under random diffusion (or random selection of senders and receivers in Model 2) the number of walkers (messages) that visit a certain node displays day to day fluctuations. For $`\mathrm{\Delta }W0`$ the fluctuations have an external component as well, as when the total number of walkers (messages) change from one day to the other, they proportionally alter the visitation of the individual nodes as well. If the magnitude of the day to day fluctuations is significant, they can overshadow the internal fluctuations $`\sigma _i^{int}`$. Indeed, if in a given time frame the total number of walkers or messages doubles, the flux on each node is expected to grow proportionally, a potentially much larger variation than the changes induced by the internal fluctuations. Therefore, for $`\mathrm{\Delta }W0`$ the external driving force, determined by the time dependent $`W(t)`$, contributes to the daily fluctuations with a dispersion $`\sigma ^{dr}(\mathrm{\Delta }W)=\sqrt{W(t)^2W(t)^2}`$. The total fluctuations for node $`i`$ are therefore given by $`\sigma _i^2=\left(\sigma _i^{int}\right)^2+\left(\sigma _i^{ext}\right)^2`$. As the effect of the driving force is felt to a different degree on each node, we can write $`\sigma _i^{ext}=A_i\sigma ^{dr}(\mathrm{\Delta }W)`$, where $`A_i`$ is a geometric factor capturing the fraction of walkers channeled to node $`i`$, and depends only on the position of node $`i`$ within the network. When $`\mathrm{\Delta }W=0`$, the external component $`\sigma ^{dr}`$ vanishes, resulting in $`\sigma _i^{int}=a_if_i^{1/2}`$, as discussed earlier, where $`a_i`$ is an empirically determined coefficient. When $`\mathrm{\Delta }W`$ is sufficiently large, so that $`A_i\sigma ^{dr}(\mathrm{\Delta }W)\sigma _i^{int}`$, then the fluctuations on each node are dominated by the changes in the external driving force. In this limit a node’s dynamical activity mimics the changes in the external driving force, allowing us to approximate the flux at node $`i`$ with $`f_i(t)=A_iW(t)`$. In this case we have $`f_i=A_iW(t)`$ and $`f_{i}^{}{}_{}{}^{2}=A_{i}^{}{}_{}{}^{2}W(t)^2`$, giving $`\sigma _i=\sqrt{f_{i}^{}{}_{}{}^{2}f_i^2}=A_i\sigma ^{dr}`$. As $`\sigma ^{dr}`$ and $`W(t)`$ are time independent characteristics of the external driving force, we find $`\sigma _i\sigma _i^{ext}=\frac{\sigma ^{dr}}{W(t)}f_i`$, providing the observed coupling (1) with $`\alpha =1`$. Note that this derivation is independent of the network topology or the transport process, predicting that any system for which the magnitude of fluctuations in the external driving force exceeds the internal fluctuations will be characterized by an $`\alpha =1`$ exponent.

These calculations imply that the fluctuations on a given node can be decomposed into an internal and an external component as

$$\sigma _i^2=a_i^2f_i+\left(\frac{\sigma ^{dr}}{W(t)}f_i\right)^2.$$
(2)
Therefore, increasing the amplitude of fluctuations $`\mathrm{\Delta }W`$ should induce a change from the $`\alpha =1/2`$ intrinsic or endogenous to the $`\alpha =1`$ driven behavior. To confirm the validity of this prediction, in Figs. 4c and f we show the average fluctuation $`\overline{\sigma _i}`$ over all nodes in function of the amplitude $`\mathrm{\Delta }W`$ of the driving force. For both models we find that for small $`\mathrm{\Delta }W`$ values $`\overline{\sigma _i}`$ remains unchanged, as in this regime $`\overline{\sigma _i}\sigma _i^{int}>\sigma _i^{ext}`$. However, after $`\mathrm{\Delta }W`$ exceeds a certain threshold, $`\overline{\sigma _i}`$ changes behavior, monotonically increasing with $`\mathrm{\Delta }W`$. In this second regime the fluctuations are driven by external forces, $`\overline{\sigma _i}\sigma _i^{ext}\overline{A_i}\sigma ^{dr}`$, and according to (2) we should observe $`\alpha =1`$. Indeed, we find that in both models the transition from the constant to the increasing $`\overline{\sigma _i}`$ regime (Figs. 4c,f) coincides with the crossover from the $`\alpha =1/2`$ to $`\alpha =1`$ (Figs. 4b,e). Note, however, that the gradual transition observed in Figs. 4b-e from $`\alpha =1/2`$ to $`\alpha =1`$ is a numerical artifact of the fitting process: in the transition regime the $`\alpha =1/2`$ and $`\alpha =1`$ scaling coexist on the same $`\sigma (f)`$ curve, giving an exponent that is different from $`1/2`$ or $`1`$. In reality the transition between the two regimes is sharp. To understand to what degree our findings depend on the specific simulation and model details we changed the topology from scale-free ba99  to random network and from undirected to directed network, as well as altering the nature of the external fluctuations by keeping $`W`$ constant in Model 1 but forcing the number of steps, $`M`$, to play the role of the stochastic external driving force. For each version we recover the transition between the $`\alpha =1/2`$ and $`\alpha =1`$ when the amplitude of the external fluctuations exceeds a certain threshold comment .

These results indicate that the $`\alpha =1/2`$ exponent captures an endogenous behavior, determined by the system’s internal collective fluctuations. In the studied model internal fluctuations are rooted in the randomness in the walkers’ arrival and diffusion; on the Internet they originate in the choices users make to where and when to send a message; for the computer chip they come from the alternating utilization of the various circuits, as required by the performed computation. In contrast, the $`\alpha =1`$ exponent describes driven systems, in which the fluctuations of individual nodes are dominated by time dependent changes in the external driving forces. Therefore, fluctuations of World Wide Web traffic, river streams and highway traffic are driven by such external factors as daily variations in the number of Web surfers, seasonal or daily changes in precipitation or daily variations in the number of drivers, respectively.

Of the two observed exponents our derivation indicates that $`\alpha =1`$ is universal, being independent of the nature of the internal dynamics or the network topology. There are no firm restrictions, however, on the scaling of the internal dynamics, raising the possibility that self-organized processes could lead to collective fluctuations that are characterized by $`\alpha `$ exponents different from $`1/2`$. Empirical evidence for potential intermediate $`\alpha `$ values comes from ecology, where (1) describes spatial and temporal variations of populations taylor . It is much debated, however, whether the observed scaling represent valid exponents, or only crossovers between $`\alpha =1/2`$ and $`1`$ anderson .

We are indebted to Jay Brockman and Steven Balensiefer for providing the data on the computer chip. This research was supported by grants from NSF, NIH and DOE.


text:
# Concrete uses of XML in software development and data analysis.

## I Introduction

Every few years there is a new technology that is touted as solving all our computing problems. In some cases these technologies turn out to be much less than they appeared to be, while other technologies do substantially advance our profession. XML is currently one of these technologies du-jour. As with all other technologies XML’s usefulness can only be determined once there has been some concrete implementations of the technology. This aim of this paper is to report concrete uses of XML within the IceCube experiment, and to demonstrate what benefits this technology actually delivers.

### I.1 IceCube Overview

The concrete examples of XML given in this paper are taken from the IceCube experiment. This experiment is a Neutrino telescope that will be deployed at the South Pole. It is the successor to the Amanda experiment which has been taking data since the mid-to-late 90’s. The experiment itself is relatively small in H.E.P. terms today with the collaboration comprising of around 150 physicists. This small size has enhanced the need for it to look to others for its software needs rather than to develop all its own software from the ground up. As will be seen throughout the rest of this paper the use XML has helped satisfy this requirement.

IceCube’s adoption of XML has also been facilitated by two other issues. The first is that there is very little legacy code or data from the Amanda experiment. What code there is is mainly restricted to the simulation portion of the software and any legacy data is in a home grown ASCII format which can be easy converted into XML. The second issue was the decision IceCube made that the core DAQ and data handling software should be written in Java. It turns out that many of the tools for, and advancements in, XML are taking place in Java first and then being migrated to other languages. This places IceCube is a good position to exploit the latest XML technologies.

## II XML

### II.1 Overview

XML (eXtensible Markup Language) is a derivative of SGML (Standard Generalized Markup Language). While SGML is, as it name suggests, very general this generality also means that it can be very time consuming to use it to implement solutions. At the other end of the spectrum a markup language such as HTML (HyperText Markup Language) is completely defined and while it is easy to use it is impossible to adapt to different needs. XML aims to take a middle road between the extremes of SGML and HTML by providing a flexible, customizable markup language that is easy to use.

The core idea of XML is that is it a text base language that defines the grammar of documents. This allows for the separation of the structure of a document from its content. Being text based also means that it can be independent of programming languages making it an ideal candidate for data transfer between programs. In fact this has lead to it becoming an industry standard for data description and exchange.

### II.2 Standards

One of the reasons XML has been successful at becoming an industry standard is that many aspects of its technology have become standards themselves, either as W3C w3c-ref  recommendations or as de-facto standards due to widespread adoption. The following is a list of the various area of XML technology and the standards which apply to them (The emphasized standards are W3C standards.)

Document Object Model (DOM) and Simple API for XML (SAX).

XML Schema and Document Type Definitions (DTD).

Java Architecture for XML Binding (JAXB).

eXtensible Stylesheet Language (XSL) and XSL Transformations (XSLT).

XPath and XQuery.

SOAP.

Many of these areas will be covered in the rest of this paper.

## III Parsing

### III.1 Initial Document Creation

The aim of parsing a document is to read its contents and make it available to a program. To start this process it is necessary to have the XML document in the first place. Once a system is up and running the creation of such documents is normally straight forward, but during the development phase of a system the documents are often hand-crafted. The approach taken on IceCube is that we start of with a complete prose description of the contents we need in a document. We then create an XML document which contains this data and is structured into elements which map onto data objects discussed in the description. Figure 1 is a simplified example of such a document. This document shows the readout of a two channel ATWD (Analog Transient Waveform Digitizer).

### III.2 Document Object Model

There are two standard approaches to parsing an XML document. The first is the Document Object Model (DOM), which is a W3C standard. In this approach the whole document is parsed all at once and a representation, using the “Composite” gamma-ref  design pattern, is created in memory. It is then up to the program to navigate around this memory structure to access the information that it needs. While this approach is very simple from the developers point of view, custom representation of the data require a second round of processing, accessing the original memory structure, to transform a document’s content in to anything other than a composite tree. Moreover for large documents this approach can end up being a resource hog as the whole document appears in memory.

### III.3 Simple API for XML

The alternate approach to parsing is the Simple API for XML (SAX). In this approach each element is handed to a call-back routine for processing. This requires that the developer of these handlers to maintain the document context as it is processed. On the other hand this allows for the processing of content as it is read and so the full document may not be held in memory after parsing has completed.

### III.4 Tools

There are a number of different tools available that can handle the parsing duties for a program. In many cases these tools have handle either approach to parsing as well as validating (see below) a document as it is read. On IceCube parsing is currently handled by the Xerces-J tool xerces-site  from Apache.

## IV Validation

An XML document is valid when its contents conform to the rules laid down in its XML Schema or DTD. As the DTD technology is being superseded in most cases by the XML Schema alternative it will not be considered any further in this paper.

Validation should not be confused with a “well formed” document. This simply means that all “begin”” elements in the document have matching “end” elements in the right place in the element hierarchy and that attributes are correctly delimited. A “well formed” document say nothing about the contents of any elements.

An XML schema describes the allowed contents of a document. It lays down the rules about how and which elements can be contained inside another element, as well as defining what attributes an element can have.

Apart from defining the legal structure of a document, an XML Schema can also specify which values of element and attributes are considered legal. Figure 2 is the XML schema against which the file shown in Figure 1 should be validated. The enumeration elements with the bitsPerSample attribute are an example of such a restriction. Full details about XML Schema can be found elsewhere xsd-primer .

All this means that it is possible to use XML validation to guarantee that an input document’s content is correct before this content is passed on to a program for processing. Thus a developer is not required to write their own validation routines.

## V Java Data Structures

As an XML Schema can define data types, Sun has developed the Java Architecture for XML Binding (JAXB) interface which enables Java classes to be created by running an XML Schema through a processor.

As well as simple class creation the reference implementation of JAXB from Sun, which is part of the Java Web Services Developer Pack Java-WSDP , provides default implementations of these classes that can easily be read from and written to XML documents. Figure 3 is an example of how easy it is to read and access the contents of the example file shown in Figure 1.

## VI Transforms

One of the powers of XML is that as the structure of a document is independent of its content it is possible to transform the contents of one document onto a different document that has an alternate structure. In IceCube we made use of this feature when we designed the “Simple Test Framework” (STF).

### VI.1 Simple Test Framework

The STF is designed to run a set of modules that test the functionality of our Digital Optical Module (DOM) hardware. These modules take as input a set of parameters and provide as a result a set out output parameters. This whole interface can be described by an XML document. Moreover the input parameters for a module, when it is to be executed, can be specified by means of a setup XML file and the output parameters returned in the form of a results XML file. Figure 4 is an example of a description file, while Figure 5 and 6 are examples out setup and results files respectively.

The description file not only associates a type with both input and output parameters, but can also specify a legal ranges and default value for input parameters. This sort of information can be used for at least two different purposes.

* The creation of a C header file that correctly declares the test signatures, any appropriate limits, and default values.
* An XSL Schema file that can be used to validate the input and output files.

In IceCube we have XSL files to handle both of these tasks. Both of the above tasks can be accomplished by using an “eXtensible Stylesheet Language” (XSL) file and a XSL transform engine, e.g. Xalan from Apache xalan-site .

### VI.2 XSLT

XSL files are a description of how the contents of an existing XML file should be transformed by an XSLT to create a new file. This new file is not required to be XML. It can be a simple text file or a file of “format objects” which can be used to generate PDF files.

In IceCube we use XSL files to generate both text files (C header files) and XML files (XML Schema files). Figure 7 contains the C header file that contains the signatures for the ExampleOne module. This file was generated by the XSL in Figure 8.

A similar XSL file can be used to create the XML Schema for setup and result XML files for any given module. This XML Schema can then be used to verify that the setup file is valid before it is passed to the STF to be executed, thus relieving the STF of the need to have custom code to check its input values.

It should also be noted that the XML description has an associated XML Schema, stfDefs.xsd (not show here), against which it should validate. By requiring any description document to validate against that schema we can then be confident that both the resultant XML schema file and C header files will be correct.

## VII Searching

### VII.1 XPath

The initial form of searching developed for XML documents is XPath. The basics from of an XPath statement is to specify a hierarchy of elements, with wildcarding allow. Using this system it is possible to select one or more elements in a document. This selection can be based on the type of elements, an elements context or even the value of its attributes or contents.

XPath is used in XSL documents to identify what transforms should be applied to which elements. The match and select attributes in Figure 8 are examples of this mechanism.

### VII.2 XQuery

Currently development is underway on an extension of the XPath mechanism which will span more than one document. This is called XQuery and is planned to form the basis of a query language for XML databases.

The format of a query in XQuery follows the same lines as XPath and the resulting elements are collected into a single XML document which acts as the results of the XQuery.

XML databases are still in their infancy, though there are some examples in existence xindice-site . At present IceCube is not using one so no examples of XQuery are currently available.

## VIII Data Transfers

The simple Object Access Protocol (SOAP) soap-primer  has become one of the main engines that has driven the adoption of XML for data transfers. This system wraps up an XML document in a “envelope” with has a “header” and a “body” such that the document can be moved from one process to another. The other process may even be on a different machine. Non-XML documents can be handled by “SOAP with Attachments” soap-attachments .

At this time IceCube has not had the opportunity to pursue the use of SOAP within its software system, but there are plans to use it as the main system for handling inter-process communications within the DAQ and production farms. This will allow components of those sub-systems to be written in either Java or C++. At present those systems are being developed in Java and RMI is being used for the inter-process communications.

## IX Conclusions

IceCube has been very fortunate in that its circumstances have made it a good platform in which to assess the promises of the XML technology. XML is quickly becoming an industry standard for both data and configuration files. This has mean that there are plenty of tools available to handle XML. Many of these tools are open source and thus are suitable for use on IceCube. A lot of work has gone into these tools and their adoption by IceCube means that we can concentrate on our own issues rather than the development of tools that are not directly related to our purpose.

Our conclusion is that XML is a good approach for text based files and small data files. It is easy to understand and customize. There is a good array of tools available that handle many of the tasks needs to read and write such files and therefore its adoption can reduce development time and, because of the tools wide adoption, improve softwares reliability.


text:
# References

June 12, 2003

Theory of strongly phase fluctuating d-wave superconductors and the spin response in underdoped cuprates

Igor F. Herbut

Physics Department, Simon Fraser University, Burnaby, British Columbia, Canada V5A 1S6

Abstract: General theory of d-wave quasiparticles coupled to phase fluctuations of the superconducting order parameter is discussed. In the charge sector the superfluid density is found to conform to the Uemura scaling. The spin susceptibility exhibits four distinct regimes with increasing frequency, and scales with the superconducting $`T_c`$, as observed.

Keywords: vortex fluctuations, underdoped cuprates, magnetic response

Fax: 1-604-291-3592, Email: iherbut@sfu.ca.

The intimate relationship between antiferromagnetism and d-wave superconductivity (dSC) has been one of the central themes in the physics of cuprates. Recently a new connection between the two has been put forward , , by which the antiferromagnetism is the result of the dynamical breaking of an approximate ”chiral” symmetry of the standard dSC. The relevant interaction between the low-energy quasiparticles is provided by the quantum vortex fluctuations, which are assumed to increase with underdoping.

The low energy theory of d-wave quasiparticles and the phase fluctuations of the superconducting order parameter at $`T=0`$ and in two dimensions may be written as  $`S=d^3x`$, $`=_\psi +_{gauge}+_\varphi `$ , where

$$_\psi =\underset{i=1}{\overset{2}{}}\overline{\mathrm{\Psi }}_i\gamma _\mu (_\mu ia_\mu )\mathrm{\Psi }_i,$$
(1)
$`_{gauge}={\displaystyle \frac{1}{\pi }}𝐚(\times 𝐀_{})+{\displaystyle \frac{1}{\pi }}𝐯(\times 𝐀_+)+`$ (2)
$`2K(𝐯+𝐀)^2+i𝐉_c(𝐯+𝐀),`$
$$_\varphi =\underset{n=1}{\overset{2}{}}[|(i𝐀_n)\mathrm{\Phi }_n|^2+\mu ^2|\mathrm{\Phi }_n|^2+\frac{b}{2}|\mathrm{\Phi }_n|^4].$$
(3)
$`\mathrm{\Psi }_i`$ describe the neutral spin-1/2 quasiparticles (spinons) near the four nodes at $`\pm 𝐊_{1,2}`$ , complex $`\mathrm{\Phi }_{1,2}`$ create the fluctuating vortex loops associated with spin up and down, and $`𝐚`$, $`𝐯`$  and $`𝐀_{1,2}=𝐀_+\pm 𝐀_{}`$ are the auxiliary gauge fields that facilitate the spinon-vortex coupling. $`𝐉_c`$ is the charge current, and $`𝐀`$ an electromagnetic gauge potential. The two quasiparticle velocities have been set to $`v_F=v_\mathrm{\Delta }=1`$ in (1). $`K`$ is the bare superfluid density.

When $`\mathrm{\Phi }_1=\mathrm{\Phi }_20`$, the above theory reduces to the $`QED_3`$ , , with the insulating spin density wave (SDW) as the ground state . Here I focus on the superconducting state, in which $`\mathrm{\Phi }_n=0`$. Integrating out $`\stackrel{}{v}`$ and the virtual vortex fluctuations yields

$$_\psi +\frac{6m}{\pi }𝐚^2+\frac{\rho _{sf}}{2}𝐀^2+iZ𝐉_c𝐀+O(𝐉_c^2),$$
(4)
with $`\rho _{sf}=48mK/(12m+4\pi K)`$. $`Z=12m/(12m+4\pi K)`$ is the charge renormalization factor. $`m^2=\mu ^2+O(b)`$ is the uniform vortex susceptibility, assumed to be proportional to doping. Since on dimensional grounds $`mT_c`$, for $`mK`$, $`\rho _{sf}T_c`$, in qualitative agreement with the Uemura scaling. Furthermore, spin may be considered separated from the charge in Eq. (4). Upon integration over the massive $`\stackrel{}{a}`$, the spin sector becomes the 2+1 dimensional Thirring model. In the RPA approximation the spin susceptibility becomes

$$\chi ^{\prime \prime }(2𝐊_i\pm 𝐪,\omega )=(\frac{12m}{\pi })^2\frac{N\mathrm{\Theta }(\omega ^2q^2)\sqrt{\omega ^2q^2}}{(96m/\pi )^2+\omega ^2q^2}.$$
(5)
This describes the soft spin mode in the dSC, the condensation of which on the insulating side yields the SDW , .

Eq. (5) implies that for $`\omega T_c`$, the maximum values of $`\chi ^{\prime \prime }`$ are located at four diagonally incommensurate wave vectors $`\pm 2𝐊_{1,2}`$. Due to their low intensity ($`\omega `$) and narrowness, these ‘mother’ peaks should be very hard to observe, and indeed have not yet been detected. As the frequency increases, however, both the intensity and the width of the peaks grow. Restoring $`v_Fv_\mathrm{\Delta }`$ in (5), and shifting the peaks to the vicinity of $`(\pi ,\pi )`$, one finds the ‘mother’ peaks overlapping first at four parallel incommensurate positions, independent of frequency. With the further increase in frequency, the initial peaks at some point begin also to overlap at $`(\pi ,\pi )`$, and for a while the commensurate response dominates. Assuming $`2𝐊_1=0.465(2\pi ,2\pi )`$ , and $`v_F=1.2`$ $`eVA`$ the energy of the ‘resonance’ may be estimated to $`60`$ $`meV`$, not far from $`40`$ $`meV`$ seen in slightly underdoped cuprates. Finally, at even larger $`\omega `$ the maximum in Eq. (5) shifts to a $`q0`$, which implies a weak redistribution of the commensurate peak to four parallel incommensurate positions, but this time with upward dispersion. The evolution of the spin response  is in agreement with the observations in underdoped YBCO .

This work has been supported by NSERC of Canada and the Research Corporation.