text:
# Nonlinear QCD Evolution: Saturation without Unitarization

## Abstract

We consider the perturbative description of saturation based on the nonlinear QCD evolution equation of Balitsky and Kovchegov (BK). Although the nonlinear corrections lead to saturation of the scattering amplitude locally in impact parameter space, we show that they do not unitarize the total cross section. The total cross section for the scattering of a strongly interacting probe on a hadronic target is found to grow exponentially with rapidity $`t=ln(s/s_0)`$, $`\sigma \mathrm{exp}\{\frac{\alpha _sN_c}{2\pi }ϵt\}`$ where $`ϵ`$ is a number of order unity. The origin of this violation of unitarity is the presence of long range Coulomb fields away from the saturation region. The growth of these fields with rapidity is not tempered by the nonlinearity of the BK equation.

preprint: nucl-th/0003???

Understanding the growth of total scattering cross sections with energy $`\sqrt{s}`$ is a longstanding problem. The unitarity, or Froissart bound states that the total inelastic cross section for the scattering of a hadronic projectile on a hadronic target can not grow faster than

$$\sigma <\pi d^2\mathrm{ln}^2(s/s_0),$$
(1)
where $`d`$ is some typical hadronic scale and $`t=ln(s/s_0)`$ the rapidity. While QCD, the theory of hadronic interactions, is a unitary theory and therefore satisfies this unitarity bound, there is no guarantee that perturbative calculations preserve this property. In fact, the linear perturbative evolution equation due to Balitsky, Fadin, Kuraeev and Lipatov (BFKL) implies an exponential growth of $`\sigma `$ with $`t`$, thus violating unitarity.

Following the pioneering works of , there has been recent progress in high energy hadronic scattering in a derivation of a nonlinear evolution equation  which tames the BFKL-type growth. These equations resum the nonlinear corrections to the QCD evolution with rapidity (energy) to all orders in partonic density and to first order in the QCD coupling. While previous studies of these equations assumed translational invariance in the impact parameter plane, here we explore for the first time their impact parameter dependence. For the equations first derived by Balitsky , we show that the total cross section does not unitarize but grows exponentially with $`t`$.

The first in the hierarchy of nonlinear BK evolution equations which govern the evolution of correlation functions of the gluon fields in the target with $`t`$, reads

$`{\displaystyle \frac{d}{dt}}\mathrm{Tr}1U^{}(x)U(y)={\displaystyle \frac{\alpha _s}{2\pi ^2}}{\displaystyle d^2z\frac{(xy)^2}{(xz)^2(yz)^2}}`$ (2)
$`N_c\mathrm{Tr}[U^{}(x)U(y)]\mathrm{Tr}[U^{}(x)U(z)]\mathrm{Tr}[U^{}(z)U(y)],`$ (3)
where $`U(x)`$ is the eikonal scattering amplitude of a fundamental probe on the target characterized by some distribution of gluon fields $`A_\mu `$. In the gauge used in

$$U(x)=P\mathrm{exp}\left[i𝑑x^{}T^aA_a^+(x)\right].$$
(4)
The averaging in (3) is taken over the ensemble characterizing the target. In the large $`N_c`$ limit, (3) simplifies to a closed equation for the scattering probability $`N(x,y)`$ of a colour singlet dipole with charges at points $`x,y`$

$$N(x,y)=\frac{1}{N_c}\mathrm{Tr}1U^{}(x)U(y).$$
(5)
This equation was independently derived by Kovchegov  in the colour dipole approach of

$`{\displaystyle \frac{d}{dt}}N(x,y)={\displaystyle \frac{\alpha _sN_c}{2\pi ^2}}{\displaystyle d^2z\frac{(xy)^2}{(xz)^2(yz)^2}}`$ (6)
$`[N(x,z)+N(y,z)N(x,y)N(x,z)N(z,y)].`$ (7)
Weigert  suceeded to reformulate Balitsky’s hierarchy in terms of a nonlinear stochastic process,

$`{\displaystyle \frac{dU(x)}{dt}}`$ $`=`$ $`gU(x)iT^a{\displaystyle \frac{d^2z}{\sqrt{4\pi ^3}}\frac{(xz)_i}{(xz)^2}\left[1\stackrel{~}{U}^{}(x)\stackrel{~}{U}(z)\right]^{ab}\xi _i^b(z)}`$ (9)
$`{\displaystyle \frac{i\alpha _s}{2\pi ^2}}{\displaystyle d^2z\frac{1}{(xz)^2}\mathrm{Tr}[T^a\stackrel{~}{U}^{}(x)\stackrel{~}{U}(z)]}.`$
Here $`U(x)`$ and $`\stackrel{~}{U}(x)`$ are the unitary matrices (4) in the fundamental and adjoint representations, respectively. The noise is characterised by Gaussian local correlations

$$\xi _i^a(t^{},z^{})\xi _j^b(t^{\prime \prime },z^{\prime \prime })=\delta ^{ab}\delta _{ij}\delta (t^{}t^{\prime \prime })\delta (z^{}z^{\prime \prime }).$$
(10)
This simple Langevin equation gives rise to an infinite number of equations for correlators of $`U`$ which coincide with those derived in .

Nonlinear evolution equations were also derived in the Wilson renormalization group approach . While results from this approach coincide with (3,7) in a certain limit, the question whether they are generally equivalent or incorporate different physics is still open. Here, we focus entirely on the BK equations (3,7,9).

From the first numerical  and analytical  studies of the BK eqs. (3,7) the following consistent picture emerges: Suppose one starts the evolution from the initial condition of small target fields (or $`N(x,y)1`$ for all $`x,y`$). Then initially the evolution follows the BFKL equation, since the nonlinear term in (7) is negligible. As the scattering probability approaches unity, the nonlinear term kicks in and eventually the growth stops as the RHS of (7) vanishes for $`N(x,y)=1`$. The larger dipoles \[$`(xy)^21/Q_s^2(t)`$\] saturate earlier, the smaller dipoles follow at later ”time” $`t`$. These features are contained in the simple parametrization

$$N(x,y)=1\mathrm{exp}\left[(xy)^2Q_s^2(t)\right].$$
(11)
The saturation momentum $`Q_s(t)`$ is a growing function of rapidity. Its exact dependence on rapidity is not known, but both, the numerical results  and simple theoretical estimates  are consistent with the exponential growth of the form

$$Q_s(t)=\mathrm{\Lambda }\mathrm{exp}\left[\alpha _sct\right]$$
(12)
with $`c`$ of order unity. This physical picture has been anticipated several years ago in .

While the BFKL equation leads to an unphysical exponential growth of the scattering probability $`N(x,y)`$ with $`t`$, the nonlinearities of eqs. (3,7) tame this growth such that $`N(x,y)<1`$, as required for a probability. This is commonly refered to as “unitarization”. However, the “saturation” of the scattering probability at fixed impact parameter does not insure that the total scattering cross section is unitary. We hence refer to the above phenomenon more accurately as ”saturation”.

To calculate the total inelastic cross section one has to integrate the scattering probability over the impact parameter. Thus in the saturation regime

$$\sigma =\pi R^2(t),$$
(13)
where $`R(t)`$ is the size of the region in the transverse plain for which the scattering probability for hadronic size ”dipoles” is unity. To satisfy the Froissart bound the radius $`R(t)`$ should grow at most linearly with $`t`$. We now present two simple calculations which establish that within the BK evolution the growth of the radius with rapidity is exponential.

First consider the Langevin equation (9). Assume that initially, at rapidity $`t_0`$ the target is black within some radius $`R_0`$. This means that for $`|z|<R_0`$ the matrix $`U(z)`$ fluctuates very strongly so that it covers the whole group space. We concentrate on a point $`x`$ which is initially outside of this black region. The matrix $`U(x)`$ then is close to unity. Thus there is no correlation between $`U(x)`$ and $`U(z)`$, and the second term on the right hand side of (9) can be set to zero. This is the random phase approximation introduced in  and used later in . As the target field ensemble evolves in rapidity, the radius of the black region grows. As long as the point $`x`$ stays outside the black region we can approximate the Langevin equation by (we drop colour indices which are inessential to our argument)

$$\frac{d}{dt}U(x)=\sqrt{\frac{\alpha _sN_c}{\pi ^2}}_{|z|<R}d^2z\frac{(xz)_i}{(xz)^2}\xi _i(z).$$
(14)
This equation neglects contributions to the derivative of $`U`$ that come from gluons originating from the sources outside the black region. Those contributions would enhance the growth of $`U`$, and so by omitting them we underestimate the rate of growth of the radius of the black region. The formal solution of eq. (14) is

$$1U(x,t)=\sqrt{\frac{\alpha _sN_c}{\pi ^2}}_{t_0}^t𝑑\tau _{|z|<R(\tau )}d^2z\frac{(xz)_i}{(xz)^2}\xi _i(z).$$
(15)
Squaring it and averaging over the noise term gives

$$(1U(x,t))^2=\frac{\alpha _sN_c}{\pi ^2}_{t_0}^t𝑑\tau _{|z|<R(\tau )}\frac{d^2z}{(xz)^2}.$$
(16)
As long as $`x`$ is outside the black region we can approximate the integral on the right hand side by

$$_{|z|<R(\tau )}d^2z\frac{1}{(xz)^2}=\pi \frac{R^2(\tau )}{x^2},$$
(17)
and eq.(16) becomes

$$(1U(x,t))^2=\frac{\alpha _sN_c}{\pi }\frac{1}{x^2}_{t_0}^t𝑑\tau R^2(\tau ).$$
(18)
Now as the black region grows, eventually it reaches the point $`x`$. At this rapidity the matrix $`U(x)`$ will start fluctuating with the amplitude of order one. Thus when $`R(t)=|x|`$, the left hand side of eq. (18) becomes a number of order one, which we call $`1/ϵ`$. We thus have an approximate equation for $`R(t)`$

$$\frac{1}{ϵ}R^2(t)=\frac{\alpha _sN_c}{\pi }_{t_0}^t𝑑\tau R^2(\tau ).$$
(19)
At large rapidities therefore the radius of the black region is exponentially large

$$R(t)=R(t_0)\mathrm{exp}\left[\frac{\alpha _sN_c}{2\pi }ϵ(tt_0)\right].$$
(20)
This is our main result.

We note that while the approximations leading to eq. (19) cease to be valid when the point $`x`$ is on the boundary of the black region, this does not affect our main conclusion. First, eq. (17) is an underestimate of the integral, thus underestimating the growth of $`R`$. Second, when $`x`$ is on the boundary of the black region and $`z`$ in the black region, although the factors $`(1U(x)U^{}(z))`$ and $`U(x)`$ in eq.(9) are not strictly unity, they are still of order one for allmost all points $`z`$. Thus, although we can not determine the exact numerical value of $`ϵ`$, the functional form of the solution as well as its parametric dependence is given correctly by eq.(20).

Note that eq. (9) refers to the evolution of matrices $`U(x)`$. This can be thought of as evolution of the scattering amplitude of a coloured probe. The preceding derivation thus refers to the growth with rapidity of the cross section for scattering of a coloured probe. In a confining theory this is not a physical quantity. However the physics of BK equation does not incorporate effects of confinement, and therefore the cross section for a colourless dipole within the BK framework must grow in the same way. To establish this point, and to make more explicit the relation between $`ϵ`$ and the BFKL dynamics, we now present an alternative derivation of eq. (20).

To this end we consider the BK evolution as the evolution of the projectile . Suppose at the initial energy the projectile is a colour dipole of size $`x_0`$. It scatters on a hadronic target of some size $`R_{\mathrm{target}}`$. As is explicit in , as the energy is increased the projectile wave function evolves according to the BFKL equation. Thus at rapidity $`t`$ the density of dipoles of size $`x`$ at transverse distance $`r`$ from the original dipole is given by the BFKL expression (see for example ):

$$n(x_0,x,r,t)=\frac{32}{x^2}\frac{\mathrm{ln}\frac{16r^2}{x_0x}}{(\pi a^2t)^{3/2}}\mathrm{exp}\left[\omega t\mathrm{ln}\frac{16r^2}{x_0x}\frac{\mathrm{ln}^2\frac{16r^2}{x_0x}}{a^2t}\right]$$
(21)
with $`\omega =4\mathrm{ln}2N_c\alpha _s/\pi `$ and $`a^2=14\zeta (3)N_c\alpha _s/\pi `$. When the density of dipoles at a given impact parameter is greater than one, multiple scatterings become important. Thus the scattering probability is not proportional to $`n`$, but is an infinite series containing all multiple scattering terms .

For our argument, it is only important that once the density of dipoles at some impact parameter $`r`$ becomes larger than some fixed critical number, the scattering amplitude at this impact parameter saturates. The exact value of this number depends on the target, but importantly it does not depend on rapidity. Thus the total cross section is given by the square of the largest impact parameter at which the dipole density in the projectile wave function is of order unity. In order to estimate this directly from eq. (21), we must choose the dipole size $`x`$ in (21) to be the smallest size which is saturated on the target at initial rapidity. Within the ansatz of eq. (11) this would be $`x=Q_s^1(t_0)`$. We then find

$`R^2(t)`$ $`=`$ $`{\displaystyle \frac{1}{16}}{\displaystyle \frac{x_0}{Q_s(t_0)}}\mathrm{exp}[{\displaystyle \frac{\alpha _sN_c}{\pi }}ϵt],`$ (22)
$`ϵ`$ $`=`$ $`7\zeta (3)\left[1+\sqrt{1+8\mathrm{ln}2/7\zeta (3)}\right].`$ (23)
Thus, also for a coulour singlet projectile, we arrive again at the exponential growth of the cross section.

The exact value of $`ϵ`$ given in eq. (23) should not be taken too seriously. The explicit form of the dipole density eq. (21) was derived by a saddle point integration, and as such is valid only for $`\mathrm{ln}\frac{16r^2}{x_0x}<\alpha _st`$. This condition is not satisfied by eq. (22). However, even beyond the saddle point approximation the density has the form

$$n(x_0,x,r,t)=\frac{1}{x^2}\mathrm{ln}\frac{16r^2}{x_0x}\mathrm{exp}\left[\alpha _stF(\frac{\mathrm{ln}\frac{16r^2}{x_0x}}{\alpha _st})\right].$$
(24)
The relevant condition is $`F=0`$. Thus, while our calculation does not specify the numerical value of $`ϵ`$, the correct solution parametrically is the same as (22).

In the target rest frame, this violation of unitarity by the BK evolution can be understood as follows: Start with a single dipole scattering on the hadronic target of transverse size $`R_{\mathrm{target}}`$. With increasing energy the projectile dipole emits additional dipoles strictly according to the BFKL evolution. The density as well as the transverse size of the projectile state thus grows. The increase in density leads to increasing importance of multiple scatterings which are properly accounted for in the BK derivation. This ensures that the scattering probability saturates locally. In the saturation regime, as long as the size of the projectile state $`R(t)`$ is smaller than the target size $`R_{\mathrm{target}}`$, the cross section grows essentially due to surface effects,

$$\sigma =\pi R_{\mathrm{target}}^2+2\pi R_{\mathrm{target}}x_0\mathrm{exp}\left[\frac{\alpha _sN_c}{2\pi }ϵt\right].$$
(25)
Thus as long as $`\alpha _sϵt<\mathrm{ln}\frac{R_{\mathrm{target}}}{x_0}`$, the cross section is practically geometrical. However once the energy is high enough so that the projectile size is larger than that of the target, the total cross section is determined by the former and grows exponentially with the logarithm of energy according to eq. (22).

This also illustrates that the applicability of the BK evolution crucially depends on the nature of the target. If the target is thick enough, so that the multiple scatterings become important before the growth of the projectile radius does, and if the target is wide enough, so that saturation occurs before the projectile radius swells beyond that of the target, then there is an intermediate regime in which the inelastic cross section remains practically constant and equal to $`\pi R_{\mathrm{target}}^2`$. Then BK applies. However, if the target is a nucleon, neither one of these conditions is satisfied. Thus the tainted infrared behaviour of the BFKL evolution of the projectile will show up right away and will invalidate the application of the BK equation.

In order to discuss the violation of unitary from the point of view of the evolution of target fields, we now go back to the stochastic process (9). The RHS of (9) describes the total Coulomb (Weizsäcker-Williams) field at point $`x`$ due to the colour charge sources at points $`z`$. Since the noise is stochastic, the colour sources are completely uncorrelated both in the transverse plain and in rapidity. For this random source, the square of the total colour charge is proportional to the area, and this is precisely the factor $`R^2`$ in eq. (18). The incoming dipole thus scatters on the Coulomb field created by the large incoherent colour charge. Because the Coulomb field is long range, the whole bulk of the region populated by the sources contributes to the evolution and leads to rapid growth of $`R`$. If the field created by the sources was screened by some mass, the evolution would be unitary. To illustrate this point, we substitute the Coulomb field $`(xz)_i/(xz)^2`$ in (eq.9) by an exponentially decaying field $`m\mathrm{exp}\{m|xz|\}`$. It is straightforward to perform now the same analysis as before. Eq. (17) is replaced by

$$_{|z|<R(\tau )}d^2zm^2\mathrm{exp}\{m|xz|\}=\mathrm{exp}\{m|xR|\}.$$
(26)
This leads to the substitution $`R^2\mathrm{exp}\{mR\}`$ in all subsequent equations with the end result that

$$R(t)=\alpha _s\frac{ϵ}{m}t,$$
(27)
which in fact saturates the Froissart bound. Thus the reason for the violation of unitarity is that the evolution is driven by the emission of the long range Coulomb field from a large number of incoherent colour sources in the target.

Cutting off the Coulomb field is not the only possibility to cure this problem. Another option is that the sources of the colour charge in the high density regime cease to be incoherent. If they have correlations ensuring that the total colour charge in a region of fixed size $`L`$ is zero, then the incoming dipole would feel the Coulomb field only within the fixed distance $`L`$ from the black region. Thus the new charges produced by the evolution would only ”split off” the edges of the black region rather than from its bulk. This scenario is equivalent to exponential decay of the field, and will lead to a unitary evolution. In a confining theory like QCD, it is likely to be materialized. We note that the desirability of such colour charge correlations was stressed in a somewhat different context in .

Although such charge correlations do not arise in the BK evolution, it is not a priori clear that they are not present in a more complete semiperturbative framework which still does not take into account the physics of confinement at low energies. In fact, the BK framework is incomplete inasmuch as it takes the evolution of the projectile wave function to be pure BFKL. One expects that once the density of gluons in this wave function becomes large, interactions should lead to saturation effects on the wave function level, i.e., the density of the dipoles should grow slower than eq. (21). Such corrections should still be semiperturbative, in the sense that they are present at small $`\alpha _s`$. However, for scattering on ”small” targets, they will become important at the same energy as the multiple scattering terms resummed in eq. (3,7). These wave function saturation effects may lead to charge correlations of the type necessary to unitarize the total cross section.

We finally note that the unitarity bound for Deeply Inelastic Scattering is different. In this case, the projectile is a virtual photon without fixed hadronic size. For transverse polarization, its perturbatively known wave function is $`\mathrm{\Phi }^2(r)\alpha _{em}\frac{1}{r^2}`$ for $`r^2Q^2`$. In such a projectile not all dipoles saturate at the same energy. The main contribution to the scattering probability comes from dipoles of size $`r>Q_s^1(t)`$ which are saturated. At high energy (i.e. for $`Q_s^1Q^1`$),

$$N(\gamma ^{})=_{Q_s^2<r^2<Q^2}d^2r\mathrm{\Phi }^2(r)\alpha _{em}\mathrm{ln}Q_s/Q.$$
(28)
With the exponential dependence (12) of $`Q_s`$ on rapidity this translates into $`N(\gamma ^{})\alpha _{em}\alpha _s\mathrm{ln}s/s_0`$, and therefore

$$\sigma _{\mathrm{DIS}}\alpha _{em}\alpha _s\pi R^2(t)t.$$
(29)
The DIS cross section has an extra power of $`t`$ relative to the cross section of a purely hadronic process. This extra power of $`t`$ is consistent with the numerical results of .

Acknowledgements This work has been supported in part by PPARC. We thank G. Milhano and H. Weigert for helpful discussions. U.A.W. thanks the Department of Mathematics and Statistics, University of Plymouth for hospitality while part of this work was done.


format: text

text:
# Sub-Planck Structure, Decoherence, and Many-Body Environments

## Abstract

In a recent paper \[Nature 412, 712 (2001)\], Zurek has argued that (1) time evolution typically causes chaotic quantum systems to generate structure that varies on the scale of phase-space volume elements of size $`(\mathrm{}^2/A)^d`$, where $`A`$ is a classical action characteristic of the state and $`d`$ is the number of degrees of freedom, and that (2) this structure implies that a small change in a phase-space coordinate $`X`$ by an amount $`\delta X\mathrm{}X/A`$ generically results in an orthogonal state. While we agree with (1), we argue that (2) is not correct if the number of degrees of freedom is small. Our arguments are based on the Berry-Voros ansatz for the structure of energy eigenstates in chaotic systems. We find, however, that (2) becomes valid if the number of degrees of freedom is large. This implies that many-body environments may be crucial for the phenomenon of quantum decoherence.

preprint:   quant-ph/0112139 NSF-ITP-01-184 December 2001
Zurek  has recently pointed out that, contrary to prior expectations , a quantum system with $`d`$ degrees of freedom will generically develop structure on the scale of phase-space volume elements of size $`(\mathrm{}^2/A)^d`$, where $`A`$ is a typical classical action associated with the state. This is diagnosed by considering the Wigner function for the quantum state $`|\psi `$,

$$W(x,p)=\frac{d^ds}{(2\pi \mathrm{})^d}\mathrm{exp}(ips/\mathrm{})\psi ^{}(x+\frac{1}{2}s)\psi (x\frac{1}{2}s),$$
(1)
where $`\psi (x)=x|\psi `$ is the position-space wave function. (We have assumed a pure state; the generalization to mixed states is given in .) By considering some simple superpositions of gaussian wave packets (coherent states), Zurek shows that there can be fluctuations of $`W(x,p)`$ in the position $`x`$ on scales $`\lambda _{x,\mathrm{min}}\mathrm{}/P`$ (where $`P`$ is a typical classical value of the momentum $`p`$), and in $`p`$ on scales $`\lambda _{p,\mathrm{min}}\mathrm{}/L`$ (where $`L`$ is a typical classical value of $`x`$). The corresponding phase-space volume is $`(\lambda _{x,\mathrm{min}}\lambda _{p\mathrm{min}})^d(\mathrm{}^2/A)^d`$, where $`APL`$.

Zurek goes on to consider the possibility that this sub-Planck structure may have important physical significance if the system in question serves as an environment for a qubit (or other two- or few-state system). In particular, he argues that if the two possible states of the qubit cause slightly different displacements in the phase space of the environment, then these two possible environmental states will be nearly orthogonal. This is then tantamount to decoherence of the quantum state of the qubit by the environment. To be specific, consider the action of a unitary displacement operator

$$D(\delta p,\delta x)=\mathrm{exp}[i(\delta p\widehat{x}+\delta x\widehat{p})/\mathrm{}],$$
(2)
where $`\widehat{x}`$ and $`\widehat{p}`$ are the position and momentum operators, on a state $`|\psi `$ of the environment. The claim of  is that

$$\psi |D(\delta p,\delta x)|\psi 0\text{for}\delta x>\mathrm{}/P\text{or}\delta p>\mathrm{}/L.$$
(3)
We will argue, however, that this is not correct in general for systems with few degrees of freedom. On the other hand, we find that eq. (3) is generically correct if the system (which is playing the role of the environment) has many degrees of freedom.

Much of the analysis of  is based on the formula

$$\left|\psi |D(\delta p,\delta x)|\psi \right|^2=(2\pi \mathrm{})^dd^dxd^dpW(x,p)W(x+\delta x,p+\delta p).$$
(4)
However, a somewhat simpler and equally valid formula is

$$\psi |D(\delta p,\delta x)|\psi =d^dxd^dpe^{i(\delta px+\delta xp)/\mathrm{}}W(x,p).$$
(5)
Here we have made use of one of the basic properties of the Wigner function: expectation values of symmetrically ordered operators can be written as simple integrals. Note that all $`x`$’s and $`p`$’s in eq. (5) are now classical integration variables. Also, in the special case $`\delta x=0`$, eq. (5) reduces to eq. (16) of .

From eq. (5), we see immediately that we do not need to have any information about the small-scale structure in $`W(x,p)`$ in order to draw conclusions about the behavior of $`D(\delta p,\delta x)`$ for small $`\delta p`$ and $`\delta x`$. This follows simply from the mathematics of the Fourier transform, which relates small-scale structure in the original variables ($`x`$ and $`p`$) to large-scale structure in the transformed variables ($`\delta p`$ and $`\delta x`$). In fact, according to eq. (5), the implication of the small-scale structure in $`W(x,p)`$ is that $`D(\delta x,\delta p)`$ will be generically nearly zero only for $`\delta x>\mathrm{}/\lambda _{p,\mathrm{min}}L`$ or $`\delta p>\mathrm{}/\lambda _{x,\mathrm{min}}P`$. That is, to get a generically orthogonal state, the shift in coordinates or momenta must be comparable to the classical scales that characterize the original state.

On the other hand, what happens to $`D(\delta p,\delta x)`$ at small $`\delta p`$ or $`\delta x`$ is governed in turn by the behavior of $`W(x,p)`$ for large $`x`$ or $`p`$. If we suppose that the corresponding state has spread out over the available phase space (which is characteristic of chaotic systems, and is explicitly assumed in ), then we can treat the large-scale structure of the Wigner function via the Berry-Voros ansatz

$$W(x,p)\delta (H(x,p)E),$$
(6)
where $`\delta (E)`$ is the Dirac delta function, $`H(x,p)`$ is the hamiltonian, and $`E`$ is the energy (assumed to be classically sharp). This ansatz simply corresponds to the ergodic distribution that covers the constant-energy surface uniformly, according to the Liouville probability measure. It is devoid of the small-scale structure elucidated in , but correctly describes the large-scale structure. When eq. (6) is used in eq. (5), the restriction of the range of integration to a thin shell of phase space will produce ringing in $`D(\delta p,\delta x)`$. The first zeroes will occur at $`\delta x\mathrm{}/P`$ and $`\delta p\mathrm{}/L`$; these are the scales where, according to , $`D(\delta p,\delta x)`$ should become small. This is correct, but then the ringing effect implies that there are recurrent maxima and minima in $`D(\delta p,\delta x)`$, with only a power-law fall-off of the peak heights. (The exponents in these power laws depend on the details of the system; more on this below.) This is contrary to the supposition of  that any recurrences will be strongly suppressed in general (even though they occur in the example of four superposed gaussians that is worked out explicitly). Finally, for $`\delta xL`$ and $`\delta pP`$, the small-scale structure of $`W(x,p)`$ becomes important, and $`D(\delta p,\delta x)`$ goes rapidly to zero.

To illustrate this, we consider a paradigmatic chaotic system, a two-dimensional billiard such as the stadium . In this case $`H(x,p)=\frac{1}{2m}p^2`$ when $`x`$ is inside, and $`H(x,p)=\mathrm{}`$ when $`x`$ is outside. If (for pedagogical simplicity) we take the billiard to be nearly circular with radius $`L`$ (e.g, a stadium with a short straight segment , or a circle with rough walls ), and the bouncing ball to have energy $`E=\frac{1}{2m}P^2`$, then from eqs. (5) and (6) we get

$$D(\delta p,\delta x)=J_0(P|\delta x|/\mathrm{})\frac{J_1(L|\delta p|/\mathrm{})}{L|\delta p|/2\mathrm{}},$$
(7)
where $`J_n(\xi )`$ is an ordinary Bessel function. Here the oscillations and slow decay of the Bessel functions correspond to the ringing phenomenon that we described above. Also, because the Berry-Voros ansatz does not have the correct small-scale structure, we do not see the cutoffs at $`\delta pP`$ and $`\delta xL`$ that would otherwise be there.

We can also go a little further than this for a chaotic system. Some time ago, Berry  conjectured that, in such a system, the energy eigenfunctions (in position space) have the statistical properties of a gaussian-random superposition of plane waves, all with the wavenumber $`k`$ given locally by classical energy conservation and the fundamental quantum relation $`k=P/\mathrm{}`$, but with wave-vectors pointing in random directions. Since then, this conjecture has been supported by a number of numerical experiments , and it underlies the successful characterization of the distribution of conductance peak heights in quantum dots in the Coulomb-blockade regime . If we use an overbar to denote the statistical random-wave average, then Berry’s conjecture can be neatly summarized via

$$\overline{W(x,p)}\delta (H(x,p)E),$$
(8)
plus the statistical statement that the wave-function $`\psi (x)`$ is a gaussian random variable. This in turn can be shown to imply

$$\overline{W(x,p)^2}=\overline{W(x,p)}^{\mathrm{\hspace{0.17em}2}}+O(\mathrm{}^{d1})$$
(9)
for $`d2`$. In this formula, $`\overline{W(x,p)}`$ itself is to be regarded as $`O(\mathrm{}^0)`$. Eq. (9) means that each particular energy eigenfunction will have only small deviations from eq. (8) on large scales. Furthermore, eqs. (8) and (9) also apply for generic time-dependent superpositions of energy eigenstates,

$$|\psi (t)=\underset{\alpha }{}c_\alpha e^{iE_\alpha t/\mathrm{}}|\alpha ,$$
(10)
and for generic mixed states,

$$\rho (t)=\underset{\alpha \beta }{}c_{\alpha \beta }e^{i(E_\alpha E_\beta )t/\mathrm{}}|\alpha \beta |,$$
(11)
provided we assume that different eigenfunctions are statistically independent superpositions of random waves. (Also, we assume that the quantum energy uncertainty is classically small.) Thus, in order to construct a state that is, for example, localized in a particular small region of phase space, we must carefully choose the magnitudes and phases of the $`c_\alpha `$’s in eq. (10) or the $`c_{\alpha \beta }`$’s in eq. (11). Once we have done this, we expect that subsequent time evolution will eventually spread the packet over the corresponding energy surface, its distribution on large scales governed by the classical ergodic distribution. In quantum mechanics, this time evolution is given entirely by the phase changes in eq. (10), and so these must be capable of producing this spreading. This is guaranteed by eqs. (8) and (9). In a many-body system, this same effect corresponds to the approach to thermal equilibrium .

We have shown that the small-scale structure in phase space discovered in  does not generically result in small values of the overlap between a displaced state and the original unless the displacements are classically large. However, this conclusion does not take into account new effects that can arise when $`d`$, the number of degrees of freedom, is large. Essentially, what happens for $`d1`$ is that the power-law exponent for the fall-off of the recurrent peaks in $`D(\delta p,\delta x)`$ becomes large with $`d`$, suppressing all the recurrences, and restoring the conclusions of . We now turn to our exposition of this phenomenon.

For definiteness, we will consider a dilute gas of $`N`$ hard spheres in a three-dimensional box of volume $`L^3`$. This system can also be viewed as a single bouncing ball in a complicated $`3N`$-dimensional chaotic billiard, and so we will assume that Berry’s random-wave conjecture holds for it at sufficiently high total energy $`E`$ . We take $`E=\frac{N}{2m}P^2`$, so that $`P`$ is the typical momentum of a single gas particle. The appropriate generalization of eq. (7) is then

$$D(\delta p,\delta x)=\frac{J_\nu (N^{1/2}P|\delta x|/\mathrm{})}{(N^{1/2}P|\delta x|/2\mathrm{})^\nu \nu !}\underset{i=1}{\overset{3N}{}}\frac{\mathrm{sin}(L\delta p_i/2\mathrm{})}{L\delta p_i/2\mathrm{}},$$
(12)
where $`\nu =(3N2)/2`$. To perform the position-space integral that results in the multiple product, we ignored the spatial regions that are excluded by the requirement that the hard-sphere centers not come closer than the hard-sphere diameter; this should not be an important effect for a dilute gas. In the large-$`N`$ limit, the $`\delta x`$ dependent factor in eq. (12) becomes $`\mathrm{exp}(P^2|\delta x|^2/6\mathrm{}^2)`$. If we assume that the $`\delta p_i`$’s are all of comparable magnitude, and that each is less than $`\mathrm{}/L`$, then the $`\delta p`$ dependent factor becomes $`\mathrm{exp}(L^2|\delta p|^2/24\mathrm{}^2)`$. Thus we have

$$D(\delta p,\delta x)=\mathrm{exp}(P^2|\delta x|^2/6\mathrm{}^2)\mathrm{exp}(L^2|\delta p|^2/24\mathrm{}^2).$$
(13)
From this formula we see that $`D(\delta p,\delta x)`$ goes rapidly to zero for either $`|\delta x|>\mathrm{}/P`$ or $`|\delta p|>\mathrm{}/L`$, which is now in accord with .

Before concluding, we briefly return to systems with few degrees of freedom, and consider the reconciliation of eqs. (4) and (5). Assuming that the state spreads out over the available phase space, we have seen that eq. (5) implies that $`D(\delta p,\delta x)`$ rings, with recurrent zeroes at multiples of $`\delta x\mathrm{}/P`$ and $`\delta p\mathrm{}/L`$. How can this behavior be reproduced by eq. (4)? It is clear that this is possible if and only if $`W(x,p)`$ in fact has the small-scale structure identified in . Thus there is an intimate relation between the large- and small-scale structure of the Wigner function.

To summarize, we find that small displacements $`\delta x>\mathrm{}/P`$ in position and/or $`\delta p>\mathrm{}/L`$ of the quantum state of a many-body chaotic system generically produce a nearly orthogonal state, provided the system has many degrees of freedom. For a chaotic system with few degrees of freedom, we find that the overlap between the displaced and the original state has a series of zeroes spaced by $`\delta x\mathrm{}/P`$ and $`\delta P\mathrm{}/L`$, separated by broad peaks whose heights fall off as a power of $`\delta x`$ or $`\delta p`$, with a final exponential cutoff at $`\delta xL`$ and $`\delta pP`$. These results imply that a many-body system, thought of as an environment, will (not unexpectedly) be more effective in causing a simple system (such as a qubit) to undergo decoherence.

###### Acknowledgements.

We thank Wojciech Zurek for discussions of his results. M.S. thanks the Institute for Theoretical Physics for hospitality. This work was supported in part by the National Science Foundation under Grant Nos. PHY00–98395 and PHY99-07949.


text:
# Comparison of variational real-space representations of the kinetic energy operator

thanks: Corresponding author. Fax: +44 (0)1223 337356

## Abstract

We present a comparison of real-space methods based on regular grids for electronic structure calculations that are designed to have basis set variational properties, using as a reference the conventional method of finite differences (a real-space method that is not variational) and the reciprocal-space plane-wave method which is fully variational. We find that a new definition of the finite differences method \[P. Maragakis, J. Soler, and E. Kaxiras, Phys. Rev. B 64, 193101 (2001)\] satisfies one of the two properties of variational behaviour at the cost of larger errors than the conventional finite differences method. On the other hand, a technique which represents functions in a number of plane-waves which is independent of system size, follows closely the plane-wave method and therefore also the criteria for variational behaviour. Its application is only limited by the requirement of having functions strictly localised in regions of real-space but this is a characteristic of most modern real-space methods as they are designed to have a computational cost that scales linearly with system size.

Electronic structure methods usually require the solution of Schrödinger’s equation

$$\widehat{H}\psi _0=\epsilon _0\psi _0$$
(1)
for the ground state energy $`\epsilon _0`$ and wavefunction $`\psi _0`$ of the Hamiltonian operator $`\widehat{H}`$ of the system. Exact solutions exist only in very simple cases. A common practical way of approximating the solution of (1) is by expanding $`\psi _0`$ in terms of a set of basis functions. It is straightforward to show that this *basis set approach* has two desirable properties:

1. The approximations to $`\epsilon _0`$ are always upper bounds to the exact $`\epsilon _0`$, so any algorithm that minimises $`\epsilon _0`$ with respect to $`\psi _0`$ will yield the optimal solution for a given basis set.
2. The leading term of the error $`\delta \epsilon _0`$ in the energy is proportional to the square of the error $`\delta \psi _0`$ in the wavefunction and thus the energy displays quadratic convergence with increasing basis set size.

These two properties are referred to as the *variational principle* and are characteristic of the basis set approximation. Basis sets in common use for calculations on molecules and solids include Gaussian functions Boys (1950); Davidson and Feller (1986), Slater functions Slater (1930), plane-waves Ihm et al. (1979), spherical Bessel functions Haynes and Payne (1997), pseudo-atomic orbitals Sankey and Niklewski (1989); Soler et al. (2001), wavelets Arias (1999), “blip” cubic splines Hernández et al. (1997), and finite elements Pask et al. (1999).

In recent years, the ability to perform calculations using only quantities that are *localised in real space*, such as the density matrix or Wannier-type orbitals, has led to linear-scaling electronic structure methods that can deal with thousands of atoms Goedecker (1999). As a consequence, the importance of localised basis sets has grown because they are required for the efficient representation of localised functions. Delocalised basis sets such as plane-waves are not suitable as they make the cost of the electronic structure calculations scale as the cube of the system size.

Basis set expansion is by no means the only practical approximation in electronic structure theory. A set of important alternatives are the “real-space” methods Beck (2000); Chelikowsky et al. (1994) which use a real-space grid to express their solutions. The representation of functions directly on a real-space grid, either regular Beck (2000) or adapted to the positions of the atoms Modine et al. (1997), simplifies the application of localisation constraints and for this reason the importance of real-space methods in recent years has grown in parallel with that of the localised basis set methods. Real-space methods based on a regular grid bear some similarity to the plane-wave pseudopotential approach as the spacing of the real-space grid defines a plane-wave kinetic energy cutoff and is the parameter with respect to which the solutions are converged. Most real-space methods are based on the finite differences (FD) approach, so they are not variational, contrary to the plane-wave basis set method. This is due to the fact that the Laplacian operator for the kinetic energy in the Hamiltonian is represented, or rather approximated, as a FD expansion

$$\frac{^2\varphi }{x^2}(x_i,y_j,z_k)\frac{1}{h_x^2}\underset{n=A/2}{\overset{A/2}{}}C_n^{(A)}\varphi (x_i+nh_x,y_j,z_k)$$
(2)
on every grid point $`(x_i,y_j,z_k)`$, where $`h_x`$ is the grid spacing in the x-direction. Using (2) is equivalent to applying the exact Laplacian operator to an approximation of the solution by a polynomial of degree $`A`$ at each point of the grid Fornberg and Sloan (1994). Therefore, the Hamiltonian operator changes as the grid spacing is varied. One consequence of this is the lack of variational behaviour in the solution process, a well known feature of real-space methods Beck (2000).

Maragakis *et al.* Maragakis et al. (2001) correctly recognised that the FD coefficients $`C_n^{(A)}`$ of equation (2) systematically underestimate the kinetic energy. They have redetermined them so that they always slightly overestimate it, provided one assumes that the solution on the grid is the real-space representation of a plane-wave expansion. This in practice means that with their new FD coefficients energies always converge from above as the grid is made finer. So at least one of the desirable features of the variational principle is restored with this amended FD approach.

A completely different approach which attempts to combine the merits of the plane-wave basis set method with the localisation properties of the real-space methods is the FFT box technique Skylaris et al. (2001a); Mostofi et al. (2001); Skylaris et al. (2001b). Here a plane-wave basis is assumed but the functions are kept localised as they are represented on a real-space grid from the outset. A subset of the plane-waves which is proportional to the number of grid points in the regions of the localised functions and independent of the size of the system is used whenever the functions are represented in reciprocal-space. Functions localised in real-space are smooth and delocalised in reciprocal-space, therefore *coarse sampling* in reciprocal-space is a very good approximation Skylaris et al. (2001a). There is no modification of the Hamiltonian operator but an intermediate truncation of the basis set is performed, so the variational behaviour can not be guaranteed but can be expected as a consequence of the close resemblance in practice of the FFT box method to the full plane-wave approach.

To assess the variational properties of the new methods, we will first use the exactly solvable model of the harmonic oscillator that Maragakis *et al*. used for their tests.

As they established in their paper that their modified FD method converges from above, here we examine if the quadratic convergence criterion of the variational principle is also satisfied. We compare their method for $`A`$=12 with the conventional FD $`A`$=12 method and with the plane-wave expansion method which is variational by definition. In Figure 1 we plot the absolute error in the ground state energy as a function of the rms error in the ground state wavefunction. The behaviour of the two FD methods is very similar, their curves in Figure 1 almost coincide and have a slope close to unity or in other words the error in the energy depends only linearly on the error in the wavefunction. On the other hand, the plane-wave curve has a slope close to 2 and becomes parallel with the parabolic line (slope 2) as the error in the wavefunction decreases, so it is truly variational as expected.

The fact that none of the FD methods exhibit the quadratic convergence of variational behaviour may not be a serious drawback for practical calculations since computational limitations always preclude high convergence with respect to grid spacing. However, the speed of convergence is important as it will determine the speed at which relative energy differences converge. To examine this aspect we need to consider an example on a real system. For this purpose, we have performed LDA density-functional calculations on a silane molecule inside a large cubic simulation cell of (40 Å)3. The atomic cores are represented by norm-conserving pseudopotentials Troullier and Martins (1991) and the charge density is expressed in terms of pseudo-atomic orbitals Sankey and Niklewski (1989) centred on the atoms and optimised for each pseudopotential in spherical regions of 8.0 a0 radius. We compare the conventional FD methods of order $`A=4,8`$ and 12 with the Maragakis *et al*. variant of order $`A=12`$, with the FFT box technique, and with the conventional plane-wave basis set approach. In all FD calculations we solve the Poisson equation for the Hartree potential using the conventional FD coefficients, as Maragakis *et al*. do, so the only difference between the two forms of FD is in the calculation of the kinetic energy. Figure 2 shows the convergence of the total energy as a function of the number of grid points in every lattice vector direction.

The convergence of the energy with the conventional FD methods is not monotonic although the effect is minimised as the order increases. At coarse grid spacings, FD converges from above since there are not enough grid points to accurately represent the charge density near the ions which is the major contribution to the total energy. When the grid spacing becomes finer than some threshold that depends on the order of the FD, the pseudopotential contribution lowers the energy, but to a value which is less than the grid-converged result beacuse of the systematic underestimation of the kinetic energy. In the case of the Maragakis *et al*. FD, the convergence with grid spacing is monotonic and always from above, as expected and as we have also confirmed in tests with other systems. In terms of the speed of convergence and deviation from the plane-wave curve which is variational by definition, the conventional FDs of order $`A=8`$ and 12 do quite well starting with errors in the order of m$`E_h`$ which decrease fast with the number of points. The convergence of the variational FD is not as good on the other hand as the size of its errors is larger and its rate of convergence is slower. The FFT box technique closely follows the plane-wave curve as it is expected to do and has errors that are at least an order of magnitude smaller than the FD $`A=12`$ curve throughout its range, so it apears to coincide with the plane-wave curve on the scale of this graph.

Some more conclusions on the behaviour of the three methods can be derived from Figure 3. Here we plot the deviation of the three best cases of Figure 2 (conventional and variational FD of $`A=12`$ and the FFT box technique) from the plane-wave curve which is variational by definition.

The FFT box curve is closest to the basis set behaviour. The conventional FD converges satisfactorily as long as the grid is dense enough but it is probably safe to say Skylaris et al. (2001a) that no matter how high an order we use, it will always be worse than the FFT box. The FD of Maragakis *et al*. has the desirable property of always converging from above but its errors are an order of magnitute larger than the conventional FD of the same order.

In summary, we have carried out a comparison of methods for electronic structure calculations directly in real-space in order to examine if and to what extent they can display the variational behaviour of basis sets. The FD representation of the Laplacian for the kinetic energy operator systematically underestimates the kinetic energy and is the source of the problems as correctly identified by Maragakis *et al*. Their proposed modification of FD coefficients systematically overestimates the kinetic energy and causes the energy to converge from above albeit with errors up to an order of magnitude larger than conventional FD. In practice this means that finer grids may be necessary and this could be a high price to pay in terms of computational cost. None of the FD methods display quadratic convergence of the energy with respect to the wavefunction which is the second characteristic of the variational principle. For cases where the functions are strictly localised in regions of the real-space grid, as with most real-space methods intended for linear-scaling calculations Fattebert and Bernholc (2000); Soler et al. (2001), the FFT box technique could be the method of choice as it is designed to closely follow the behaviour of the fully variational plane-wave method.

C.K.S. would like to thank the EPSRC (grant number GR/M75525) for postdoctoral research funding. O.D. would like to thank the European Commission for a Research Training Network fellowship (grant number HPRN-CT-2000-00154). P.D.H. would like to thank Magdalene College, Cambridge, for a research fellowship.


text:
# Temporal fractal structures: Origin of power-laws in the world-wide Web

## I Introduction

Complex evolving networks differ from static random graphs in that their size increases in time, thus impacting the linking process in a nontrivial manner. Hence the emergent structure of links is related to salient features of the growth processes in the network, which is governed by the individual linking rules . In the case of the scale-free structures emerging at large evolution times in networks with preferential linking, the universality classes characterized by the same set of scaling exponents can be distinguished, that are based on several relevant details of the microscopic linking properties. The theoretical background of the universality classes of dynamic networks is still missing. By study of many particular networks it has been recognized that certain dynamic constraints on the linking processes can change the emergent scale-free behavior .

The world-wide Web differs from the generic scale-free evolving networks in two ways: (1) it is represented by a directed graph and, (2) most importantly, it has a variable wiring diagram. Frequent updates of the out-going links, that are peculiar for conduct of the agents in the real Web, makes the wiring diagram of the Web graph changing at the same paste as the graph grows. Whereas, wiring diagram of some other networks changes on much slower scale or not at all . The intimate relationship between structural and growth properties leads to a specific architecture of links in the Web. Recent measurements in the real Web have shown  that both out- and in-degree distributions are power-law with different exponents, as well as the size of the connected clusters out of the giant component.

Occurrence of power-laws is a remarkable feature in large number of complex evolving networks, that indicates presence of underlying self-organization while the network grows. The emergent hierarchical organization of node ranks is highly relevant for the stability of the network under attacks , and for the character of other dynamic processes on the network, such as random walk processes . Therefore understanding the mechanisms of self-organization that lead to power-laws in the world-wide Web is crucial both for its functional stability and for designing efficient search algorithms  and transport processes  on the Web graph.

As a step towards realistic modeling of the dynamics of world-wide Web we proposed recently the model  which takes into account the basic relevant features of the Web growth: directed linking, rewiring of preexisting links, and bias activity of agents and bias attachment of links. It was shown  that, when the degree of rewiring in the graph is adjusted to $`\beta 3`$ (i.e., to each new added link in the graph there are in the average three updated links among preexisting nodes), the model reproduces fairly well the emergent power-law distributions of the out- and in-degree, and the scaling exponent of the connected components. In the present work we use this model to study the details of the growth process that precedes the emergence of link structure in the Web graph. Particularly, we demonstrate that a spatio-temporal fractal structure of linking activity occurs on the growth time scale by successive addition of nodes and the average increase of the number of links by $`M=\beta +1`$. The fractal properties of the structure are measured by scaling behavior of the distribution of time intervals of the first-return activity to a given node. We show that this activity pattern results in the algebraic increase of the average local connectivity $`<q_\kappa (s,t)>`$ at node $`s`$ with time ($`\kappa `$ refers to “out” and “in” links), implying scaling behavior in the underlying local probability distribution. In view of the scaling theory the local connectivity is then related to the emergent degree distribution in the limit of large evolution times  $`t\mathrm{}`$. We demonstrate these steps by directly simulating the appropriate quantities, both for out- and in-links. In addition, we show that dynamic processes, such as random walks on the Web graph , that use the information on the local connectivity may also result in the power-law distributions. Here we compute (for the same graph parameters) the distributions of distances on node hierarchy made by an ensemble of random walkers which utilize parts of locally available information on the Web structure.

## II Linking Rules and Fractal Growth Patterns

The growth model is defined by the dynamic rules  that can be summarized as follows: At each time step $`t>M_0M`$ add a node $`i=t`$ and create $`M`$ links. A link is first attempted from the new added node with probability $`\alpha `$ to a target node $`k`$ that is selected with the probability $`p_{in}(k,t)`$, specified below. Else, a link is created between a pair of preexisting nodes (updated link) as follows: A link from node $`n<i`$ to target node $`k<i`$ at time $`t=i`$ occurs with the probability

$$C(n,k,t)=(1\alpha )p_{out}(n,t)\times p_{in}(k,t),$$
(1)
where both probability to select an origin of the link $`p_{out}(n,t)`$ and to select a target node $`p_{in}(k,t)`$ depend on the current connectivity of these nodes $`q_{out}(n,t)`$ and $`q_{in}(k,t)`$, respectively,

$$p_{out}(n,t)=\frac{\alpha +q_{out}(n,t)/M}{(1+\alpha )t},p_{in}(k,t)=\frac{\alpha +q_{in}(k,t)/M}{(1+\alpha )t}.$$
(2)
At the moment of addition $`q_{out}(i,i)=q_{in}(i,i)=0`$ and increasing in time. Therefore, a ratio of the number of added and updated links $`\beta (1\alpha )/\alpha `$, which is independent on the actual number of links $`M`$, is the control parameter in the model. For simplicity we keep $`M`$ fixed, assuming that the number of links fluctuate in time around the average value $`M`$. Motivation for the above linking rules are discussed in detail . For $`\beta 3`$ the results of numerical simulations within this model (see Ref. ) agree satisfactorily with the empirical data on the real Web . We would like to stress that the property of rewiring while the graph grows, which is enabled by $`C(n,k,t)>0`$ in Eq. (1), yields qualitatively new scaling features, compared with the graphs with frozen links ($`\beta 0`$). For instance, one of the important consequences of rewiring is the appearance of the scale-free structure of the out-degree distribution. The entire class of graphs generated with the above rules for varying $`\beta `$ in the range $`0<\beta <\mathrm{}`$ was studied recently .

Here we concentrate to the growth phase of the Web graph. We fix $`\beta =3`$ and chose $`M=4`$. In principle, share between added and rewired links is statistical, thus the universal degree-distributions are independent on $`M`$ in the scaling region . However, some local properties can depend on the actual growth rate of the number of links. To show these dependences is another goal of this work.

In Fig. 1 we show the activity pattern of each node during the growth time $`t=1\mathrm{}N`$, with $`N=10^4`$ nodes-steps. A node $`n`$ can be linked only in the time after its appearance $`tn`$. Value on the vertical axis of a point in the plain $`(k,t)`$ represents the moment of time $`t`$ when the node $`k`$ got an in-link (upper-left part), or a node $`n`$ that created an out-link at time $`t`$ (lower-right part). Quantitative characterization of these patterns can be done in two ways.

First, we notice that the time intervals between two consecutive linking activities at a given node are irregular, resembling a fractal set. For this reason we measure the distribution of successive time intervals $`\mathrm{\Delta }t`$ for in-linking and out-linking separately. The results are given in Fig. 2 (left), where the algebraic decay of the distributions $`P_\kappa (\mathrm{\Delta }t)(\mathrm{\Delta }t)^{\mu _\kappa }`$, confirms the fractal character of these sets. The slopes of the two curves are $`\mu _{out}=0.82\pm 0.01`$ and $`\mu _{in}=0.87\pm 0.01`$. Second, fixing a node $`s`$ we watch how the number of links accumulates at that node with time. Averaging over a large ensembles of networks, we find that for $`ts`$

$$<q_\kappa (s,t)>t^{\gamma _\kappa },$$
(3)
where $`\kappa `$ “out” or “in”, as shown in Fig. 2 (right), where the scaling exponents are $`\gamma _{out}=0.66\pm 0.03`$ and $`\gamma _{in}=0.87\pm 0.03`$.

## III Local Connectivity and Emergent Degree Distributions

The average connectivity at a node increases as a power of evolution time for times $`ts`$, which is compatible with scaling behavior of the local probability distribution $`\rho _\kappa (q,s,t)`$ that the node $`s`$ collected $`q`$ links in time up to the step $`t`$. It was shown  analytically that preference linking with the probability $`p_{in}(k,t)`$ given in Eq. (2) leads to the power-law behavior of the $`\rho _{in}(q,s,t)`$ both in $`q`$ and $`t/s`$ arguments. Our results in Fig. 2 suggest that, due to rewiring with the probability $`C(n,k,t)`$ in Eq. (1), the distribution $`\rho _{out}(q,s,t)`$ also exhibits scaling behavior but with different exponent ($`\gamma _{out}\gamma _{in}`$). The emergent degree distributions $`P_\kappa (q)`$ are defined as

$$P_\kappa (q)=\underset{t\mathrm{}}{lim}\underset{s<t}{}\rho _\kappa (q,s,t)q^{\tau _\kappa },$$
(4)
which we extend to both out- and in-links. In addition, the exact scaling relation that applies to in-link distributions  can be easily extended to out-links , i.e.,

$$\tau _\kappa =1/\gamma _\kappa +1.$$
(5)
Here we assume that the general scaling form applies both for in- and out-links

$$\rho _\kappa (q,s,t)(s/t)^\mu f(q^x(s/t)^\mathrm{\Delta }),$$
(6)
with conserved number of links of both kind, i.e., $`_q\rho _\kappa (q,s,t)=1`$ where $`\kappa =`$ “in”, “out”. Then together with Eqs. (3)-(4) we find $`\mu =\mathrm{\Delta }/x=\gamma `$ and $`\tau =(1+\gamma )x/\mathrm{\Delta }`$, leading to Eq. (5). The measured distributions of emergent node ranks $`P_{in}(q)`$ and $`P_{out}(q)`$ after $`N=10^5`$ added nodes are shown in Fig. 3. The slopes $`\tau _{out}1=1.70\pm 0.03`$ and $`\tau _{in}1=1.26\pm 0.02`$ obey the scaling relation (5) with the respective values for $`\gamma _{out}`$ and $`\gamma _{in}`$ taken from average connectivity in Fig. 2.

## IV Walker on the Web: Local Structures

We have demonstrated that evolution of local organization both of out- and in-links at individual nodes is responsible for the global scaling of the emergent node distributions for large times and number of nodes. Here we would like to show that dynamic processes on grown network (i.e., at different time scale) which use the information on these local properties also obey certain scaling laws. Such processes on the Web are different kinds of random walks related, for example, to search algorithms. We define two types of random walks : The adaptive random walk (ARW) that selects the target node $`k_{\mathrm{}}`$ with the weight which is proportional to in-linking probability of the visited node, and a naive random walk (NRW) selecting one of the out links with equal probability. The corresponding weights are

$$w_{ARW}(n,k_{\mathrm{}})=p_{in}(k_{\mathrm{}})/\underset{\mathrm{}=1}{\overset{q_{out}(n)}{}}p_{in}(n,k_{\mathrm{}}),w_{NRW}(n,k_{\mathrm{}})=1/q_{out}(n).$$
(7)
In Fig. 4 we show the distributions of distances in hierarchy levels $`\mathrm{\Delta }q_\kappa `$ inside the clusters of connected nodes which are visited in cumulative time by an ensemble of walkers. As the Fig. 4 shows, these local clusters are organized scale-free structures with $`W(\mathrm{\Delta }q_\kappa )(\mathrm{\Delta }q_\kappa )^{\delta _\kappa }`$. Notably the scaling exponents $`\delta _{out}\delta _{in}`$ expressing the correlations due to normalization of weights in Eq. (7). For the case of ARW the exponents are close to $`\tau _{in}`$ of the global structure of in-links, whereas for the clusters scanned by NRW they are reduced by approximately unity, i.e., $`2.07\pm 0.04`$ and $`1.10\pm 0.03`$ for ARW and NRW, respectively.

## V Conclusions

The probabilistic character of linking rules for $`\alpha <1`$ with self-consistently varying linking probabilities leads to power-law behavior of local and global (emergent) link structures, both for out- and in-links. This numerical results (see also references ) are in agreement with the analytical results obtained by rate equation approach in the same model . We have demonstrated here that the basis of these scaling laws lies in the occurrence of dynamic fractals and hence the algebraic growth of the local connectivity in time, from which then the hierarchical global structure emerges at large times.

In addition, such local connectivity affects random-walk processes on grown networks. The connected clusters (subgraphs) scanned by the random-walk ensembles also exhibit scaling behavior of distances in node degrees. The scaling properties of these subgraphs on the Web strictly depend on the applied random walk strategy, i.e., in the degree of information about local connectivity that the walkers use. Moreover, the scaling exponents of the distributions of distances in these subgraphs on the Web decrease with increasing rate $`M`$, in contrast to the global structure of the graph, which is universal for large network size $`N`$.

## Acknowledgments

This work was supported by the Ministry of Education, Science and Sports of the Republic of Slovenia.


text:
# X-ray Spectrum and Pulsations of the Vela Pulsar

## 1. Introduction

The Vela pulsar (B0833-45), with a period of 89 millisecond and a characteristic age of 11,000 years, it is the prototype of the young “Vela-like” pulsars. The excellent angular resolution of Chandra allows us to separate the pulsar from its synchrotron nebula and study its X-ray properties.

## 2. Spectrum

We extracted a total of 16000 counts in the range 0.25–2.0 keV from two HRC-S/LETG observations of 25 ks each. No deviations of counts in the individual bins from the mean value of neighboring 10–20 bins were found with statistical significance higher than 2.7$`\sigma `$. This may indicate that there are no elements heavier than H present on the NS surface. The spectral fits with thermal models show excess counts at the high energy end.

To get a handle on the harder component, we used the publicly available ACIS/HETG data obtained in the Continuous Clocking (CC) mode. We used only the zero order image, which is equivalent to about 3 ks observation without the grating. We applied manual corrections for the dither and the SIM motion for both timing and the event position in sky plane (see Zavlin et al. 2000 for details). The background was estimated by interpolating between the neighboring regions in the 1D image of the pulsar and its nebula. Since the ACIS response is poorly known at E $`<`$ 0.5 keV, and the background dominates above 8.0 keV in the CC mode, we use the events in the energy range 0.5–8.0 keV.

The combined spectral fit to the HRC-S/LETG and ACIS/HETG-CC data clearly shows two components – a soft thermal component, which fits equally well with a blackbody or a magnetic neutron star atmosphere model, and a harder component, which we interpret as a power-law (PL). The parameters of the fit depend on the model chosen for the soft thermal component. Table 1 shows the parameters of fits to the combined data (see Pavlov et al. 2001 for details) . Figure 1 shows the broadband spectrum of the Vela pulsar. The soft X-ray spectral parameters are for the atmosphere model. The extrapolation of the PL component matches both the optical and the hard X-ray/$`\gamma `$-ray flux. This is the first conclusive observation of both the thermal and the non-thermal radiation from the Vela pulsar.

## 3. Timing

The Vela pulsar was observed twice with the HRC-I and twice with the HRC-S/LETG. Both of the HRC-I observations and the first LETG observation suffered from the HRC timing problem. We corrected the times by shifting to the time of the next event in the level 1 event list taking only the events with time correction less than 4 ms. The second LETG observation was in the timing mode, so we could recover all the event times.

We used the radio ephemeris of the Vela pulsar based on the observations at the Hobart and Parkes radio observatories. The phase of each event was determined and folded to get the pulse profiles. The left panel of Figure 2 shows the observed pulse profiles of the HRC observations. The pulse profiles show three peaks separated by about 1/3 of phase, with a total pulsed fraction of $`9\%`$ (see also Helfand et al. 2001). The pulsed fluxes in the three peaks are shown in Figure 1.

For the ACIS/HETG-CC observation, the recorded times are the event readout times. We corrected the times for the dither and SIM motion, and the delay from the event to its readout. The correction for absolute time depends on the actual position of the source on the CCD, which is poorly known at present. Therefore, we are unable to get the absolute times for the CC observation. The relative times, on the other hand, are accurate. The right panel of Figure 2 shows the ACIS pulse profiles in two energy bands ($`E<1.8`$ keV and $`E>1.8`$ keV). The difference in the profile shapes and the pulsed fraction is striking. We estimate the intrinsic pulsed fraction of the pulsar to be $`8\pm 2\%`$ and $`62\pm 20\%`$ in the $`E<1.8`$ keV and $`E>1.8`$ keV bands, respectively.

## 4. Summary

* The HRC-S/LETG spectrum shows no conclusive evidence for lines in the 0.25–2.0 keV range. The spectrum shows two components (thermal + non-thermal), similar to middle-aged pulsars. The thermal component is consistent with the emission from a magnetic Hydrogen atmosphere. Extrapolation of the PL component matches the optical and the hard X-ray points.
* The thermal luminosity of the Vela pulsar is lower than that predicted by the “standard” neutron star cooling models (Tsuruta 1998).
* HRC data show three peaks with total pulsed fraction of about $`9\%`$. Energy-resolved ACIS pulse profiles show energy-dependent shape and pulsed fraction. The pulse profile at $`E<1.8`$ keV is similar to that observed with HRC. At $`E>1.8`$ keV, the pulse profile shows only two peaks with estimated intrinsic pulsed fraction of $`62\pm 20\%`$.

### Acknowledgments.

This work was supproted by NASA grants NAG5-7017, NAG5-10865, NAS8-3852, and SAO grant GO1-2071X.

## References

Helfand, D. J., Gotthelf, E. V., Halpern, J. P. 2001, ApJ, 556, 380

Kanbach, G. et al. 1994, A&A, 289, 855

Nasuti F. P. et al. 1997, A&A, 323, 839

Pavlov, G. G. et al. 2001, ApJ, 552, L129

Schönfelder, V. et al. 2000, A&AS, 143, 145

Strickman, M. S. et al. 1996, ApJ, 460, 735

Tsuruta, S. 1998, Phys. Rev., 292, 1

Zavlin, V. E., Pavlov, G. G., Sanwal, D. & Trümper, J. 2000, ApJ, 540 L25


text:
# Abstract

## Abstract

I discuss the remit of Starlink’s software strategy groups and a particular item on the agenda of a meeting of the image processing software strategy group held on 26th January 2001: ‘Why don’t people use Starlink software?’. The answer to this question was speculated to be primarily a sociological effect: those people supervising UK astronomy Ph.D. students are largely people who had learnt their trade at a time when Starlink had a less than perfect reputation. I report on the recommendations made to Starlink to counter this effect.

## Introduction

In the United Kingdom, the Starlink Project (see http://www.starlink.rl.ac.uk) has provided interactive data processing facilities for astronomers for over 18 years. Starlink has provided both hardware and software on behalf of PPARC; ‘the primary purpose …\[being\] …to maximize the return on PPARC’s expenditure on astronomical computing’ (http://www.starlink.rl.ac.uk). In 1999, Starlink’s budget stood at approximately $`2.5`$ million pounds sterling per annum. Out of this sum, approximately $`0.6`$ million pounds is ear-marked for the provision of software. These provisions include the development, distribution and documentation of the software.

The Starlink software collection is a set of software written both by members of the Starlink Project and imported from outside the Project for the use of the UK astronomical community. Many utilities and packages exist, mostly focussing on data reduction but there are also some (limited) theoretical tools.

## Software Strategy Groups

A mainstay of Starlink’s software development is the existence of Software Strategy Groups (SSGs). There are six such SSGs whose remit covers different aspects of astronomical research: spectroscopy, image processing111The Author of this article has been a member of the Image Processing SSG since 1999., information services & databases, radio, mm, & sub-mm astronomy, theory & infrastructure and X-ray astronomy. The SSGs were set up after the 1991/1992 review of Starlink to provide feedback, re-assess software priorities and identify future objectives. Membership of the SSGs is composed of both Starlink’s programmers and expert users drawn from across a broad spectrum of personal experience and opinion. The SSGs are annually asked to advise on strategies to be adopted over in the coming years in order to ‘satisfy (and continue to satisfy) the needs identified by users in the \[Starlink Users\] questionnaire. …They will also be asked to identify in more detail those software projects that should be carried out over the next year in order to implement the strategy’ (http://star-www.rl.ac.uk/star/docs/sgp44.htx/node3.html).

## Sociological Issues

Third on the agenda of the image processing SSG for $`26^{\mathrm{th}}`$ January 2001 was the loaded question ‘Why don’t people use Starlink software?’ This question came as a surprise for a few members of the SSG as many of the users used Starlink’s software on a daily basis in the research.

Starlink software, however, is by no means the only available source of software for astronomers. The major competition for users is the successful Image Reduction and Analysis Facility (IRAF). IRAF is a publicly available ‘general purpose software system for the reduction and analysis of scientific data’ (http://star-www.rl.ac.uk/iraf/web/faq/FAQsec01.html)222IRAF written and supported by the IRAF programming group at the National Optical Astronomy Observatories (NOAO) in Tucson, Arizona.. IRAF overlaps many of the areas that the Starlink software does (e.g. general purpose image reduction software). One may expect that there would be substantial competition for users between the two, however, this is not the case. UK astronomers have *tended* in the past to migrate towards using the IRAF option in preference to Starlink’s software. Contributory factors to this may be that in many observatories (e.g. the Isaac Newton Group) IRAF has historically been much more readily available than Starlink and UK astronomers have used IRAF to be more compatible with their American collaborators. So why do many UK astronomers use IRAF when much of the equivalent Starlink software is arguably better, more stable and more reliable? In answering this question, the image processing SSG expediently pointed out that this was not always the case. During the 1980’s, Starlink software had a less than perfect reputation. But why in 2001 is there still such a trend for astronomers to be preferentially using IRAF packages over Starlink ones? The image processing SSG speculated that the major reason was a sociological one. Post-graduate students who are presently commencing work on their Ph.D. theses will primarily be introduced to the tools of their trade by their Ph.D. supervisors. Their supervisors, on the whole, in turn learnt their trade during, or lived through the aforementioned time when Starlink had a less than perfect reputation. Their supervisors prejudices and perspectives naturally influence their students: in guiding their post-graduate students they are likely to recommend IRAF in preference to Starlink software to their charges. Thus, Starlink software has become less used over time, potentially creating a further problem: a lack of local expert users. Even if new post-graduate students chose to use Starlink software, they would probably have little to no personal guidance available to them whereas if they were to chose IRAF, many people within the given astronomy department would potentially be able to assist.

## Discussion

To remedy this situation, the image processing SSG has recommended a number of changes be implemented. Firstly, to give the software a higher visibility, the Starlink newsletter would be re-introduced to publicize new software developments. More ideas were to implement a method for disseminating what software is available on-line; the introduction of a ‘man’ page; and an education scheme to provide novice and infrequent users reach rapid decisions about what Starlink software could be helpful. Although these recommendations are a start, it is probably going to take some time to wash away the memories of the imperfect reputation that Starlink had in the 1980s and counter the present sociological prejudice against Starlink software. In the immediate future, many new software products will be released. Of particular note is the Interactive Data Language (IDL). The recent rise and increasing use of IDL has shown that quality, new software will always have a place within the astronomical community and Starlink must re-double its efforts to compete.

## Acknowledgements

I sincerely acknowledge the data analysis facilities provided by the Starlink Project which is run by CCLRC on behalf of PPARC, without which my Ph.D. would not be what it is. I thank Peter Draper and Nigel Metcalfe for useful discussions on the factual content of this article and the image processing SSG that inspired this.