text:
# Linking Vibrational Dynamics to Local Electronic Structure: Local Analysis of Dynamics of the relaxed Si87 Cluster

## I Introduction

In recent years, total-energy-vibrational-spectrum studies have been carried out for extended periodic systems as well as complex systems with reduced symmetry because of the development of efficient and reliable methods of calculating the total energyKim95 ; Sti95 ; Sti96 . These studies start with the calculation of the total energy using semi-empirical tight-binding (TB) approachesKim95  or first principles methodsSti95 ; Sti96  such as those based on the density functional theory (DFT). The force constant matrix between a pair of atoms are then calculated as the appropriate second partial derivatives of the total energy with respect to the displacements of the atoms about their equilibrium positions. While it is indeed more satisfying to calculate the vibrational properties using force constant matrices determined from the total energy based on quantum mechanical calculations rather than from empirically-fitted classical potential, a key advantage of the total energy approaches over the approaches based on the classical potentials is that the former allows the possibility of linking the interesting dynamical features, in particular vibrational modes associated with surfaces, interfaces, defects, impuritites, etc., directly to the local electronic structure. However, this issue has not yet been systematically addressed. We present here a scheme that provides a framework through which dynamic features of a vibrating system can be related to the electronic structure. We also demonstrate the application of this scheme using the vibration dynamics of the relaxed Si87 cluster as the working example. We will show how all the interesting dynamical features of the relaxed Si87 clusters can be traced to their electronic origins, including the local electronic structure such as bond charges, dangling bonds, etc.

## II Methods

The dynamics of a system is described by the equation of motion

$$\underset{j\beta }{}H_{i\alpha ,j\beta }c_{j\beta }^\lambda =\omega _\lambda ^2c_{i\alpha }^\lambda $$
(1)
where $`H_{i\alpha ,j\beta }=\varphi _{i\alpha ,j\beta }/\sqrt{m_im_j}`$, $`c_{i\alpha }^\lambda =\sqrt{m_i}u_{i\alpha }^\lambda `$, with $`\varphi _{i\alpha ,j\beta }`$ being the force constant matrix element between atom $`i`$ of mass $`m_i`$ along the direction $`\alpha `$ ($`x,y,`$ or $`z`$) and atom $`j`$ of mass $`m_j`$ along the direction $`\beta `$, and $`u_{i\alpha }^\lambda `$ the eigen-displacement of the $`\lambda `$-mode (with the eigen-frequency $`\omega _\lambda `$) at site $`i`$ along the $`\alpha `$-direction. Traditionally, the vibrational density of states (VDOS) projected on sites (individual atoms) is the tool to analyze the nature and properties of a defect or impurity mode. This is usually accomplished by introducing the sum rule $`_{i\alpha }(c_{i\alpha }^\lambda )^2=1`$ into the defining equation of the VDOS such that

$$\rho (\omega )=\underset{\lambda }{}\delta (\omega \omega _\lambda )=\underset{\lambda i\alpha }{}(c_{i\alpha }^\lambda )^2\delta (\omega \omega _\lambda )$$
(2)
In this way,

$$\rho (\omega )=\underset{i\alpha }{}\rho _{i\alpha }(\omega )$$
(3)
where

$$\rho _{i\alpha }(\omega )=\underset{\lambda }{}(c_{i\alpha }^\lambda )^2\delta (\omega \omega _\lambda )$$
(4)
is the local VDOS at site $`i`$ along the $`\alpha `$ polarization. However, since the force field characterizing a particular vibrational mode is determined by the interaction of the atom at a given site with its neighbors, it is more illuminating if, in addition, the VDOS projected on pairs of atoms is used in combination with the site VDOS to analyze the mode. The interactions of a certain atom with its neighbors are described by the force constant matrices. Thus a sum rule that involves the force constant matrices would be an appropriate vehicle to construct the VDOS projected onto pairs of atoms. Such a scheme can be set up as follows.

Multiplying the equation of motion with $`c_{i\alpha }^\lambda `$ (Eq. (1)), we haveJay97a ; Jay97b

$`(c_{i\alpha }^\lambda )^2`$ $`=`$ $`{\displaystyle \underset{j\beta }{}}{\displaystyle \frac{H_{i\alpha ,j\beta }c_{i\alpha }^\lambda c_{j\beta }^\lambda }{\omega _\lambda ^2}}`$ (5)
$`=`$ $`{\displaystyle \underset{ji,\beta }{}}{\displaystyle \frac{(u_{i\alpha }^\lambda u_{j\beta }^\lambda u_{i\alpha }^\lambda u_{i\beta }^\lambda )\varphi _{i\alpha ,j\beta }}{\omega _\lambda ^2}}`$
In arriving at Eq. (5), the relation $`\varphi _{i\alpha ,i\beta }=_{ji}\varphi _{i\alpha ,j\beta }`$ has been used. Equation (5) can serve two useful functions. First, the summation of Eq. (5) over both $`i`$ and $`\alpha `$ leads to

$$\underset{i\alpha ,ji,\beta }{}\frac{(u_{i\alpha }^\lambda u_{j\beta }^\lambda u_{i\alpha }^\lambda u_{i\beta }^\lambda )\varphi _{i\alpha ,j\beta }}{\omega _\lambda ^2}=\underset{i\alpha }{}(c_{i\alpha }^\lambda )^2=1$$
(6)
This is precisely the sum rule needed to project the VDOS onto pairs of atoms $`i`$ and $`j`$ as Eq. (6) also involves the force constant matrix elements $`\varphi _{i\alpha ,j\beta }`$. Thus, when Eq. (6) is inserted into the defining equation for the VDOS (Eq. (2)), we obtain

$`\rho (\omega )`$ $`=`$ $`{\displaystyle \underset{\lambda ,i\alpha ,ji,\beta }{}}{\displaystyle \frac{(u_{i\alpha }^\lambda u_{j\beta }^\lambda u_{i\alpha }^\lambda u_{i\beta }^\lambda )\varphi _{i\alpha ,j\beta }}{\omega _\lambda ^2}}\delta (\omega \omega _\lambda )`$ (7)
$`=`$ $`{\displaystyle \underset{i\alpha ,ji,\beta }{}}\rho _{i\alpha ,j\beta }(\omega )`$
where $`\rho _{i\alpha ,j\beta }(\omega )`$, the pair VDOS of atom $`i`$ along the $`\alpha `$-direction and atom $`j`$ along the $`\beta `$-direction linked by the force constant matrix element $`\varphi _{i\alpha ,j\beta }`$, is given by

$$\rho _{i\alpha ,j\beta }(\omega )=\underset{\lambda }{}\frac{(u_{i\alpha }^\lambda u_{j\beta }^\lambda u_{i\alpha }^\lambda u_{i\beta }^\lambda )\varphi _{i\alpha ,j\beta }}{\omega _\lambda ^2}\delta (\omega \omega _\lambda )$$
(8)
The second useful function of Eq. (5) is that it provides the framework to link the local VDOS at a particular site to the interactions with its neighbors. In this way, the vibrational mode localized in the vicinity of a given site can be correlated with the vibrational motions of its interacting neighbors. Thus, from Eqs. (4) and (8), we have

$$\rho _{i\alpha }(\omega )=\underset{ji,\beta }{}\rho _{i\alpha ,j\beta }(\omega )$$
(9)
Equation (5) can also be rewritten as

$$(c_{i\alpha }^\lambda )^2=\frac{_{\beta \alpha }\varphi _{i\alpha ,i\beta }u_{i\alpha }^\lambda u_{i\beta }^\lambda +_{ji,\beta }\varphi _{i\alpha ,j\beta }u_{i\alpha }^\lambda u_{j\beta }^\lambda }{\omega _\lambda ^2H_{i\alpha ,i\alpha }}$$
(10)
$`_{i\alpha }(c_{i\alpha }^\lambda )^2=1`$ then leads to

$$\rho (\omega )=\underset{i\alpha ,j\beta }{}\overline{\rho }_{i\alpha ,j\beta }(\omega )$$
(11)
with

$$\overline{\rho }_{i\alpha ,j\beta }(\omega )=\underset{\lambda }{}\frac{\varphi _{i\alpha ,j\beta }u_{i\alpha }^\lambda u_{j\beta }^\lambda \varphi _{i\alpha ,i\alpha }(u_{i\alpha }^\lambda )^2\delta _{ij}\delta _{\alpha \beta }}{\omega _\lambda ^2H_{i\alpha ,i\alpha }}\delta (\omega \omega _\lambda )$$
(12)
Equations (10) and (12) also yield

$$\rho _{i\alpha }(\omega )=\underset{j\beta }{}\overline{\rho }_{i\alpha ,j\beta }(\omega )$$
(13)
Equation (13) allows the decomposition of the site VDOS in terms of on-site as well as off-site contributions.

While the scheme based on Eq. (5) provides the vehicle to decompose the VDOS in terms of contributions from pairs of atoms directly interacting with each other via $`\varphi _{i\alpha ,j\beta }`$, one may also devise a scheme to project the VDOS onto pairs of atoms interacting directly as well as those not interacting directly. Such a scheme is useful for examining collective (extended) as well as more localized vibrational modes. To develop such a scheme, we note that

$$\underset{j}{}c_{j\beta }^\lambda =0$$
(14)
Multiplying Eq. (14) by $`c_{i\alpha }`$ and summing over $`i,\alpha `$ and $`\beta `$ lead to

$$\underset{i\alpha ,j\beta }{}c_{i\alpha }^\lambda c_{j\beta }^\lambda =\underset{i\alpha }{}(c_{i\alpha }^\lambda )^2+\underset{i\alpha ,j\beta }{}\{c_{i\alpha }^\lambda c_{j\beta }^\lambda (c_{i\alpha }^\lambda )^2\delta _{ij}\delta _{\alpha \beta }\}=0$$
(15)
Hence,

$$\underset{i\alpha ,j\beta }{}\{c_{i\alpha }^\lambda c_{j\beta }^\lambda (c_{i\alpha }^\lambda )^2\delta _{ij}\delta _{\alpha \beta }\}=\underset{i\alpha }{}(c_{i\alpha }^\lambda )^2=1$$
(16)
Equation (16) then leads to

$$\rho (\omega )=\underset{i\alpha ,j\beta }{}\widehat{\rho }_{i\alpha ,j\beta }(\omega )$$
(17)
with

$$\widehat{\rho }_{i\alpha ,j\beta }(\omega )=\underset{\lambda }{}\{c_{i\alpha }^\lambda c_{j\beta }^\lambda (c_{i\alpha }^\lambda )^2\delta _{ij}\delta _{\alpha \beta }\}\delta (\omega \omega _\lambda )$$
(18)
According to Eq.(16), we can also write

$$\rho _{i\alpha }(\omega )=\underset{j\beta }{}\widehat{\rho }_{i\alpha ,j\beta }(\omega )$$
(19)
This equation then provides the link of the site VDOS to the pair VDOS $`\widehat{\rho }_{i\alpha ,j\beta }`$. Similarly $`_jc_{j\beta }^\lambda =0`$ leads to

$$\underset{i\alpha ,j}{}c_{i\alpha }^\lambda c_{j\alpha }^\lambda =\underset{i\alpha }{}(c_{i\alpha }^\lambda )^2+\underset{i\alpha ,ji}{}c_{i\alpha }^\lambda c_{j\alpha }^\lambda =0$$
(20)
or

$$\rho (\omega )=\underset{i\alpha ,ji}{}\stackrel{~}{\rho }_{i\alpha ,j\alpha }(\omega )$$
(21)
with

$$\stackrel{~}{\rho }_{i\alpha ,j\alpha }(\omega )=\underset{\lambda }{}c_{i\alpha }^\lambda c_{j\alpha }^\lambda \delta (\omega \omega _\lambda )$$
(22)
Equation (20) also leads to

$$\rho _{i\alpha }(\omega )=\underset{j}{}\stackrel{~}{\rho }_{i\alpha ,j\alpha }(\omega )$$
(23)
Equations (5) through (23) provide some examples of projecting the VDOS onto pairs of atom and relating site VDOS to pair VDOS for the purpose of studying the correlation between the vibrations at different localities for localized as well as more extended vibrational modes.

Within the framework of a total-energy-vibrational-spectrum approach, the bonding configurations of atoms in the system, such as the bond charges, the nature of the bonds, and the number of bonds associated with each atom, are already determinedAlf99 . Therefore, it is straightforward to relate a particular vibrational mode to its electronic origin via either the site VDOS or the pair VDOS. For example, the VDOS can be partitioned according to the contributions from different ranges of the total bond charge $`q_i=(1/2)_j\sigma _{ij}`$ associated with the atom at site $`i`$ in the system where $`\sigma _{ij}`$ is the bond charge between the pair of atoms $`i`$ and $`j`$. This can be accomplished by summing the site VDOS $`\rho _i=_\alpha \rho _{i\alpha }`$ over those sites $`i`$ where $`q_i`$ falls within the range $`[a,b]`$, i.e.,

$$\rho (\omega ;aq_ib)=\underset{aq_ib}{}\rho _i(\omega )$$
(24)
Since the bond charge is closely related to the strength of the bond, such a partitioning can shed light on whether a particular mode is associated with the strengthening or softening of the bonds. Similarly, the VDOS can also be grouped in terms of contributions from atoms possessing a certain number of bonds, namely

$$\rho (\omega ;n_i=k)=\underset{n_i=k}{}\rho _i(\omega )$$
(25)
where $`n_i`$ denotes the number of bonds possessed by the particle $`i`$. This presentation of the vibrational spectrum provides the means to identify the modes associated with dangling bonds and those associated with distorted bonding configurations. To examine how an individual bond affects a particular mode, one may group the pair VDOSs $`\rho _{ij}=_{\alpha \beta }\rho _{i\alpha ,j\beta }`$ according to the range of the bond charge $`\sigma _{ij}`$ such that

$$\rho (\omega ;a\sigma _{ij}b)=\underset{a\sigma _{ij}b}{}\rho _{ij}(\omega )$$
(26)
The characteristics of the bonds is usually described by molecular orbitals such as $`\tau =ss\sigma `$, $`sp\sigma `$, $`pp\sigma `$, or $`pp\pi `$. To determine the nature of the bonds responsible for the vibrational modes, one may classify the contribution to a particular mode by a certain bonding type $`\tau `$ using

$$\rho (\omega ;\tau )=\underset{i,ji}{}\frac{\sigma _{ij}^\tau }{\sigma _{ij}}\rho _{ij}(\omega )$$
(27)
where

$$\sigma _{ij}=\underset{\tau }{}\sigma _{ij}^\tau $$
(28)
The VDOS can also be decomposed in terms of quantities such as the distance, the bond charge, etc. via $`\overline{\rho }_{i\alpha ,j\beta }`$, or $`\widehat{\rho }_{i\alpha ,j\beta }`$ or $`\stackrel{~}{\rho }_{i\alpha ,j\alpha }`$, following similar procedure as described in Eq. (26). Thus, a very flexible framework that allows the convenient decomposition of the VDOS to pair VDOS via $`\rho _{i\alpha ,j\beta }`$, or $`\overline{\rho }_{i\alpha ,j\beta }`$, or $`\widehat{\rho }_{i\alpha ,j\beta }`$, or $`\stackrel{~}{\rho }_{i\alpha ,j\alpha }`$, and the linking of the site VDOS to these pair VDOS (see Eqs. (9), (13), (19), and (23)) can be set up. This framework can then be profitably used to analyze interesting dynamical features of vibrations, in particular how these features are correlated with the local electronic structure, of complex systems with reduced or no symmetry.

## III Results and Discussion

To demonstrate the application of the scheme outlined in Eqs. (5) to (27) to examine the properties of the vibrational modes and to trace these modes to their electronic origins, we use the dynamics of a relaxed Si87 as the working example. The initial configuration of the cluster is an 87-atom tetrahedral network with an atom at the central site. The $`x`$-, $`y`$-, and $`z`$-axis are along the cube edges of the cube-unit cell with the central site at a corner of the cell. The equilibrium configuration of this cluster (see Fig. 1) was first determined by the MD scheme using the non-orthogonal tight-binding Hamiltonian developed by Mennon and SubbaswamyMen97 . The force constant matrices of the relaxed cluster were then numerically calculated from the total energy.

In Fig. 2, the VDOS of the relaxed cluster, together with contributions to the VDOS from those sites with the total bond charge in a certain range (see Eq. (24)), is shown. We use the bonding configuration of the bulk Si as the reference to choose the ranges of the total bond charge because Si-based systems favor the tetrahedral $`sp^3`$ bonding configuration. In bulk Si, each Si atom contributes 0.25$`e`$ to one of the four equivalent bonds, giving rise to a total bond charge of 1$`e`$ for every atom. Thus, we have selected the following four ranges, namely $`q_i>1e`$, $`1e>q_i>0.8e`$, $`0.8e>q_i>0.6e`$, and $`0.6e>q_i`$, to examine the effect of relaxation on the bonding configuration, and hence on the vibrational spectrum of Si87. From Fig. 2, it can be seen that the majority of the contributions to the VDOS originate from those sites with $`1e>q_i>0.8e`$ or $`0.8e>q_i>0.6e`$, indicating that the relaxation has probably caused a general softening of the bonds and thus resulting in an overall shift of the VDOS towards lower frequencies compared to the bulk VDOS. Figure 2 also shows that the mode with the highest frequency is a defect mode originating from the range $`q_i>1e`$, thus a strengthened bonding configuration, while the mode with the lowest frequency originates from the range $`0.6e>q_i`$, a softened bonding configuration or a configuration with dangling bonds. In Fig. 3, the contribution to the VDOS is decomposed according to the number of bonds associated with each atom (see Eq. (25)). The existence of a bond between a pair of atoms is established by the criterion $`\sigma _{ij}>0.04e`$ rather than relying on the distanceAlf99 . Figure. 3 shows that the majority of the contributions still originate from atoms with 4 bonds, with a significant minority contribution from atoms with 3 bonds. The defect mode at the highest frequency can be seen as associated with an atom with 4 bonds. Together with the evidence in Fig. 2, it can be concluded that this mode results from the strengthened bonding configuration. The lowest (non-zero) frequency mode, on the other hand, is seen to originate from an atom with 2 bonds. The evidences provided by Figs. 2 and 3 thus indicate that this mode is a dangling-bond mode, but with remaining bonds strengthened. To further explore this issue, we present in Fig. 4 the decomposition of the VDOS in terms of the ranges of bond charge $`\sigma _{ij}`$ (see Eq. (26)). It can be seen that indeed both the highest- and the lowest-frequency mode are associated with $`\sigma _{ij}>0.5e`$, indicating that the bonds involved in these vibrations are strengthened bonds (compared to the bulk bond). Figure 4 also shows that the majority of the contribution to the VDOS come from the range $`0.4e<\sigma _{ij}<0.5e`$, the range with bond charge not too much deviated from the bulk bonding configuration, thus reinforcing the observation provided by Figs. 2 and 3. Another interesting feature in Fig. 4 is the negative contribution in the low frequency VDOS from the range $`\sigma _{ij}<10^3e`$. The negative contribution is an indication that for contributions to a certain vibrational mode by pairs of atoms $`i`$ and $`j`$ interacting via $`\varphi _{i\alpha ,j\beta }`$, some of those with no chemical bond formed between them may have the force constant too “soft” to support the vibration at that frequency. In Fig. 5, we show the decomposition of the VDOS in terms of bonding types (see Eq. (27)). It can be clearly seen that the main contribution to all the vibrational modes comes from $`pp\sigma `$-bonds, with a substantial minority contribution from $`sp\sigma `$-bonds. The directional characteristics of the $`pp\sigma `$-bonds indicates the dominance of the stretching vibrations as is expected for a Si-based structure.

Figure 6 presents the decomposition of the VDOS according to the ranges of the distance between atoms via $`\widehat{\rho }_{i\alpha ,j\beta }`$ (see Eqs. (18)). Since $`\widehat{\rho }_{i\alpha ,j\beta }`$ only involves $`c_{i\alpha }^\lambda `$ and $`c_{j\beta }^\lambda `$, it provides an useful vehicle to analyze contributions to a certain vibrational mode from the on-site environment, pair of atoms directly interacting (via $`\varphi _{i\alpha ,j\beta }`$) with each other, and pairs of atoms not directly interacting. The three ranges of distances chosen for the decomposition in Fig. 6 correspond roughly to 1st nearest-neighbor distances ($`r_{ij}<3.10\AA `$), 2nd nearest-neighbor distances ($`3.10\AA <r_{ij}<4.15\AA `$), and 3rd nearest-neighbor distances and beyond ($`r_{ij}>4.15\AA `$) respectively. Figure 6 shows that there are contributions to the VDOS from all three ranges to most of the low- to intermediate-frequency modes and the contributions are consistent with the number of pairs of atoms in each range. This is an indication that these modes are extended modes. On the other hand, the dominant contribution to the high-frequency modes mainly originates from the range $`r_{ij}<3.10\AA `$. Considering the fact that for any given atom in the cluster, there are many more neighbors at farther distances than close distances, the dominance of the contribution to VDOS from $`r_{ij}<3.10\AA `$ for high-frequency modes suggests that these modes are more localized. For example, the highest-frequency mode appears to be a localized mode as contributions to this mode from second nearest neighbors and beyond are out-weighted by the contribution from the nearest neighbors, even though there are many more farther neighbors than nearest neighbors for every atom in the cluster. The inset of Fig. 6 shows in more detail the contributions to the non-zero lowest-frequency mode. The dominant contribution to this mode is apparently from the on-site environment. To gain more insight to these two interesting modes, one may first plot $`_\alpha (c_{i\alpha }^\lambda )^2`$ for the mode $`\omega _\lambda `$ as a function of site index $`i`$ to determine where the mode is localized and then decompose the site VDOS for that site according to distances (1st-, 2nd-, etc., based on Eqs. (8) and (9)) to shed light on the correlation between the vibration of the site and those of its neighbors. The upper inset of Fig. 7 is the plot of $`_\alpha (c_{i\alpha }^\lambda )^2`$ for the highest-frequency mode ( 108 THz) as a function of the site index. The mode is obviously localized in the immediate neighborhood of site 1 (the central site). The main panel of Fig. 7 shows the site VDOS $`\rho _1=_\alpha \rho _{1\alpha }(\omega )`$ decomposed in terms of 1st-, 2nd-, etc., nearest-neighbor distance using Eqs. (8) and (9). It can be seen that the contribution to this mode is almost entirely from the coupling of the central site (site 1) with its 4 nearest neighbors. In the lower inset of Fig. 7, the nature of the coupling is examined by decomposing $`\rho _1`$ in terms of $`_{j1}\rho _{1\alpha ,j\beta }`$. It is seen that the majority of the contribution are from $`xx`$-, $`yy`$-, and $`zz`$-nearest neighbor coupling. It is also seen that the mode is actually a composite of three almost degenerate modes (107.2, 108.0, and 107.5 THz), an indication of the distortion of the immediate environment of the central site from the perfect tetrahedral configuration.

The analysis of the non-zero lowest-frequency mode is shown in Fig. 8. The upper inset shows that the mode is an on-site mode at site 82 located on the surface of the cluster (see Fig. 1). The site VDOS $`\rho _{82}`$, decomposed according to Eq. (13), confirms this picture. The lower inset of Fig. 8 is the decomposition of $`\rho _{82}`$ according to Eq. (13) so that the on-site coupling can be examined. It shows that two almost equivalent on-site $`xy`$\- and $`yx`$-coupling are the main contributors to this mode, indicating that this is mainly a wagging mode in the $`xy`$ plane anchored by the two strengthened bonds of the atom at site 82 (a site with two dangling bonds, see Fig. 1).

The relaxed Si87 cluster is a system with no symmetry. The demonstration discussed in this section shows the flexibility of the methodologies developed in Sec. II to carry out detailed analysis of the dynamics of such a system and to trace interesting dynamical features to their electronic origin. Hence the framework described in Sec. II is expected to be extremely useful for the study and the analysis of vibrational dynamics of complex systems with low or no symmetry.

###### Acknowledgements.

This work was supported by the NSF (DMR-011284) and the DOE (DE-FG02-00ER45832) grants.

FIGURES

FIG. 1 The equilibrium configuration of the relaxed Si87 cluster. The $`z`$-axis is the vertical axis and the $`x`$-axis is the horizontal axis. Atom 82 (see the result and discussion in the text) is highlighted by the dark ball. FIG. 2 The decomposition of the VDOS in terms of the range of the bond charge $`q_i`$ associated with the atoms in the cluster. FIG.3 The decomposition of the VDOS in terms of the number of bonds associated with the atoms in the cluster. FIG.4 The decomposition of the VDOS in terms of the bond charge $`\sigma _{ij}`$ between pairs of atoms in the cluster. FIG.5 The decomposition of the VDOS according to the bond types. FIG.6 The decomposition of the VDOS in terms of the distance between pairs of atoms according to $`\widehat{\rho }_{i\alpha ,j\beta }`$ (Eq. (17)). The inset is the enlargement of the curves in the immediate neighborhood of the mode with the lowest non-zero frequency. Fig. 7 Upper inset shows the plot of $`_\alpha (c_{i\alpha }^\lambda )^2`$ vs. the site index $`i`$ for the highest frequency mode ( 108 THz). The main panel shows the site VDOS for site 1 (the central site), $`\rho _1`$, decomposed according to 1st nearest neighbor, 2nd nearest neighbor, etc., distance using Eqs. (8) and (9). The lower inset shows the decomposition of $`\rho _1`$ in terms of $`_{j1}\rho _{1\alpha ,j\beta }`$ in the immediate neighborhood of the mode with the highest frequency. Fig. 8 Upper inset shows the plot of $`_\alpha (c_{i\alpha }^\lambda )^2`$ vs. the site index for the lowest non-zero frequency mode. It can be seen that this mode is confined at the site 82 (a surface site, see Fig. 1). The main panel shows the site VDOS $`\rho _{82}`$ and its decomposition according to Eq. (13). The lower inset shows the decomposition of $`\rho _{82}`$ in terms of $`_j\overline{\rho }_{82\alpha ,j\beta }`$ in the immediate neighborhood of the mode with the lowest non-zero frequency.


format: text

text:
# 1 Introduction

## 1 Introduction

Nanoscopic magnetic arrays have been proposed recently as devices of ultrahigh-density information storage . On the other hand, it was announced recently  that an array of bistable magnetic wires can display a complex behavior. Namely, the authors discussed the effect of self-organized criticality (commonly abbreviated as SOC) . The problem is that a system in the critical state cannot be used to store any information, because a flip of one magnetic moment can change the magnetic structure of the whole array. The aim of this work is to check by means of a computer simulation if the array is indeed in the critical state. The values of the simulation parameters are selected as close as possible to the experimental system described in .

The dynamics of the total magnetization, i.e. its time dependence in the presence of the oscillating magnetic field, is governed by two agents: the long-range magnetostaticinteraction between each two elements of the array, and the switching field(coercive field), which varies from one element to another. The nonzeroswitching field enables to preserve information in the array at positive temperatures. The interaction leads to flips of magnetic moments and theinformation — otherwise frozen in the array — is destroyed, at leastpartially.

Let us recall the characteristic features of SOC. According to an illuminating paperby Flyvbjerg , excitations wandering in a system can lead it to a self-organized state, i.e. a spontaneously formed state far from thermalequilibrium. By criticality we mean that in this state, excitations can be of any order of magnitude. In fact there is no characteristic scale in a criticalstate: the array is scale-free in the same way as a cluster of spins “up” in a ferromagnet at its Curie temperature.

Basically, we search for excitations which spread over the lattice. The technique applied is known as damage spreading . It is believed, that if SOC is present in a system, a local variation of its structure can lead to a global change . Then, existence of such a spreading is an argument for the presence of criticality. On the other hand, one of the arguments of Refs.  is based on the calculation of the box-counting fractal dimension  of the cluster of elements with a given orientation of magnetic moments. That is why we investigate also this fractal dimension.

## 2 The system and the Pardavi-Horvath algorithm

Our simulations concentrate on the nanoscopic array described in Ref. . Monodomain magnetic cylinders of nickel, 57 nm of diameter and 115 nm length, form a square array with the period $`100`$ nm. The magnetization is $`M_s=370`$ emu/cm3. The mean switching field of one element is $`H_s=710`$ Oe, with the standard deviation $`105`$Oe. The switching field is constant in time for each element. The magnetostatic interaction is calculated by means of the rectangular prism approximation formula . The only modification introduced to the investigation of SOC are the periodic boundary conditions, to avoid boundaries which could alter the results.

The dynamics of the system is simulated by means of the Pardavi-Horvath (PH) algorithm . Initial parameters of this algorithm are an external magnetic field $`H_{ex}`$, the magnetic states of each element of the array and their switching fields $`H_s`$. At the begining usually $`H_{ex}=0`$ and the magnetic states of the elements are random. Then for every element of the array, the interaction field $`H_i`$ can be calculated which comes from other elements. The numbers which represent this interaction field are preserved in the computer memory. After that an element is sought which is most wiling to flip. This element has to fulfil three conditions:

1. Magnetization of this element must have opposite direction to the total field $`H_{tot}`$ exerted on this element ($`H_{tot}=H_i+H_{ex}`$).
2. Total field must have greater value than the switching field of the element ($`H_{tot}>H_s`$).
3. This element has maximal value of difference between the total field exerted on the element and its switching field.

The selected element is fliped, and all numbers which represent the interaction fields are corrected by the change of the interaction field coming from the flipped element.
This procedure is repeated until no one element can be flipped. Then, the external field is changed. Subsequent changes of this field are of the values which make possible a flip of at least one element. The whole procedure is repeated. Finally the hysteresis loop of the whole array is obtained.

The algorithm is checked by the calculation of thehysteresis loop (see Fig. 2) and a comparison with the experimental one  — the accordance is quantitative and good.

## 3 Results of simulation

### 3.1 The damage spreading

The damage spreading technique is applied for an array of $`100\times 100`$ elements. Two arrays are stored simultanously in the computer memory, in the same randomly selected initial state. The PH procedure is applied to lead these arrays to stable (stationary) states, where no flips occur. Then, one magnetic moment is flipped in one array, and wecheck if the flip is stable. Subsequently we apply a periodic magnetic field of amplitude $`H_m`$. We consider the case of small frequency, where the whole arrays go to stable states each time before the field is varied. In this case a particular value of the frequency is not relevant. The damage is defined as the Hamming distance between two arrays, defined as the number of elements of the arrays with opposite directions of the magnetic moments.

The result is that the Hamming distance increases only during some transient time. Then, the system reaches a limit cycle, with the length equal to some multiple of the period of the applied field. A typical result is shown in Fig. 3. We note that if a system is critical, the size of damages should increase until the system boundaries are reached . In our case damages remain limited to a closed area at the lattice center.

The size $`A`$ of this maximal damage has a maximum for the amplitude of the applied field $`H_m`$ close to 1300 Oe, and it vanishes in most cases for $`H_m>1400`$ Oe. It is obvious that any damage must disappearat strong fields, when the system is saturated. Besides, $`A`$ varies strongly from a sample to asample. As a rule, the damage is localized as a shapeless and seemingly random pattern of the array elements, formed of several clusters, separated but close to each other.

### 3.2 The fractal dimension

The fractal dimension $`D`$ of the cluster of the “up” (or “down”) oriented elements is calculated for the $`1024\times 1024`$ array.Obviously, the investigated magnetic state of the system should not be arbitrary, but it should be a consequence of the physical properties and processes. For a random initial state the fractal dimension $`D`$ is equal to 2.0, what means that the system has no fractal character.On the other hand, such a state is an artifact of the (pseudo)random numbers generator,and therefore it cannot be treated as realistic.The only physical condition of such a state is its stability, but in our opinion this condition leaves too much freedom. Therefore we have calculated the time evolution of the fractal dimension $`D`$ during the virtually performed hysteresis experiment. It is obvious, that $`D=2.0`$ or $`D=0.0`$ in the saturated states. We found that this bistable character of $`D`$ is preserved also for small amplitude $`H_m`$ of the applied external field. Between these two values, the scaling is not proper. An example is shown in Fig. 4(a). Typically, a part of the curve for small size of the box shows the inclination close to zero, and another part — close to two. Similar plots have been presented also in Refs. . As the applied field changes, the bent part of the curve is shifted left and finally disappears. In fact, plots like Fig. 4(a) appear when we try to assign the fractal character to a random configuration of spins. We demonstrate it in Fig. 4(b), where the magnetic configuration of the array produces the curve shown in Fig. 4(a). Surely, this cluster does not exhibit any self-similarity.

## 4 Conclusions

The fractal dimension $`D`$ obtained in Refs.  is $`1.97`$, which is close to non-fractal value 2.0. On the other hand, in the Ising ferromagnet at the Curie point $`D=187/961.95`$ . Then, the value of $`D`$ itself does not allow to state if the system is in the critical state.

However, the criterion of criticality can be drawn from the obtained curve. Its bent character means that there is a characteristic length in the investigated structure. The value of this length is just at X coordinate of the bend point. We demonstrate it in Fig. 5. The two curves are obtained with the same box-counting technique for the conventional Ising ferromagnet. The curve for temperature close to the $`T_c`$ is a fairly straight line, except its upper part where a continuous deviation from linearity is observed. This deviation can be assigned to a short waiting time (equivalent to $`N_{run}`$), which should be very large in the critical region, and to a finiteness of the system. The curve obtained for $`T=0.44T_c`$ displays the same characteristic length, as in Fig. 4(a).This is an indication, that the system is not scale-free, and therefore not in a critical state.

This conclusion is confirmed by the results obtained by means of the damage-spreading technique. In all cases investigated, we have not observed a damage which reaches the boundaries of the system.

However, this conclusion is valid only for the system investigated experimentally in Ref. . We stress this point because there are serious differences between the this system and the model system discussed in Refs. . In particular, in the latter the lattice is triangular, what makes the system frustrated because of the antiferromagnetic-like character of the magnetostatic interaction. Frustration enhances both the ground state energy and its degeneration. The system can wander in the phase space, and therefore damages can spread more easily. Then, the existence of SOC in an array with the triangular structure remains to be checked.

## Acknowledgements

The numerical calculations were partially carried out in ACK-CYFRONET-AGH. The machine time is financed by the Ministry of Scientific Research and Information Technology in Poland, grant No. KBN/HP-K460-XP/AGH/032/2002.


text:
# Renormalization of Anisotropy and Glueball Masses on Tadpole Improved Lattice Gauge Action

## I Introduction

Compact $`U(1)`$ gauge theory in (2+1) dimensions is one of the simplest models with dynamical gauge degrees of freedom and possesses some important similarities with QCD ham94 . The model has two essential features in common with QCD, confinement pol77 ; pol78  and chiral symmetry breaking fie90 . The theory is interesting in its own right, for it has analytically been shown to confine electrically charged particles even in the weak-coupling regime (at zero temperature)pol77 ; pol78 ; ban77 ; banks78 ; dre79 ; gof82 . The confinement is understood as a result of the dynamics of the monopoles which emerge due to the compactness of the gauge field. The string tension as a function of the coupling behaves in a similar fashion to that of the 4-dimensional $`SU(N)`$ lattice gauge theory. This model also allows us to work with large lattices with reasonable statistics. Other common features of compact $`U(1)_{(2+1)}`$ and QCD are the existence of a mass gap and of a confinement-deconfinement phase transition at some non-zero temperature. Thus, being reasonably simple and theoretically well understood in the weak-coupling limit, the $`U(1)`$ model provides a good testing ground for the development of new methods and new algorithmic approaches. In a recent paper loan02  we have obtained the first clear picture of the static quark potential, showing very clear evidence of the linear confining behaviour at large distances. The evidence of the the scaling behaviour of the string tension and the mass gap has also been observed in this model.

In the present paper we want to extend the analysis of Ref. loan02  in various respects. Since the measured ratio of spatial to temporal lattice spacings is not the same as the input parameter in the action, it becomes important to determine the true or renormalized anisotropy $`\xi _{phys.}`$ as a function of the bare anisotropy $`\xi _0`$. An important advantage of using anisotropic lattices has been the need to measure the renormalization of anisotropy in the simulation. The existing theoretical snipp97 ; garcia97  and numerical studies klas98  with Wilson action for $`SU(3)`$ lattice gauge theory have shown that at finite coupling $`g`$, the renormalized anisotropy $`\xi _{phys.}`$ differs appreciably form bare anisotropy $`\xi _0`$. However, it has recently been shown that the use of improved actions, supplemented by tadpole improvement, besides providing a better discretization scheme for QCD, offer the advantage of a significant reduction of renormalization of $`\xi _0`$ to a few percent morn97 ; mor96 . The anisotropy parameter $`\eta `$ (ratio of the renormalized and bare anisotropy) and anisotropic coefficients have been calculated to one-loop order for improved actions in various recent studies garcia97 ; sakai00 ; eng00 ; alford01 ; drummond02 . These calculations have provided very reliable results and the observed behaviour is confirmed non-perturbatively by large scale simulations on fine lattices morn97 ; shakes98 ; alford98 . We investigate the influence of tadpole improvement on isotropic and anisotropic lattices for the $`U(1)`$ model in (2+1) dimensions in reducing the renormalization of the bare anisotropy at weak and strong couplings. We also apply tadpole improved $`U(1)`$ lattice gauge theory to calculations of the static quark potential, the string tension and the scalar glueball masses and compare the results with simulations of the Wilson action.

The rest of this paper is organized as follows. After outlining the tadpole improved $`U(1)`$ gauge model in (2+1) dimensions in Sect. II, we describe the method for determination of the renormalized anisotropy in Sect. III. We present results from our simulations on anisotropic lattices using both standard and tadpole improved Wilson gauge action in Sect. IV. We compare bare and renormalized anisotropies, static quark potential and scalar glueball masses from these actions. We conclude in Sect. V with a summary and outlook on future work.

## II Compact $`U(1)`$ model in (2+1) dimensions

The tadpole-improved $`U(1)`$ gauge action on an anisotropic lattice can be written in the following form klas98 :

$$S=\beta \left[\underset{r,i>j}{}\frac{\xi _0}{u_s^4}(1P_{ij}(r))+\underset{r,i}{}\frac{1}{\xi _0u_s^2u_t^2}(1P_{it}(r))\right]$$
(1)
where $`P_{\mu \nu }`$ is the plaquette operator, $`\xi _0=\mathrm{\Delta }\tau =a_t/a_s`$ is the bare anisotropy at the classical level and $`u_s`$ and $`u_t`$ are the mean fields for the tadpole improvement. The notation used in Eq.(1) differ slightly from that used in Refs. garcia97 ; drummond02 , where the spatial and temporal mean-field improvement factors, $`u_s`$ and $`u_t`$ were absorbed into definition of $`\beta `$ and $`\xi _0`$. This, however, follows the notation introduced in Ref. morn97 .

On the anisotropic lattice, the mean fields are determined using the measured values of the average plaquettes mor97 . We first compute $`u_s`$ from spatial plaquettes, $`u_s^4=P_{ij}`$, and then we compute $`u_t`$ from temporal plaquettes, $`u_t^2u_s^2=P_{it}`$. Another way to determine the mean fields is to use the mean links in Landau gauge lep98

$$u_t=\text{ReTr}U_t,u_s=\text{ReTr}U_i$$
(2)
where the lattice version of the gauge condition is obtained by maximizing the quantity,

$$\underset{r,\mu }{}\frac{1}{u_\mu a_\mu ^2}\text{Tr}U_\mu (r)$$
(3)
Since the temporal lattice spacing in our simulations is very small, we adopt the following convention mor97 ; mor96 ; morn97  for the mean fields in tadpole improvement

$$u_t1,u_s=P_{ij}^{1/4}.$$
(4)
This prescription eliminates the need for gauge fixing and the results yield values for $`u_s`$ which differ from those using Landau gauge by only a few percent.

## III Renormalization of Anisotropy

Following the procedure of Klassen klas98  and Shakespeare and Trottier shakes98 , we measure the static quark potential extracted from Wilson loops in the spatial and temporal directions. Accordingly on an anisotropic lattice there are two potentials, $`V_{xt}(R)`$ and $`V_{xy}(R)`$. The two potentials differ by a factor of $`\xi _{phys.}`$ and by an additive constant, since the self-energy corrections to the static potential are different if the quark and anti-quark propagate along the temporal or a spatial direction. Thus $`\xi _{phys.}`$ can be determined by comparing the static quark potential computed from the logarithmic ratio of time-like Wilson loops $`R(x,\tau )`$, where

$$R(x,\tau )\frac{W_{xt}(x,\tau +1)}{W_{xt}(x,\tau )},$$
(5)
with the potential computed from that of the space-like Wilson loops $`R(x,y)`$, where

$$R(x,y)\frac{W_{xy}(x,y+1)}{W_{xt}(x,y)}.$$
(6)
Asymptotically, for large $`\tau `$ and $`y`$, the ratios $`R(x,\tau )`$ and $`R(x,y)`$ approach

$`R(x,\tau )`$ $`=`$ $`Z_{x\tau }\text{e}^{\tau V_{xt}}+(\text{excited state contr.})`$ (7)
$`R(x,y)`$ $`=`$ $`Z_{xy}\text{e}^{yV_{xy}}+(\text{excited state contr.}).`$ (8)
To suppress the excited state contributions, a simple APE smearing technique alb87 ; teper86 ; ishi83  was used. In this technique an iterative smearing procedure is used to construct Wilson loop (and glueball) operators with a very high degree of overlap with the lowest-lying state. In our single-link smoothing procedure, we replace every space-like link variable by

$$U_iP\left[\alpha U_i+\frac{(1\alpha )}{2}\underset{s}{}U_s\right]$$
(9)
where the sum over ‘$`s`$’ refers to the “staples”, or 3-link paths bracketing the given link on either side in the spatial plane, and P denotes a projection onto the group $`U(1)`$, achieved by renormalizing the magnitude to unity. We used a smearing parameter $`\alpha =0.7`$ and up to ten iterations of the smearing process. To reduce the statistical errors, the time-like Wilson loops were constructed from “thermally averaged” time-like links teper86 ; ishi83 .

The links making up the space-like and time-like Wilson loops are smeared by the same amount so that the ratios $`R(x,\tau )`$ and $`R(x,y)`$ have the same excited-state contribution. Similarly, finite-volume corrections to the $`R(x,\tau )`$ and $`R(x,y)`$ are the same if the temporal and spatial extents are equal in physical units, i.e. $`N_s=\xi _{phys.}N_t`$ in lattice units. These statements are expected to hold only for large $`x`$, $`y`$ and $`\tau `$; otherwise there can be large $`O(a_s^2,a_t^2)`$ lattice errors.

The physical anisotropy is determined from the ratio of the potentials $`V_{xt}(R)`$ and $`V_{xy}(R)`$ estimated from $`R_{xt}`$ and $`R_{xy}`$ respectively. The unphysical constant in the potentials is removed by subtraction of the simulation results at two different radii

$$\xi _{phys.}=\frac{V_{xt}(R_2)V_{xt}(R_1)}{V_{xy}(R_2)V_{xy}(R_1)}.$$
(10)
The measured renormalization of the anisotropy, $`\eta `$ is then determined from

$$\eta \frac{\xi _{phys.}}{\xi _0}.$$
(11)
## IV Simulation and Results

Simulations were performed on four lattices of $`N_s^2\times N_t`$ sites, with $`N_s=16`$ and $`N_t`$ ranging from 32 to 48 with mean-link improvement, and four lattices with the Wilson action. Configurations were generated by using the Metropolis algorithm. The details of the algorithm are discussed elsewhere loan02 . 50000 sweeps were performed for thermalization of the configurations and self-consistent determination of the tadpole factors. Configurations are stored every 250 sweeps thereafter. Ensembles of about 1000 configurations were used to measure the static quark potential, while 1,400 configurations, at coupling values from $`\beta =1.0`$ to 2.5, were generated for the glueball mass. We fixed $`\xi _0=16/N_t`$ in the first pass, so that the lattice size remains fixed at $`16a_s`$ in all directions. The simulation parameters of the lattices analyzed here are given in Table 1.

After measuring the Wilson loops at fixed values of $`\beta `$, we compute the ratios $`R_{xt}`$ and $`R_{xy}`$. We find that the individual ratios reach their plateaus for $`\tau 3`$ and $`y3`$ for fixed $`x`$ as shown in Figures 1 and 2. These ratios are are expected to be independent of $`\tau `$ and $`y`$ for $`\tau `$, $`y`$ $`3`$ respectively. The estimates of the potentials $`V_{xt}(R)`$ and $`V_{xy}(R)`$ can now be found from these ratios.

Figure 3 shows a graph of the static quark potentials, computed from spatial and temporal Wilson loops, as a function of radius $`R`$ at $`\beta =1.306`$ and $`\mathrm{\Delta }\tau =0.5`$. The potential in the lattice units obtained from the ratio of the time-like Wilson loops has been rescaled by the input anisotropy. To extract the string tension, the time-like potential is well fitted by a form

$$V_{xt}(R)=a+b\mathrm{ln}R+\sigma R,$$
(12)
including a logarithmic Coulomb term as expected for classical QED in (2+1) dimensions which dominates the behaviour at small distances, and a linear term as predicted by Polyakov pol78  and Göpfert and Mack gof82  dominating the behaviour at large distances and showing a clear evidence of the linear confining behaviour at large distances.

We measured each anisotropy twice, using two different radii $`R_1`$ for subtraction, with fixed $`R_2`$. Setting $`R_2=4`$, we computed the anisotropy with $`R_1=2`$ and $`\sqrt{2}`$. The two determinations of anisotropy, shown in Table 2, are in excellent agreement. The numerical values of the renormalization of the anisotropy parameter $`\eta `$ appears to be equal to unity even at large $`\beta `$. It is seen that for with mean-field improvement the input anisotropy is renormalized by few percent over the range of lattices analyzed here, whereas the measured value of anisotropy is about $`1520\%`$ lower than the bare anisotropy with the standard Wilson action. This can be seen from Figure 4 where the renormalization of the anisotropy is plainly visible as a difference in slope of the potentials computed from $`R_{xt}`$ and $`R_{xy}`$.

Glueball correlation functions $`C(\tau )`$ were also calculated

$$C(\tau )=\underset{\tau _0}{}0\overline{\mathrm{\Phi }}_i(\tau +\tau _0)\overline{\mathrm{\Phi }}_i(\tau _0)0$$
(13)
where $`\overline{\mathrm{\Phi }}_i(\tau )`$ is the optimized glueball operator found by a variational technique, following Morningstar and Peardon morn97  and Teper tep99 , from a linear combination of the basic operators $`\varphi _i`$,

$$\mathrm{\Phi }(\tau )=\underset{\alpha }{}v_{i\alpha }\varphi _{i\alpha }(\tau )$$
(14)
where the index $`\alpha `$ runs over the rectangular Wilson loops with dimensions $`l_x=[n1,n+1]`$, $`l_y=[n1,n+1]`$ and smearing $`n_s=[m1,m+1]`$, making 27 operators in all.

The optimized correlation function was fitted with the simple form

$$C_i=c_0+c_1\mathrm{cosh}m_i(T/2\tau )$$
(15)
to determine the glueball mass estimates.

The results for the symmetric and the anti-symmetric glueball masses over the square root of the string tension are shown in Table 3, along with the mean plaquette values at different $`\beta `$ at $`\mathrm{\Delta }\tau =1.0`$. Figure 5 shows the behaviour of the logarithm of the antisymmetric mass gap over the square root of the string tension as a function of $`\beta `$. It can be seen that that the ratio scales exponentially to zero in the weak-coupling limit as should be in three-dimensional confining theories. The solid line is a fit to the data over the range $`0.916\beta 2.12`$. The slope of the data matches the predicted form gof82 , however, the intercept of the scaling curve is large by a factor of 2 (our previous estimates of constant coefficient for standard Wilson action are large by a factor of 5.2 loan02 ). It would be interesting to test the sensitivity of the slope and intercept of the scaling curve by including the radiative corrections.

A plot of mass ratio against the effective lattice spacing $`a_{eff}`$ loan02  is shown in Figure 6. At weak coupling, the theory is expected to approach a theory of free bosons gof82  so that symmetric state will be composed of two $`0^{}`$ bosons and the mass ration should approach two. The mass ratio of the lowest glueball states scale against the effective lattice spacing towards a value close to 2.0, as expected for a theory of free scalar bosons. We note that mass ratio exhibits scaling behaviour, even with the Wilson action loan02 , however, in contrast with the Wilson action, a significant reduction in the errors with tadpole improved action is apparent in the mass gap and the mass ratio.

## V Summary and Outlook

Mean field improved $`U(1)`$ lattice gauge theory in (2+1) dimensions was applied to calculations of the static quark potential, the renormalized lattice anisotropy and the scalar glueball masses. We analyzed the mean-link improved action on isotropic and anisotropic lattices and comparisons were made with the simulations of the Wilson action. By comparing the static quark potentials computed from space-like and time-like Wilson loops, we determined the physical anisotropy of the tadpole improved Wilson action. We found that with mean-link improved Wilson action, the bare anisotropy is renormalized by less than a few percent, in contrast with the standard Wilson action, where the measured value of anisotropy is found to be about $`1520\%`$ lower than bare anisotropy on the lattices analyzed here. We found that tadpole improvement significantly reduces discretization errors in the static quark potential and the glueball masses. The mass ratio of the two lowest glueball states scales against the effective lattice spacing towards a value close to 2.0, as expected for a theory of free scalar bosons.

We intend to extend PIMC techniques to Symanzik improved $`U(1)`$ lattice gauge theory. The intention is to study the effects of improvement on the scaling slope, the constant coefficients and scaling behaviour observed in the weak-coupling regime of the theory. We also plan the study the one-loop correction to the anisotropy factor for Symanzik improved $`U(1)`$ gauge action in three dimensions. We shall report on this work in the near future.

###### Acknowledgements.

This work was supported by the Australian Research Council. We are grateful for access to the computing facilities of the Australian Centre for Advanced Computing and Communications (ac3) and the Australian Partnership for Advanced Computing (APAC).


text:
# Core Lexicon and Contagious Words

## Abstract

We present the new empirical parameter $`f_c`$, the most probable usage frequency of a word in a language, computed via the distribution of documents over frequency $`x`$ of the word. This parameter allows for filtering the core lexicon of a language from the content words, which tend to be extremely frequent in some texts written in specific genres or by certain authors. Distributions of documents over frequencies for such words display long tails as $`x>f_c`$ representing a bunch of documents in which such words are used in abundance. Collections of such documents exhibit a percolation like phase transition as the coarse grain of frequency $`\mathrm{\Delta }f`$ (flattening out the strongly irregular frequency data series) approaches the critical value $`f_c`$.

Studies of lists of words arranged in terms of their frequencies belong to the most important domains of quantitative linguistics Zipf . Detection of most frequent words constituting the core lexicon of a language is important not only for foreign language learners, but also for various practical applications, including text compression (this was recognized as early as in S ), speech recognition Lau , information retrieval SJ , etc. It is easy to order words with respect to their mean frequency of uses, which is typically measured as the number of instances of a word normalized by the sample of one million words (ipm, instances per million words); though the notion of a word should cover all word forms, like goes, went, gone for go.

However, the mean frequency is not a sufficient selection criterion, because of the large relative dispersion of the word frequencies which vary very much from one text to the next especially in ample and diverse collections of documents. Some words (like prepositions) occur in many texts with predictable rates, others (like pronouns or mental verbs) are significantly more frequent for certain writers or genres, while some are ”contagious”: these words (such as proper names, technical terms, abbreviations, etc.) appear in just a few documents, but when they appear, they are often found in abundance CG . The variability of rates of words can be characterized in a variety of ways, including the Poisson K-mixtures CG . It can be measured by the coefficient of variation (the standard deviation divided by the mean), as in Kilgariff . However, the coefficient of variation as a measure of relative dispersion is not very useful when the average frequency is close to zero, which occurs quite often for the semantically loaded words. Another way to measure the variability of rates for contagious words is to compute the document frequency (or inverse document frequency, SJ ) by counting the number of documents the given word is mentioned in, the burstiness parameter that is the mean frequency, except that it ignores documents with no instances of the word (see references in CG ), etc. As is evident, each of these parameters does not capture much of the heterogeneous structure of word rates series and none of them provides any general approach to describe it.

In this paper, we approach the problem of detection of contagious words and selection for the core lexicon of a language from a probabilistic point of view. In accordance to it, the frequency $`x`$ of a word (counted as the number of its instances per million words observed in any document of a given language) is a random positive variable distributed with some (unknown) probability density function $`\rho (x)`$. Strongly irregular frequency data series can be flattened out by introducing the frequency coarse grain $`\mathrm{\Delta }f.`$ The statistics of any word $`w`$ can be characterized by the number of documents $`N_w(n,\mathrm{\Delta }f)=N\left[(n1)\mathrm{\Delta }f<xn\mathrm{\Delta }f\right]`$ for which the rate of the word $`w`$ drops into the $`n`$th frequency interval for various $`\mathrm{\Delta }f`$.

We have found that for any coarse grain $`\mathrm{\Delta }f`$ the distributions of documents $`N_w(n,\mathrm{\Delta }f)`$ over $`n`$ are the asymmetric curves having one maximum $`f_c`$ that is the most probable frequency at which the word $`w`$ would appear in a randomly chosen document written in the given language. The value of $`f_c`$ is independent of genres, authors, and topics of documents and is an intrinsic characteristic of the word in a contemporary language. Obviously, $`f_c`$ varies as the language evolves approaching zero as the word becomes obsolete. For the frequencies $`x`$ close to the distribution maximum $`f_c`$, distributions of documents $`N_w`$ over $`n`$ are bell shaped, but have anomalous tails as $`|xf_c|`$ is large enough. We can estimate the most probable frequencies $`f_c`$ of words independently by two methods: first, from the distributions of documents $`N_w(n,\mathrm{\Delta }f)`$ as $`\mathrm{\Delta }f1`$ and, second, from the distributions of authors using the same word in their texts. For any word, both methods gave identical values of $`f_c`$.

The most probable frequency $`f_c`$ helps to detect the content words and to select words for the core lexicon of a language. Common words which appear uniformly in most documents (like prepositions, conjunctions, relational verbs, some size adjectives, etc.) have usually relatively high mean rates $`\overline{f}`$ (ipm), and $`f_c\overline{f}`$. They obviously belong to the core lexicon of a language.

The use of semantically loaded words depends essentially on authors, genres, and topics. Despite they are not found in a bunch of documents, their mean rates $`\overline{f}`$ are still very high because of their excessive popularity in certain collections of texts, but their most probable rates get down $`f_c\overline{f}`$ indicating the presence of long tails in distributions of documents $`N_w`$ over $`n`$ as $`n>f_c/\mathrm{\Delta }f.`$ Eventually, for contagious words found in abundance in just a few documents, the distributions $`N_w`$ over $`n`$ have long tails, however their mean rates $`\overline{f}`$ are very low since $`x=0`$ for almost all texts, and $`f_c\overline{f}`$ or even $`f_c>\overline{f}.`$

Our study of word rates is based on the reference corpus of Russian corpus , which includes more than $`410^7`$ words in $`1566`$ texts, which are balanced in their coverage of various genres: fiction, newspapers, various informative texts originally written in Russian from 1980 to 2002. Unlike earlier corpora (of about 1 million words), reference corpora of this size are close to saturation, namely, any collection of new documents added to the corpus does not cause statistically significant changes to the frequency and patterns of uses of its words.

The corpus informs us about the list of 5000 words most frequently used in modern Russian 5000 . The study of $`f_c`$ shows that for about 60$`\%`$ of them their most probable frequencies is $`f_c50`$ ipm, and just 2.51$`\%`$ have $`f_c600`$ ipm. Among the words having typically very high most probable frequencies, one can mention conjunctions and prepositions, pronouns and relational verbs, some motion verbs, and size adjectives. The proper nouns, acronyms, technical terms and other semantically loaded words typically have comparably small values of $`f_c`$.

Document counts in the $`n`$th frequency interval obviously decrease with $`n`$ as $`n>f_c/\mathrm{\Delta }f`$. We have found that it decays exponentially with $`n`$ for rather small coarse grains $`\mathrm{\Delta }ff_c,`$

$$N_w(n,\mathrm{\Delta }f)\mathrm{exp}\left[\frac{n\mathrm{\Delta }ff_c}{\xi _w}\right],\mathrm{\Delta }ff_c,$$
(1)
where $`\xi _w`$ plays the role of a ”correlation length” of the word $`w`$ and diverges as the scale $`\mathrm{\Delta }f`$ approaches the critical value $`f_c`$ as $`\xi _w\left|f_c\mathrm{\Delta }f\right|^{\alpha _w}`$ with the positive index $`\alpha _w`$ which is close to unity for the majority of words (see the data of Table 1). The parameter $`\xi _w`$ casts the characteristic excess of the word frequency $`x`$ over $`f_c`$.

We have observed that the values of $`\alpha _w`$ for words belonging to the same semantic group (such as relational verbs, motion verbs, perception verbs, some size adjectives, pronouns, conjunctions, and prepositions) are very close even if other values of empirical parameters measuring the variability of their rates are rather diverse.

For larger scales $`\mathrm{\Delta }ff_c,`$ for the frequency intervals with $`n>f_c/\mathrm{\Delta }f`$, the distributions $`N_w`$ are scale free,

$$N_w(n,\mathrm{\Delta }f)(n\mathrm{\Delta }ff_c)^{\beta _w},$$
(2)
where the index $`\beta _w>1`$ (see the data of Table 1).

The data of Table 1 show that the value of $`\beta _w`$ grows up with $`f_c`$ almost linearly (the coefficient of linear correlation between $`f_c`$ and $`\beta _w`$ in Table 1 is 0.93). For the supercritical phase $`\mathrm{\Delta }ff_c,`$ the tail of the probability distribution $`N_w`$ forms a stretched exponential.

Let us note that the asymptotic behaviors (1) and (2) are typical for the subcritical phase and the critical regime of percolation systems percol . Herewith, $`\mathrm{\Delta }f`$ plays the role of the order parameter, and $`f_c`$ is its critical value. A percolation-like phase transition observed with respect to a word in the collections of documents in which the frequency $`x`$ of this word exceeds its most probable rate $`f_c`$ gives us an evidence of existence of the genres of literature. Content words are of particular interest since their usage features usually the content of a text, so that the corpus of texts in which such words pile up can be interpreted as the literature genre possessing a special lexicon.

In this brief report we have studied the empirical distributions of documents over the frequencies of Russian words computed on the linguistic corpus of 1566 texts (of $`410^7`$ words). The approach to the word frequency analysis which we have proposed is very general and can be applied for any other human (or artificial) language. We have introduced the new empirical parameter $`f_c,`$ the most probable usage frequency of a word in a contemporary language. The value of $`f_c`$ is independent on authors, genres, and topics, but obviously varies in time as the language evolves. The most probable frequency of a word could be useful in studies devoted to the evolution of languages. This parameter helps us to handle the heterogeneous structure of word rate series and to determine whether they represent core lexicon. Distributions of documents over frequencies for the semantically loaded words which are found in abundance in a few documents have remarkably long tails. The typical excess of the word rate over $`f_c`$ which plays the role of the correlation length in (1) can be used in the automatic recognition of grammatical functions of words in a language. We have shown that the collections of documents accumulating content words exhibit a percolation like phase transition uprising the certain genres of literature appropriate of specialized lexicons or terminologies.

One of the authors (S.S.) is supported by the Alexander von Humboldt Foundation (Germany). D.V. benefits from the BiBoS Research Center (Bielefeld University) as a guest researcher. Authors thank T. Krüger for the fruitful discussion.


text:
# The surface density of Extremely Red Objects in high-z quasar fields

## 1 Introduction

EROs are galaxies with very red optical to near-infrared colours. Several definitions of EROs are in use, but the among the most common are $`RK5`$ and $`IK4`$. The colour criterion selects $`z`$ 1–2 early-type galaxies with passively evolving stellar populations characterized by strong 4000 Å breaks in their spectral energy distributions, but also dusty starbursting galaxies which have their UV-flux attenuated by internal dust. The first EROs were discovered in the deep near-infrared surveys by Elston, Rieke & Rieke elston88 , and since then there have been several surveys aimed at studying their spectroscopic cimatti02 , photometric ss00 ; smail02 ; mannucci02  and clustering properties daddi00 ; roche02 .

Here, we have searched for EROs in the fields of $`z=1.8`$–3.0 quasars. It has been suggested that the fields of high-redshift radio galaxies and quasars often have an excess of EROs mccarthy92 ; hr94 ; cimatti00 ; thompson00 ; hall01 . E.g. Cimatti et al. cimatti00  find an excess of $`RK6`$ EROs in 14 radio galaxy and quasars fields at $`1.5<z<2.0`$ and Hall et al. hall01  find that EROs are about 3 times more numerous in the fields of $`z1.5`$ quasars. Several possibilities for the ERO excess exist, but the favored explanation seems to be that the EROs are early-type galaxies associated with the quasars, perhaps signifying a forming cluster or a galaxy group at the quasar redshift. Another possibility is that the ERO excess traces foreground clusters or overdensities in the ERO population which gravitationally magnify the quasars, boosting them into flux-limited samples. In either case, the ERO excess may represent clusters in formation at cosmologically significant epochs. Care must be taken, however, because the angular distribution of EROs on the sky is characterized by overdensities and large voids daddi00 ; roche02 . This causes a strong field-to-field variation in the ERO surface density, especially if the fields are small.

## 2 The quasar sample

The quasar fields presented here are part of a larger ongoing survey to study the host galaxies (with the NIRC111Near Infrared Camera on Keck I) and environments of radio-loud and radio-quiet quasars at $`1.8<z<3.0`$. The full sample consists of 40 quasars from the catalogs by Barthel et al. barthel88  and Hewitt & Burbidge hb93 , and are nearly evenly distributed between radio-loud and radio-quiet systems with comparable $`V`$ magnitudes and redshifts. In order to study the environments of these quasars, we have obtained deep $`R`$ and $`K_s`$ images with the COSMIC222Carnegie Observatories Spectroscopic Multislit and Imaging Camera and the PFIRCAM333Prime Focus Infrared Camera instruments on the Palomar 200 inch telescope. The $`R`$-band imaging is completed, and a subset of 13 radio-loud quasar fields have been imaged in $`K_s`$. The ERO samples are drawn from an area of 2.1 $`\times `$2.1 arcmin2 around each quasar (the region of overlap of the $`R`$ and $`K_s`$ imaging). The 5$`\sigma `$ detection limits in the $`R`$\- and the $`K_s`$-filters are $`25`$ and $`20.5`$, respectively, and this has allowed us to search for EROs in the quasar fields using the $`RK5`$ ERO definition.

## 3 Results

The colour-magnitude diagram of objects detected in the quasar fields is shown in Fig. 1. In total, there are 160 objects with $`RK5.0`$. Using the sextractor bertin96  star-galaxy classifier, and treating all objects at $`K_s<18`$ with a star-galaxy class greater than 0.85 as stars, we find that six of the 160 EROs are possible stars. Excluding these yields a surface density of 1.26 EROs per arcmin2 at $`K_s19`$, and 2.35 arcmin-2 at $`K_s20`$. This surface density is two–three times higher than what is found in wide-field random surveys, typically 0.5 at $`K_s19`$ and 1.4 arcmin-2 at $`K20`$ daddi00 ; roche02 . Although there is a large variation in the ERO surface density from field to field, approximately 30–40 % of the quasar fields have notably large overdensities, a relatively large fraction of the sample in studies like this.

At the resolution of our images ($`1`$ arcsec), it is difficult to distinguish between stars and extended objects, and sextractor has problems with star-galaxy separation at faint magnitudes. We therefore performed the star-galaxy separation only at $`K_s<18`$. The contamination of faint red stars among EROs is uncertain, but in our case it is likely to be small since the quasar fields lie at high galactic latitudes ($`>35`$ deg, except three which lie between 21 and 27 deg). For completeness, we also carried out the analysis using star-galaxy separation regardless of $`K_s`$ flux. This resulted in $`20`$ % (as compared to 4 %) of the EROs being classified as stars. This of course lowers the ERO counts, but our conclusions remain unaltered. The surface densities using this approach are 1.08 at $`K_s19`$ and 1.93 arcmin-2 at $`K_s20`$. We note that Mannucci et al. mannucci02  find $`10`$ % stars at $`K^{}<20`$ based on $`R`$, $`J`$ and $`K^{}`$ imaging.

In Fig. 2, we compare the average cumulative surface density of EROs in our quasar fields with random field surveys from the literature. The field counts are represented by the surveys of Daddi et al. daddi00 , Roche et al. roche02  Cimatti et al. cimatti02  and Scodeggio & Silva ss00 , and they are seen to agree well with each other, except the counts in the Chandra Deep Field (CDF) by Scodeggio & Silva. This discrepancy is most likely caused by the inhomogeneous distribution of EROs on the sky, given that the area of the CDF is only 43 arcmin2. The counts in our quasar fields are seen to be consistently higher than the random field counts at all flux levels. This is not what we expect if we are sampling the random field ERO population with our 13 widely separated fields. The average ERO surface density from several small, widely separated fields should in principle be similar to the surface density in a large, contiguous field. Hence, it appears that the quasar fields are biassed toward regions of high ERO density.

Using the average of Daddi et al’s and Roche et al’s counts, we predict that our quasar fields should contain $`27`$ EROs at $`K_s19`$, whereas the observed number is 79. There is thus a clear excess of EROs in the quasar fields, but given that the EROs cluster strongly and therefore that the field-to-field variation is large, is the excess we observe significant?

In order to answer this question, we calculate (a) the expected number of EROs over the total area of our fields from the literature surveys, $`N_{field}`$, and (b), the expected variance, $`\sigma ^2`$, of the counts given that the ERO distribution is non-random, i.e. $`\sigma ^2=N_{field}\left(1+N_{field}A_\omega C\right)`$ daddi00 , where $`C`$ is a factor that depends on the area of the field and $`A_\omega `$ is the amplitude of the ERO angular two-point correlation function as found by Daddi et al. daddi00  and Roche et al. roche02 . The significance of the excess in the quasar fields is thereafter evaluated as $`(N_{qso}N_{field})/\sigma `$, where $`N_{qso}`$ is the observed number of EROs in the quasar fields. The result of this calculation shows that the significance of the excess at $`K_s18`$, $`19`$ and $`20`$ is 3.1,4.6 and 2.3$`\sigma `$, respectively. The ERO excess in the quasar fields therefore appears to be significant, even after taking into account the large fluctuations in the ERO surface density.

## 4 Discussion

Our analysis strongly suggests that the quasar lines of sight are biassed toward regions of high ERO density, and our result is consistent with that of Cimatti et al. cimatti00  and Hall et al. hall01 . At the moment it is impossible to tell whether the EROs are at the same redshift as the quasars or if they are in the foreground. In order to address this, either photometric or spectroscopic redshifts are needed. We therefore offer two different explanations for the observed excess and discuss each of them below:

The excess might be caused by galaxies, both early types and dusty starbursts, at the quasar redshifts. In this case, a major part of the EROs could be the early-type population of forming(?) clusters or groups of galaxies physically associated with the quasars. The evidence in favor of this consists of observations finding that luminous quasars and radio galaxies often reside in regions of high galaxy density. Another possibility is that the quasars are markers of large-scale structure traced by the EROs, such as walls or sheets of galaxies, as suggested by Thompson, Aftreth & Soifer thompson00 . Our current $`K_s`$-band images are too small to determine whether the EROs cluster around the quasars or if they are just part of a larger-scale overdensity, but we are currently awaiting wide-field (9$`\times `$9 arcmin2) near-infrared data to properly address this issue.

There is no evidence for a significant excess in the $`R`$\- and $`K_s`$-band counts in our fields, which we might expect if there are relatively rich clusters in the fields. However, this does not exclude the possibility that the excess is caused by clusters or groups of galaxies because most cluster members would simply be swamped in the large numbers of foreground sources. By selection on $`RK_s`$ colour we can find the ellipticals at $`z`$ 1–2 preferentially against much of the confusing foreground galaxies.

Another interesting possibility is that the excess is caused by EROs foreground to the quasars, either associated with clusters or groups of galaxies, or just overdensities (‘clumps’) in the ERO population. In this case, lensing magnification bias is likely to be important for the quasars. By pointing the telescope to a distant quasar, one is more likely to pick out lines of sight where there are mass concentrations along the line of sight which gravitationally amplify the quasar light.

## 5 Future and ongoing work

In order to address the question of how extended the ERO excess is, we have started an imaging programme in $`J`$ and $`K_s`$ using the new WIRC444Wide-field IR Camera instrument on the Palomar 200 inch. The WIRC has a field of view of 9$`\times `$9 arcmin2 which matches the size of our wide-field $`R`$-band images. We are also planning to obtain spectra and morphologies of the EROs. This will allow us to determine redshifts, spectral types and structural parameters of the EROs, thereby opening up for an investigation of the Fundamental Plane of early-type galaxies at cosmologically significant epochs.


text:
# Screening and Anti-Screening Effects in 𝐽/𝜓 Production on Nuclei

## Abstract

Nuclear effects in $`J/\psi `$ hadro- and electroproduction on nuclei are considered in framework of reggeon approach. It is shown that screening regime which holds for electroproduction at $`x_F0.7`$ and for hadroproduction at $`x_F(0.3÷0.4)`$ is changed with anti-screening regime for smaller $`x_F`$ values.

Heavy-quark production on nuclei provides an important information on strong-interaction mechanisms and is intensively discussed (see e.g. refs  for review of experimental data, discussion of some phenomenological models and more references). In this note we will discuss a phenomenon of changing screening regime to anti-screening one in $`J/\psi `$ hadro- and electroproduction when $`x_\psi x_F`$ decreases in the framework of the BCKT model  based on the reggeon approach. For $`J/\psi `$ hadroproduction it happens in the region of negative $`x_\psi `$. An investigation of this region can provide an important information on dynamics of charmonium production and may discriminate different dynamical models .

Nuclear effects are usually discussed in terms of conventional power-law parameterization $`F_A(x_\psi )F_N(x_\psi )A^{\alpha (x_\psi )}`$, where $`F_N(x_\psi )`$ ( $`F_A(x_\psi )`$ ) is inclusive $`J/\psi `$ production cross section on a nucleon (on a nucleus). The function $`\alpha (x_\psi )`$ characterizes nuclear effects at different longitudinal momentum fraction $`x_\psi `$.

Experimental data on $`J/\psi `$ hadroproduction reveal a striking contradiction with simplest theoretical expectations. Experimentally  function $`\alpha (x_\psi )`$ decreases from $`0.93÷0.95`$ at $`x_\psi 0`$ to values $`0.75`$ at $`x_\psi 0.8`$ thus indicating an increase of the absorption as $`x_\psi `$ increases. Formation-time mechanism predicts an opposite behavior. Instead of expected scaling with $`p_{J/\psi }^{lab}`$ experiment shows an approximate Feynman scaling with $`x_\psi `$ . Comparison of $`J/\psi `$ production data at different energies exhibits an explicit breakdown of QCD factorization theorem for this process .

On the other hand when comparing $`x`$ dependence of $`\alpha `$ for $`J/\psi `$ hadroproduction with one for light-quark particle hadroproduction (inclusive production of pions, nucleons, lambdas, etc.), one observes the same trend – a small absorption at $`x0`$, large absorption at $`x1`$ and approximate Feynman scaling. The difference is only quantitative one – $`\alpha `$ is about $`0.4÷0.5`$ at $`x1`$ for light quarks instead of $`0.7`$ for heavy quarks. This behavior of $`\alpha (x)`$ for light hadrons allows a natural explanation in Regge theory. Small absorption at $`x0`$ is due to Abramovskii–Gribov–Kancheli (AGK) cancellation for inclusive spectra , and increasing absorption at high $`x`$ is due to violation of AGK rules because of momentum conservation requirement .

In ref. the model for heavy quark and lepton pair production was constructed taking into account the most essential aspects of reggeon approach. In the spirit of parton picture  it was suggested that a fast projectile looks due to quantum-mechanical fluctuations like a cloud of virtual particles consisting of light partons (quarks and gluons, or light-quark hadrons, mainly pions) and (with some small probability) of heavy partons (say $`c\overline{c}`$ pair).333 In reggeon theory a mechanism of emitting these partons is supposed to be a multiperipheral one. Creation of heavy quarks in this model can be considered as an example of ”intrinsic charm” mechanism which was discussed by Brodsky et al. . In reggeon approach one can consider on equal footing any heavy component (e.g. “an intrinsic lepton pair” and so on). Different constituents of this fluctuation interact with nuclear matter and determine a dependence upon atomic number $`A`$. There are both elastic and inelastic interactions which from the viewpoint of reggeon diagrams are the different discontinuities of the same reggeon diagram. As a result there exist definite numerical relations between the discontinuities of different types (AGK rules ).

As it was stressed in ref. one has to distinguish two types of reggeon diagram cuttings. For diagrams of the first class the registered particle is contained inside a pomeron and appears in intermediate state only due to cutting of this pomeron. In this case the AGK rules are always valid and all contributions of many-pomeron diagrams to inclusive spectra cancel. It happens because each additional pomeron can be both cut and uncut, and due to the opposite signs of inelastic and elastic rescatterings they cancel. As a result, only one-pomeron diagram contribution survives giving a linear A dependence of spectra. This situation is typical for low $`x`$ particles.

Another situation occurs if a particle to be registered is contained inside “vertex” of reggeon diagram (we use a term “vertex” for part of diagram common for several pomerons attached to it). In this case AGK cancellation is not valid in general. The striking manifestation of this mechanism comes from $`A`$ dependence of inclusive spectra at large $`x`$ close to 1 (momentum conservation mechanism ). The fastest particles belong to the “vertex” but not to the pomeron, and corresponding discontinuities give $`\alpha 1/3`$ for large $`A`$ compared to $`\alpha 1`$ for small $`x`$.

Similar effect exists if one observes the “vertex” particle of some particular type, e.g. the $`J/\psi `$ meson. There are two sources of nuclear effects in this case . The first one is connected with rescatterings of light partons in fluctuation on different target nucleons or, alternatively in the anti-laboratory frame, with fusion of fluctuations originated from different nucleons of fast nucleus. This mechanism plays a small role at present energies $`p_{lab}10^3GeV`$  and its contribution is inessential for anti-screening effects we are interested in.

Another contribution is connected to rescatterings of the charmed state itself. One of the important ingredient of the model of ref. is an assumption that it is not the $`J/\psi `$ meson that propagates along nucleus but some primary colorless system containing both $`c\overline{c}`$ quarks (in color state) and light quarks (to screen this color charge). This state denoted as $`X`$ can be $`D\overline{D}`$ mesons, or $`D^{}\overline{D}^{}`$ state or something else, but the crucial point is that this system is of large size and therefore strongly interacts with nucleons (with cross section of order of $`20÷30mb`$). Only at last stage this system $`X`$ converts into $`J/\psi `$ meson (or in $`\psi ^{}`$) which is registered experimentally.

After inelastic rescattering the state $`X`$ can be lost for registration (e.g. due to conversion to a state with small projection into $`J/\psi `$). One can describe this effect quantitatively introducing a probability $`ϵ`$ to find the state $`X`$ after its inelastic rescattering. Then for $`ϵ=1`$ one has exact AGK cancellation, for $`ϵ=0`$ the situation is similar to momentum conservation mechanism for light mesons at $`x`$ close to 1 (no AGK cancellation) with $`\alpha 1/3`$, and for $`0<ϵ<1`$ the situation is an intermediate one with incomplete AGK cancellation. Moreover, the momentum of the state $`X`$ decreases after each inelastic rescattering.This redistribution of the longitudinal momentum of the state $`X`$ is essential for anti-screening effects discussed below.

Let us denote the momentum fraction distribution of system $`X`$ immediately after its production as $`F_0(x)`$. The state $`X`$ has some projection into state $`J/\psi `$ (as well as into $`\psi ^{}`$, etc.)444 For simplicity we suggest that only one state $`X`$ has a noticeable projection into $`J/\psi `$. In order to describe properly the observed difference of nuclear dependencies for $`J/\psi `$ and $`\psi ^{}`$, one has to consider at least two primary states, $`X_1`$ and $`X_2`$, with different projections to $`J/\psi `$ and $`\psi ^{}`$ and slightly different interactions with nuclear matter. . The $`x_\psi `$ distribution can be obtained from $`F_0(x)`$ by convolution with some projection function $`G_\psi `$:

$$F_N^{(\psi )}(x_\psi )=F_1G_\psi $$
(1)
where the notation was used

$$(fg)(x)=_x^1\frac{dz}{z}f(z)g(x/z),$$
(2)
or, in the rapidity variables,

$$(fg)(\overline{y})=_0^{\overline{y}}𝑑\xi f(\xi )g(\overline{y}\xi ),$$
(3)
where $`\overline{y}=Yy`$ is a rapidity in the anti-lab frame.

The state $`X`$ travelling through nuclear matter suffers both elastic and inelastic rescatterings. As in ref.  the change of $`x`$ distribution due to inelastic rescattering will be described as a convolution with some function $`ϵG(z)`$, where $`G(z)`$ is normalized $`(\underset{0}{\overset{1}{}}G(z)𝑑z/z=1)`$, and the parameter $`ϵ`$ determines a probability to have again the state $`X`$ after inelastic interaction (though with different momentum). After $`k`$ inelastic rescatterings of the state $`X`$ we have $`k`$-fold convolution:

$$F_k^{(\psi )}=F_1\underset{k\text{ times}}{\underset{}{ϵG\mathrm{}ϵG}}G_\psi ,$$
(4)
It is essential that the operation of convolution is commutative and associative, so we can first convolute $`F_1`$ with $`G_\psi `$ in the Eq.(4) and thus get the $`k`$-fold convolution of the $`J/\psi `$ spectrum on the nucleon, $`F_N^\psi (x)`$, with rescattering functions $`G`$:

$$F_k^\psi =F_N^\psi \underset{k\text{ times}}{\underset{}{ϵG\mathrm{}ϵG}}.$$
(5)
In order to obtain the $`J/\psi `$ inclusive spectrum on a nucleus, the functions $`F_k^{(\psi )}`$ should be weighted with cross sections $`\sigma _A^{(k+1)}`$ for $`k`$ inelastic rescattering after $`X`$ production:

$$F_A^\psi =\underset{k=0}{\overset{\mathrm{}}{}}F_k(x_\psi )\sigma _A^{(k+1)}.$$
(6)
The explicit form of the weights $`\sigma _A^{(k+1)}`$ depends upon energy region . However the change $`\sigma _A^{(k+1)}`$ is rather smooth and for simplicity we will use the same low-energy expressions as in ref.:

$$\sigma _A^{(k+1)}=\sigma _{aN}^{(\psi )}d^2b_0^{T_A(b)}𝑑v\mathrm{exp}(v\sigma _X)(v\sigma _X)^k/k!,$$
(7)
where $`\sigma _{aN}^{(\psi )}`$ is the cross section of $`J/\psi `$ production on a nucleon, $`\sigma _X`$ is the cross section of interaction of the state $`X`$ with a nucleon, and $`T_A(b)`$ is two-dimensional nuclear density as a function of impact parameter $`b`$ (nuclear profile). These formulas correspond to rescatterings ordered in longitudinal direction (low-energy regime) and after summation over $`k`$ give well known optical-type formula for the production cross section. Thus we come to Eq.(6) for $`J/\psi `$ nuclear inclusive cross section where $`F_k(x)`$ and $`\sigma _A^{(k+1)}`$ are defined by Eqs.(5) and (7).

Let us discuss parameterizations of distribution functions entering the model. The function $`F_{aN}^\psi (x_\psi )`$ can be determined from experimental data on $`J/\psi `$ production on nucleon. Its form strongly depends on a type of a projectile $`a`$. For hadron beam it can be parameterized as

$$F_{aN}^\psi (x)=C_a(1x)^{\beta _a},$$
(8)
where $`\beta _\pi 2`$ for pion and $`\beta _p3`$ for proton beam. For $`J/\psi `$ production by photons (and electrons) the form of this function is strongly different: produced $`\psi `$ particles are concentrated near $`x_\psi 1`$, and $`x_\psi `$ distribution can be parameterized in a power-like way

$$F_{\gamma N}^{(\psi )}(x)=\beta _\gamma x^{\beta _\gamma },$$
(9)
with $`\beta _\gamma 6`$. We choose the value $`\beta _\gamma =6`$.

As to the function $`G(x)`$, it is natural to suggest its similarity to the function $`F_{\gamma N}^{(\psi )}(x)`$. For simplicity we choose exactly the same form

$$G(x)=\beta x^\beta ,$$
(10)
with $`\beta `$ considered as a free parameter.

Terms of the series (6) in number of rescatterings $`k`$ decrease fastly with $`k`$, and it is enough to take into account several first terms ($`4÷5`$). The corresponding convolutions at fixed impact parameter $`b`$ are estimated analytically, and the $`b`$ integration was performed numerically. We use for calculations the Woods-Saxon parameterization of nuclear density with standard values of parameters taken from ref..

Before comparison of theoretical predictions with experimental data let us discuss firstly quantitative character of absorptive corrections. As we discussed above at small $`x0`$ we have AGK cancellation, and the extent of its violation is proportional to $`1ϵ`$, a probability to lose the state $`X`$ after inelastic rescattering. The particle can be lost for registration not only if it disappears after interaction but also if it loses part of its momentum and is shifted into different region in $`x`$. As a result at large $`x`$ values the AGK rules are violated and absorption effects can be observed (momentum conservation effect). However it is only redistribution (not a loss) of particles for the whole $`x`$ region (see Fig.1 for illustration of this fact). For cross section integrated over $`x`$ screening (absorption) effects are connected only with $`1ϵ`$ value but not with momentum conservation. For $`J/\psi `$ photoproduction the absorption effects were usually analyzed only for integrated cross section and this is a reason for rather small values of effective cross sections extracted for $`J/\psi `$ absorption ($`3÷5mb`$). We want to demonstrate here 555 First this phenomenon was discussed in the talk  based on unpublished results by K. Boreskov and A. Kaidalov. that as a function of $`x`$ nuclear effects can be rather big (with effective cross section of order of $`20÷30mb`$) but integrally they cancelled because of cancelling contributions from screening and anti-screening regions.

This is demonstrated at Fig.2 where we show description of the NMC data  on ratio $`R(Sn/Be)`$ of $`J/\psi `$ electroproduction cross sections as a function of $`xx_F`$. Since in this case the $`J/\psi `$ spectrum is concentrated at $`x`$ close to 1, the anti-screening regime takes place already at $`x0.8`$. Note that theoretical predictions at small $`x`$ value ($`x<0.6`$) are less reliable due to their sensitivity to small-$`x`$ parameterization of functions $`F_{\gamma N}^{(\psi )}`$ and $`G(x)`$.

For $`J/\psi `$ hadroproduction the anti-screening regime can be expected (see Fig.1b) only for negative $`x`$. Fig.3 and Fig.6 show that for $`J/\psi `$ and $`\psi ^{}`$ hadroproduction the change of regime happens at $`x(0.3÷0.5)`$. This value depends on parameters of the model – the less are momentum losses (the larger is $`\beta `$) and the less is $`ϵ`$ the more negative $`x`$’s are necessary for the anti-screening onset (see Fig.3). For example, for the value of $`\beta =4`$ the crossover happens already at $`x0.2`$.

It is important to note that it is not entirely adequate to analyze such delicate effects as variation of dependence on atomic number with change of $`x`$ or $`p_T`$ in terms of the standard parameterizations $`A^\alpha `$. In general, $`\alpha `$ is not only a function of $`x`$ and $`p_T`$ but depends also on atomic number $`A`$. To illustrate this fact we show at Fig.4 the ratios of spectra $`R(A/N)`$ versus $`A`$ at fixed $`x_\psi `$ which by no means are described with a power dependence. The same effect is demonstrated in Fig.3: there is a noticeable difference of solid and dashed theoretical curves which correspond to different choices of $`A`$ ranges. Therefore the $`\alpha `$ value extracted from experiments depends significantly on a chosen range of $`A`$ and on a way of analysis. To exclude this uncertainty it is preferable to compare with theory ratios of cross sections for different nuclei. Description of such data from E772 experiment is shown in fig.5.

A difference in nuclear effects in $`J/\psi `$ and $`\psi ^{}`$ hadroproduction should be connected in the model with presence of several primary states, say $`X`$ and $`X^{}`$, which interact with nuclear matter in different ways and have different branchings into $`J/\psi `$ and $`\psi ^{}`$. In Fig.6 we give as an example a comparison of data on $`\psi ^{}`$ hadroproduction from refs  with model calculation with a single $`X`$ state at the same parameter values as in fig.2. The two-channel generalization is evident.

Let us mention that in some theoretical models  used for description of $`J/\psi `$, $`\psi ^{}`$ hadroproduction in the negative $`x`$ region the function $`\alpha (x)`$ decreases (screening effects increase) as $`|x|`$ increases in this region.

Thus, an experimental investigation of nuclear effects for charmonium production in the whole region of $`x`$ can provide a valuable information for testing different dynamical models. In the reggeon approach one can distinguish three different regions in rapidity or in Feynman variable for $`J/\psi `$ hadroproduction – in central region one has strong cancellation of screening diagrams and $`\alpha `$ is close to one (small absorption cross section); at $`x>0.2`$ due to violation of cancellation rule the absorption effects increase and effective cross section is about $`20÷30mb`$; at negative $`x`$ less than $`(0.3÷0.4)`$ the anti-screening regime should be observed. Data on $`A`$ dependence of charmonia production in the negative $`x`$ region will be available soon from HERA-B experiment .
Acknowledgements. We are grateful to A. Capella and M. Danilov for stimulating discussions. This work has been partially supported by RFBR grants 01-02-17383 and 00-15-96786, INTAS grant 00-00366 and DFG grant 436 RUS 113/721/0-1.