text:
# A spatial network explanation for a hierarchy of urban power laws

## Abstract

In this Letter we report an application of spatial complex networks to simulate the evolution of the spatial economic system. The model depicts the spatial trade network with land areas as nodes and trade relations between activities as connections. The network grows by addition of new pairs of activities using simple market economic assumptions. The definition of nodes and connections allows us to compare model and reality on the level of market land values. The model provides an explanation to several empirically observed statistical phenomena in spatial economics. Furthermore, market demand is an important driving force for most other spatially extended systems of human origin. Thus, a model that simulates the geographic distribution of economic activity may shed light on a range of other systems that have been investigated separately using complex networks.

Standard models of urban economics have difficulties describing economic growth and as a consequence they cannot reproduce important empirical observations, maybe most notably Zipf’s LawZipf\_1949 ; Ioannides\_Overman\_2002 ; Fujita\_Krugman\_Venables\_1999 . One of the basic tenets of market economics is that economic growth arises as a result of mutually beneficial trade between specialized agents. We use this mechanism as the basis of an abstract network representation of economic activity and trade on a surface.

Casting the urban economy as a complex network puts a range of new tools at our disposal for analyzing this important systemAlbert\_Barabasi\_2002 ; Dorogovtsev\_Mendes\_2002 ; Newman\_2003 . The benefit goes beyond direct observations about land uses and land values since many other growing spatial systems of human origin grow as a direct response to demand from within the economy. Because of this, a model that simulates the geographic distribution of economic activity can also shed light on the workings of many other social, technical and economical systems – many of which has been studies as complex networks in their own right such as the Internet, phone call graphs and human sexual contacts Albert\_Jeong\_Barabasi\_1999 ; Pastor-Satorras\_Vazquez\_Vespignani\_2001 ; Aiello\_Chung\_Lu\_2000 ; Liljeros\_Edling\_Amaral\_Stanley\_Aberg\_2001 ; Albert\_Barabasi\_2002 ; Dorogovtsev\_Mendes\_2002 ; Newman\_2002a . It is quite conceivable that reported power law distributions in these and other such systems may stem more or less directly from near-optimal distribution of services to meet demand within the geographic distribution of the urban system.

We take fixed-size non-overlapping land areas to be network nodes, and the trade relations between them to be connections, in the network model. The activity within a node is defined as the number of connections connecting it to the rest of the network. We hereby obtain an abstract representation that can capture heterogeneity in producer specialization and spatiality at the same time. This allows us to introduce a basic economic unit, the Economically Correlated Pair (ECP), as an abstract expected future stream of trade between two activities. Per definition, an ECP has unit economic value to both its end-points – a reasonable simplification given the very large number of trade relations in an economy.

The simulations we use for the results in this Letter start from a lattice of undeveloped land areas and proceeds to simulate the allocation of new urban activities and their trade relations. We compare statistics of simulated node degrees with empirically observed land values. The comparison assumes a linear relationship between the amount of money generated and the value of land. This relationship follows from our definition of ECPs and from i) market pricing of goods and services and ii) the connection between trade benefits and land value.

i) Market pricing of commodities provides an adaptive measure that allows us to compare the activities that generate them. Hence, on average, an edge contributes identically to the value of both nodes to which it connects. This contribution is exactly our definition of activity, which implies that the degree of a node is proportional to its benefits due to trade.

ii) This connection consists of two proportionalities. For a node $`i`$ we have

$$v_ir_ix_i,$$
(1)
where $`v_i`$ is the value of the corresponding land area, $`r_i`$ is the bid-rentAlonso\_1964 , and $`x_i`$ is the total trade benefits as outlined above. Capitalizing periodic rent income from the site $`i`$ gives land value $`v_i=\frac{r_i}{i}`$ where $`i`$ is the interest rateOSullivan\_2002 . The second proportionality is a weak form of the leftover principle from urban economics, which states that, in a competitive land market, rent equals the amount of money left after all expenses (except rent) are paid. This amount of money equals the sum of all trade benefits at the site. For our results it is sufficient that, on average, a certain proportion of each new unit of trade benefit goes to the landowner.

We formulate the network model by dividing a geographic area into $`N`$ non-overlapping fixed-size cells $`\{1,2,\mathrm{},N\}`$ and taking these to be network nodes. The undirected connections between nodes are trade streams (ECPs) between activities in the cells, these can also connect a cell to itself and two nodes can be connected by any number of mutual ECPs. The amount of profit (before rent), generated in a node $`i`$, is per the definition of ECPs equal to $`mx_i`$, where $`m`$ is the value of an ECP and $`x_i`$ is the degree.

Network evolution is the process by which a network changes over time. New activity occurs in the form of additions of ECPs to the system. The growth mechanisms include balanced node addition and preferential attachment which gives a strong similarity to the Barabási-Albert model, and a scale-free node distribution can be expectedBarabasi\_Albert\_1999 .

The network is initialized by the pairwise connection of $`n_0`$ (even) spatially uniformly distributed nodes. Addition of ECPs are separated into primary and secondary growth, corresponding to the placing of first and second edge end-points in the network. The placement of an end-point is done either preferentially (per-activity) or uniformly (per-node); corresponding to *in situ* expansion or establishment of new activity, respectively. Expansion of node activity can correspond to existing activities becoming more space efficient as well as them being out-bid and replaced by more efficient activities.

Preferential growth corresponds to trade with an existing activity. Activities are considered to be average activities and thus they all interact to the same extent with each other. If we assume constant fractions $`q_1`$ and $`q_2`$ for uniform and preferential growth respectively, with $`q_1+q_2=1`$, the probability of selecting a node $`i`$ preferentially in primary growth is

$$\mathrm{\Pi }_i^{1,pref}=q_2\frac{x_i}{_jx_j}.$$
(2)
Uniform growth corresponds to average activities selecting sites in a competitive land market. The probability of selecting a node $`i`$ uniformly as a primary effect is

$$\mathrm{\Pi }_i^{1,uni}=q_1\frac{\delta _i^{(D)}+b\delta _i^{(P)}+bϵ\frac{n_t^{(P)}}{n_t^{(E)}}\delta _i^{(E)}}{n_t^{(D)}+b(1+ϵ)n_t^{(P)}},$$
(3)
where $`n_t^{(D)}`$ is the number of developed nodes at time $`t`$ and $`\delta _j^{(D)}=1`$ if node $`j`$ is developed and $`\delta _j^{(D)}=0`$ otherwise. The meanings of $`n_t^{(P)}`$, $`n_t^{(E)}`$, $`\delta _j^{(P)}`$ and $`\delta _j^{(E)}`$ are analogous to $`n_t^{(D)}`$ and $`\delta _j^{(D)}`$, with $`P`$ referring to perimeter nodes and $`E`$ referring to external nodes. Perimeter nodes are all undeveloped nodes that are adjacent to a developed node and external nodes are undeveloped nodes without any developed neighbors. The parameter $`b`$ controls the fraction of perimeter nodes that are open for development and $`ϵ`$ controls how many external nodes that are open for development compared to the number of perimeter nodes. The interpretation of this is that external growth is identical to perimeter growth but takes place on ambient infrastructure not explicitly represented in the model. Infrastructure tends to grow as the urban system grows which makes it reasonable to assume that $`ϵ`$ is constant over time.

The probability of secondary preferential growth at site $`i`$ as a consequence of primary growth at site $`j`$ is

$$\mathrm{\Pi }_i^{2,pref}=q_2\frac{D_{ij}x_i}{_kD_{kj}x_k},$$
(4)
and for secondary uniform growth, it is

$$\mathrm{\Pi }_i^{2,uni}=q_1\frac{D_{ij}\left(\delta _i^{(D)}+b\delta _i^{(P)}+bϵ\frac{n_t^{(P)}}{n_t^{(E)}}\delta _i^{(E)}\right)}{_kD_{kj}\left(\delta _k^{(D)}+b\delta _k^{(P)}+bϵ\frac{n_t^{(P)}}{n_t^{(E)}}\delta _k^{(E)}\right)}.$$
(5)
$`D_{ij}`$ is a matrix of site-to-site interaction strengths, which decay with transportation costs. We have used $`D_{ij}=(1+cd(i,j))^\alpha `$, where $`d(i,j)`$ is the Euclidean distance between sites $`i`$ and $`j`$. The non-negative parameters $`c`$ and $`\alpha `$ controls the impact of spatiality.

To simplify the analysis of the model it is useful to assume that development of new land (addition of an active node) takes place at a constant rate $`q_A`$ (compared to other types of growth, not as a function of physical time). To motivate the assumption, let us consider the growth of developed clusters and their perimeters. It is true in simulations of the model, and it can be empirically verified (see Figures 1c and 1d), that the cluster area distribution is close to a simple power-law with density function $`f(A)A^\beta `$, and that the relation between cluster perimeter $`P`$ and cluster area $`A`$ has the form $`PA^\lambda `$, with $`\lambda <1`$. From this we observe that for the entire system of clusters we have

$$\frac{n_t^{(P)}}{n_t^{(D)}}\frac{_1^{\mathrm{}}A^\beta A^\lambda 𝑑A}{_1^{\mathrm{}}A^\beta A𝑑A}=\frac{\beta 2}{\beta \lambda 1},$$
(6)
assuming that $`\beta >2`$ and $`\lambda <1`$. We define $`q_1^{}`$ as the fraction of activity increments that occur uniformly on developed nodes, and Eq. (3) gives

$$q_1^{}=\underset{i}{}\delta _i^{(D)}\mathrm{\Pi }_i^{1,uni}=q_1\left(1+b(1+ϵ)\frac{n_t^{(P)}}{n_t^{(D)}}\right)^1,$$
(7)
which, because of Eq. (6), is approximately constant. This means that the rate of node activation, $`q_A=1q_1^{}q_2`$, also can be considered constant.

If we assume that the spatial bias between sites is smallAndersson\_Hellervik\_Lindgren\_2003 , the time evolution of expected activity on a developed site $`i`$ follows the equation

$$x_i(t+1)=x_i(t)+2q_1^{}\frac{1}{n_t^{(D)}}+2q_2\frac{x_i(t)}{_jx_j(t)},$$
(8)
which is solved by the continuous-time method introduced by Barabási et al Albert\_Barabasi\_Jeong\_1999 . After sufficiently long time the degree distribution approaches the form $`P[x_i=x](x+B)^\gamma `$, with

$$B=\frac{q_1^{}}{q_2(1q_1^{}q_2)}$$
(9)
and

$$\gamma =1+\frac{1}{q_2}.$$
(10)
Most of the parameters in the model can be readily estimated from empirical data. If $`n_t^{(C)}`$ is the number of clusters, we have $`ϵ=n_t^{(C)}/n_t^{(D)}`$. The exponent $`\gamma `$ in the distribution of activity gives information about the size of $`q_2`$ (and $`q_1`$) via Eq. (10). The perimeter parameter $`b`$ can be determined by $`q_A=\frac{n_t^{(D)}}{_ix_i(t)}`$ and $`b=q_A\left((q_1q_A)(1+ϵ)\frac{n_t^{(P)}}{n_t^{(D)}}\right)^1`$.

In principle, these parameters are not empirical in nature, they are aggregated measures of fundamental properties of the current economic system (and region) of interest. The spatial parameters $`c`$ and $`\alpha `$ reflect the statistics of transport characteristics of the economic configuration. They are not as easily estimated from data as the other parameters but results seem robust to their exact values and the functional form of $`D_{ij}`$.

The data used for all empirical results is based on a database delivered by Sweden Statistics that covers estimations of the market value of all land in Sweden (2.9 million data points). Cluster measurements was obtained by identifying clusters using land value as the threshold parameter (75 kSEK/Ha) and then comparing accumulated price and population within the obtained areas. To identify clusters in simulated and empirical data we use a computer program that masks away all data points below a threshold value and then treats all contiguous (8-cell neighborhood) areas as clusters.

Once the frequency distribution of land value per unit area is correct, the challenge in realizing the cluster level is to faithfully capture the spatial sorting of the node degrees. Using the model presented in this Letter we obtain excellent agreement between simulated and empirical statistics for a range of non-trivial higher-order structures in successive orders of upward causation (cells to clusters), see Figures (1b) through (1b). These results are not trivial consequences of the land price distribution; it is perfectly possible to arrange the developed cells into any system of clusters. The same is true for the relation between cluster area and perimeter. Also, in Figure (1f) we demonstrate that geographical distribution of population and land value can be observed interchangeably. This is a likely example of how an earlier studied spatially growing system of human origin likely is just a response to the growth of the economic network. The growth of the latter is defined by activity inter-dependency and transportation costs – a more fundamental candidate for a generating process than more vague postulations about why and when people re-locateMarsili\_Zhang\_1998 .

All figures reflect data from the same Monte Carlo simulation run. Figure (1e) shows a part of the spatial configuration from this run. The parameters for the reported simulation, $`m=75`$kSEK ($`10`$SEK$`1`$USD), $`b=0.15`$ and $`ϵ=0.25`$ have been estimated from empirical data. By studying the exponent in the empirical distribution of land values, an initial estimate of $`q_1=0.1`$ was obtained. For the reported results $`q_1`$ was adjusted to $`0.2`$. The spatial parameters were $`c=0.1`$ and $`\alpha =2`$. Investigation of sensitivity shows that exponents and proportions change slowly and smoothly with all parameters. A square grid of $`2000\times 2000`$ cells was used and the number of iterations was $`10^6`$.

Macroscopic models, such as Simon’s Model and later derivations, are based on the rather unrealistic assumption that individual cities would react in unison to external eventsSimon\_1955 ; Gabaix\_1999 ; Schweitzer\_Steinbrink\_1998 . Although this may in some sense be trueIoannides\_Overman\_2002 , such models hold limited explanatory and predictive power since they can not be mapped to a realistic situation where the internal structure of cities is important or when clusters interfere and coagulate. By making a similar assumption for small fixes-size land areas our representation remains valid regardless of the dynamics of higher-order structures such as cities.

The relation between pattern (probability distributions etc.) and process (microdynamics) is not one-to-one. Rather, for each observed macroscopic pattern there is likely a whole class of processes capable of generating it. Apart from arguing for the network model ontology on the micro level there is a number of implications that remain to test against empirical data. Among these are the time evolution of land values, explicit measurement of the trade network and making studies similar to the one presented here for other countries.

We thank Anders Hagson and Jonas Tornberg, City and Mobility, Chalmers University of Technology, for providing preprocessed empirical data. We also thank Martin Nilsson and Steen Rasmussen for valuable discussions and input.


format: text

text:
# 1. Introduction

## 1. Introduction

The nature of the vibrational excitations of silicate glasses and melts represents a challenging problem of condensed matter physics, glass and earth sciences. Significant experimental and theoretical efforts have been made in order to understand the vibrational properties of these materials, because they provide valuable information about the underlying microscopic structure and other anomalous physical properties of silicate glasses . However, the lack of long-range translational order considerably complicates this task and requires new approaches for measurement, analysis and simulations of the vibrational spectra of glasses and melts as well as a careful investigation on how the vibrational spectra depend on the details of the considered glass (chemical composition, thermal history, etc.).

Therefore we compare in the present paper the vibrational density of states (VDOS) and several vibrational characteristics of an ab initio model of sodium tetrasilicate glass Na2Si4O9 (denoted hereafter NS4) calculated by two different methods. The first calculation is based on an ab initio treatment of the interatomic forces and is performed in the framework of the density functional theory, while the second one is performed in the approximation of a classical valence force fields (VFF) potential which is especially suited for describing the dynamics of partially covalent materials like silicate glasses.

The structure of the NS4 glass has been intensively investigated by different experimental methods , reverse Monte Carlo (RMC)  as well as classical molecular dynamics (MD) simulations . Recently, the structural and electronic properties of the NS4 glass have been modeled by a combined Car-Parrinello and classical MD simulation . The attention given to this system is justified by the fact that it can be used as a prototype for more complicated aluminosilicate and hydrous silicate glasses.

The aims of the present study are: (i) to test the quality of the ab initio NS4 vibrational dynamics by comparing the calculated and experimental polarized Raman spectra; (ii) to test the quality and the transferability of a previously proposed  valence force fields potential for sodium silicate glasses by comparison with ab initio calculations; and (iii) to analyse the effects of adding Na2O to a-SiO2 on the character of the vibrational characteristics of the glass. We would like to note that the present ab initio calculations using the Car-Parrinello (CP) method involve various approximations (pseudopotentials, truncated plane-wave basis, etc.) Nevertheless we may expect that the CP approach leads to a better estimation of the real interatomic forces by comparison to the VFF approach. On the other hand, the computation time needed to perform a VFF type calculation is generally very small compared to the one used in first principles calculations.

## 2. Calculation procedures

### 2.1. The glass models

The amorphous NS4 model used in this study contains 90 atoms (24 silicon atoms, 54 oxygen atoms and 12 sodium atoms) confined in a cubic box of edge length 10.81 Å corresponding to the experimental NS4 mass density of 2.38 g/cm3 . The glass model has been generated using combined classical and ab initio molecular dynamics (MD) simulations. For this we have firstly performed the liquid equilibration, the quench, and the initial part of the low temperature relaxation of the system by using a classical force field. Subsequently we have refined the obtained structure within the framework of the ab initio MD. We note that this approach was successfully employed for the study of the structural, electronic and vibrational properties of vitreous SiO2 , as well as for a structural and electronic study of vitreous NS4 .

The classical MD simulations were performed with an interatomic potential which is a modification of the potential proposed by van Beest et al. , derived in order to study sodium silicates . The NS4 model was generated by quenching a well equilibrated liquid at 3500 K to 300 K, with a quench rate equal to $`510^{13}`$ K/s. The glass model obtained in this way was annealed at 300 K for 70 ps. The final atomic coordinates and velocities after the classical relaxation, were used as initial coordinates and velocities for a short ($`0.5`$ ps) ab initio MD simulation. The ab initio simulations were performed in the framework of the Car-Parrinello (CP) method  using the CPMD software .

In the first principles calculations, the electronic structure calculations were treated via the Kohn-Sham formulation  of density functional theory, within the generalized gradient approximation employing the B-LYP functional . The valence Kohn-Sham orbitals were expanded in a plane-wave basis set defined by an energy cutoff of 70 Ry, at the $`\mathrm{\Gamma }`$-point of the supercell. The core-valence interactions were described by norm-conserving pseudopotentials of the Bachelet-Hamann-Schlüter type for Si atoms  and of the Trouiller-Martins type for O atoms . Further details are given in Ref. . The CP dynamics was performed using a fictitious electronic mass of 800 a.u. and a time step of 0.097 fs. In the present CP calculations, we used for the sodium atoms a semi-core Trouiller-Martins norm-conserving pseudopotential  instead of the Goedecker semi-core type used in Ref. . This change allowed a better estimation of the internal pressure of the model. Nevertheless, the same structural modifications as reported in  occurred immediately after the CP dynamics was switched on, i.e. the shortening of the Si-NBO bond and the lengthening of the Si-BO and Na-NBO bonds (in the remainder of the text, ’BO’ denotes the bridging oxygens, while ’NBO’ denotes the non-bridging oxygens). The CP run leads to structural relaxation without changes in the network topology, the ring and Q-species distributions being determined in the classical stage.

In order to produce the first NS4 glass model (called the ab initio NS4 model) used further for the vibrational calculations, we relaxed to $`0`$ K the glass model obtained at the end of the CP simulation described above. The second glass model was produced by performing a relaxation of the ab initio NS4 model with a harmonic valence force field (VFF) potential of the Kirkwood-type (see Ref.  for further details as well as the next subsection where we give the values of the VFF potential parameters used in this work). At the end of the VFF relaxation, the NS4 model had practically the same average Si-O and Na-NBO bond lengths as the starting ab initio model. However we found a small decrease ($`3\%`$) of the Na-BO bond length as well as a decrease of the standard deviation of the average Si-O bond length distribution, from 0.033 Å to 0.019 Å. In addition we observed a decrease of the standard deviation of the intratetrahedral O-Si-O angle distribution (from 3.7 degrees to 2.2 degrees).

### 2.2. VDOS and Raman spectra calculations

Within the first principles approach, the dynamical matrix was computed numerically by evaluating the second derivatives of the total energy with respect to atomic displacements ($`5.310^3`$ Å for each atom), which were calculated using the finite differences of the atomic forces. (We note that, for these latter calculations as well as for the ab initio relaxation to $`0`$ K mentioned in the previous subsection, we used an energy cutoff of 90 Ry for the plane-wave basis expansions of the Kohn-Sham orbitals.) The direct diagonalization of the dynamical matrix provides the eigenvalues $`\{\omega _p^{\mathrm{CP}}\}`$ and the corresponding normalized eigenvectors $`\{𝐞^{\mathrm{CP}}(\omega _p^{\mathrm{CP}})\}`$, $`p=1,\mathrm{},3N`$, where $`N=90`$ is the number of atoms in the model. Each eigenvector is a $`3N`$-dimensional vector which components are proportional to the displacements of the atoms in mode $`p`$.

Within the VFF approach, the dynamical matrix was calculated analytically and was diagonalized using the Householder method. The detailed expressions of the dynamical matrix elements are given elsewhere . The obtained vibrational frequencies will be denoted $`\{\omega _p^{\mathrm{VFF}}\}`$ and the corresponding eigenvectors $`\{𝐞^{\mathrm{VFF}}(\omega _p^{\mathrm{VFF}})\}`$, $`p=1,\mathrm{},3N`$. In the present work we used the following values for the stretching force constants $`\alpha _{ij}`$ : $`\alpha _{\mathrm{Si}\mathrm{BO}}=465\mathrm{N}/\mathrm{m}`$, $`\alpha _{\mathrm{Si}\mathrm{NBO}}=655\mathrm{N}/\mathrm{m}`$, $`\alpha _{\mathrm{Na}\mathrm{BO}}=25\mathrm{N}/\mathrm{m}`$, and $`\alpha _{\mathrm{Na}\mathrm{NBO}}=30\mathrm{N}/\mathrm{m}`$, while the bending force constants $`\beta _{ijk}`$ were $`\beta _{\mathrm{BO}\mathrm{Si}\mathrm{BO}}=\beta _{\mathrm{BO}\mathrm{Si}\mathrm{NBO}}=\beta _{\mathrm{NBO}\mathrm{Si}\mathrm{NBO}}=35\mathrm{N}/\mathrm{m}`$ and $`\beta _{\mathrm{Si}\mathrm{BO}\mathrm{Si}}=14\mathrm{N}/\mathrm{m}`$. We note that we have used larger Na-O stretching and slightly larger Si-BO-Si bending force constants than the values reported in Ref.  since these values led to a better agreement between the calculated and experimental heat capacities for sodium silicate glasses .

The knowledge of all the eigenmodes allows us to calculate the reduced polarized Raman spectra of the NS4 glass model using the bond polarizability approximation (BPA)  presented in details in Refs. . The BPA parameters used in the present work are : $`A_{\mathrm{Si}\mathrm{BO}}^{}=0.5\mathrm{\AA }^3`$, $`\gamma _{\mathrm{Si}\mathrm{BO}}=0.05\mathrm{\AA }^3`$, $`A_{\mathrm{Si}\mathrm{NBO}}^{}=1\mathrm{\AA }^3`$, $`\gamma _{\mathrm{Si}\mathrm{NBO}}=0.05\mathrm{\AA }^3`$, $`A_{\mathrm{Na}\mathrm{O}}^{}=0.05\mathrm{\AA }^3`$, $`\gamma _{\mathrm{Na}\mathrm{O}}=0\mathrm{\AA }^3`$, where $`A_{\alpha \beta }^{}`$ is the derivative of the parallel bond polarizability parameter with respect to the length of the bond between atoms $`\alpha `$ and $`\beta `$, while $`\gamma _{\alpha \beta }`$ is the perpendicular bond polarizability parameter for the $`\alpha \beta `$ bond ($`\alpha ,\beta `$=Si,Na,BO,NBO,O). Recently, a comparison of the Raman scattering mechanism in $`\alpha `$-quartz using both ab initio and BPA methods has shown that the BPA reproduces the ab initio Raman intensities within $`15\%`$ . For amorphous silica models generated by ab initio molecular dynamics simulations, the use of the BPA gives rise to Raman spectra in good agreement with experimental data . This gives strong support for the application of the BPA for the calculation of the Raman spectra of large models of silicate glasses. The comparison of the Raman spectra calculated independently from the CP and the VFF eigenmodes will allow us to investigate separately the effects of the BPA parameters and the eigenmodes on the accuracy of the Raman spectra calculations.

## 3. Results

In Fig. 1 we compare the CP and VFF vibrational densities of states. Both VDOS are normalized to one, and we have used the same uniform Gaussian broadening of full width at half maximum $`2\sigma =40\mathrm{cm}^1`$. In both spectra we can distinguish three ranges: a low-frequency range ($`0500\mathrm{cm}^1`$), an intermediate frequency range ($`500900\mathrm{cm}^1`$), and a high-frequency range ($`9001200\mathrm{cm}^1`$ ). The low-frequency bands arise generally from Na-O stretching and bending as well as from Si-O bending and rocking motions. The strong mid-frequency band near $`800\mathrm{cm}^1`$ arises from complex motion with large displacements of the Si atoms against the BO atoms . The high-frequency modes correspond mainly to Si-BO and Si-NBO stretching motions .

In spite of the fact that the intensities of the CP low-frequency bands are slightly higher than the VFF ones, both approaches yield almost identical band positions and shapes. Larger discrepancies between the two VDOS are found in the mid- and the high-frequency ranges since the VFF mid-frequency bands are shifted to higher frequencies compared to the CP bands. Nevertheless the two approaches yield practically the same intensities (same number of vibrational modes). In the high-frequency range the shape of the vibrational bands is generally the same but the VFF bands are shifted to higher frequencies, have larger intensities and are more narrow. The latter effect may be attributed to the decrease of the standard deviation of the Si-O bond length distribution.

Unfortunately, to the best of our knowledge, there are no VDOS data for the NS4 glass measured by inelastic neutron scattering to compare with our calculations. Therefore we compare in Fig. 2 the calculated and experimentally reduced polarized (VV and VH) Raman spectra. The polarization is VV if the electric fields of the incident and the scattered light are parallel and is VH when they are perpendicular. We recall that the so-called reduced spectrum is obtained by multiplying the experimental measured spectrum with the correction factor $`\omega /[n(\omega )(\omega \omega _0)^4]`$ in order to have a temperature independent quantity ($`n(\omega )`$ is the Bose factor). The polarized Raman spectra were measured on doubly-polished glass plates in backscattering geometry using the $`514.5\mathrm{nm}`$ line of an Ar+ laser and a XY triple spectrometer equipped with a liquid-nitrogen cooled CCD detector and confocal entrance optics with $`300\mathrm{s}`$ accumulation time per spectral window. The CCD spectra were binned with a step of $`4\mathrm{cm}^1`$. For the calculations of the theoretical spectra, the final spectra have been obtained by applying a uniform Gaussian broadening with full width at half maximum equal to $`30\mathrm{cm}^1`$ to the calculated Raman intensities. It should be noted that generally it is not possible to determine experimentally absolute Raman intensities. Therefore the VV spectra (Fig. 2a) were normalized so that the strongest peaks (at about $`1100\mathrm{cm}^1`$, and characterized by a small depolarization ratio $`I_{VH}/I_{VV}0.07`$) have the same maximum unit intensity, while the VH spectra (Fig. 2b), were normalized so that the strongest depolarized peaks at $`800\mathrm{cm}^1`$ have the same intensity (with a depolarization ratio $`I_{VH}/I_{VV}0.35`$).

In the VV Raman spectra (Fig. 2a) we distinguish also three bands. The previous calculation  of the VV Raman spectra on NS4 glass models have assigned the low-frequency band at about $`520\mathrm{cm}^1`$ to symmetric bending of the Si-BO-Si linkages, the mid-frequency band at about $`780\mathrm{cm}^1`$ to the Si motions against the BO atoms and the high-frequency band at about $`1100\mathrm{cm}^1`$ with a shoulder at about $`920\mathrm{cm}^1`$ to the Si-NBO stretching motions.

For the position and the relative intensity of the low-frequency band, a good overall agreement between the experimental and the CP calculated curve is observed although, e.g., the width of the CP band is smaller and the shoulder in the experimental data at about $`600\mathrm{cm}^1`$ is not reproduced. The agreement between the low-frequency VFF band and the experimental spectrum is less good - the calculated VFF intensity is stronger and the maximum peak position is shifted to slightly higher frequencies. In other words, the VFF potential overestimates the amplitudes of the Si-BO-Si symmetric bending vibrations.

The intensity of the CP calculated mid-frequency band is in good agreement with the experiment but its position is shifted approximately $`80\mathrm{cm}^1`$ to lower frequencies. The shape of the VFF corresponding band is similar to the experimental one, but we note a smaller intensity and a slight shift to lower frequencies (see the inset in Fig. 2a).

Concerning the high-frequency band, we note that the position of the main peak is shifted approximately $`20\mathrm{cm}^1`$ downwards in the CP case and approximately $`20\mathrm{cm}^1`$ upwards in the VFF case indicating slightly weaker and slightly stronger Si-NBO interactions, respectively. Both approaches yield much stronger intensity for the $`920\mathrm{cm}^1`$ band which arises mainly from Si-NBO stretching vibrations in SiO4 tetrahedra with 2 NBO per Si ($`Q^2`$ species) .

In the depolarized (VH) spectra, the calculated low-frequency bands seem to be correctly positioned (see Fig. 2b, where we have lopped off the top of the VH spectra in order to see the details of the spectra in the low- and mid-frequency ranges). Regarding the relative intensities, we note a good agreement for the CP band. In the mid-frequency range, the CP band is broadened and shifted to slightly lower frequencies, while the VFF band is shifted to slightly higher frequencies. In a previous calculation of the VH Raman spectra using a NS4 model generated by RMC simulations , the agreement was better in this frequency range. Since the mid-frequency band arises mainly from Si vibrations in $`Q^4`$species , the differences between the present VFF calculation and the results presented in Ref.  are to be attributed to the different degree of polymerization of the corresponding models .

In the high-frequency range (see the inset in Fig.2b) the calculated depolarized Raman intensities are much weaker in both calculated spectra which indicates that the corresponding Si-O bending BPA parameters are too small.

## 4. Discussion

In order to understand the differences in the VDOS and Raman spectra computed within the CP and VFF approaches, we have analysed first the corresponding partial VDOS. The partial VDOS curves $`g_\alpha (\omega )`$ in Fig. 3 ($`\alpha =`$ Si, BO, NBO, Na) have been defined as :

$`g_\alpha (\omega _p)=g(\omega _p){\displaystyle \underset{i\alpha }{}}|𝐞_i(\omega _p)|^2`$ (1)
where $`𝐞_i(\omega _p)`$, $`i=1,\mathrm{},N`$, are 3-component real space eigenvectors.

For the low-frequency bands, we have an overall good agreement between the CP and VFF partial VDOS. The band maxima positions for the Si (Fig. 3a), BO (Fig. 3b) and NBO (Fig. 3c) atoms are very close ($`\pm 30\mathrm{cm}^1`$). In this frequency range the NBO atoms participate mainly in O-Si-O bending vibrations . The biggest intensity differences are observed in the BO (Fig. 3b) and Na (Fig. 3d) partial VDOS, which can explain the differences in the total VDOS in this frequency range.

As can be seen in Fig. 3d, both the CP and VFF calculations predict a dominant Na contribution only in the low-frequency range, but the VFF peak is shifted ($`40\mathrm{cm}^1`$) as a whole to lower frequencies. From a previous VFF analysis  of the vibrational characteristics of two NS4 structural models obtained by classical MD and RMC simulations, it followed that the Na motion dominated below $`100\mathrm{cm}^1`$. In the present VFF calculation, the Na contribution reaches a maximum at $`170\mathrm{cm}^1`$. Thus the Na-O stretching force constants used in the present work describe better the Na-O motion than the force constants reported in Ref. .

At intermediate frequencies, the CP and the VFF partial VDOS for the Si, BO and NBO atoms have qualitatively the same shape but the ones for Si and BO are shifted to lower frequencies in the CP calculation (Fig. 3a and b). These shifts are similar to the observed shift of the mid-frequency band in the VV calculated using the CP eigenmodes. So it seems that the CP approach does not describe properly the vibrational dynamics of the Si and BO atoms in the mid-frequency range. Interestingly, the Si motions are enhanced in the CP calculation while the BO motions are enhanced in the VFF calculation.

In contrast to this, the partial VDOS show that in the high-frequency range all of the VFF partial VDOS bands are stronger and more narrow (especially for the BO and Si atoms) than the CP ones. This feature may be related to the stronger localization of the CP modes in this frequency range, as it will be discussed below.

As a first measure for the localization of the modes we consider the participation ratio $`p_c`$ which, for a given eigenmode $`p`$, is defined as follows  :

$`p_c(\omega _p)={\displaystyle \frac{(_{i=1}^N|𝐮_i(\omega _p)|^2)^2}{N_{i=1}^N|𝐮_i(\omega _p)|^4}},`$ (2)
where $`𝐮_i(\omega _p)=𝐞_i(\omega _p)/\sqrt{M_i}`$ is the atomic displacement of atom $`i`$ in mode $`p`$ and $`M_i`$ is the mass of atom $`i`$. The frequency dependences of the CP and VFF participation ratios are shown in Figs. 4a and 4c. As can be seen from Fig. 4a and 4c, the CP modes are generally more localized than the VFF ones, especially around 200 and above $`900\mathrm{cm}^1`$. This is an interesting result which merits to be studied in more details. The localization is also enhanced near the edges of the three main VDOS bands for both the CP and the VFF participation ratios. These are the so-called band tails containing localized modes, which have been already observed for a larger NS4 glass models in Ref. , as well as for silica glass models . In the CP case, these tails are more pronounced, in particular at the top of the lower band and on both sides of the high-frequency bands. The analysis of the frequency dependences of the vibrational correlation lengths $`L_c(\omega _p)`$ which are plotted in Fig. 4b and 4d, confirms the localization of the high-frequency modes. The vibrational correlation length is defined as:

$`L_c(\omega _p)=\sqrt{{\displaystyle \frac{_{i=1}^N|𝐫_i𝐫_p|^2|𝐮_i(\omega _p)|^2}{_{i=1}^N|𝐮_i(\omega _p)|^2}}}`$ (3)
where $`𝐫_i`$ is the position of atom $`i`$, $`𝐫_p`$ is the center of mode $`p`$ displacements, defined as $`{\displaystyle \underset{i=1}{\overset{N}{}}}𝐫_i|𝐮_i(\omega _p)|^2`$. In this sense, the vibrational correlation length  is the standard deviation of $`𝐫_p`$, and gives a measure of the spatial localization of the vibrational modes by indicating that the amplitude of the atomic vibrations decreases significantly beyond this correlation length.

## 5. Conclusions

In this work, we have investigated some aspects of the vibrational dynamics of a model of sodium tetrasilicate glass, constructed by combined classical and Car-Parrinello molecular dynamics simulations. Its vibrational dynamics has been then studied using two different approaches - a first principles one and a classical valence force fields one. The valence force field serves as a natural approach to analyze the character of the vibrational modes in terms of simple motions like stretching, bending and rocking. Within the bond polarization approximation, we have calculated the polarized Raman spectra and compared with experimental Raman data.

The VDOS calculated by the two methods have very similar shape, especially for the low-frequency bands below $`500\mathrm{cm}^1`$. In the high-frequency range, the VFF partial VDOS for the Si, BO and NBO atoms are shifted to higher frequencies and they are narrower. This may indicate that the potential VFF parameters have to be adjusted in order to better reproduce the CP vibrational characteristics of our NS4 model (an increasing number of low frequency modes, the position as well as the localization of the high frequency modes).

The comparison of the calculated and experimental VV Raman spectra shows that a good agreement between the relative intensities of the low- and high-frequency bands is achieved using the CP eigenmodes. This indicates that the corresponding BPA parameters are approximately correct and can be used for calculation of the VV Raman spectra of structural models of silicate glasses that are still too large to be fully treated ab initio. However, despite the fact that the CP eigenvectors are calculated from first principles they do not reproduce well the position of the mid-frequency VV and VH Raman bands at around $`780\mathrm{cm}^1`$. Additional improvement of the BPA parameters describing the high-frequency depolarized Raman bands is also necessary.

The analysis of the character of the vibrational modes shows that the localization of the CP vibrational modes is stronger than the one of the VFF modes, in particular on the high $`\omega `$ side of the low-frequency band ($`200\mathrm{cm}^1`$) and on both sides of the high-frequency bands.

## Acknowledgments

We thank Magali Benoit for very interesting and stimulating discussions. Part of this was supported by the European Community’s Human Potential Program under contract HPRN-CT-2002-00307, DYGLAGEMEM. The financial support from the German Science Foundation under project SFB 408 for N. Zotov is kindly acknowledged. The ab initio simulations have been performed on the IBM/SP3 at CINES, Montpellier, FRANCE.

FIGURE CAPTIONS

Figure 1 : Vibrational densities of states of the ab initio NS4 glass model obtained within the Car-Parrinello ab initio approach (solid line) and within a classical valence forcefield approximation (dashed line).

Figure 2 : Frequency dependence of the VV (a) and the VH (b) Raman intensities calculated within the bond-polarizability approximation using the CP eigenmodes (dotted lines) and the classical VFF eigenmodes (dashed lines). The solid lines are the experimental spectra.

Figure 3 : Partial VDOS for Si (a), BO (b), NBO (c), and Na (d). The solid and dashed lines correspond to the CP and VFF calculations, respectively.

Figure 4 : Frequency dependences of the participation $`p_c`$ and of the vibrational correlation length $`L_c`$ in the CP case (panels a) and b) ) and in the VFF case (panels c) and d)). The horizontal lines at a height equal to 5.405 Å (i.e. half of our simulation box length) are drawn as a guide for the eye in order to estimate the number of the normal modes with correlation length smaller or bigger than this value.


text:
# The PROOF Distributed Parallel Analysis Framework based on ROOT

## I INTRODUCTION

The Parallel ROOT Facility, PROOF, is an extension of the well known ROOT system (http://root.cern.ch/) ROOT96  that allows the easy and transparent analysis of large sets of ROOT files in parallel on remote computer clusters. The main design goals for the PROOF system are transparency, scalability and adaptability. With transparency we mean that there should be as little difference as possible between a local ROOT based analysis session and a remote parallel PROOF session, both being interactive and giving the same results. With scalability we mean that the basic architecture should not put any implicit limitations on the number of computers that can be used in parallel. And with adaptability we mean that the system should be able to adapt itself to variations in the remote environment (changing load on the cluster nodes, network interruptions, etc.).

Being an extension of the ROOT system, PROOF is designed to work on objects in ROOT data stores. These objects can be individually keyed objects as well as TTree based object collections. By logically grouping many ROOT files into a single object very large data sets can be created. In a local cluster environment these data files can be distributed over the disks of the cluster nodes or made available via a NAS or SAN solution.

In the near future, by employing Grid technologies, we plan to extend PROOF from single clusters to virtual global clusters. In such an environment the processing may take longer (not interactive), but the user will still be presented with a single result, like the processing was done locally.

The PROOF development is a joint effort between CERN and MIT.

## II THE PROOF SYSTEM

### II.1 System Architecture

PROOF consists of a 3-tier architecture, the ROOT client session, the PROOF master server and the PROOF slave servers. The user connects from his ROOT session to a master server on a remote cluster and the master server in turn creates slave servers on all the nodes in the cluster. Queries are processed in parallel by all the slave servers. Using a pull protocol the slave servers ask the master for work packets, which allows the master to distribute customized packets for each slave server. Slower slaves get smaller work packets than faster ones and faster ones process more packets. In this scheme the parallel processing performance is a function of the duration of each small job, packet, and the networking bandwidth and latency. Since the bandwidth and latency of a networked cluster are fixed the main tunable parameter in this scheme is the packet size. If the packet size is chosen too small the parallelism will suffer due to the communication overhead caused by the many packets sent over the network between the master and the slave servers. If the packet size is too large the effect of the difference in performance of each node is not evened out sufficiently. This allows the PROOF system to adapt itself to the performance and load on each individual cluster node and to optimize the job execution time.

### II.2 TSelector Framework

The TSelector framework plays a central role in the PROOF system. It allows users to write their analysis code in a way that enables them to process small data samples locally on a workstation and using the same code to process large data sets in parallel on a cluster or the grid using PROOF.

To use the framework a user derives a class from TSelector and implements the member functions that are part of the protocol. The framework also specifies how input objects are made available to the selector object and how output objects are returned to the user or client session.

The TTree::MakeSelector() function generates a skeleton version of a TSelector derived class that should be used as the basis for the analysis code. The user can either edit the generated files directly or create a class that inherits from the generated class. The second option makes sense if the definition of the tree is expected to change. It allows for the class to be easily regenerated without modifying the user class. Table 1 lists the most important functions and members which we will describe in some more detail.

The input objects are made available to the selector in the list pointed to by fInput. The list is available throughout the life of the selector. The objects in the list are owned by the system. Each of the slaves has a full copy of the list.

The Begin() function is called before any object or event is processed. It is a good place to create things like histograms and initialize data that is needed to process each event. The Begin() function is called once in each slave.

The Init() function is called each time a new file is opened by the system. The generated code takes care of setting up the branches in case of a tree analysis. The user can extend this routine for example to read a calibration object that pertains to the objects or tree in that file.

The Process() function is called for each object or event in the tree. When processing a tree it is up to the user to decided which branches should be read and when. Generally it is most efficient to read branches as they are used. E.g. read the branch(es) that are used to make the selection and only if the event is to be processed read the rest of the branches that are needed. When processing keyed objects the pointer to the object is stored in the fObject member.

The Terminate() function is called after all objects or events have been processed. It allows the user to do a cleanup, like deleting temporary objects and closing auxiliary files. It is also the place to do final calculations. The Terminate() function is called once in each slave.

Finally the contents of the fOutput lists of all the slaves are merged in the master and send back to the client. The merging uses the special Merge(Tlist\*) API. The objects are fist grouped by name. Using the introspection features provided by CINT CINT , the C++ interpreter used by ROOT, it is then determined if the Merge() function is implemented for those objects. For histograms this is available by default. User classes can implement the API if merging makes sense. If the merge function is not available the individual objects are returned to the client.

### II.3 Specifying Data Sets

The TDSet class is used to specify the collection of objects or trees that are to be processed with PROOF. When analyzing trees the user must pass "TTree" to the TDSet constructor. The TDSet::Add() function is then used to add a file, directory within that file and tree name to the set.

To process a collection of objects, the user must pass the name of the class of the objects to the constructor. In that case TDSet::Add() is used to add files and the directories within them that contain the objects.

The contents of one TDSet can also be added to another TDSet. The TDSet implements the Print() function which allows the contents to be inspected. Users can create TDSets ”by hand” but it is also foreseen that (Grid) catalogs will provide query interfaces that return TDSets. It will be possible to use logical filenames rather then physical filenames to specify files. The translation will be done using the API defined by the abstract class TGrid in combination with plugins implementing the communication with available catalogs.

### II.4 The PROOF Package Manager

In complex analysis environments the analysis scripts very likely depend on one or more external libraries containing common algorithms. To be able to run these scripts successfully on PROOF it is required that these external libraries are available on each slave node.

The PROOF package manager has been designed to distribute and install these libraries on a PROOF cluster. Packages are compressed tar files containing the library source or binaries with one additional PROOF-INF directory (like Java JAR files). This PROOF-INF directory contains a SETUP.C script and, optionally, a BUILD.sh shell script. The BUILD.sh script is used to compile the library source on the PROOF cluster, which might have a different type of CPU or OS than the user’s local system. The SETUP.C script is used to load the libraries into the PROOF server processes. Package files are called PAR (PROOF ARchive) files and must have a .par extension.

For example an Event.par file which provides a library called libEvent.so has the following BUILD.sh:

```
  make libEvent.so

```

and SETUP.C:

```
  Int_t SETUP()
  {
    gSystem->Load("libEvent");
    return 1;
  }

```

In a two step process packages are first uploaded to PROOF and then enabled as required. The PROOF package API (see table 2) gives the user full control. To avoid unnecessary transfer of par files the upload command first asks the MD5 checksum of the remote version of the package and, if it exists, checks it with the one of the local version. Only when the checksums are different will the par file be transfered to the PROOF cluster.

## III A REAL LIFE EXAMPLE

In this section we will use a simple but full fledged example to illustrate the use of PROOF and to describe some of the features of the implementation.

### III.1 Local Development

The data model underlying the tree in this example is typical for a HEP experiment (see fig 2). The tree contains a TClonesArray for a hit, track and vertex class and a single event class. Relations between hits, tracks and vertexes are expressed using TRef and TRefArray.

For this example we use a very simple example script. The relevant parts are included below. All error checking has been omitted for brevity. In the Init() function a single histogram is created which is stored in the data member fVtx\_x of the selector class.

```
  void antsel::Begin(TTree *tree)
  {
    // Initialize the tree branches.
    Init(tree);

    fVtx_x = new TH1F("v","v",100,-10.,10.);
  }

```

The Process() function reads the data, makes a simple selection and fills the histogram. It uses the fRMSSelVtx reference to get the proper vertex. In a more elaborate version it would only read the event branch for the test and read the vertex branch as needed.

```
  void antsel::Process(Int_t entry)
  {
    fChain->GetTree()->GetEntry(entry);

    if (eventInfo->fPdlMean > 1500) {
      TPhAnTVertex *vtx = (TPhAnTVertex*)
        eventInfo->fRMSSelVtx->GetObject();
      fVtx_x->Fill(vtx->fPos.X());
    }
  }

```

The Terminate() function adds the histogram to the output list fOutput such that it will be returned to the client by PROOF.

```
  void antsel::Terminate()
  {
    fOutput->Add(fVtx_x);
  }

```

Running the script on a typical laptop we can analyze 2000 events, or 8 Mbyte, in about 6 seconds (ROOT I/O has compressed the data in the tree by a factor 5.) Once the script works as required we are ready to run it using PROOF.

### III.2 Session Setup

We start by creating a PROOF session. The argument specifies the PROOF cluster we want to connect to.

```
  root[1] gROOT->Proof("pgate.lns.mit.edu")

```

We first have to authenticate ourselves to the PROOF server. All authentication methods implemented in the standard ROOT TAuthenticate class are available. The PROOF master server is then created. The master server reads the proof.conf configuration file and uses the information to start a number of slave servers. The configuration file is provided by the administrator of the PROOF cluster. It contains information about the nodes that make up the PROOF cluster, the relative performance of the nodes and which nodes share a file system.

When all the slaves are started we are ready to run queries. But before that we can configure packages. For this example we’ll upload and enable a single package

```
  root[2] gProof->UploadPackage("ant.par")
  root[3] gProof->EnablePackage("ant")

```

Normally we only need to upload a package if it has changed. The EnablePackage() function will cause the package to be build in the case of a source package. Finally the SETUP.C script is run on all slaves, loading our library and making available its classes.

We will use a previously loaded script to create the data set

```
  root[4] TDSet *d = get_dataset()

```

Instead of creating the TDSet it might be returned by a database or catalog query.

### III.3 PROOF Run

We are now ready to use PROOF to process our query. Using the same script that was used previously and the data set we just created we issue the command

```
  root[5] d->Process("antsel.C", "", 60000)

```

The system goes through a number of steps to implement this command. First the selector script and possibly its corresponding include file are send to the master and from the master to the slaves. The system will optimize this step if an up to date version of the script is already available in the cluster. It also makes use of the shared file system(s) if available.

Then the client sends the query message, including the input objects, option and parameters, to the master. The master determines the total number of events in the dataset by having the slaves open the files in the dataset in parallel. This information could in the future also be obtained from a database, avoiding this step in he process. The master server creates a list of nodes with the files to be processed and their sizes. This list will be used to optimize the distribution of work over the slaves. The files have to be specified as rootd URLs for this optimization to be available, e.g. "root://proof.mit.edu/data/file.root".

The master now forwards the query to the slaves. Each slave starts the script and enters a loop in which it asks the master for work packets. The master will allocate work to each slave using the information about the location of the files. A slave will first be assigned files which are on its local disk. When a slave has exhausted all local files the rootd protocol is used to process files on other nodes. If no location information was available, .e.g if the files are stored on a central NFS server, the algorithm reduces to a round robin assignment of the files to the slaves. An heuristic is used to determine the packet size. It takes into account the number of slaves and their relative performance as well as the total number of events to be processed. The master monitors the real-time performance of each slave allowing this heuristic to be refined. The master also records which packet is processed by which slave allowing error recovery to be implemented in the case of slave failure.

When all events are processed the slaves send the partial results to the master. Using the previously described Merge() algorithm the master combines the partial results and sends them to the client.

### III.4 The Results

The processing of the above query ran in about 12 seconds using 8 slaves on 4 dual Athlon 1.4 GHz machines. The amount of data processed was 240 Mbyte. This shows the large improvement of even a small dedicated cluster over a typical desktop workstation. A more precise measurement of the PROOF performance will be presented in the next section.

## IV SCALABILITY MEASUREMENTS

First performance measurements show a very good and efficient scalability (see fig 4). For the tests we used a Linux cluster of 32 nodes. Each node had two Itanium 2 1 GHz CPU’s, 2x75 GB 15K SCSI disk, 2 GB RAM and fast ethernet. The data set to be analyzed consisted of 128 files totaling 8.8 GB of data (9 million events). Each cluster node had 4 files (277 MB). Processing these 128 files using one node took 325 seconds and 32 nodes in parallel only 12 seconds. One node processing only its 4 local files took 9 seconds. This shows that the efficiency is about 88%. We also ran tests using all 64 CPU’s in the cluster (two slaves per node). This showed the same linearity but less efficiency, due to overhead in the Linux SMP implementation and resource contention in each node. All in all we expect an efficient scalability to at least a 100 node cluster.

## V PROOF IN THE GRID ENVIRONMENT

To be able to build a global virtual PROOF cluster we need to use the Grid services that are currently being developed. The interfacing of PROOF to the Grid can be done at several levels. The following levels have been identified:

* Interface to the grid file catalog allowing a user to select a data set based on tags or logical file names (using wildcards etc).
* Interface to the grid resource broker to find the best location(s), based on the data set, where to run the query. This could trigger the replication of some missing files to a cluster (only when the amount of data in the files is relatively small).
* Interface to the grid job queue manager to start proof master and slave daemons on the remote cluster. The ROOT client will then connect to these pre-started daemons to create a PROOF session. This will require the grid queuing system to support interactive high priority jobs.

We are currently working with the AliEn ALIEN  Grid developers on a prototype that implements step wise the above scenarios. AliEn (http://alien.cern.ch) is a grid solution tuned for typical HEP data processing. It has been developed by the ALICE collaboration but is experiment independent. It provides all commonly understood Grid services, like: file catalog, replication service, job queue manager, resource broker, authentication service, monitoring, etc. It is entirely implemented in Perl and has a C client API. AliEn is simple to install and works reliably. In ALICE it is routinely used for simulation and reconstruction. AliEn has also been chosen as the grid component for the EU MammoGrid project.

The ROOT and PROOF interface to AliEn, and other Grid middleware, will be via a TGrid abstract interface.

## VI FUTURE WORK

Currently we are fine tuning PROOF for a single cluster and are working with some ”early adopter sites” that have setup dedicated clusters.

Several smaller and larger additions and refinements of the system are already foreseen. General infrastructure for dynamic startup of the session will allow PROOF to co-exist with traditional batch system as well as Grid based environments. This could be extended to dynamic allocation and release of slaves during the session, e.g. based on the data set to be processed.

We are also working at combining multiple clusters at geographically separated sites which will require a hierarchy of master servers. Further integration of PROOF and ROOT will drive the implementation of event lists and friend trees. We are also working on TTree::Draw() style functionality.

The Grid middle-ware is in wild development at the moment and we are following closely these developments to make sure it will support the features we need, especially facilities for interactive high priority jobs and data services.

## VII CONCLUSIONS

We demonstrated the PROOF system, part of the widely adopted ROOT analysis environment, that allows users to harness the power of a large cluster of workstations from their desktops. PROOF greatly extends the amount and range of data that can be interactively analyzed.

###### Acknowledgements.

This material is based upon work supported by the National Science Foundation under Grant No. 0219063.


text:
# Vector meson production and inclusive 𝐾⁰_𝑠𝐾⁰_𝑠 final state at HERA

## 1 Vector Meson Production at HERA

Production of light and heavy vector mesons (V) in exclusive and in proton dissociation processes has been studied in $`ep`$ reactions at HERA in a wide range of the $`\gamma p`$ centre-of-mass energy $`W`$ and photon virtuality $`Q^2`$. For very low $`Q^2`$ and $`W>10`$ GeV, these reactions display features characteristic of a soft diffractive process. The cross sections rises weakly with the energy $`W`$, $`\sigma _V(W)W^\delta `$ ($`\delta 0.22`$), and has a steep exponential $`t`$ dependence, where $`t`$ is the squared four-momentum transfer at the proton vertex. Such processes are well described within the framework of Regge phenomenology  and the Vector-Meson Dominance model (VDM) , where the photon is assumed to fluctuate into a vector meson before scattering from the proton via Pomeron exchange. However, this approach fails at high values of $`Q^2`$ or high vector meson masses. In the hard regime, models based on perturbative QCD (pQCD)  can be used to describe the vector meson production. The photon fluctuates into a $`q\overline{q}`$ state and the quark dipole interacts with the proton in the lowest order via two gluons exchange . The exchange of the gluon ladder has also been calculated in the leading logarithm approximation (LLA) . The cross section is related to the rise of the gluon density in the proton as $`x`$ decreases, where $`x`$ is the Bjorken scale variable, and has a strong W dependence, $`\sigma _V(W)W^\delta `$ ($`\delta 0.8`$.)

In this contribution, perturbative QCD models are compared with results from vector meson production in photoproduction and deep inelastic scattering at HERA.

### 1.1 Exclusive Vector Meson Production

Figure 1 shows $`\sigma _{tot}^{\gamma ^{}pJ/\psi p}`$ from ZEUS  and H1  as a function of W for fixed $`Q^2`$, and also in photoproduction. The data is fitted with $`W^\delta `$. The ZEUS measurements are shown in table 1. The results are compatible with results in photoproduction from H1  $`\delta =0.83\pm 0.07(statsyst)`$ and ZEUS  $`\delta =0.69\pm 0.02(stat)\pm 0.03(syst)`$. The cross sections are also compared to pQCD predictions from Frankfurt, Koepf and Strikman (FKS)  using CTEQ4M PDF , and from Martin, Ryskin and Teubner (MRT)  using CTEQ5M PDF . The results are sensitive to the PDF used and have large theoretical uncertainties, but they qualitatively describe the strong W-dependence of the cross section. They also indicate that at high vector meson masses $`M_V`$, a hard regime already exist at very low values of $`Q^2`$, and that $`M_V^2`$ may set a hard scale.

The $`Q^2`$-dependence of the $`\delta `$ parameter is shown in Fig. 2 for the light vector meson $`\rho ^0`$ from ZEUS  and H1 , and for $`J/\psi `$ in photoproduction from H1 . At very low values of $`Q^2`$, Regge phenomenology and VDM models describe the $`\rho ^0`$ meson production. However, the energy W-dependence becomes steeper at high $`Q^2`$, with a smooth transition from soft to hard regime. On the other hand, the $`J/\psi `$ cross section has a strong energy dependence already at $`Q^2=0`$. These results indicate that $`Q^2`$ may also set a hard scale.

### 1.2 Vector Meson Production with Proton Dissociation

It was predicted  that hard processes may occur at high $`t`$. Typical elastic vector meson production has $`t<1`$ GeV2, and cannot be used to verify these predictions. On the other hand, proton-dissociation processes have higher values of $`t`$. Figure 3 shows the differential cross sections of $`\rho `$, $`\varphi `$ and $`J/\psi `$ vector mesons in the photoproduction proton-dissociative reaction $`\gamma pVY`$ in the range $`80<W<120`$ GeV and $`x>0.01`$ from ZEUS , where $`Y`$ is the dissociated hadronic system. The data is compared to predictions from Forshaw and Poludniowski  using LLA BFKL . The results are in good agreement, indicating that pQCD can describe vector meson production at high $`t`$ regime.

### 1.3 Summary

The results presented in this contribution indicates a smooth transition from soft to hard regime. They also show that pQCD describes the vector production in the hard regime, and that $`Q^2`$, $`M_V^2`$ and $`t`$ can set a hard scale.

## 2 Inclusive $`K_s^0K_s^0`$ Final State in Deep Inelastic Scattering at HERA

The $`K_s^0K_s^0`$ system is expected to couple to scalar and tensor glueballs. This has motivated intense experimental and theoretical study during the past few years . Lattice QCD calculations  predict the existence of a scalar glueball with a mass of $`1730\pm 100`$ MeV and a tensor glueball at $`2400\pm 120`$ MeV. The scalar glueball can mix with $`q\overline{q}`$ states with $`I=0`$ from the scalar meson nonet, leading to three $`J^{PC}=0^{++}`$ states whereas only two can fit into the nonet. Experimentally, four states with $`J^{PC}=0^{++}`$ and $`I=0`$ have been established : $`f_0(980)`$, $`f_0(1370)`$, $`f_0(1500)`$ and $`f_0(1710)`$.

The state most frequently considered to be a glueball candidate is $`f_0(1710)`$  , but its gluon content has not yet been established. This state was first observed in radiative $`J/\psi `$ decays  and its angular momentum $`J=0`$ was established by the WA102 experiment using a partial-wave analysis in the $`K^+K^{}`$ and $`K_s^0K_s^0`$ final states . Observation of $`f_0(1710)`$ in $`\gamma \gamma `$ collisions may indicate a large quark content. A recent publication from L3  reports the observation of two states in $`\gamma \gamma `$ collisions above $`1500`$ MeV, the well-established $`f_2^{}(1525)`$  and a broad resonance at 1760 MeV. It is not clear if the latter state is the $`f_0(1710)`$.

The $`ep`$ collisions at HERA provide an opportunity to study resonance production in a new environment. Production of $`K_s^0`$ particles has been studied previously at HERA . In this contribution, the first observation of resonances in the $`K_s^0K_s^0`$ final state in inclusive deep inelastic $`ep`$ scattering (DIS) is reported .

### 2.1 Results

The data used for this study correspond to a total integrated luminosity of 120 pb-1 collected in ZEUS during the 1996-2000 running period.

Oppositely charged track pairs reconstructed by the ZEUS central tracking detector (CTD) and assigned to a secondary vertex were selected and combined to form $`K_s^0`$ candidates. Both tracks were assigned the mass of a charged pion and the invariant-mass $`M(\pi ^+\pi ^{})`$ was calculated. Only events with at least one pair of $`K_s^0`$ candidates were selected. The invariant mass of the $`K_s^0`$ pair candidate $`M(K_s^0,K_s^0)`$ was reconstructed in the range $`0.995<M(K_s^0K_s^0)<2.795`$ GeV.

Figure 4 shows the distribution in $`x`$ and $`Q^2`$ of selected events containing at least one pair of $`K_s^0`$ candidates. The virtual photon-proton center of mass energy is in the range $`50<W<250`$ GeV.

Figure 5 shows the $`M(\pi ^+\pi ^{})`$ distribution in the range $`0.45<M(\pi ^+\pi ^{})<0.55`$ GeV after the $`K_s^0`$ pair candidate selection.

A cut $`cos\theta _{K_s^0K_s^0}<0.92`$ is applied to the selected $`K_s^0`$ pair candidates, where $`\theta _{K_s^0K_s^0}`$ is the opening angle between the two $`K_s^0`$ candidates in the laboratory frame. This cut removes phase-space effects and the presence of $`f_0(980)`$/$`a_0(980)`$ at the threshold, which decays to collinear $`K_s^0`$ pairs affecting the measurement in the 1300 MeV mass region.

Figure 6 shows the $`K_s^0K_s^0`$ invariant-mass spectrum in the range $`0.995<M_{K_s^0K_s^0}<2.795`$ GeV for data with $`cos\theta _{K_s^0K_s^0}<0.92`$ (filled circles with error bar). Two clear peaks are seen, one around 1500 MeV, consistent with the well-established $`f_2^{}(1525)`$ , and the other around 1700 MeV, close to $`f_0(1710)`$ . There is also an enhancement around 1300 MeV, consistent with the $`f_2(1270)/a_2(1320)`$ interference.

It was found that most of the $`K_s^0`$ pair candidates selected after all cuts are in the “gluon-rich” region $`x_p=2p_B/Q>1`$ of the target region of the Breit frame, which corresponds to the remnant of the proton  , where $`p_B`$ is the absolute momentum of the $`K_s^0K_s^0`$ in the Breit frame.

### 2.2 Summary

The first observation in deep inelastic $`ep`$ scattering of a state near $`1525`$ MeV, consistent with $`f_2^{}(1525)`$, and another close to $`f_0(1710)`$ is reported.

An enhancement which can be interpreted as the $`f_2(1270)`$/$`a_2^0(1320)`$ interference is observed, but its measurement is affected by the presence of the $`f_0(980)`$/$`a_0(980)`$ at the $`K_s^0K_s^0`$ threshold.